Zonal Statistics & Raster Sampling in Python

Zonal statistics is the operation that turns a continuous surface into a table: mean elevation per catchment, dominant land-cover class per census tract, 90th-percentile temperature per neighbourhood. It is where the vector and raster halves of Spatial Analysis & Advanced Query Techniques meet, and it is deceptively easy to get subtly wrong — every implementation has to answer the same awkward question about partially covered pixels, and they answer it differently. This topic covers the pixel-in-polygon decision, nodata semantics, categorical versus continuous summaries, point sampling, and the memory strategy that keeps a per-zone loop from swallowing a multi-gigabyte scene. It consumes the datasets prepared in Raster Data Handling with Rasterio and xarray & rioxarray Raster Cubes, and sits beside Geometric Intersections & Overlays, which answers the same "how much of A is inside B" question for two vector layers.

Architecture & Data Structures

Strip away the libraries and a zonal statistic is a weighted group-by. The grouping key is a zone identifier taken from a GeoDataFrame; the values are pixels drawn from a 2D array; and the weight is how much of each pixel belongs to that zone. Four pieces of state control the whole computation:

• The affine transform — a six-parameter mapping from (row, col) to map coordinates. It, not the pixel array, defines where a pixel is.
• The grid CRS — the coordinate system those map coordinates live in. Pixels have no opinion about the zones' CRS, and nothing will warn you if they disagree.
• The nodata value — the sentinel that must be excluded before any arithmetic. If it leaks into a sum, the mean is silently wrong rather than obviously wrong.
• The coverage rule — how a pixel straddling the zone boundary is counted. This is the only genuinely ambiguous part, and it is the one that changes your answer.

import geopandas as gpd
import rasterio

zones = gpd.read_file("parcels.gpkg", layer="parcels")

with rasterio.open("ndvi_utm32n.tif") as src:
    grid_crs = src.crs               # EPSG:32632 — what the pixel coordinates mean
    transform = src.transform        # affine: (row, col) -> (easting, northing)
    fill = src.nodata                # the sentinel that must never enter a mean
    pixel_area = abs(transform.a * transform.e)   # m² per pixel in a metric CRS

# A zonal statistic is a group-by: key = zone id, values = pixels, weight = coverage
assert zones.crs is not None and zones.crs.equals(grid_crs), "align the CRSs first"
print(f"{pixel_area:.1f} m² per pixel")

The coverage rule has three practical settings. The centroid rule — GDAL's default, and what rasterstats does with all_touched=False — includes a pixel if and only if its centre falls inside the polygon. All-touched includes any pixel the polygon boundary so much as clips, which inflates the sample and drags the statistic toward whatever surrounds the zone. Exact coverage gives each pixel a weight between 0 and 1 equal to the fraction of its area inside the polygon, which is what exactextract computes analytically. On a zone that spans hundreds of pixels the three agree to a rounding error; on a zone a few pixels across they can disagree badly.

Notice the two extremes in the exact-coverage panel: cells at 0.01 and 0.00 are counted in full under all_touched and dropped entirely under the centroid rule, despite contributing almost nothing. That single observation explains most of the disagreement you will ever see between two zonal-statistics implementations.

Environment Configuration & Dependency Resolution

Four libraries cover the whole problem space, and they layer on the same GDAL/GEOS/PROJ foundation described in Raster Data Handling with Rasterio. Pin them together so a GDAL upgrade in one wheel cannot desynchronise the others.

python -m pip install \
  "rasterio>=1.3,<2" "geopandas>=1.0" "shapely>=2.0" \
  "rasterstats>=0.19" "exactextract>=0.2" \
  "rioxarray>=0.15" "xarray>=2024.3" "regionmask>=0.12" "numpy>=1.26"

