Computing Zonal Index Statistics per Plot Polygon

A variety trial has 420 plots. The index raster is 24,000 pixels square. The obvious implementation — rasterise each plot polygon over the raster extent, mask, summarise — does 420 full-extent rasterisations of a grid that is 99.9 % irrelevant to each one, and takes most of an hour.

The same job done windowed takes about twelve seconds. This page covers that, along with the parts that matter more than speed: which statistics to compute, how to handle plots that are partly outside the flight, and how to produce a table that actually joins to the previous flight’s. It is the summarising step of computing vegetation indices with rasterio.

Why the windowed form is not just faster

Reading only each plot’s own window changes the scaling from plots × raster area to total plot area, which on a trial is a factor of hundreds. That alone would justify it, but two other properties matter more.

Memory stays bounded. A full-extent boolean mask for a 24,000-square raster is 576 MB per plot; the windowed version allocates kilobytes. On a machine also holding the raster, the difference decides whether the job runs.

Partial plots are handled honestly. A plot extending past the flight boundary has real, partial coverage. The windowed read with boundless=True returns that partial coverage with NaN outside, which is the truth; a naive implementation either raises or silently clips.

Full-extent against windowed zonal statistics Two panels over the same raster. On the left, a full-extent approach rasterises a mask the size of the whole raster for each of four hundred and twenty plots, reading five hundred and seventy-six megabytes per plot and taking about fifty minutes in total. On the right, a windowed approach reads only each plot's own bounding window, a few tens of kilobytes each, and completes in about twelve seconds. A note records that the windowed version also handles plots partly outside the flight boundary correctly, returning their partial coverage rather than raising. full extent per plot windowed per plot 576 MB read · 420 times ≈ 50 minutes a few tens of kB each ≈ 12 seconds Two hundred times faster, and correct on plots that cross the flight boundary. Memory stays bounded, which is what decides whether the job runs at all.

Figure 1 — The windowed form, and the three things it buys.

Minimal reproducible solution

import numpy as np
import rasterio
from rasterio.features import geometry_mask
from rasterio.windows import from_bounds


def zonal_index_stats(index_path: str, plots: list[dict], *,
                      nir_path: str | None = None,
                      min_pixels: int = 50) -> list[dict]:
    """Per-plot index statistics, reading one window per plot.

    Each plot carries an id and a shapely geometry in the raster's CRS.
    Plots partly outside the raster return their partial coverage with a
    flag, because a partly covered plot is a real observation with a caveat
    rather than an error.
    """
    results = []
    with rasterio.open(index_path) as src:
        raster_bounds = src.bounds
        for plot in plots:
            geom = plot["geometry"]
            minx, miny, maxx, maxy = geom.bounds
            outside = (maxx < raster_bounds.left or minx > raster_bounds.right
                       or maxy < raster_bounds.bottom or miny > raster_bounds.top)
            if outside:
                results.append({"id": plot["id"], "usable": False,
                                "note": "plot lies entirely outside the raster"})
                continue

            win = from_bounds(minx, miny, maxx, maxy, transform=src.transform)
            data = src.read(1, window=win, masked=True, boundless=True,
                            fill_value=np.nan).filled(np.nan)
            wt = src.window_transform(win)
            inside = ~geometry_mask([geom], out_shape=data.shape, transform=wt,
                                    invert=False)

            values = data[inside & np.isfinite(data)]
            coverage = float(values.size / max(int(inside.sum()), 1))
            if values.size < min_pixels:
                results.append({"id": plot["id"], "usable": False,
                                "pixels": int(values.size), "coverage": coverage,
                                "note": "too few valid pixels"})
                continue

            q10, q25, q50, q75, q90 = np.percentile(values, [10, 25, 50, 75, 90])
            results.append({
                "id": plot["id"], "usable": True,
                "pixels": int(values.size), "coverage": coverage,
                "median": float(q50), "iqr": float(q75 - q25),
                "p10": float(q10), "p90": float(q90),
                "partial": coverage < 0.95,
            })
    return results

The coverage figure — valid pixels over plot pixels — is the field that turns a table of medians into something a scientist can filter on. A plot at 0.4 coverage has a median computed from less than half its area, and whether that is acceptable depends on the analysis rather than on the pipeline.

Which statistics to compute

Four numbers per plot cover almost every downstream use, and each answers a different question.

The median is the central tendency, robust to a poor patch that a mean would follow. The interquartile range describes within-plot variability, which in a trial is often the finding rather than the noise — a plot with a high median and a wide spread is establishing unevenly. The tenth percentile captures the worst part of the plot, which is what determines yield in a stress trial. And the pixel count with coverage is what lets somebody decide whether to trust the other three.

A mean and a standard deviation are the conventional choices and are worse for this data, because index distributions within a plot are routinely skewed and both statistics follow the tail.

import numpy as np


def compare_statistics(values: np.ndarray) -> dict:
    """Show why a median is preferred, on a plot with a poor corner."""
    v = values[np.isfinite(values)]
    return {"mean": float(np.mean(v)), "median": float(np.median(v)),
            "sd": float(np.std(v, ddof=1)),
            "iqr": float(np.percentile(v, 75) - np.percentile(v, 25)),
            "skew_indicator": float(np.mean(v) - np.median(v))}

