3

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?

1 Answer 1

5

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())
5
  • Thanks, this works. Though I am a bit confused with how to handle the output. I tried findings just the "mean", and I see that the output stats is a list, showing {'mean': 0.0} for each row. How can I append the stats output to my geopandas dataframe so that there is a new mean column, showing the mean value for each polygon/feature/row? Jan 23, 2022 at 0:02
  • 1
    gis.stackexchange.com/a/418278/2856
    – user2856
    Jan 23, 2022 at 0:32
  • Got it, thanks, that source showed me how to get that column of means for each polygon appended. Though all zonal stats mean values shown for each polygon are 0.0, which must be another issue I need to work out. Jan 24, 2022 at 22:15
  • @LostinSpatialAnalysis Do you resolved your issue with 0 for each polygon?? Sep 1, 2022 at 4:24
  • 1
    @GeorgeBoldeanu ensure your vector data is in the same CRS as your raster. If you still get the issue, post a new question with self-contained code attempt that demonstrates the issue (and enough details about both the raster and vector data that we can replicate).
    – user2856
    Sep 1, 2022 at 4:30

Your Answer

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

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