Zonal Statistics for Polygon Aggregation

Tabular machine learning models — gradient boosting, random forests, linear regressors — consume one row per sample with a fixed set of numeric columns. But a huge share of geospatial signal lives in rasters: elevation, spectral indices, temperature, land cover. To feed that signal into a tabular model you have to collapse the many pixels that fall inside each analysis unit — an agricultural field, a census tract, a watershed — into a small vector of summary numbers. That reduction is zonal statistics: for every polygon, compute the mean, standard deviation, percentiles, or dominant class of the raster values it covers, and emit one feature row per polygon.

This is part of Spatial Feature Engineering for Machine Learning, which covers the full path from raw geodata to model-ready feature tables. Where raster band math and index calculation produces continuous surfaces such as NDVI, zonal statistics is the step that turns those surfaces into per-polygon predictors — for example, mean NDVI, NDVI standard deviation, and the 90th-percentile NDVI for each field in a crop-yield dataset. This guide covers the complete workflow: matching CRS, running rasterstats.zonal_stats with the right statistics and nodata handling, adding categorical majority for land-cover layers, assembling the feature table, verifying correctness, and tuning performance at scale.

Problem Framing

The core operation is a spatial reduction. A raster is a dense regular grid; a polygon layer is a sparse set of irregular shapes. Zonal statistics overlays the two, selects the pixels that fall within each polygon, and applies an aggregation function to their values. The output has exactly as many rows as there are polygons, which is precisely the shape a tabular estimator expects.

The choice of statistics is a feature-engineering decision, not a formality:

• Mean captures the central tendency — average greenness of a field, average slope of a parcel.
• Standard deviation captures within-polygon heterogeneity. A field with high NDVI variance may be partially stressed; a uniform field will have low variance. This is often more predictive than the mean alone.
• Percentiles (10th, 50th, 90th) describe the distribution’s shape and are robust to a handful of outlier pixels in a way the mean is not.
• Min / max / range flag extremes — the coldest pixel in a frost-risk analysis, the steepest pixel in a landslide model.
• Majority (mode) is the only sensible summary for categorical rasters such as land cover, where averaging integer class codes is meaningless.
• Count records how many valid pixels contributed, which doubles as a data-quality and confidence signal.

A single continuous raster can yield five or six columns per polygon, and stacking several rasters (elevation, NDVI, temperature, precipitation) quickly builds a rich tabular feature matrix. Zonal statistics is therefore one of the highest-leverage steps in the whole pipeline: it is where raster signal becomes model input. Its close relatives are spatial lag and neighborhood statistics, which summarise values from neighbouring polygons rather than underlying pixels, and vector proximity and buffer generation, which measures distances rather than aggregating grids.

Prerequisites & Environment Setup

Zonal statistics sits on top of the standard raster/vector stack. Pin these versions so nodata and all_touched behaviour is reproducible across machines.

# Pinned requirements for zonal statistics
rasterstats==0.19.0
rasterio==1.3.10
geopandas==0.14.4
numpy==1.26.4
shapely==2.0.5
pyproj==3.6.1

Install with:

pip install "rasterstats==0.19.0" "rasterio==1.3.10" "geopandas==0.14.4" \
            "numpy==1.26.4" "shapely==2.0.5" "pyproj==3.6.1"

GDAL/PROJ system dependencies (Ubuntu/Debian):

sudo apt-get install -y gdal-bin libgdal-dev libproj-dev

CRS match is a hard prerequisite

rasterstats does not reproject on your behalf. It overlays polygon coordinates directly onto the raster’s coordinate space. If the polygon layer and the raster are in different coordinate reference systems, the geometries land in the wrong place — usually completely off the raster — and every polygon silently returns None. The single most important setup step is to read the raster’s CRS and reproject the polygons to match it before doing anything else. Full reprojection patterns live in CRS alignment and projection handling.

import rasterio
import geopandas as gpd

RASTER_PATH = "ndvi.tif"

polygons = gpd.read_file("fields.geojson")
with rasterio.open(RASTER_PATH) as src:
    raster_crs = src.crs
    raster_nodata = src.nodata