• rasterstats is pure Python over rasterio and numpy, so it installs anywhere rasterio does. It offers zonal_stats, point_query, categorical=True, and an add_stats hook for arbitrary callables.
• exactextract is a compiled C++ core with Python bindings. Binary wheels exist for mainstream CPython versions on Linux, macOS and Windows; if pip starts invoking CMake, you have fallen off the wheel matrix — install from conda-forge instead of building.
• rioxarray plus regionmask is the route for multi-band or time-series cubes, covered in depth under xarray & rioxarray Raster Cubes.
• rasterio.mask with plain NumPy is always available and is what you drop to when you need a statistic no library exposes.

On Windows, keep the whole stack from one source. A pip install rasterio alongside a conda install gdal puts two GDAL DLLs on the search path and produces an ImportError: DLL load failed at the first import rasterio — never a message about GDAL. If you already use conda, take rasterio, geopandas, rasterstats and exactextract from conda-forge in a single solve.

Vectorized Operations & Core Workflow

The canonical pipeline is four steps: load the zones, reproject them into the raster's CRS, run one summarising call over the whole layer, and join the resulting table back onto the GeoDataFrame by index. Everything else is a variation on it.

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

zones = gpd.read_file("census_tracts.gpkg")

with rasterio.open("ndvi_utm32n.tif") as src:
    if zones.crs != src.crs:
        zones = zones.to_crs(src.crs)   # move the vector; never resample the raster
    fill = src.nodata

# Passing the PATH (not an array) lets rasterstats read one window per feature
stats = zonal_stats(
    zones,                       # consumed via __geo_interface__ — the CRS is NOT carried
    "ndvi_utm32n.tif",
    stats=["count", "mean", "std", "median", "percentile_90"],
    nodata=fill,                 # explicit, in case the file carries no nodata tag
    all_touched=False,           # centroid rule
)

zones = zones.join(pd.DataFrame(stats, index=zones.index))
print(zones[["tract_id", "count", "mean", "percentile_90"]].head())

Two details in that block matter more than they look. zonal_stats consumes the zones through __geo_interface__, the GeoJSON-like mapping that a GeoDataFrame exposes — and __geo_interface__ carries geometry but no CRS. That is the mechanical reason a mismatch is silent rather than fatal, and the reason the to_crs line comes first. Second, results come back as a plain list of dicts in input order, so pd.DataFrame(stats, index=zones.index) is the correct join: never sort or filter the zones between the call and the join.

The exactextract equivalent is a single call that returns a DataFrame directly and carries the identifier column through for you:

import geopandas as gpd
from exactextract import exact_extract

zones = gpd.read_file("census_tracts.gpkg").to_crs("EPSG:32632")

summary = exact_extract(
    "ndvi_utm32n.tif",                      # path, rasterio dataset, or xarray DataArray
    zones,                                  # path or GeoDataFrame
    ["count", "mean", "stdev", "min", "max"],
    include_cols=["tract_id"],
    output="pandas",
)

# count is the SUM OF COVERAGE FRACTIONS — 9.14, not 9. It is a float on purpose.
print(summary.head())

A head-to-head of the two — API shape, weighting behaviour, and how each scales as zone count rises — is in Zonal Statistics with rasterstats vs exactextract. The short version: rasterstats for breadth and hackability, exactextract when partial pixels carry real weight or when there are a great many zones.

Zonal Aggregation & Sampling Details

Continuous versus categorical

A continuous raster (elevation, NDVI, temperature) supports the whole arithmetic family: mean, standard deviation, median, percentiles. A categorical raster (land cover, soil class, zoning code) supports none of them — the mean of class codes 10 and 50 is 30, which is a different class or no class at all. For categories you want tallies and fractions.

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

CLASS_MAP = {10: "tree_cover", 20: "shrubland", 40: "cropland", 50: "built_up", 80: "water"}

zones = gpd.read_file("census_tracts.gpkg").to_crs("EPSG:32632")

cover = zonal_stats(
    zones, "worldcover_utm32n.tif",
    categorical=True,
    category_map=CLASS_MAP,
    nodata=0,
    all_touched=False,
)
# -> [{'tree_cover': 812, 'cropland': 4410, 'built_up': 96}, ...]  raw pixel tallies

