70 zonal stats

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

Uncomment the following line to install leafmap if needed.

In [1]:

# %pip install -U leafmap

# %pip install -U leafmap

In [2]:

# %pip install -U rasterstats geopandas

# %pip install -U rasterstats geopandas

In [3]:

import leafmap
import geopandas as gpd

import leafmap import geopandas as gpd

In [4]:

dsm = "https://open.gishub.org/data/elevation/dsm.tif"
hag = "https://open.gishub.org/data/elevation/hag.tif"
buildings = "https://open.gishub.org/data/elevation/buildings.geojson"

dsm = "https://open.gishub.org/data/elevation/dsm.tif" hag = "https://open.gishub.org/data/elevation/hag.tif" buildings = "https://open.gishub.org/data/elevation/buildings.geojson"

In [5]:

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

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

Out[5]:

In [6]:

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

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

Out[6]:

84

In [7]:

gdf.head()

gdf.head()

Out[7]:

	HEIGHT	SQMETERS	geometry
0	17.530001	3525.113281	MULTIPOLYGON (((-77.05121 38.81215, -77.05132 ...
1	6.620000	51.499035	MULTIPOLYGON (((-77.05339 38.80930, -77.05341 ...
2	5.970000	233.687225	MULTIPOLYGON (((-77.05272 38.81064, -77.05282 ...
3	7.190000	320.350342	MULTIPOLYGON (((-77.05211 38.81024, -77.05213 ...
4	6.340000	85.387062	MULTIPOLYGON (((-77.05285 38.81006, -77.05303 ...

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 [8]:

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

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

Out[8]:

84

In [9]:

stats[:5]

stats[:5]

Out[9]:

[{'min': 10.539999961853027,
  'max': 24.950000762939453,
  'mean': 19.82754578143535,
  'count': 843},
 {'min': 6.880000114440918,
  'max': 10.819999694824219,
  'mean': 8.442221747504341,
  'count': 18},
 {'min': 6.079999923706055,
  'max': 19.959999084472656,
  'mean': 9.612962510850695,
  'count': 54},
 {'min': 7.300000190734863,
  'max': 18.920000076293945,
  'mean': 11.190888129340278,
  'count': 90},
 {'min': 5.929999828338623,
  'max': 18.399999618530273,
  'mean': 12.659524100167411,
  'count': 21}]

In [10]:

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

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

Out[10]:

84

In [11]:

stats_geojson[0]

stats_geojson[0]

Out[11]:

{'id': '0',
 'type': 'Feature',
 'properties': {'HEIGHT': 17.530000686645508,
  'SQMETERS': 3525.11328125,
  'mean': 19.82754578143535,
  'count': 843},
 'geometry': {'type': 'MultiPolygon',
  'coordinates': [(((321904.9249656575, 4297929.359149771),
     (321894.1639320927, 4297864.087817658),
     (321872.0912812228, 4297867.69244786),
     (321879.23566416965, 4297911.170307515),
     (321851.01798099, 4297915.801282184),
     (321845.1891851835, 4297880.621811273),
     (321824.5034280424, 4297884.084395529),
     (321831.04482289474, 4297923.911486235),
     (321828.3629885514, 4297924.415869579),
     (321831.17563400406, 4297941.341652784),
     (321904.9249656575, 4297929.359149771)),)]},
 'bbox': (321824.5034280424,
  4297864.087817658,
  321904.9249656575,
  4297941.341652784)}

In [12]:

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

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

Out[12]:

84

In [13]:

stats_gdf.head()

stats_gdf.head()

Out[13]:

	geometry	HEIGHT	SQMETERS	mean	count
0	MULTIPOLYGON (((-77.05121 38.81215, -77.05132 ...	17.530001	3525.113281	19.827546	843
1	MULTIPOLYGON (((-77.05339 38.80930, -77.05341 ...	6.620000	51.499035	8.442222	18
2	MULTIPOLYGON (((-77.05272 38.81064, -77.05282 ...	5.970000	233.687225	9.612963	54
3	MULTIPOLYGON (((-77.05211 38.81024, -77.05213 ...	7.190000	320.350342	11.190888	90
4	MULTIPOLYGON (((-77.05285 38.81006, -77.05303 ...	6.340000	85.387062	12.659524	21

In [14]:

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

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

Out[14]: