Forestry GIS Python · Forestry · Ecology

Zonal Statistics for Forest Inventory with rasterstats

Forest inventory constantly asks the same question of a raster: what does this surface look like inside each stand, plot, or ownership parcel? Mean canopy height per stand, total biomass per parcel, the 95th-percentile height per plot, or the dominant land-cover class per management unit — all of these are zonal statistics, a per-polygon reduction of raster cells. The rasterstats library computes them in one call, reading the raster window-by-window so a national mosaic never lands in RAM at once. This guide covers that narrow task for inventory work, including a rioxarray/numpy alternative when you need full control. It sits within the Raster-Vector Overlay Techniques workflow, part of the broader Ecological GIS Data Foundations in Python framework.

When to Use rasterstats vs a rioxarray/numpy Approach

Both approaches answer “summarise a raster inside these polygons,” but they trade convenience against control. rasterstats.zonal_stats is the fast path: it handles rasterisation, nodata, and windowed reads internally and returns a list of dicts. A hand-rolled rioxarray + numpy approach is worth it when you need a reduction rasterstats does not expose (a weighted percentile, a custom biomass allometry per cell), or when the raster is already an in-memory xarray cube from an earlier step.

Requirement Use rasterstats Use rioxarray + numpy
Standard reducers (mean, sum, min, max, std, count) ✔ one call ✔ but you write the reduction
Percentiles per zone ✔ percentiles=[...] ✔ via np.nanpercentile
Categorical majority / class counts ✔ categorical=True ✔ via np.bincount
Raster too large for RAM ✔ windowed reads by path needs manual windowing
Custom / weighted reducer limited (add_stats) ✔ full control
Raster already an xarray cube reopen from disk ✔ operate in place
Exact area-weighted (fractional cells) ✗ (centroid rule) use exactextract instead

For routine per-stand inventory metrics, reach for rasterstats. Drop to rioxarray only when the reducer is non-standard or the data already lives in memory. If you need true fractional-cell weighting on small polygons, neither is exact — use a coverage-fraction engine such as exactextract, as discussed in the parent Raster-Vector Overlay Techniques guide.

Minimal Reproducible Example

This computes mean, sum, standard deviation, count, and the 90th and 95th percentiles of a canopy height model for every stand polygon, then joins the results back onto the original attribute table so stand identifiers travel with their statistics. Note the two non-negotiables: the vector must be in the raster’s CRS, and the raster must carry an explicit nodata value.

import geopandas as gpd
import rasterio
from rasterstats import zonal_stats

RASTER = "chm.tif"          # 1 m canopy height model, metres
VECTOR = "stands.gpkg"      # forest stand polygons

stands = gpd.read_file(VECTOR)

# 1. Align the vector to the raster CRS (never warp the raster for one overlay).
with rasterio.open(RASTER) as src:
    raster_crs = src.crs
    nodata = src.nodata
if nodata is None:
    raise ValueError("Raster nodata is None; set it before zonal_stats.")
stands = stands.to_crs(raster_crs)

# 2. Reduce the raster inside every polygon.
stats = zonal_stats(
    vectors=stands,
    raster=RASTER,                       # pass the PATH so reads are windowed
    stats=["mean", "sum", "std", "count"],
    percentiles=[90, 95],                # adds percentile_90, percentile_95
    nodata=nodata,
    all_touched=False,                   # centroid rule — unbiased for large plots
    geojson_out=False,
)

# 3. Join the statistics back onto the stand attributes.
for key in ("mean", "sum", "std", "count", "percentile_90", "percentile_95"):
    stands[f"chm_{key}"] = [s.get(key) for s in stats]

stands.to_file("stands_with_chm.gpkg", driver="GPKG")

For a categorical raster — a land-cover, forest-type, or fuel-class grid — switch to class tallies. categorical=True returns, per zone, a dict of {class_code: pixel_count}; a small helper turns that into the majority class and its purity fraction, which is the usual inventory need.

from rasterstats import zonal_stats

cat = zonal_stats(
    vectors=stands, raster="landcover.tif",
    categorical=True, nodata=0, all_touched=False,
)

def majority(counts: dict) -> tuple[int | None, float]:
    """Dominant class code and its share of valid pixels in the zone."""
    if not counts:
        return None, 0.0
    total = sum(counts.values())
    code, n = max(counts.items(), key=lambda kv: kv[1])
    return code, n / total

stands["lc_majority"], stands["lc_purity"] = zip(*(majority(c) for c in cat))

