#!/usr/bin/env python
# coding: utf-8

# [![image](https://jupyterlite.rtfd.io/en/latest/_static/badge.svg)](https://demo.leafmap.org/lab/index.html?path=notebooks/70_zonal_stats.ipynb)
# [![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/notebooks/70_zonal_stats.ipynb)
# [![image](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/opengeos/leafmap/HEAD)
# 
# **Calculating zonal statistics - summarizing geospatial raster datasets based on vector geometries**
# 
# Uncomment the following line to install [leafmap](https://leafmap.org) if needed.

# In[ ]:


# %pip install -U leafmap


# In[ ]:


# %pip install -U rasterstats geopandas


# In[ ]:


import leafmap
import geopandas as gpd


# In[ ]:


dsm = "https://opengeos.org/data/elevation/dsm.tif"
hag = "https://opengeos.org/data/elevation/hag.tif"
buildings = "https://raw.githubusercontent.com/opengeos/data/refs/heads/main/elevation/buildings.geojson"


# In[ ]:


m = leafmap.Map()
m.add_cog_layer(dsm, name="DSM", palette="terrain")
m.add_cog_layer(hag, name="Height Above Ground", palette="magma")
m.add_geojson(buildings, layer_name="Buildings")
m


# In[ ]:


gdf = gpd.read_file(buildings)
len(gdf)


# In[ ]:


gdf.head()


# The `leafmap.zonal_stats()` function wraps the [`rasterstats.zonal_stats()`](https://pythonhosted.org/rasterstats/index.html) function and performs reprojection if necessary. 
# 
# By default, the zonal_stats function will return the following [statistics](https://pythonhosted.org/rasterstats/manual.html#statistics):
# 
# * min
# * max
# * mean
# * count
# * 
# Optionally, these statistics are also available.
# 
# * sum
# * std
# * median
# * majority
# * minority
# * unique
# * range
# * nodata

# In[ ]:


stats = leafmap.zonal_stats(gdf, hag, stats=["min", "max", "mean", "count"])
len(stats)


# In[ ]:


stats[:5]


# In[ ]:


stats_geojson = leafmap.zonal_stats(gdf, hag, stats=["mean", "count"], geojson_out=True)
len(stats_geojson)


# In[ ]:


stats_geojson[0]


# In[ ]:


stats_gdf = leafmap.zonal_stats(gdf, hag, stats=["mean", "count"], gdf_out=True)
len(stats_gdf)


# In[ ]:


stats_gdf.head()


# In[ ]:


m = leafmap.Map()
m.add_gdf(stats_gdf, layer_name="Zonal Stats")
m