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.

What a zonal summary reduces Three stand polygons of different sizes lie over a raster grid. The pixels whose centres fall inside each polygon are shaded, and the counts differ sharply: the largest stand contributes ninety-four pixels, the middle one twenty-two, and the smallest only three. A table beside them lists the mean and standard deviation each polygon yields, with a note that the three-pixel stand's statistics carry almost no information. stand A 94 px · mean 21.4 · sd 4.8 stand B 22 px · mean 18.9 · sd 5.6 stand C 3 px · mean 24.1 · sd 9.2 too few pixels to report Always return the pixel count — a mean over three cells and a mean over ninety look identical without it.

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,
    })
Three ways to summarise a polygon smaller than the pixel A long narrow riparian polygon crosses a coarse raster grid diagonally. Under centre-based selection no pixel centre falls inside it, so the summary returns nothing. Under all-touched selection every pixel the polygon grazes is included at full weight, so the summary describes mostly upland. Under fractional weighting each pixel contributes in proportion to the area of overlap, which is the only option that returns a defensible number. a 30 m riparian strip on a 50 m grid centre-based 0 pixels — returns nothing all touched 12 pixels at full weight — mostly upland fractional weighting 12 pixels weighted 0.04–0.38 — defensible When polygons approach pixel size, the weighting scheme is the analysis.

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.

Reopening the raster per polygon against one windowed pass Runtime is plotted against the number of polygons summarised. The naive approach, which reopens and reads the raster for every polygon, rises steeply and reaches several minutes by five thousand polygons. The windowed approach, which opens once and reads only the window each polygon needs, stays nearly flat. The gap widens continuously, so the choice matters more as inventories grow. 0 1000 2500 4000 5000 polygons summarised runtime reopen per polygon one open, windowed reads Sort polygons by position first and consecutive windows will hit the same cached blocks.

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.