Computing Zonal Statistics with rasterstats

To summarise raster values inside each polygon of a vector layer, reproject the geometries into the raster’s CRS and call zonal_stats:

import geopandas as gpd
import rasterio
from rasterstats import zonal_stats

with rasterio.open("ndvi_20230615.tif") as src:
    zones = gpd.read_file("fields.gpkg").to_crs(src.crs)

stats = zonal_stats(
    zones.geometry,
    "ndvi_20230615.tif",
    stats=["mean", "median", "std", "count"],
    all_touched=False,          # centre-based inclusion; state it explicitly
    nodata=None,                # use the file's declared value
)
zones = zones.join(gpd.pd.DataFrame(stats))

This is the library route through the workflow described in Zonal Statistics and Vector–Raster Integration.

Why This Arises in Remote Sensing Workflows

Almost every operational use of satellite imagery ends in a table. A crop monitoring service reports mean NDVI per parcel; a fire agency reports burned hectares per district; a water authority reports surface extent per reservoir. The imagery is an intermediate, and the deliverable is a number attached to a place.

rasterstats exists because that operation is common enough to standardise and fiddly enough to get wrong. It handles the window derivation, the geometry rasterisation and the nodata exclusion in one call, and it returns plain dictionaries that join cleanly onto a dataframe.

What it does not do is guess your conventions. It will not reproject your geometries, it will not decide the boundary rule for you, and it will not tell you that a zone returned statistics from three pixels rather than three thousand unless you ask for count. Those are exactly the decisions covered here.

Environment & Setup

pip install "rasterstats>=0.19" "geopandas>=0.13" "rasterio>=1.3.0"

Complete Working Example

This function computes statistics for one raster and returns a tidy dataframe, with the identifier, the date and the pixel count preserved — the shape a multi-date run needs.

from pathlib import Path

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

STATS = ["mean", "median", "std", "min", "max", "count"]


def zonal_table(
    raster_path: str,
    zones: gpd.GeoDataFrame,
    *,
    id_field: str = "field_id",
    date: str | None = None,
    all_touched: bool = False,
    band: int = 1,
) -> pd.DataFrame:
    """Zonal statistics for one raster as a tidy table, one row per zone."""
    with rasterio.open(raster_path) as src:
        if zones.crs != src.crs:
            # Move the vectors, never the pixels
            zones = zones.to_crs(src.crs)

    records = zonal_stats(
        zones.geometry,
        raster_path,
        stats=STATS,
        all_touched=all_touched,
        band=band,
        geojson_out=False,
        nodata=None,          # honour the file's declared nodata
    )

    table = pd.DataFrame(records)
    table.insert(0, id_field, zones[id_field].to_numpy())
    table["date"] = date or Path(raster_path).stem.split("_")[-1]
    table["all_touched"] = all_touched       # record the convention in the data
    return table


if __name__ == "__main__":
    fields = gpd.read_file("fields.gpkg")
    df = zonal_table("ndvi_20230615.tif", fields, date="2023-06-15")
    print(df.head())
    print("zones with no observed pixels:", int((df["count"] == 0).sum()))

Two lines are load-bearing beyond the call itself. Reprojecting inside the function means callers cannot forget, and the failure mode it prevents — every zone returning None — looks like a data problem rather than a CRS problem. Recording all_touched in the output means a dataset assembled by several people over several months can be audited rather than argued about.

Variant Patterns

1. Categorical rasters: class fractions per zone

For a land-cover map or a scene classification layer, the useful summary is a fraction per class, not a mean.

import pandas as pd
from rasterstats import zonal_stats

counts = zonal_stats(zones.geometry, "landcover.tif", categorical=True, all_touched=False)

frac = pd.DataFrame(counts).fillna(0)
frac = frac.div(frac.sum(axis=1), axis=0)      # counts → fractions per zone
frac.columns = [f"class_{int(c)}_frac" for c in frac.columns]

Because the classes are codes, all_touched changes both the numerator and the denominator, so the fractions stay well-formed either way — but the numbers differ between conventions, which is why the earlier advice to record it applies with extra force here. The class semantics for Sentinel-2 are covered in Masking Clouds with the Sentinel-2 SCL Band.

2. Custom statistics: percentiles, thresholds, coverage

add_stats accepts a callable that receives the masked array for each zone, which covers everything the built-in list does not.