# Turn tallies into shares; absent classes become 0, not NaN
shares = pd.DataFrame(cover, index=zones.index).fillna(0)
shares = shares.div(shares.sum(axis=1), axis=0)

exactextract expresses the same idea with area weights rather than pixel counts, which is the more defensible number when zones are small relative to the cell size:

import geopandas as gpd
from exactextract import exact_extract

zones = gpd.read_file("census_tracts.gpkg").to_crs("EPSG:32632")

frac = exact_extract(
    "worldcover_utm32n.tif", zones,
    ["frac", "majority", "variety"],
    include_cols=["tract_id"],
    output="pandas",
)
# `frac` expands into one frac_<value> column per class present, summing to 1 per zone

majority gives the modal class, variety the number of distinct classes — a cheap heterogeneity index. Whatever you do, never resample a categorical raster with bilinear or cubic interpolation on the way in; nearest neighbour is the only correct choice, as covered under Resampling & Overviews When Writing COGs.

Nodata, and what count really means

Nodata handling is where quiet errors live. If the GeoTIFF declares a nodata value, rasterio and both statistics libraries honour it. If it does not — common with older DEMs that use -9999 as a convention — every one of those sentinels is treated as a real measurement and your catchment mean comes back around -3000. Pass nodata= explicitly whenever you did not write the file yourself. NaN is a special case: it never compares equal to itself, so a nodata=nan tag only works through masked reads, not through equality tests you write by hand.

Once nodata is excluded, count acquires three different meanings that are easy to conflate. In rasterstats it is the number of valid pixels selected by the coverage rule. In exactextract it is the sum of coverage fractions over valid pixels — a float. And neither is the same as the number of cells in the zone's bounding box. When a zone returns count of 0 and every other statistic as None, it means no pixel satisfied the rule: the zone is smaller than a cell, sits entirely on nodata, or — most often — is in the wrong CRS.

Reading only what a zone needs

Reading a whole scene into memory to summarise a handful of parcels is the single most common way a zonal job dies. A 40 000 × 40 000 float32 band is 6.4 GB before you have computed anything. The fix is a windowed read per zone: convert the zone's bounds to a pixel window, read only that window, build the mask at that window's transform, and discard both before moving on. Peak memory then tracks the largest zone, not the raster. rasterstats already does this internally when you hand it a path — one of the better reasons to pass a path rather than a pre-read array.

When a library's built-in statistics are not enough — a trimmed mean, a bespoke index, a histogram you want to keep — write the loop yourself. rasterio.features.geometry_mask burns the geometry onto exactly the window you read, so the mask and the data always share a transform:

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import geometry_mask
from rasterio.windows import Window, from_bounds


def zone_pixels(src, geom):
    """Return the valid pixel values of one zone, reading only its window."""
    win = from_bounds(*geom.bounds, transform=src.transform)
    win = win.round_offsets().round_lengths()
    win = win.intersection(Window(0, 0, src.width, src.height))
    if win.width < 1 or win.height < 1:
        return np.array([], dtype="float32")

    data = src.read(1, window=win, masked=True)          # masked array honours nodata
    inside = geometry_mask(
        [geom],
        out_shape=data.shape,
        transform=src.window_transform(win),
        invert=True,          # True where the pixel centre falls INSIDE the geometry
        all_touched=False,
    )
    return data[inside & ~data.mask].astype("float32")


catchments = gpd.read_file("catchments.gpkg")

with rasterio.open("dem_utm32n.tif") as src:
    catchments = catchments.to_crs(src.crs)
    for row in catchments.itertuples():
        values = zone_pixels(src, row.geometry)
        if values.size == 0:
            continue                                      # zone off-grid or all nodata
        p5, p50, p95 = np.percentile(values, [5, 50, 95])
        print(f"{row.catchment_id}: relief {p95 - p5:.1f} m, median {p50:.1f} m")

