Calculating zonal statistics - summarizing geospatial raster datasets based on vector geometries

Uncomment the following line to install leafmap 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() function and performs reprojection if necessary.

By default, the zonal_stats function will return the following 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