The rioxarray/numpy alternative, useful when the reducer is custom, rasterises each geometry to a mask and reduces the masked array directly:

import numpy as np
import rioxarray
from rasterio.features import geometry_mask

chm = rioxarray.open_rasterio("chm.tif", masked=True).squeeze("band", drop=True)
transform = chm.rio.transform()
records = []
for _, row in stands.to_crs(chm.rio.crs).iterrows():
    mask = geometry_mask([row.geometry], out_shape=chm.shape,
                         transform=transform, invert=True)
    cells = chm.values[mask]
    cells = cells[np.isfinite(cells)]
    records.append({
        "mean": float(np.mean(cells)) if cells.size else None,
        "p95": float(np.percentile(cells, 95)) if cells.size else None,
    })

Parameter Reference

These zonal_stats arguments materially change the result; pin them deliberately rather than relying on defaults.

Argument Type Default Recommended Ecological rationale
vectors GeoDataFrame / path — GeoDataFrame in raster CRS Geometries must match the raster grid or every zone reads nodata
raster path / ndarray — file path A path enables windowed reads; an array forces the whole raster into RAM
stats list of str ["count","min","max","mean"] name reducers explicitly Documents intent; add sum for biomass, std for heterogeneity
percentiles list of number none [95] for canopy Percentiles resist outliers better than max for dominant height
nodata number raster’s nodata must equal raster fill If wrong, fill values are averaged as real data
all_touched bool False False for stands Centroid rule; True over-includes edge cells and inflates thin polygons
categorical bool False True for class rasters Returns per-class pixel counts instead of summary statistics
add_stats dict none for custom reducers Register a named callable (e.g. a weighted mean) without leaving the library

Expected Output and Verification

zonal_stats returns a list the same length and order as the input vectors, one dict per feature. Three checks confirm the numbers are real before you report them: every zone has a positive pixel count, no statistic exceeds the raster’s physical ceiling, and the joined frame did not silently reorder.

import numpy as np

chm_count = np.array([s["count"] for s in stats])
chm_max = np.array([s["max"] if s["max"] is not None else np.nan for s in stats])

# 1. No empty zones (a zero count means off-grid or all-nodata).
assert (chm_count > 0).all(), "Some stands returned zero valid pixels — check CRS/nodata."

# 2. Physical plausibility (a 1 m CHM should never report 400 m).
assert np.nanmax(chm_max) < 120.0, "Implausible height — nodata likely leaked in."

# 3. Row alignment preserved after the join.
assert len(stats) == len(stands), "Stats length != feature count."

A zero count almost always means a CRS mismatch left the polygon off-grid, or the stand fell entirely on nodata; either way, fix it before reporting rather than treating the None as a real value. Once verified, these per-stand heights commonly pair with a Canopy Height Model Creation product upstream and feed forest-health checks via Vegetation Index Calculation in Python.

Common Pitfalls

  • nodata left as None. rasterstats then counts the fill value (often -9999 or 0) as data, dragging means toward the fill. Read src.nodata, assert it is set, and pass it explicitly.
  • CRS mismatch. If the vector and raster CRS differ, polygons land off-grid and return all-nodata zones. Reproject the vector to the raster CRS before the call, never the reverse for a single overlay.
  • all_touched=True on plots. Setting it inflates small or thin polygons — narrow riparian buffers, fixed-radius plots a few cells wide — by pulling in every edge cell. Use the centroid rule (False) unless omission is worse than over-inclusion.
  • Passing an in-memory array for a huge raster. Handing zonal_stats a numpy array forces the whole raster into RAM. Pass the file path so it reads only the window each polygon needs.

Frequently Asked Questions

How do I compute a statistic rasterstats does not provide, like a weighted mean?

Use the add_stats argument to register a named callable that receives the masked pixel array for each zone and returns a scalar. For example, add_stats={"p10": lambda a: float(np.nanpercentile(a.compressed(), 10))} adds a p10 field. For reductions that need per-cell weights (an area or biomass weighting), the rioxarray/numpy approach gives you the mask and the weights together and is usually clearer.

Why does a stand return None for its mean?

The zone contained no valid pixels — every cell was nodata, or the polygon fell off the raster after a CRS mismatch. Confirm the vector was reprojected to the raster CRS, that the raster has an explicit nodata value, and that the polygon actually overlaps the raster extent. A non-zero count in the result is the signal that a mean is trustworthy.