geometry_mask returns True outside the geometry by default, which suits NumPy's masked-array convention but is the reverse of what most people expect — invert=True flips it to "True inside". Getting that backwards produces statistics for the complement of every zone, and the numbers look plausible enough to ship. Full masking recipes, including crop=True and multipart handling, are in Masking Rasters with Polygons in Python.

Sampling at point locations

Point sampling is the degenerate case where the zone is a single coordinate. rasterio's sample() is the fastest route: it takes an iterable of (x, y) pairs in the raster's CRS and yields one array per point.

import geopandas as gpd
import numpy as np
import rasterio

sensors = gpd.read_file("air_quality_sensors.geojson")

with rasterio.open("lst_utm32n.tif") as src:
    sensors = sensors.to_crs(src.crs)                  # sample() takes raster-CRS coords
    coords = [(geom.x, geom.y) for geom in sensors.geometry]
    sensors["lst_c"] = [record[0] for record in src.sample(coords, indexes=1)]
    fill = src.nodata

# Points outside the grid come back as the nodata value — no exception is raised
sensors["lst_c"] = sensors["lst_c"].replace(fill, np.nan)
print(sensors["lst_c"].isna().sum(), "sensors fell outside the scene")

sample() is nearest-neighbour: it returns the value of the cell containing the point, with no interpolation. For a smooth surface such as elevation or temperature that produces a visible staircase along a transect; rasterstats.point_query(..., interpolate="bilinear") blends the four surrounding cells instead. For a class raster, bilinear is meaningless and interpolate="nearest" is mandatory. The trade-offs, plus how to sample a time series of bands at once, are worked through in Sampling Raster Values at Point Locations.

Cubes: one mask, every time step

When the raster is a stack — twelve monthly composites, a decade of daily temperature — you do not want twelve independent zonal passes. Build the zone mask once and reduce along the spatial dimensions, letting xarray broadcast over time and dask chunk the work:

import geopandas as gpd
import numpy as np
import regionmask
import rioxarray  # noqa: F401 — registers the .rio accessor
import xarray as xr

cube = xr.open_dataset("modis_lst_2024.nc", chunks={"time": 12})["lst"]
zones = gpd.read_file("admin2.gpkg").to_crs("EPSG:4326")

regions = regionmask.from_geopandas(zones, numbers="zone_id", name="admin2")
mask3d = regions.mask_3D(cube["lon"], cube["lat"])      # bool, dims (region, lat, lon)

# On a lon/lat grid a cell's ground area shrinks toward the poles — weight by cos(lat)
weights = np.cos(np.deg2rad(cube["lat"]))
zonal_mean = cube.weighted(mask3d * weights).mean(dim=("lat", "lon"))

zonal_mean.compute().to_dataframe(name="lst_mean").to_parquet("admin2_lst_2024.parquet")

Two regionmask details bite people. It wraps longitudes into −180…180 by default, so on a projected grid with x/y coordinates you must pass wrap_lon=False or the mask comes back empty. And mask_3D applies the centroid rule; regionmask>=0.12 adds mask_3D_frac_approx for approximate fractional coverage on regular lon/lat grids when small regions matter.

CRS Alignment & Projection Pipeline

Every zonal statistic assumes zones and pixels are indexed in the same coordinate system, and nothing in the stack enforces it. rasterstats receives bare GeoJSON geometries with the CRS already stripped; rasterio.mask compares raw numbers; a NumPy mask has no concept of projection at all. When the two disagree the coordinates simply do not overlap, and you get an empty result rather than an error.

Always reproject the vector, not the raster. Reprojecting zones is a lossless coordinate transform on a few thousand vertices. Reprojecting the raster resamples every pixel, changes values through interpolation, and costs orders of magnitude more. Establish the canonical CRS with Coordinate Systems with PyProj and let the raster define it for this operation.

import geopandas as gpd
import rasterio

