I'm truly stumped on why I'm getting incorrect results with zonal_statistics. I'm passing a geometry and the raster in the normal fashion that's worked for me before, with 'data' as the raster and 'gdf' a geopandas dataframe:

zonal_stats(gdf.geometry, data.read(1), affine=data.transform, stats="count", categorical=True, nodata=data.nodata)

My array has values from 0 to 584 (np.max(data.read(1)) returns 584). But regardless of what geometry I pass, it always returns the number 64537, which is obviously nowhere present in my array.

Where is it getting this number from?

Here I will provide steps to reproduce: The files (example_shapefiles.geojson and raster_file.tiff) can be downloaded from here. Here is the snippet to load in the files:

import geopandas as gdf
import rasterio
from rasterstats import zonal_stats


The shapefiles initially have EPSG:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

I convert to the same EPSG as the raster with:


Now they both have EPSG:

<Projected CRS: EPSG:5070>
Name: NAD83 / Conus Albers
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Albers Equal Area
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Now running zonal_statistics with

zonal_stats(gdf.geometry, data.read(1), affine=data.transform, stats="count", categorical=True, nodata=data.nodata)

The results is

[{64537: 24603, 'count': 24603},
 {64537: 9846, 'count': 9846},
 {64537: 4540, 'count': 4540},
 {64537: 3747, 'count': 3747},
 {64537: 1919, 'count': 1919},
 {64537: 1497, 'count': 1497},
 {64537: 1304, 'count': 1304},
 {64537: 1255, 'count': 1255},
 {64537: 1112, 'count': 1112},
 {64537: 1062, 'count': 1062},
 {64537: 1053, 'count': 1053}]

I don't know what else I might be doing wrong. I try to figure things out on my own I'm pretty new to GIS and I just can't tell what I'm doing wrong.


1 Answer 1


Your polygons don't intersect your raster:

import geopandas as gpd
import rasterio
from shapely.geometry import box

gdf = gpd.read_file('example_shapefiles.geojson')
data = rasterio.open('raster_file.tiff')

raster_bbox = box(*data.bounds)

for feat in gdf.geometry:

The datasets have to overlap and be provided in the same projection

  • Thanks for your answer, it looks like you're right about this. I was starting with ~50,000 parcel shapefiles from a specific county and my raster was coming from the Microsoft GIS data catalog that has separate tiles in a Cloud Optimized GeoTiff. Previously I would take the centroid of the county and pass that to the data catalog to get the right tile. For reasons I don't understand, this didn't work here - as you observed, the tile didn't intersect with my parcels. Instead I passed the county's bounding box, which returned 2 tiles instead of one, with the other one containing my parcels.
    – Samuel T
    Oct 1, 2023 at 17:59

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.