Zonal Statistics and Vector–Raster Integration

Most remote sensing questions are eventually asked about a place, not about a pixel: the mean NDVI of a field, the burned fraction of a protected area, the water extent inside a catchment. Answering them means joining vector geometries to raster values, and the join has more failure modes than it looks. This topic covers those mechanics — alignment, boundary conventions, validity handling and batch shape — inside the wider context of Core Raster Fundamentals & STAC Mapping.

The specific challenge is that a polygon and a pixel grid disagree about geometry. Pixels are fixed squares in a projected grid; polygons are arbitrary shapes in whatever CRS the surveyor used. Every zonal statistic is therefore an approximation, and the quality of the approximation depends on choices that are usually left at their defaults.

Prerequisites

pip install "rasterio>=1.3.0" "geopandas>=0.13" "rasterstats>=0.19" "shapely>=2.0" "numpy>=1.23"

Conceptually you need windowed reads from Optimizing rasterio Window Reads for Memory Efficiency, CRS handling from Mastering CRS Transformations in rasterio, and the validity rules from Extracting nodata and dtype from a GeoTIFF.

How a polygon becomes a set of pixels

Every zonal statistic passes through the same four steps, and each one is a place where the answer can change.

Step-by-step workflow

1. Put both layers in one CRS — by moving the vectors

import geopandas as gpd
import rasterio

with rasterio.open("ndvi_20230615.tif") as src:
    raster_crs = src.crs

zones = gpd.read_file("fields.gpkg")
if zones.crs != raster_crs:
    zones = zones.to_crs(raster_crs)     # exact: no pixels are touched

Reprojecting the vector layer transforms a few hundred coordinate pairs. Reprojecting the raster resamples millions of pixels and changes the values you are about to average, which is the argument made at length in Reprojecting a Raster from UTM to WGS84.

2. Read a window, not a scene

from rasterio.windows import from_bounds

def zone_window(src, geom, pad: int = 1):
    """Pixel window covering a geometry, rounded outward and clamped to the scene."""
    win = from_bounds(*geom.bounds, transform=src.transform)
    win = win.round_lengths(op="ceil").round_offsets(op="floor")
    # pad so all_touched and edge pixels are inside the read
    win = win.round_offsets(op="floor")
    return win.intersection(rasterio.windows.Window(0, 0, src.width, src.height))

3. Rasterise the geometry into the window

from rasterio.features import geometry_mask

def zone_mask(src, geom, window, *, all_touched: bool = False):
    """Boolean mask, True inside the geometry, for the pixels covered by `window`."""
    transform = src.window_transform(window)
    inside = geometry_mask(
        [geom],
        out_shape=(int(window.height), int(window.width)),
        transform=transform,
        invert=True,               # True = inside the polygon
        all_touched=all_touched,
    )
    return inside

4. Aggregate with validity handled explicitly

import numpy as np

def zone_stats(src, geom, *, band: int = 1, all_touched: bool = False) -> dict:
    window = zone_window(src, geom)
    if window.width < 1 or window.height < 1:
        return {"count": 0}

    data = src.read(band, window=window, masked=True)   # honours nodata and masks
    inside = zone_mask(src, geom, window, all_touched=all_touched)
    values = np.ma.masked_array(data, mask=data.mask | ~inside)

    count = int(values.count())
    if count == 0:
        return {"count": 0}
    return {
        "count": count,
        "mean": float(values.mean()),
        "median": float(np.ma.median(values)),
        "std": float(values.std()),
        "min": float(values.min()),
        "max": float(values.max()),
    }

masked=True is what keeps declared fill values out of the mean; combining that mask with the geometry mask is what keeps pixels outside the polygon out of it. Both are easy to omit and neither omission raises.

5. Emit a table, not a pile of rasters

One row per geometry per date, carrying the identifier, the statistics and the pixel count. The count is not optional: it is the only field that tells a consumer whether a mean rests on four pixels or four thousand, and it is what makes cloud-masked time series interpretable.

Parameter reference

The boundary convention, quantified

The single largest source of disagreement between two zonal statistics implementations is not the arithmetic — it is which pixels each one counted.

