11

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?

2

2 Answers 2

14

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'
)
0
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')

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.