zones = gpd.read_file("parcels.gpkg")

with rasterio.open("ndvi_utm32n.tif") as src:
    if zones.crs is None:
        raise ValueError("zones carry no CRS — declare it with set_crs() before anything")
    if not zones.crs.equals(src.crs):
        zones = zones.to_crs(src.crs)

    # Fail loudly when the extents do not actually meet, instead of returning None
    left, bottom, right, top = src.bounds
    zx0, zy0, zx1, zy1 = zones.total_bounds
    if zx1 < left or zx0 > right or zy1 < bottom or zy0 > top:
        raise ValueError(f"zones {zones.total_bounds} lie outside raster {src.bounds}")

The choice of which projected CRS still matters once they agree. Area-weighted statistics assume every pixel covers the same ground area, which holds in a local UTM zone and fails in a geographic CRS, where a cell's north–south extent is constant in degrees but its east–west extent shrinks with the cosine of latitude. Over a country-scale extent use an equal-area projection (Lambert Azimuthal Equal Area, Albers) so pixel area stays constant; over a lon/lat cube apply the cos(lat) weighting shown above. Never compute area-weighted statistics in EPSG:3857 — Web Mercator's scale factor grows away from the equator, so pixels near the poles claim many times the ground area they actually cover. Web Mercator belongs at the tile-render boundary and nowhere else.

One more alignment trap: if you are combining several rasters per zone — NDVI over a cloud mask, say — they must share a grid, not merely a CRS. Two scenes in EPSG:32632 with different origins or resolutions will not align cell-for-cell. Snap one to the other with rioxarray's reproject_match, or read both through a WarpedVRT onto a common grid, before any joint statistic. The worked example for the NDVI case is Computing NDVI Zonal Means per Parcel.

Production Export & Integration

Zonal statistics almost always feed something else: a dashboard, a model, a choropleth. The output is a tidy table keyed by zone id, so integration is mostly about keeping that key trustworthy.

import geopandas as gpd
from exactextract import exact_extract

parcels = gpd.read_file("parcels.gpkg").to_crs("EPSG:32632")

summary = exact_extract(
    "ndvi_utm32n.tif", parcels, ["mean", "count"],
    include_cols=["parcel_id"], output="pandas",
).rename(columns={"mean": "ndvi_mean", "count": "cell_equivalents"})

parcels = parcels.merge(summary, on="parcel_id", validate="one_to_one")

parcels.to_parquet("parcel_ndvi.parquet")                       # GeoParquet, CRS in metadata
parcels.drop(columns="geometry").to_csv("parcel_ndvi.csv", index=False)

The validate="one_to_one" argument is worth the keystrokes: it turns a duplicated parcel_id — the usual consequence of an earlier explode or overlay — into an immediate MergeError rather than a row-count explosion nobody notices until the choropleth looks odd. Beyond the file drop:

• GeoParquet keeps the CRS in column metadata and round-trips cleanly, which is why it beats CSV plus a sidecar. See Cloud-Native Geospatial Formats.
• PostGIS takes the joined frame through GeoDataFrame.to_postgis(); store the statistic columns beside the geometry so the map server does one query, not two. The connection pattern is in PostGIS Integration with Python.
• Web maps consume the result directly — a ndvi_mean column joined onto the zone geometry is exactly what a choropleth needs, per Web Mapping & Interactive Visualization.
• Remote rasters work unchanged if you point the reader at a COG URL, but each zone becomes an HTTP range request. Sort zones by a spatial key first so consecutive reads hit nearby blocks, and read the guidance in Windowed Reads from Cloud-Optimized GeoTIFF.
• Scale out when one process is not enough: the loop is embarrassingly parallel over zone chunks, so joblib across cores or the partitioning approach in Scaling with Dask-GeoPandas both apply. Give each worker its own open dataset handle — rasterio datasets are not fork-safe.

