Computing Zonal Statistics with rasterstats

Compute per-polygon zonal statistics — mean, std, median, percentiles — from a raster with rasterstats, handle nodata and all_touched, and join results to a GeoDataFrame.

TL;DR Answer

Call rasterstats.zonal_stats(vectors, raster, stats=[...]) to reduce every pixel inside each polygon to a summary value. Pass stats=["mean", "std", "median", "percentile_90"] for a continuous band, categorical=True for a land-cover raster to get a per-class pixel tally and a majority label, and always supply nodata so masked pixels never enter the aggregation. The function returns one dictionary per input feature in the original order, so you can turn it into a DataFrame and attach it to the source GeoDataFrame by position. This is the concrete, hands-on companion to Zonal Statistics for Polygon Aggregation, which covers the underlying concepts.

from rasterstats import zonal_stats

stats = zonal_stats(
    "fields.geojson",
    "ndvi.tif",
    stats=["mean", "std", "median", "percentile_90"],
    nodata=-9999.0,
    all_touched=False,
)
print(stats[0])
# {'mean': 0.61, 'std': 0.08, 'median': 0.63, 'percentile_90': 0.71}

Part of: Zonal Statistics for Polygon Aggregation

The zonal_stats Call Surface

rasterstats.zonal_stats takes two positional arguments and a handful of keywords that decide exactly which pixels are summarised and how. Getting a feature table you can trust for model training is almost entirely about setting these keywords correctly rather than writing much code.

  • vectors — a path, a GeoJSON-like mapping, or an iterable of geometries. Passing a file path is convenient, but for a production pipeline you usually already hold a GeoDataFrame in memory and want to feed its geometries directly.
  • raster — a path to a GDAL-readable file or an in-memory (array, affine) pair. When you pass a path, rasterstats reads only the windows it needs, which keeps memory flat even for a large mosaic.
  • stats — a list of reducer names. Continuous bands support count, min, max, mean, sum, std, median, majority, minority, unique, range, and percentile_<q> for any q. Ask only for what a downstream model consumes — every extra statistic is another pass over the pixel mask.
  • nodata — the sentinel value to exclude. If the raster’s own nodata tag is reliable you can omit it, but passing it explicitly removes any ambiguity and is the single most common fix for wrong means.
  • all_touched — False (default) counts pixels whose centroid is inside the polygon; True counts every pixel the boundary intersects.
  • categorical — True switches the reducer to a per-value pixel count, ideal for discrete land-cover or soil-class rasters.

A percentile_90 in the list gives you the 90th percentile of the in-polygon distribution — far more robust to a handful of anomalous pixels than the raw max, and a strong feature for models that care about the upper tail of a spectral index.

How a Polygon Becomes One Feature Row

The diagram below traces a single field polygon through the aggregation. rasterstats rasterises the polygon against the raster grid to build a boolean mask, applies the nodata mask on top, then runs each requested reducer over the surviving pixels to emit one row of features.

zonal_stats data flow: polygon plus raster to feature row A polygon geometry and a raster grid are overlaid to form a masked pixel set, nodata pixels are dropped, and the remaining pixels are reduced by mean, std, median and percentile_90 into a single feature row. 1. Inputs raster grid + polygon 2. Pixel mask green = in-polygon red = nodata (dropped) 3. Reducers mean std median percentile_90 4. One feature row joined back to the GeoDataFrame (order preserved) field_id A-17 ndvi_mean 0.61 ndvi_std 0.08 ndvi_median 0.63 ndvi_p90 0.71

A Production-Ready Typed Function

The helper below wraps zonal_stats with the guarantees a training pipeline needs: a CRS check so the vectors and raster genuinely overlap, an explicit nodata, prefixed column names so several bands can be merged without collisions, and a positional join that preserves row order. It returns a new GeoDataFrame rather than mutating the input.

from __future__ import annotations

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

DEFAULT_STATS = ["mean", "std", "median", "percentile_90"]


def add_zonal_features(
    gdf: gpd.GeoDataFrame,
    raster_path: str,
    stats: list[str] = DEFAULT_STATS,
    prefix: str = "z",
    nodata: float | None = None,
    all_touched: bool = False,
    band: int = 1,
) -> gpd.GeoDataFrame:
    """Attach per-polygon zonal statistics from a raster to a GeoDataFrame.

    Parameters
    ----------
    gdf : gpd.GeoDataFrame
        Polygons to aggregate over. Must carry a defined CRS.
    raster_path : str
        Path to a GDAL-readable raster.
    stats : list[str]
        Reducers passed to rasterstats (e.g. mean, std, percentile_90).
    prefix : str
        Column-name prefix so multiple bands can be merged safely.
    nodata : float or None
        Sentinel value to exclude. Falls back to the raster's own tag.
    all_touched : bool
        If True, count every pixel the polygon boundary intersects.
    band : int
        1-based raster band index to read.

    Returns
    -------
    gpd.GeoDataFrame
        A copy of gdf with one new column per requested statistic.
    """
    if gdf.crs is None:
        raise ValueError("Input GeoDataFrame has no CRS; set one before aggregating.")

    with rasterio.open(raster_path) as src:
        raster_crs = src.crs
        raster_nodata = src.nodata

    if raster_crs is not None and gdf.crs != raster_crs:
        raise ValueError(
            f"CRS mismatch: polygons are {gdf.crs}, raster is {raster_crs}. "
            "Reproject the polygons to the raster CRS first."
        )

    results = zonal_stats(
        gdf.geometry,
        raster_path,
        stats=stats,
        nodata=nodata if nodata is not None else raster_nodata,
        all_touched=all_touched,
        band=band,
        geojson_out=False,
    )

    stats_df = pd.DataFrame(results).add_prefix(f"{prefix}_")

    out = gdf.reset_index(drop=True).copy()
    out = pd.concat([out, stats_df.reset_index(drop=True)], axis=1)
    return gpd.GeoDataFrame(out, geometry=gdf.geometry.name, crs=gdf.crs)

