Using GeoPandas with rasterstats for zonal statistic gives error

Question

I'm trying to use GeoPandas and zonal statistic from rasterstats together.

I run this code

import geopandas as gpd
from rasterstats import zonal_stats
    
geodf = gpd.read_file("SHP/shape.shp")
zonal_stats(geodf, "raster.tif")

but I get this error:

ParseException: Unknown type: 'FID' ParseException: Unexpected EOF parsing WKB Traceback (most recent call last): File "/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py", line 2820, in run_code exec code_obj in self.user_global_ns, self.user_ns File "", line 1, in zonal_stats(geodf,"raster.tif",stats='mean') File "/usr/local/lib/python2.7/dist-packages/rasterstats/main.py", line 21, in zonal_stats return list(gen_zonal_stats(*args, **kwargs)) File "/usr/local/lib/python2.7/dist-packages/rasterstats/main.py", line 128, in gen_zonal_stats for i, feat in enumerate(features_iter): File "/usr/local/lib/python2.7/dist-packages/rasterstats/io.py", line 107, in features_iter = (parse_feature(x) for x in obj) File "/usr/local/lib/python2.7/dist-packages/rasterstats/io.py", line 70, in parse_feature raise ValueError("Can't parse %s as a geojson Feature object" % obj) ValueError: Can't parse FID as a geojson Feature object

The GeoDataFrame looks like this:

    In[43]: geodf
    Out[42]: 
         FID                                      geometry

    0     0  POLYGON ((7.662965420348624 45.05600185397847,...
    1     1  POLYGON ((7.665210146514162 45.05607590098681,...
    2     2  POLYGON ((7.663353281768154 45.05601113630467,...
    3     3  POLYGON ((7.660229137844505 45.0570570693801, ...
    4     4  POLYGON ((7.660047499767138 45.0572381357898, ...
    ....

Any suggestion how to solve this problem?

Answer 1

For anyone else that stumbles upon this you can do:

To get just one stat:

gdf['mean'] = pd.DataFrame(
    zonal_stats(
        vectors=gdf['geometry'], 
        raster='raster.tif', 
        stats='mean'
    )
)['mean']

To get all computed stats:

gdf = gdf.join(
    pd.DataFrame(
        zonal_stats(
            vectors=gdf['geometry'], 
            raster='raster.tif', 
            stats=['max','min']
        )
    ),
    how='left'
)

Answer 2

For some reason @Taras's solution did not work for me, I got the error:

TypeError: 'list' object is not callable

However, I bypass the problem with the following. Not the best solution, but I got the same result.

origShp = gpd.read_file(r'E:\...\myshapefile.shp')

zonal_stats = zonal_stats(r'E:\...\myshapefile.shp',
            r'E:\...\myRaster.tif',
            stats="min mean max")

df = pd.DataFrame(zonal_stats)
df_redux = pd.concat([df, origShp], axis=1)

cols_at_end = ['min', 'mean', 'max']
final_gf = df_redux[[c for c in df_redux if c not in cols_at_end]
        + [c for c in cols_at_end if c in df_redux]]

Out[328]: 
           min       mean        max
0    -0.035565  17.707176  34.366177
1    -0.156580  16.005182  32.490967
2    -0.256113  13.957398  31.997578
3    -0.032301   7.868770  23.284115
4    -0.172403  14.436392  32.300911

Or more compact:

import geopandas as gpd
from rasterstats import zonal_stats

def ZonalStats(shape, raster, stats):
    # shape - shapefile path
    # raster - raster path
    # stats - stats as list, f.e. 'min mean max' ; 'min'
    # the result is final_gdf as GeoDataFrame

    shape_gdf = gpd.read_file(shape)
    zonalSt = zonal_stats(shape, raster, stats = stats)
    df = pd.DataFrame(zonalSt)
    df_concat = pd.concat([df, shape_gdf], axis=1)
    final_gdf = gpd.GeoDataFrame(df_concat, geometry=df_concat.geometry)
    return final_gdf
    
new_df = ZonalStats(shpfile, myraster,'mean')