Zonal statistics with rasterstats: "Invalid out_shape, must be 2D" (Python)

Question

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.

Code

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

    # CALCULATING ZONAL STATS

    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 
     48     # CALCULATING ZONAL STATS
     49 
---> 50     municipal_pop_sum = zonal_stats(municipal_filter, pop_array, affine=affine, stats=['sum'], geojson_out=True)
     51 
     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))
     32 
     33 

~\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)
    157 
    158             # rasterized geometry
--> 159             rv_array = rasterize_geom(geom, like=fsrc, all_touched=all_touched)
    160 
    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)
    343 
    344         if len(out_shape) != 2:
--> 345             raise ValueError('Invalid out_shape, must be 2D')
    346 
    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.

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)