The pd.concat on axis=1 after a shared reset_index(drop=True) is what guarantees each statistic lands on the polygon it was computed from. Because zonal_stats returns results in input order, aligning by position is both correct and cheaper than any spatial join. Reproject the polygons up front if the CRS check fails — the same discipline covered under CRS Alignment and Projection Handling applies here.

Walkthrough: Mean NDVI per Field

1. Load the polygons and confirm the CRS matches the raster

import geopandas as gpd
import rasterio

fields = gpd.read_file("fields.geojson")

with rasterio.open("ndvi.tif") as src:
    print("raster CRS:", src.crs)          # EPSG:32633
    print("raster nodata:", src.nodata)    # -9999.0

fields = fields.to_crs("EPSG:32633")       # align to the raster

If your polygons come from a source with self-intersections or ring errors, clean them before this step — invalid geometries rasterise unpredictably. See Fixing Invalid Geometries Before Buffering for a robust repair recipe that applies equally to zonal work.

2. Compute the statistics with the typed helper

fields = add_zonal_features(
    fields,
    "ndvi.tif",
    stats=["mean", "std", "median", "percentile_90"],
    prefix="ndvi",
    nodata=-9999.0,
    all_touched=False,
)

print(fields[["field_id", "ndvi_mean", "ndvi_std", "ndvi_p90"]].head())

The NDVI band itself is produced by band math over red and near-infrared reflectance; if you still need to generate it, How to Calculate NDVI and EVI with rasterio walks through that step.

3. Add a categorical majority for land cover

For a discrete raster — a land-cover map where each integer codes a class — switch to categorical=True. rasterstats then returns a pixel count per class, and adding majority to stats gives you the dominant class label per polygon.

from rasterstats import zonal_stats

lc = zonal_stats(
    fields.geometry,
    "landcover.tif",
    categorical=True,
    stats=["majority"],
    nodata=0,
)
print(lc[0])
# {'majority': 3, 2: 41, 3: 610, 5: 12}

The integer keys are the class codes with their pixel tallies; majority (here class 3, cropland) is the single most common class inside that field — a compact categorical feature for a classifier.

4. Handle polygons smaller than a pixel

Small fields whose centroid misses every cell return None. Re-run those with all_touched=True so any pixel the boundary clips is included:

missing = fields["ndvi_mean"].isna()
if missing.any():
    patch = zonal_stats(
        fields.loc[missing, "geometry"],
        "ndvi.tif",
        stats=["mean"],
        nodata=-9999.0,
        all_touched=True,
    )
    fields.loc[missing, "ndvi_mean"] = [r["mean"] for r in patch]

Verification

Never trust an aggregation you have not spot-checked. Pick one polygon, mask the raster to it by hand with rasterio, and confirm the manual mean matches what rasterstats reported. Agreement to several decimal places proves your nodata and all_touched settings are doing what you expect.

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

def manual_masked_mean(raster_path: str, geometry, nodata: float) -> float:
    """Mean of valid pixels inside a single geometry, computed directly."""
    with rasterio.open(raster_path) as src:
        out_image, _ = mask(src, [geometry], crop=True, all_touched=False)
    band = out_image[0].astype("float64")
    valid = band[(band != nodata) & ~np.isnan(band)]
    return float(valid.mean())

row = fields.iloc[0]
manual = manual_masked_mean("ndvi.tif", row.geometry, nodata=-9999.0)
reported = row["ndvi_mean"]

print(f"manual: {manual:.6f}  rasterstats: {reported:.6f}")
assert np.isclose(manual, reported, atol=1e-6), "mismatch — check nodata/all_touched"

If the two disagree, the usual causes are a nodata value that rasterstats never received, an all_touched setting that differs between the two calls, or a lingering CRS mismatch that shifts the mask by a pixel. Once the aggregated table is verified, it is ready to feed a model alongside distance and neighbourhood features — for example the pairwise structure built in Creating Distance Matrices for Spatial Features.

FAQ

Why do my rasterstats results contain None or NaN for some polygons?

A None result means no valid raster pixels fell inside that polygon. It happens when the polygon lies entirely over nodata, when it is smaller than a single pixel and its centroid misses every cell, or when it sits outside the raster extent. Set all_touched=True to include every pixel the boundary intersects, confirm the nodata value is passed correctly, and verify both layers share the same CRS before aggregating.

Does rasterstats preserve the row order of my input polygons?

Yes. zonal_stats returns one dictionary per input feature in the exact order the features were read, so you can build a DataFrame from the result and attach it to the original GeoDataFrame by position with pd.concat on axis=1 after a shared reset_index. Never join on a spatial predicate or a mutable label index — positional alignment is what keeps each statistic matched to its polygon.

What is the difference between all_touched=True and the default pixel selection?

By default rasterstats includes only pixels whose centroid falls inside the polygon, which can miss narrow or small features. With all_touched=True, every pixel the polygon boundary touches is counted. all_touched raises the pixel count and is safer for small polygons, but it double-counts edge pixels shared between adjacent zones, so for wall-to-wall polygons that tile a surface the default centroid rule usually gives cleaner, non-overlapping aggregates.

Part of: Zonal Statistics for Polygon Aggregation — Spatial Feature Engineering for Machine Learning