I have a shapefile with overlapping polygons and a set of rasters. I would like to summarize some statistics on my polygons. However, as my polygons are overlapping, to avoid issues with zonal_stats, I loop over the polygons of my shapefile. I also would like to loop over the rasters. This is the code I am using:

import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterstats import zonal_stats
from shapely.geometry import shape

raster1=data_directory + '/path/raster1.tif'
raster2=data_directory + '/path/raster2.tif'
raster3=data_directory + '/path/raster3.tif'
RASTS = [raster1, raster2, raster3]

shapefile = gpd.read_file(data_directory + 'path/shapefile.shp')

dfs = []
for rast in RASTS:
    stats_list = []
    for index, row in shapefile.iterrows():
            #Take the geometry of the current polygon
            geometry = row['geometry']
            #Create a GeoDataFrame with the current polygon
            gdf_polygon = gpd.GeoDataFrame(geometry=[geometry], crs=shapefile.crs)
            #Perform zonal_statistis on the current polygon
            stats_results = zonal_stats(gdf_polygon, rast, stats = ['mean'])
            #Append the statistics to the list
            # Create a DataFrame for the current raster's statistics
            df = pd.DataFrame(stats_list)
            # Append the DataFrame to the list

# Concatenate the DataFrames along the columns
polygons_stats = pd.concat([shapefile] + dfs, axis=1)

The desired goal is to obtain my original shapefile with three columns containing the names of the rasters corresponding to the values of the three rasters. However, instead of the additional three columns, I obtain 6.177 columns called "mean". I have no idea of why this is happening and how I can fix it. Any suggestions?

I am using Python 3.9

  • is it possible to share sample of shapefile and few rasters that we can use to test? Commented Feb 23 at 11:04
  • Unfortunately, I can share the rasters but not the shapefile. Not sure how I can add more information to be more helpful.
    – OgeiD
    Commented Feb 23 at 11:15
  • Based on my experiences, overlapping polygons are no problem for zonal stats. Do you have other experiences?
    – Pieter
    Commented Feb 23 at 11:55
  • @Pieter I think it can have problems as documented here desktop.arcgis.com/en/arcmap/latest/tools/…. I also asked a couple of weeks ago how to deal with that gis.stackexchange.com/questions/475689/…
    – OgeiD
    Commented Feb 23 at 11:59
  • 1
    Rasterstats deals differently with this than the zonal statistics tool of ArcMap, so the shortcomings of ArcGIS zonal stats don't apply. The answer to your previous question was also that rasterstats will just process each polygon independently regardless of overlaps, so you don't need to do any looping yourself. Just call zonal_stats once for all your polygons and the necessary looping will be done for you under the hood.
    – Pieter
    Commented Feb 23 at 12:29

1 Answer 1


As also discussed in the comments, rasterstats will correctly calculate zonal statistics for overlapping polygons, so you don't need to do the looping over polygons yourself. Just call zonal_stats once for your collection of polygons.

If the loop over all polygons is removed, the many column issue will disappear as well.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.