assert polygons.crs is not None, "Polygon layer has no CRS — cannot align to raster"
if polygons.crs != raster_crs:
    polygons = polygons.to_crs(raster_crs)
assert polygons.crs == raster_crs, "CRS alignment failed"

Step-by-Step Implementation

Step 1 — Align polygon CRS to the raster

Always drive the alignment from the raster, because reprojecting a raster is lossy (it resamples pixels) while reprojecting vectors is exact. Read the raster CRS and nodata value once, then push the polygons onto it.

import rasterio
import geopandas as gpd

def load_aligned(vector_path: str, raster_path: str) -> tuple[gpd.GeoDataFrame, float | None]:
    """Load polygons reprojected to the raster CRS and return the raster nodata value.

    Args:
        vector_path: Path to the polygon vector file.
        raster_path: Path to the raster to sample.
    Returns:
        (polygons in raster CRS, raster nodata value or None).
    """
    polygons = gpd.read_file(vector_path)
    with rasterio.open(raster_path) as src:
        raster_crs = src.crs
        nodata = src.nodata

    if polygons.crs is None:
        raise ValueError("Polygon layer has no CRS; set it before running zonal stats")
    if polygons.crs != raster_crs:
        polygons = polygons.to_crs(raster_crs)
    return polygons.reset_index(drop=True), nodata

polygons, nodata = load_aligned("fields.geojson", "ndvi.tif")
print(f"Aligned {len(polygons)} polygons to raster CRS; nodata = {nodata}")

Step 2 — Run rasterstats.zonal_stats for continuous statistics

zonal_stats accepts the polygons and the raster path and returns a list of dictionaries — one per polygon, in input order. Request the continuous statistics you need and pass the nodata value explicitly so masked pixels never contaminate the aggregates.

from rasterstats import zonal_stats

stats = zonal_stats(
    polygons,                 # GeoDataFrame or list of geometries
    "ndvi.tif",
    stats=["count", "mean", "std", "min", "max"],
    nodata=nodata,            # exclude these pixel values from every statistic
    all_touched=False,        # only pixels whose centre is inside the polygon
    geojson_out=False,        # return plain dicts, not GeoJSON features
)

print(stats[0])
# {'count': 812, 'mean': 0.612, 'std': 0.041, 'min': 0.44, 'max': 0.71}

all_touched=False is the correct default for area-representative continuous features: a pixel counts toward a polygon only when its centre lies inside the boundary. This keeps edge pixels from being double-counted across adjacent polygons. Reserve all_touched=True for the small-polygon case discussed in Step 3 and in Troubleshooting.

Step 3 — Add percentiles and categorical majority

Percentiles capture the shape of the distribution and are far more robust than mean/std when a polygon contains a few anomalous pixels (a cloud shadow, a water body). Request them with the percentile_<n> syntax.

from rasterstats import zonal_stats

continuous_stats = zonal_stats(
    polygons,
    "ndvi.tif",
    stats=["mean", "std", "percentile_10", "percentile_50", "percentile_90"],
    nodata=nodata,
    all_touched=False,
)
print(continuous_stats[0])
# {'mean': 0.612, 'std': 0.041, 'percentile_10': 0.55,
#  'percentile_50': 0.62, 'percentile_90': 0.66}

For a categorical raster such as a land-cover classification, a numeric mean over integer class codes is nonsense. Use categorical=True to get a per-class pixel count, and request majority to get the single dominant class — a clean categorical feature for the model.

lc_stats = zonal_stats(
    polygons,
    "landcover.tif",
    stats=["majority", "count"],
    categorical=True,          # also return {class_code: pixel_count} maps
    nodata=0,                  # 0 is the unclassified/nodata class here
    all_touched=True,          # small parcels still capture their dominant class
)
print(lc_stats[0])
# {'majority': 3, 'count': 812, 3: 640, 2: 150, 5: 22}

The dominant-class integer can be fed directly to tree models, or one-hot / target-encoded first — see encoding categorical geographic features for the encoding choices that matter for spatial categoricals.

