Performing a zonal statistics looping over rasters and over polygons

Question

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
            stats_list.append(stats_results[0])
            # Create a DataFrame for the current raster's statistics
            df = pd.DataFrame(stats_list)
            # Append the DataFrame to the list
            dfs.append(df)

# 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

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.