Python Zonal Statistics – Solving Attribute Error in Raster Analysis

gdalpythonrasterrasterstatszonal statistics

I am trying to carry out a simple zonal statistics task in python, and am having trouble getting this to work.

I have a raster layer (in black/grayscale) and a roads polygons layer (in green, actually polygons since I drew buffers around a lines layer), shown below:

enter image description here

The roads layer has many different road "segments", each of which is a separate polygon/feature/row in the shapefile attribute table. What I am trying to do is, for each road segment, find the average pixel value of the underlying raster layer, within that specific road segment polygon, and finally append that new column of average values to the road layer's attribute table. This task works fine in QGIS, using the "Zonal statistics". However, I cannot seem to get this same process to work in python.

I am using the following zonal_stats function from the rasterstats module to try and perform my zonal statistics task with the following, referring to this zonal_stats documentation: https://pythonhosted.org/rasterstats/

This looks like such a simple function to use, with simple, straightforward arguments, and so I am simply running:

zonal_stats("roads.shp", "raster.tif", stats="count min mean max median")

And I get the following error traceback:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_17488/2186187301.py in <module>
      5 roads_buffers = gpd.read_file('roads_buffers.shp')
      6 
----> 7 zonal_stats("roads_buffers.shp", "Fresno_Rasters/Fresno_NOx_PopWeighted.tif", stats="count min mean max median")

~\AppData\Roaming\Python\Python39\site-packages\rasterstats\main.py in zonal_stats(*args, **kwargs)
     29     The only difference is that ``zonal_stats`` will
     30     return a list rather than a generator."""
---> 31     return list(gen_zonal_stats(*args, **kwargs))
     32 
     33 

~\AppData\Roaming\Python\Python39\site-packages\rasterstats\main.py in gen_zonal_stats(vectors, raster, layer, band, nodata, affine, stats, all_touched, categorical, category_map, add_stats, zone_func, raster_out, prefix, geojson_out, boundless, **kwargs)
    151         features_iter = read_features(vectors, layer)
    152         for _, feat in enumerate(features_iter):
--> 153             geom = shape(feat['geometry'])
    154 
    155             if 'Point' in geom.type:

~\AppData\Roaming\Python\Python39\site-packages\shapely\geometry\geo.py in shape(context)
    100     else:
    101         ob = context
--> 102     geom_type = ob.get("type").lower()
    103     if 'coordinates' in ob and _is_coordinates_empty(ob['coordinates']):
    104         return _empty_shape_for_no_coordinates(geom_type)

AttributeError: 'NoneType' object has no attribute 'get'

I have tried searching online, and diagnosing this, but I cannot figure out what AttributeError: 'NoneType' object has no attribute 'get' is telling me, because I do not know what 'get' means or where this function is coming from. How can I fix this issue so that I can get my zonal statistics operation running?

Best Answer

You probably have a null geometry, try filtering them out.

import geopandas as gpd
import pandas as pd
from rasterstats import zonal_stats


shp = 'vector.shp'
ras = 'raster.tif'

gdf = gpd.read_file(shp)

# filter out empty or null geometries
gdf = gdf[~(gdf['geometry'].is_empty | gdf['geometry'].isna())]

# zonal stats
stats = zonal_stats(gdf['geometry'], ras, stats="count min mean max median")

# Join stats to gdf 
df = pd.DataFrame(stats)
df = pd.concat([df, gdf], axis=1)
gdf = gpd.GeoDataFrame(df, geometry=df.geometry)
print(gdf.head())
Related Question