Step 4 — Assemble the feature table and join back

zonal_stats returns a list aligned to the input row order, so converting it to a DataFrame and joining on the shared index is safe as long as you never reorder or filter the polygons between the two operations. Reset the index first to make the contract explicit.

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

def zonal_features(
    polygons: gpd.GeoDataFrame,
    raster_path: str,
    stats: list[str],
    nodata: float | None,
    prefix: str,
    all_touched: bool = False,
) -> gpd.GeoDataFrame:
    """Compute zonal statistics and join them back to the polygons as prefixed columns.

    Args:
        polygons: Polygon layer already in the raster CRS.
        raster_path: Raster to sample.
        stats: rasterstats statistic names.
        nodata: Raster nodata value to exclude.
        prefix: Column-name prefix, e.g. 'ndvi'.
        all_touched: Whether edge-touching pixels are included.
    Returns:
        Copy of polygons with one prefixed column per requested statistic.
    """
    polygons = polygons.reset_index(drop=True)
    records = zonal_stats(
        polygons, raster_path, stats=stats,
        nodata=nodata, all_touched=all_touched,
    )
    feats = pd.DataFrame(records, index=polygons.index).add_prefix(f"{prefix}_")
    assert len(feats) == len(polygons), "Row-count mismatch between stats and polygons"
    return polygons.join(feats)

ndvi_feats = zonal_features(
    polygons, "ndvi.tif",
    stats=["count", "mean", "std", "percentile_10", "percentile_90"],
    nodata=nodata, prefix="ndvi",
)
print(ndvi_feats[["ndvi_mean", "ndvi_std", "ndvi_percentile_90"]].head())

Stacking several rasters is just repeated calls with different prefixes:

features = polygons.reset_index(drop=True)
for path, prefix, nd in [
    ("ndvi.tif", "ndvi", nodata),
    ("elevation.tif", "elev", -9999),
    ("temperature.tif", "temp", None),
]:
    features = zonal_features(
        features, path,
        stats=["mean", "std", "percentile_50"],
        nodata=nd, prefix=prefix,
    )

# Drop geometry to hand a pure feature matrix to the model
X = features.drop(columns="geometry")
print(f"Feature matrix: {X.shape[0]} polygons x {X.shape[1]} columns")

For a deeper walkthrough of the rasterstats API — including affine transforms, in-memory arrays, and streaming large vector files — see computing zonal statistics with rasterstats.

Verification & Testing

Never trust an aggregate you have not checked against an independent computation. The most reliable test is to reproduce one polygon’s mean by masking the raster manually with rasterio.mask and comparing.

import numpy as np
import rasterio
from rasterio.mask import mask
from rasterstats import zonal_stats

def verify_polygon_mean(polygons, raster_path, nodata, idx=0, tol=1e-6):
    """Assert zonal_stats mean matches a manual rasterio.mask mean for one polygon."""
    geom = [polygons.iloc[idx].geometry.__geo_interface__]
    with rasterio.open(raster_path) as src:
        out_image, _ = mask(src, geom, crop=True, all_touched=False)

    arr = out_image[0].astype("float64")
    valid = arr != (nodata if nodata is not None else np.nan)
    if nodata is None:
        valid = np.isfinite(arr)
    manual_mean = arr[valid].mean()

    zs_mean = zonal_stats(
        polygons.iloc[[idx]], raster_path,
        stats=["mean"], nodata=nodata, all_touched=False,
    )[0]["mean"]

    assert abs(manual_mean - zs_mean) < tol, (
        f"Mismatch: manual={manual_mean:.6f} vs zonal_stats={zs_mean:.6f}"
    )
    print(f"Polygon {idx} verified: mean = {zs_mean:.6f}")

verify_polygon_mean(polygons, "ndvi.tif", nodata, idx=0)

Two more cheap assertions catch the most common silent failures:

records = zonal_stats(polygons, "ndvi.tif", stats=["count", "mean"], nodata=nodata)

# 1. No polygon should have collapsed to zero pixels unexpectedly.
empty = [i for i, r in enumerate(records) if r["count"] == 0]
assert not empty, f"Polygons with zero valid pixels: {empty} — check CRS / nodata / size"