For area-like statistics — counting burned pixels, water pixels, a class fraction — the convention matters even at large sizes, because it changes the denominator as well as the numerator. State it explicitly in any published figure.

Fractional coverage, and when a binary mask is not enough

The mask produced by geometry_mask is binary: a pixel is in or out. For most statistics on reasonably sized polygons that is fine, but three situations need the fraction of each pixel that falls inside the geometry rather than a yes-or-no answer.

The first is small polygons, where boundary pixels are a large share of the total. A field five pixels across has sixteen boundary pixels out of twenty-five, so a binary decision on each of them determines most of the answer. Weighting each boundary pixel by its overlap fraction gives a mean that changes smoothly as the polygon moves, instead of jumping as pixel centres cross the edge.

The second is area accounting that must sum correctly. If you compute the burned area of every administrative unit in a region and the units tile the region exactly, binary masks will not sum to the region’s total: pixels on shared boundaries are counted twice or not at all, depending on the convention. Fractional weights make the sum exact by construction, which matters when the numbers are reported publicly and have to reconcile.

The third is statistics that are sensitive to weighting rather than to membership — an area-weighted mean, a population-weighted average, or any per-hectare figure. Here the weight is doing real work and a binary approximation biases the result in a direction that depends on the polygon’s shape.

Computing coverage fractions means rasterising at a finer resolution and averaging, or using a library that computes exact polygon-pixel intersections. The simple version is to rasterise the geometry at, say, ten times the raster resolution and take the block mean, which converges quickly and costs a hundred times the mask memory for the window — acceptable, because the window is small.

import numpy as np
from rasterio.features import geometry_mask


def coverage_fraction(src, geom, window, factor: int = 10) -> np.ndarray:
    """Approximate per-pixel coverage fraction by rasterising at `factor` times resolution."""
    h, w = int(window.height), int(window.width)
    fine_transform = src.window_transform(window) * rasterio.Affine.scale(1 / factor)
    fine = geometry_mask([geom], out_shape=(h * factor, w * factor),
                         transform=fine_transform, invert=True)
    # Block-average the fine mask back down to one weight per output pixel
    return fine.reshape(h, factor, w, factor).mean(axis=(1, 3))

The weighted mean is then np.average(values, weights=weight * valid), where valid excludes nodata. Record which method produced a published figure: a coverage-weighted mean and a centre-based mean are different statistics, and the difference is exactly the quantity plotted earlier.

Scaling the join across dates and scenes

A single date over a few hundred polygons is a loop. A five-year time series over a national field register is a data-engineering problem, and it has three levers.

The first is loop order, and it is the one people get wrong. Open each raster once and iterate the geometries that intersect it, rather than iterating geometries and reopening rasters. Every open against remote storage costs requests and a header parse; multiplying that by the polygon count is the single largest avoidable cost in the whole workflow.

The second is a spatial pre-filter. Most polygon-scene pairs do not intersect at all, and finding that out with a windowed read is the most expensive possible way to learn it. Build an index over the geometries once, test it against each scene’s footprint — which is already in the STAC item, as covered in Querying STAC Catalogs Programmatically — and only then read anything.

The third is parallelism that matches the bottleneck. Zonal statistics against remote rasters are network-bound, so threads work well and processes mostly waste memory. Against local NVMe storage with heavy arithmetic per polygon, the balance flips. Measure one tile before choosing, using the same reasoning as Scaling Raster Processing with Dask applies to array workloads.

Two smaller habits pay for themselves at this scale. Cache the rasterised mask per geometry when the raster grid is stable across dates — the mask depends only on the geometry and the transform, so recomputing it for 180 dates is pure waste. And write results incrementally, appending each scene’s rows to a partitioned table rather than accumulating them in memory, so that a run interrupted at eighty percent leaves eighty percent of the answers on disk instead of nothing.

---

Library or hand-written loop?

rasterstats implements everything above and is the right default for exploratory work and for one-off analyses: one call, sensible defaults, and support for categorical counts as well as continuous statistics. It reads the raster itself, handles the masking, and returns a list of dictionaries ready to become a dataframe.