Performance checklist. Pass a path, not a pre-read array, so reads stay windowed. Keep zones and raster in the same CRS to skip a reprojection per call. Build overviews on the raster if a coarse statistic is acceptable and read the reduced resolution instead. Dissolve zones you will only ever report together, using the patterns in Dissolving & Aggregating Features by Attribute, so you compute one statistic instead of two hundred. And prefer one call over the whole layer to a Python loop of single-feature calls — both libraries amortise setup across features.

Windows / Platform Edge Cases & Debugging

• ImportError: DLL load failed while importing _base on Windows — two GDAL builds on the path from mixing pip rasterio with conda gdal. Rebuild the environment from a single channel.
• pip install exactextract starts running CMake — no wheel matches your Python version or platform. Install from conda-forge, or move to a Python version that has wheels.
• Every zone returns count: 0 and mean: None — a CRS mismatch, an off-grid extent, or zones smaller than a cell. Run the extent assertion above, then retry with all_touched=True to distinguish the last case.
• ValueError: Input shapes do not overlap raster from rasterio.mask — the geometry is fully outside the scene. Filter zones by src.bounds first, or catch and record a null result.
• Means that look like -3271.4 — the file declares no nodata and its -9999 fill entered the arithmetic. Pass nodata=-9999 explicitly.
• Statistics computed for the complement of every zone — geometry_mask returns True outside by default; you need invert=True.
• A MergeError or silently duplicated rows after the join — multipart zones were exploded upstream and the id is no longer unique. Dissolve back, or key on a surrogate index.
• An empty regionmask mask on a projected grid — longitude wrapping. Pass wrap_lon=False when coordinates are eastings and northings rather than degrees.
• Backslash paths failing in GDAL readers — use forward slashes or raw strings; UNC paths need //server/share/... form.
• A COG zonal run that crawls — thousands of tiny range requests. Set GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR, open the dataset once outside the loop, and sort zones spatially so consecutive windows share blocks.
• Attempt to read from a closed dataset — a statistic was deferred until after the with block. Capture src.crs, src.nodata and src.transform into locals while the dataset is open.

Frequently Asked Questions

Should I use all_touched=True or the default? Default (centroid) when zones are large relative to the cell size — it matches GDAL, QGIS and PostGIS Raster, so your numbers reconcile with everyone else's. all_touched=True only when zones are small or narrow enough that the centroid rule would return nothing, and then accept that boundary values are over-represented. If the answer needs to be defensible, use exact coverage weighting instead of choosing between two approximations.

Why does exactextract report a fractional count? Because count is the sum of coverage fractions, not a pixel tally. A zone covering nine whole cells plus slivers of eight more reports something like 9.14. That is the correct denominator for an area-weighted mean, and it is also a useful diagnostic: a count far below 1 tells you the zone is smaller than a cell and the mean rests on a single pixel.

Can I reproject the raster instead of the zones? You can, but you should not. Reprojecting the vector moves a few thousand vertices exactly; reprojecting the raster resamples every pixel and changes the values you are about to average. Reserve raster reprojection for the case where several rasters must be forced onto one grid, and even then use reproject_match with nearest for categorical data.

How do I get zonal statistics without loading the whole raster? Hand the reader a file path rather than a NumPy array. rasterstats and exactextract both read one window per feature, so peak memory follows the largest zone rather than the scene. Only pass an array plus an affine when the data is already in memory for another reason.

What is the right statistic for a land-cover raster? Never the mean — class codes are labels, not quantities. Use categorical=True in rasterstats for per-class pixel tallies, or frac, majority and variety in exactextract for area shares, modal class and class diversity. If you also resample such a raster anywhere in the pipeline, nearest neighbour is the only valid method.

Do zonal means need an equal-area projection? For a single UTM zone, no — pixel area is effectively constant, so a plain mean is already area-weighted. For continental extents or lon/lat grids, yes: pixel ground area varies with latitude, so use an equal-area CRS or apply cos(lat) weights. EPSG:3857 is wrong for this in every case.