A skew_indicator far from zero flags a plot whose mean and median disagree, which is exactly the plot where the choice of statistic changes the conclusion.

Three ways a plot mean misleads Three rows. Edge pixels partially covered by the polygon contribute a mixture of plot and alley, and on small plots these can be a substantial fraction of the total, biasing the mean toward the alley's value. Soil and shadow pixels within the plot drag the mean down by an amount that depends on canopy closure rather than on plant health, so a sparse early-season plot and an unhealthy late-season one produce similar means for different reasons. A single outlying pixel from a misregistration artefact or a nodata leak can move a small plot's mean noticeably, which is why a median or a trimmed mean is usually the better statistic. partially covered edge pixels mix plot and alley — a large fraction of a small plot soil and shadow within the plot drag the mean by canopy closure, not by plant health a single outlying pixel moves a small plot noticeably — prefer a median or trimmed mean Buffer the polygon inward, mask soil, and report a robust statistic. All three are cheap.

Figure 3 — Three biases, all of which survive a visually convincing raster.

Edge-case matrix

Situation Naive result Handling
Plot outside the raster Exception or empty Return unusable with a note
Plot partly outside Silently clipped boundless=True, report coverage
Plot smaller than a few pixels Unstable statistics Minimum pixel count
Overlapping plot polygons Pixels counted twice Acceptable; note it, or resolve the geometry
Polygon in a different CRS Wrong window, empty result Assert the CRS before starting
Multipolygon plot Only the first part masked Pass the geometry as-is; rasterio handles parts
Plot with a hole Hole filled geometry_mask honours interiors by default
All pixels masked by canopy masking No statistics Report coverage zero, not an error

Verification snippet

import numpy as np


def verify_plot_table(rows: list[dict], *, expected_ids: set,
                      min_usable_fraction: float = 0.9) -> dict:
    """Checks a plot statistics table should pass before it is used."""
    ids = {r["id"] for r in rows}
    problems = []

    missing = expected_ids - ids
    if missing:
        problems.append(f"{len(missing)} expected plots are absent: {sorted(missing)[:5]}")
    if len(ids) != len(rows):
        problems.append("duplicate plot ids in the table")

    usable = [r for r in rows if r.get("usable")]
    frac = len(usable) / max(len(rows), 1)
    if frac < min_usable_fraction:
        problems.append(f"only {frac:.0%} of plots are usable")

    partial = [r["id"] for r in usable if r.get("partial")]
    medians = np.array([r["median"] for r in usable])
    return {"plots": len(rows), "usable": len(usable),
            "partial_coverage": partial,
            "median_range": (float(medians.min()), float(medians.max()))
            if medians.size else None,
            "problems": problems}

The duplicate-id check catches a mistake that is otherwise very hard to see: a plot layer with a repeated identifier produces a table that joins to the previous flight’s table and quietly doubles some rows.

Joining across flights

The output that a trial actually consumes is one row per plot per flight, and joining those correctly needs the plot identity to be stable. Three practices make it so.

Use the trial’s own plot identifiers, not row and column positions or an index into a file. Positions change when a layer is re-exported; identifiers do not.

Store the plot geometry once with the trial and reuse it for every flight. Re-digitising introduces a few percent of change in every plot, indistinguishable from an effect.

Carry the flight date and the index provenance on every row, so a long table can be filtered and audited without reference to an external log.

from datetime import date


def to_long_table(rows: list[dict], *, flight: date, index_name: str,
                  index_version: str, masks: dict) -> list[dict]:
    """One row per plot per flight, carrying everything needed to interpret it."""
    return [{**r, "flight_date": flight.isoformat(),
             "index": index_name, "index_version": index_version,
             "mask_policy": ";".join(f"{k}={v}" for k, v in sorted(masks.items()))}
            for r in rows]
A long plot table that joins across flights A table structure with one row per plot per flight. Columns are plot identifier, flight date, index name and version, median, interquartile range, pixel count, coverage and mask policy. Three rows are shown for plot A-14 across three flight dates, with the median rising from zero point five one to zero point seven two while coverage falls from zero point ninety-eight to zero point ninety-one. A note states that the mask policy column is what makes a comparison between two rows defensible. plot flight index@ver median iqr pixels coverage mask policy A-142026-04-02ndvi@2 0.510.094 812 0.98soil=0.2 A-142026-05-11ndvi@2 0.660.074 640 0.95soil=0.2 A-142026-06-18ndvi@2 0.720.064 402 0.91soil=0.2 Every row carries what is needed to interpret it, without an external log. The mask policy column is what makes a comparison between two rows defensible. Coverage falling while the median rises is a real signal, and only visible because both are recorded.

Figure 2 — The table shape that survives a season. Every column earns its place at analysis time.

When to escalate

  • Plot identifiers differ between the field layout and the vector layer. That is a data management problem and it will silently misattribute results. Resolve it before computing anything.
  • A large share of plots are partial. The flight boundary was drawn too tight. Re-fly with a wider margin; a trial where the edge plots are partial has lost its replicates at the edges.
  • The statistics look right and the trial analysis disagrees with field observation. Check the masking composition before the index. A change in canopy fraction is the most common explanation.

← Computing Vegetation Indices with Rasterio