A hand-written loop earns its keep in three situations. The first is when the read pattern matters — when you need one open per scene with many windows inside it, or when the raster lives behind a signed URL that must be refreshed. The second is when the statistic is not in the library’s vocabulary: a percentile of a masked, weighted sample, a fraction above a threshold that varies per geometry, or anything that needs the raw values rather than a summary. The third is when the same mask is reused across dates, because a library call recomputes it every time and a loop can cache it.

The two are not exclusive. A common shape is to use rasterstats for the exploratory pass that decides what to compute, then reimplement the chosen statistic in a loop for the production run, asserting on one tile that both agree to within floating-point tolerance. That assertion is worth writing down: it is the cheapest available check that a hand-written mask is not inverted, and it catches the boundary-convention mismatch immediately if the two are configured differently.

Whichever you choose, keep the conventions in one place. A small dataclass holding the boundary rule, the band, the nodata override and the statistic list — passed into every call — makes the configuration visible and reviewable, and stops one analyst’s all_touched=True quietly entering a dataset built by everyone else with the default. The same argument applies to the index definitions in Building a YAML-Driven Multi-Index Pipeline: configuration that lives in a shared object is configuration that can be audited.

---

Verification and testing

Three checks catch the errors that matter:

import numpy as np

# 1. A synthetic raster of known constant value must return that value exactly
#    for any polygon fully inside it, under either convention.
assert abs(zone_stats(src, big_interior_polygon)["mean"] - 42.0) < 1e-6

# 2. The pixel count must scale with area: doubling a polygon's linear size
#    should roughly quadruple the count.
assert 3.5 < stats_large["count"] / stats_small["count"] < 4.5

# 3. A polygon outside the raster must return count 0 rather than raising.
assert zone_stats(src, far_away_polygon)["count"] == 0

The first check is the one that catches a mask that is inverted — the single most common bug in hand-written implementations, because an inverted mask still produces plausible numbers when the surrounding area resembles the interior.

Troubleshooting

Every zone returns count: 0

The vector layer and the raster are in different CRSs, so the geometries land outside the scene. Compare zones.total_bounds against src.bounds before anything else; the magnitudes alone usually reveal it, as described in Fixing EPSG Mismatches in rasterio.open.

The mean is suspiciously close to the nodata value

masked=True was omitted, so fill pixels are being averaged in. Read masked, or apply src.dataset_mask() explicitly.

Statistics disagree with a desktop GIS by a few percent

Boundary convention. Set all_touched to match the other tool before concluding that either is wrong.

A run over thousands of polygons is impossibly slow

The loop is opening the raster once per polygon, or reading the whole scene per polygon. Open once, iterate geometries inside that context, and read per-geometry windows — the ordering argument developed in Automated Image Clipping and Cropping.

Memory grows steadily through a long run

Masked arrays are being retained — often accidentally, by appending the array rather than the summary to a list. Keep only the statistics.

Frequently Asked Questions

Q: Why do my zonal means differ from a desktop GIS? Almost always the boundary convention. Desktop tools often include every pixel the polygon touches, while the default in rasterio and rasterstats includes only pixels whose centre falls inside. On small polygons the two differ by several percent.

Q: Should I reproject the raster or the vector layer? The vector layer, in almost every case. Reprojecting geometries is exact and cheap; reprojecting a raster resamples every pixel, costs far more, and introduces interpolation error into the values you are about to summarise.

Q: How do I handle polygons smaller than a pixel? Decide explicitly rather than letting the default decide. With centre-based inclusion such a polygon can return no pixels at all; with all_touched it returns one to four pixels whose values represent a much larger area. Record the pixel count so consumers can filter on it.

Q: Can I compute zonal statistics without loading the raster? You still need the pixels inside each geometry, but not the rest of the scene. A windowed read per geometry keeps the transfer proportional to the polygons rather than to the imagery, which is what makes the operation viable against remote archives.

---

Related

• Computing Zonal Statistics with rasterstats — the library implementation of this workflow, and where its defaults differ.
• Rasterizing Vector Polygons onto a Raster Grid — turning geometries into aligned label rasters.
• Extracting Pixel Values at Point Locations — the point-sampling case, which has its own alignment traps.
• Automated Image Clipping and Cropping — the same geometry-to-window machinery, used to produce rasters instead of tables.
• Mastering CRS Transformations in rasterio — the alignment layer every zonal join depends on.