import numpy as np
from rasterstats import zonal_stats


def p90(masked):
    """90th percentile over the valid pixels of a zone."""
    return float(np.percentile(masked.compressed(), 90)) if masked.count() else None


def frac_above_05(masked):
    """Fraction of valid pixels above 0.5 — e.g. vegetated fraction from NDVI."""
    if not masked.count():
        return None
    return float((masked.compressed() > 0.5).mean())


stats = zonal_stats(
    zones.geometry, "ndvi_20230615.tif",
    stats=["count"], add_stats={"p90": p90, "veg_frac": frac_above_05},
)

The array arrives already masked, so compressed() gives only the valid pixels and the guard on masked.count() is what keeps a fully clouded zone from raising.

3. Many dates, one geometry set

import pandas as pd

fields = gpd.read_file("fields.gpkg")
with rasterio.open(rasters[0]) as src:
    fields = fields.to_crs(src.crs)          # once, outside the loop

frames = [zonal_table(path, fields, date=date) for date, path in dated_rasters]
tidy = pd.concat(frames, ignore_index=True)
tidy.to_parquet("field_ndvi_2023.parquet", partition_cols=["date"])

Writing partitioned parquet means a later query for one date reads one partition, and an interrupted run leaves the completed dates on disk.

Reading the Output Critically

A zonal table is easy to produce and easy to over-trust. Three columns turn it from a list of numbers into something defensible.

The pixel count is the first, and it should be plotted before anything else is concluded. Zones whose count collapses on a particular date are telling you about cloud cover, not about vegetation, and a time series that ignores the count will show a “dip” that is really an absence. A simple rule that survives review: discard any zone-date whose count falls below a fixed fraction — often a quarter — of that zone’s typical count.

The standard deviation is the second. A zone with a high internal spread is a zone where the mean describes little: a parcel that is half harvested, a district containing both a lake and a city. Reporting the mean alone in those cases is not wrong so much as uninformative, and the standard deviation is what lets a reader see it.

The convention flag is the third, and it only matters when datasets are combined. Two analysts running the same code with different all_touched settings produce tables that look identical and disagree systematically on small zones. Carrying the flag in the data makes that visible in a join instead of invisible in a report — the same argument for recording provenance that applies to Writing and Validating Cloud-Optimized GeoTIFFs.

One last habit is worth adopting: keep the zone identifiers as they arrived, including their type. Silently converting a string identifier to an integer, or reordering the geometry list without carrying the identifier alongside, produces a table that joins cleanly onto the wrong rows. Because zonal_stats returns results positionally, that failure is entirely silent, and it is the single most damaging mistake available in this workflow.

Common Errors

Every zone returns None for every statistic

The geometries are in a different CRS from the raster, so their windows fall outside the scene. Compare zones.total_bounds with src.bounds; the magnitude mismatch is usually obvious, and the diagnostic is the same one used in Fixing EPSG Mismatches in rasterio.open.

The mean is dragged toward a large negative number

The file declares no nodata, so the sentinel is being treated as data. Pass nodata=-9999 explicitly for that raster, and fix the declaration at the source if you own the write step.

zonal_stats is slow on a large geometry set

It reopens the raster per geometry when handed a path and a long list. For big runs, pass the raster once per scene and iterate geometries inside, or use the windowed loop from the parent topic.

Frequently Asked Questions

Q: Does rasterstats reproject geometries for me? No. It assumes the geometries are already in the raster’s CRS and will silently return empty statistics if they are not, because the geometries fall outside the raster’s extent.

Q: How do I get the fraction of each land-cover class per zone? Pass categorical=True. rasterstats then returns a dictionary of class code to pixel count per zone, which you divide by the total count to get fractions.

Q: Why is the count lower than the polygon area suggests? Nodata pixels are excluded from the count. In a cloud-masked scene that is the intended behaviour, and the count is exactly the field that tells you how much of the zone was observed.

Related

• Zonal Statistics and Vector–Raster Integration — the parent topic, including the hand-written loop and weighting.
• Rasterizing Vector Polygons onto a Raster Grid — the inverse operation, and the source of the masks used here.
• Extracting Pixel Values at Point Locations — the point-sampling case.
• Temporal Aggregation and Time-Series Analysis — what to do with the multi-date table this produces.