# 2. Every polygon should have produced a numeric mean.
none_means = [i for i, r in enumerate(records) if r["mean"] is None]
assert not none_means, f"Polygons returning None mean: {none_means}"

print(f"All {len(records)} polygons produced valid statistics")

If the empty list is non-empty, the failure is almost always one of the four issues in the next section. Because these features feed a model, undetected None or zero-count rows will propagate as NaN into training and either crash the fit or silently bias it, so gate them here.

Troubleshooting & Common Errors

Every polygon returns None — CRS mismatch

Cause: The polygons and raster are in different CRSs, so the geometries fall outside the raster extent and intersect no pixels. rasterstats does not warn; it just returns None for mean, std, and friends.

Fix: Reproject the polygons to the raster CRS before calling zonal_stats (Step 1). Confirm the overlap by checking that polygons.total_bounds lies inside the raster bounds:

with rasterio.open("ndvi.tif") as src:
    rb = src.bounds
pb = polygons.total_bounds  # (minx, miny, maxx, maxy)
assert pb[0] >= rb.left and pb[2] <= rb.right and \
       pb[1] >= rb.bottom and pb[3] <= rb.top, "Polygons fall outside raster extent"

all_touched over- or under-counts edge pixels

Cause: With all_touched=False, a pixel is included only if its centre is inside the polygon, so thin or small polygons capture too few pixels. With all_touched=True, every pixel the boundary merely grazes is included, inflating counts and letting edge pixels be shared between neighbouring polygons.

Fix: Use all_touched=False for area-weighted continuous means on reasonably large polygons; switch to all_touched=True only for small or sliver geometries that would otherwise return zero pixels. For strict area-weighting, use exactextract (see Performance Optimisation), which weights partial pixels by their true overlap fraction instead of an all-or-nothing rule.

Tiny polygons smaller than one pixel return None

Cause: If no pixel centre falls inside the polygon, all_touched=False yields zero pixels and a None mean. This is common for point-buffer parcels or fine cadastral slivers against a coarse raster.

Fix: Set all_touched=True so any grazed pixel counts, guaranteeing at least one value. When even that is unreliable, resample the raster to a finer resolution before aggregating, or sample the raster at each polygon centroid as a fallback:

import rasterio

def centroid_fallback(polygons, raster_path):
    """Sample the raster at each polygon centroid as a last-resort value."""
    coords = [(p.centroid.x, p.centroid.y) for p in polygons.geometry]
    with rasterio.open(raster_path) as src:
        return [v[0] for v in src.sample(coords)]

Nodata pixels contaminate the mean

Cause: The nodata argument was omitted or set incorrectly, so sentinel values such as -9999 or 0 are treated as real data and drag the mean toward the sentinel.

Fix: Always read the nodata value from the raster (src.nodata) and pass it to zonal_stats. When the raster header lacks a nodata tag, supply it explicitly. A mean far outside the raster’s physical range (for example, a negative NDVI mean where NDVI should be 0–1) is the tell-tale sign that nodata leaked in.

rasterio / GDAL cannot open the raster

Cause: A missing GDAL system library or an unsupported driver, surfacing as RasterioIOError.

Fix: Confirm the GDAL apt packages above are installed and that the file path and driver are valid. gdalinfo ndvi.tif should print a header; if it fails, the problem is the environment, not the Python code.

Feature columns are misaligned after the join

Cause: The polygons were filtered or sorted between calling zonal_stats and building the DataFrame, so positional order no longer matches.

Fix: Call reset_index(drop=True) immediately before zonal_stats and build the stats DataFrame with index=polygons.index, exactly as in Step 4. Never reorder the layer in between.

Performance Optimisation

Pre-clip the raster to the polygon extent

zonal_stats reads raster windows per polygon, but if your polygons cover a small part of a large mosaic, first clip the raster to their combined bounding box so the library scans far fewer pixels.

import rasterio
from rasterio.mask import mask
from shapely.geometry import box

