Python – Zonal Statistics Spatial Join Raster to Polygon

pythonrasterzonal statistics

  • I have a data frame that consists of 3 polygons
  • I have a raster with a single band

Seems like Zonal statistics as a table in ArcGIS does the job but I am using pure python script on this one; with no access to ArcPy.

how do I join the raster value (min, max, average, etc.) to the 3 polygons?

import geopandas as gpd

gdf = gpd.read_file("polygons.shp")
print("length of gdf: ", len(gdf))
# length of gdf: 3

def zonstats(gdf, raster_path):
    ...
    ...
    ...
    return gdf_joined

Best Answer

You can use rasterstats.zonal_stats

import geopandas as gpd
from rasterstats import zonal_stats

zones = "zones.shp"
values = "values.tif"

gdf = gpd.read_file(zones)
print(gdf.head())

stats = gpd.GeoDataFrame(zonal_stats(gdf, values, stats=["min", "max", "mean"]))
gdf = gdf.join(stats)

print(gdf.head())

Example output:

   id                                           geometry
0   1  POLYGON ((440819.096 3751134.570, 440983.634 3...
1   2  POLYGON ((441147.373 3751255.178, 441265.585 3...
2   3  POLYGON ((441438.110 3751077.860, 441789.549 3...
3   4  POLYGON ((440808.713 3750651.340, 440925.327 3...
4   5  POLYGON ((441174.530 3750450.061, 441371.017 3...
   id                                           geometry  ...    max        mean
0   1  POLYGON ((440819.096 3751134.570, 440983.634 3...  ...  173.0  133.909091
1   2  POLYGON ((441147.373 3751255.178, 441265.585 3...  ...  140.0  120.900000
2   3  POLYGON ((441438.110 3751077.860, 441789.549 3...  ...  165.0  122.714286
3   4  POLYGON ((440808.713 3750651.340, 440925.327 3...  ...  140.0  118.800000
4   5  POLYGON ((441174.530 3750450.061, 441371.017 3...  ...  173.0  122.555556
Related Question