I'm trying to produce zonal statistics (summing pixel values representing population estimates by polygons), roughly following GeoDelta Labs (19:20 onwards. A ValueError "Invalid out_shape, must be 2D" is raised for the zonal_stats function, though.

Below my code and the full error log.


import numpy as np
import geopandas as gpd
import cartopy.crs as ccrs
from cartopy.feature import ShapelyFeature
import rasterio as rio
from rasterio.plot import show
from rasterstats import zonal_stats, gen_zonal_stats
import matplotlib.pyplot as plt

# Define study area to be applied to raster window
states = gpd.read_file('my.shp')
study_area = states[states['FIELD'] == 'Value X']
xmin, ymin, xmax, ymax = study_area.total_bounds

# Load municipal boundaries and filter those falling within study area
municipal_all = gpd.read_file('....shp')
municipal_filter = municipal_all[municipal_all['FIELD'] == 'Value X']

# Load raster to extent of study area
# Source: CIESIN HRSL: https://ciesin.columbia.edu/data/hrsl/#data

with rio.open('...tif') as dataset:

    top, lft = dataset.index(xmin, ymax) 
    bot, rgt = dataset.index(xmax, ymin)
    pop_array = dataset.read(window=((top, bot), (lft, rgt)))

    affine = dataset.transform


    municipal_pop_sum = zonal_stats(municipal_filter, pop_array, affine=affine, stats=['sum'], geojson_out=True)

Error log

import numpy as np...
C:\Users\...\envs\xxx\lib\site-packages\rasterstats\io.py:301: UserWarning: Setting nodata to -999; specify nodata explicitly
  warnings.warn("Setting nodata to -999; specify nodata explicitly")
ValueError                                Traceback (most recent call last)
c:\Users\...\script.py in 
---> 50     municipal_pop_sum = zonal_stats(municipal_filter, pop_array, affine=affine, stats=['sum'], geojson_out=True)
     52 """

~\anaconda3\envs\xxx\lib\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))

~\anaconda3\envs\xxx\lib\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, **kwargs)
    158             # rasterized geometry
--> 159             rv_array = rasterize_geom(geom, like=fsrc, all_touched=all_touched)
    161             # nodata mask

~\anaconda3\envs\xxx\lib\site-packages\rasterstats\utils.py in rasterize_geom(geom, like, all_touched)
     39     """
     40     geoms = [(geom, 1)]
---> 41     rv_array = features.rasterize(
     42         geoms,
     43         out_shape=like.shape,

~\anaconda3\envs\xxx\lib\site-packages\rasterio\env.py in wrapper(*args, **kwds)
    382     def wrapper(*args, **kwds):
    383         if local._env:
--> 384             return f(*args, **kwds)
    385         else:
    386             with Env.from_defaults():

~\anaconda3\envs\xxx\lib\site-packages\rasterio\features.py in rasterize(shapes, out_shape, fill, out, transform, all_touched, merge_alg, default_value, dtype)
    344         if len(out_shape) != 2:
--> 345             raise ValueError('Invalid out_shape, must be 2D')
    347         out = np.empty(out_shape, dtype=dtype)

ValueError: Invalid out_shape, must be 2D

What worked

The types for the inputs to zonal_stats are what I think is expected according to the manual, i.e. first a GeoDataFrame (municipal_filtered, with a shape of (20, 17)), then an nd array (pop_array, with a shape of (1, 1179, 5935)), and the affine. Is the issue the third dimension for the array here, meaning I should somehow remove it?

Changing zonal_stats to gen_zonal_stats works, though I'm not sure what to do with the resulting generator object.

Displaying both the pop_array (via imshow()) and the polygons (via add_feature) on the same figure in matplotlib also works fine.

My system

I'm running Python 3.8.8. Installed rasterstats today via Anaconda so I presume all the latest dependencies are there.

1 Answer 1


So I figured it out with some kind help :)

I needed to specifically define the band when reading the array from the source raster (even when there's only 1 band)

array = dataset.read(1)

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.