def clip_to_polygons(raster_path, polygons, out_path):
    """Clip a raster to the bounding box of the polygon layer to shrink I/O."""
    bbox = [box(*polygons.total_bounds).__geo_interface__]
    with rasterio.open(raster_path) as src:
        clipped, transform = mask(src, bbox, crop=True)
        meta = src.meta.copy()
    meta.update(height=clipped.shape[1], width=clipped.shape[2], transform=transform)
    with rasterio.open(out_path, "w", **meta) as dst:
        dst.write(clipped)
    return out_path

Parallelise over geometry chunks

zonal_stats is CPU-bound and releases work per geometry, so splitting the polygons into chunks and processing them across a process pool scales nearly linearly with cores. Each worker opens its own raster handle.

import numpy as np
import pandas as pd
from concurrent.futures import ProcessPoolExecutor
from rasterstats import zonal_stats

def _chunk_stats(args):
    geoms, raster_path, stats, nodata = args
    return zonal_stats(geoms, raster_path, stats=stats, nodata=nodata)

def parallel_zonal(polygons, raster_path, stats, nodata, n_workers=4):
    """Compute zonal statistics across a process pool, preserving input order."""
    chunks = np.array_split(polygons.reset_index(drop=True), n_workers)
    tasks = [(c, raster_path, stats, nodata) for c in chunks]
    with ProcessPoolExecutor(max_workers=n_workers) as ex:
        results = list(ex.map(_chunk_stats, tasks))
    flat = [rec for chunk in results for rec in chunk]
    return pd.DataFrame(flat)

Because np.array_split preserves order and the results are concatenated in the same order, the row-to-polygon contract from Step 4 still holds.

Use exactextract for large jobs and exact weighting

For large polygon sets, or whenever you need true area-weighted statistics rather than the centre-in/centre-out rule, the exactextract library computes exact pixel-polygon intersection fractions and runs several times faster than rasterstats on big workloads. It is the right tool when edge pixels materially affect the aggregate — coarse rasters, small polygons, or high-precision agronomic features. The rasterstats workflow above remains the most convenient choice for moderate polygon counts and quick iteration. For downstream normalisation of the resulting columns before modelling, see feature scaling for geospatial inputs.

FAQ

Why does zonal_stats return None for some polygons?

A None result means no valid pixels intersected the polygon. The usual causes are a CRS mismatch that places the polygon outside the raster extent, a polygon smaller than a single pixel, or every overlapping pixel being nodata. Reproject the polygons to the raster CRS first, and set all_touched=True for tiny geometries so any grazed pixel is counted.

Should I use all_touched=True or all_touched=False?

Use all_touched=False (the default) for area-representative continuous statistics, so only pixels whose centre falls inside the polygon are counted and edge pixels are not shared between neighbours. Use all_touched=True for small polygons or line-like features that would otherwise capture zero pixels, accepting that boundary pixels get over-counted. When precise partial-pixel weighting matters, prefer exactextract.

How do I compute zonal statistics for a categorical land-cover raster?

Pass categorical=True to zonal_stats, or request the majority statistic. This returns the dominant class per polygon plus a per-class pixel count, instead of a meaningless numeric mean over class codes. Feed the majority class to tree models directly, or encode it first following encoding categorical geographic features.

What is faster than rasterstats for very large polygon sets?

The exactextract library computes area-weighted zonal statistics using exact pixel-polygon intersection and is several times faster than rasterstats on large jobs. Pre-clipping the raster to the polygon extent and parallelising over geometry chunks — both shown above — give further large speedups without changing the output.

Related

• Raster Band Math and Index Calculation — produce the NDVI and spectral-index surfaces that zonal statistics aggregates
• Vector Proximity and Buffer Generation — distance-based features that complement per-polygon aggregates
• Spatial Lag and Neighborhood Statistics — summarise neighbouring polygons rather than underlying pixels
• CRS Alignment and Projection Handling — the reprojection step that makes zonal overlays valid
• Computing Zonal Statistics with rasterstats — a focused deep dive on the rasterstats API

Part of: Spatial Feature Engineering for Machine Learning