How to calculate mean value of a raster for each polygon in a shapefile

Question

I wrote a Python script who calculate, for rasters, mean value for each polygon of a shapefile. But I think that it's a "war machine" (it takes many treatment time) and I would like to simplify it.

Is somebody know a simple way to do it?

I'm using geopandas for shapefile manipulations.

Answer 1

You can use the rasterstats package for this.

For example (assuming you have a geopandas.GeoDataFrame called gdf):

from rasterstats import zonal_stats

with rasterio.open("/path/to/raster.tif") as src:
    affine = src.transform
    array = src.read(1)
    df_zonal_stats = pd.DataFrame(zonal_stats(gdf, array, affine=affine))

# adding statistics back to original GeoDataFrame
gdf2 = pd.concat([gdf, df_zonal_stats], axis=1)