Zonal Statistics of Rainfall Rasters over Field Boundaries

One-sentence answer: convert the grid to millimetres of depth, subset it to the farm with the latitude slice in the product’s own order, and compute an exact area-weighted mean per boundary — then check the season total against a gauge before trusting any of it.

Context

Weather aggregation looks like the simplest step in a farm platform, and it produces more silently wrong numbers than any other. The reasons are mundane: precipitation stored as a flux rather than a depth, latitude stored descending in one product and ascending in the next, and centroid sampling that assigns a whole cell’s value to a field that only partly overlaps it. None of them raise an exception. All of them produce a plausible-looking daily series that is wrong by a factor of 86,400, or empty, or discontinuous.

This guide is the aggregation half of weather and agronomic data integration; its output is what computing growing degree days accumulates.

Three ways to aggregate a grid to a field, and what each returns The same field polygon over the same four grid cells, aggregated three ways. Exact area weighting returns 18.6 millimetres. All-touched averaging weights every intersecting cell equally and returns 18.4. Centroid sampling takes the single cell containing the centroid and returns 21.0. Exact area weighting 14 mm 21 mm 16 mm 20 mm 18.6 mm weights 14 / 39 / 11 / 36 % All-touched, equal weights 14 mm 21 mm 16 mm 20 mm 18.4 mm a 2% sliver counts as much as 39% Centroid sampling 14 mm 21 mm 16 mm 20 mm 21.0 mm jumps to 20 if the boundary shifts

Prerequisites

Beyond the parent topic’s stack: exactextract 0.2.* (or rasterstats 0.19.* as a fallback), a gridded product opened with xarray, and field boundaries in EPSG:4326 from PostGIS schema design for farm data.

Step-by-step

1. Establish the units from the variable’s own attributes, and convert once.

The order that keeps a unit error from becoming a drought Five steps: read the declared unit from the variable's attributes, convert once to millimetres of depth, subset the grid respecting its latitude ordering, aggregate with exact area weights, and assert that every daily value falls between zero and three hundred millimetres. Read the unit from the attribute never from habit Convert once to mm of depth Subset latitude order matters Weight by area exact cell fractions Assert 0–300 mm per day The range assertion is what turns a silent factor-of-86,400 error into a failure on the first day of data.

2. Subset to the farm, checking the latitude order rather than assuming it.

3. Aggregate with exact area weights.

4. Fail on missing values — a null here becomes a null season total three steps later.

5. Compare with a gauge.

PYTHON 
import geopandas as gpd
import numpy as np
import pandas as pd
import xarray as xr
from exactextract import exact_extract

SECONDS_PER_DAY = 86_400.0


def to_mm_per_day(da: xr.DataArray) -> xr.DataArray:
    """Normalise precipitation to millimetres per day from whatever the product stores."""
    units = str(da.attrs.get("units", "")).strip().lower()
    if units in ("mm", "mm/day", "mm d-1"):
        out = da
    elif units in ("m", "metre", "meter"):
        out = da * 1000.0
    elif units in ("kg m-2 s-1", "kg/m2/s", "kg m**-2 s**-1"):
        out = da * SECONDS_PER_DAY          # 1 kg m-2 == 1 mm of depth
    elif units in ("kg m-2", "kg/m2"):
        out = da
    else:
        raise ValueError(
            f"unrecognised precipitation units {units!r} — add the conversion explicitly "
            "rather than assuming millimetres")
    out.attrs["units"] = "mm"
    return out


def subset_to_farm(ds: xr.Dataset, fields: gpd.GeoDataFrame, pad_deg: float = 0.2) -> xr.Dataset:
    """Subset a lat/lon grid to the farm, handling either latitude ordering."""
    minx, miny, maxx, maxy = fields.to_crs(4326).total_bounds
    lat = ds["latitude"].values
    descending = bool(lat[0] > lat[-1])
    lat_slice = (slice(maxy + pad_deg, miny - pad_deg) if descending
                 else slice(miny - pad_deg, maxy + pad_deg))
    out = ds.sel(latitude=lat_slice, longitude=slice(minx - pad_deg, maxx + pad_deg))
    assert out.sizes["latitude"] > 0 and out.sizes["longitude"] > 0, (
        "empty subset — latitude ordering or longitude convention (0–360 vs −180–180) is wrong")
    return out


def field_daily(ds: xr.Dataset, fields: gpd.GeoDataFrame, var: str = "precip") -> pd.DataFrame:
    """Area-weighted daily mean of `var` per field. One row per field-day."""
    da = to_mm_per_day(ds[var])
    frames = []
    for day in da["time"].values:
        layer = da.sel(time=day).rio.write_crs("EPSG:4326")
        stats = exact_extract(layer, fields, ["mean"], output="pandas")
        frames.append(pd.DataFrame({
            "field_id": fields["field_id"].to_numpy(),
            "date": pd.Timestamp(day).normalize(),
            "precip_mm": stats["mean"].to_numpy(),
        }))
    out = pd.concat(frames, ignore_index=True)

    missing = out["precip_mm"].isna()
    assert not missing.any(), (
        f"{int(missing.sum())} field-day(s) with no value — "
        f"fields {sorted(out.loc[missing, 'field_id'].unique())[:5]} fall outside the grid")
    assert out["precip_mm"].between(0, 300).all(), (
        f"daily rainfall outside 0–300 mm (max {out['precip_mm'].max():.1f}) — unit conversion")
    return out

Inline verification — the gauge comparison, which is the only external evidence available:

PYTHON 
season = field_daily(ds, fields).query("field_id == @gauge_field_id")
gridded_total = float(season["precip_mm"].sum())
gauge_total = float(gauge_records["precip_mm"].sum())

print(f"gridded {gridded_total:.0f} mm vs gauge {gauge_total:.0f} mm "
      f"({(gridded_total - gauge_total) / gauge_total:+.1%})")
assert abs(gridded_total - gauge_total) < max(25.0, 0.15 * gauge_total), (
    "gridded and gauge season totals disagree by more than 15% — check units, "
    "the accumulation window, and whether the grid blends forecast with analysis")

Settings worth being deliberate about

Setting Default here Why it matters
Aggregation exact area weighting Makes a field’s value a smooth function of its geometry. Centroid sampling jumps whenever a boundary edit moves the centroid across a cell line
Unit handling read the attribute, convert once Flux in kg m⁻² s⁻¹ becomes millimetres per day by ×86,400. Guessing here produces a season total near zero and a reported drought
Bounding-box pad 0.2° Enough that a field on the subset edge still has whole cells around it; larger pads read data no field touches
Plausible daily range 0–300 mm A cheap assertion that catches unit errors, sign errors and a corrupted variable in one line
Missing-value policy fail, never fill A NaN field-day propagates into a null season total three steps later, at which point its cause is unrecoverable
Gauge tolerance 15% or 25 mm Gridded and gauge season totals rarely agree exactly; a persistent divergence beyond this is a unit or window problem, not weather

Temperature and rainfall deserve different confidence. The spatial gradient in daily temperature across a farm is far smaller than any grid’s own error, so gridded temperature is trustworthy; convective rainfall is not resolvable at 4 km and should be treated as an estimate to be checked.

Gotchas and edge cases

  • Longitude conventions differ as much as latitude ordering. Products on a 0–360 grid return nothing for a farm at −93° until the longitudes are shifted. Detect it by checking whether the grid’s longitude maximum exceeds 180, and convert rather than assuming.
Where each rainfall source can and cannot be trusted Two panels contrasting a four-kilometre gridded rainfall product with an on-farm gauge under convective rain. The grid averages a storm across the cell so both halves of the farm read the same, which is adequate for seasonal accounting; the gauge measures its own point exactly and says nothing about a field several kilometres away. Gridded rainfall, 4 km cell A summer storm crosses half the farm One cell averages the wet half and the dry half. Both fields are reported at 15 mm. The seasonal total stays roughly right. Good for accounting, wrong for one event. An on-farm gauge Measures one point exactly Reads 30 mm; the far field had nothing. Says nothing about the field eight km away. Drift from the grid usually means a blocked gauge. Good for one event, and for validation.

  • Convective rainfall is not resolvable at any usable grid resolution. A summer thunderstorm drops 30 mm on one half of a farm and nothing on the other; a 4 km cell reports the average of both. Use gridded rainfall for seasonal accounting and a gauge for anything where a single event matters.

  • Accumulation windows are not always midnight to midnight. Some products accumulate to the timestamp, others from it, and some are UTC on a farm eight hours away. A one-day systematic offset between rainfall and the operations it should explain is usually this.

  • Forecast blended with analysis changes the past. Products that splice a forecast onto the end of the analysis produce covariates that change retrospectively as the forecast is replaced. Read the product’s status flag and keep forecast rows in a separate column.

  • exactextract needs a CRS on the raster. An xarray layer without rio.write_crs silently mismatches the polygons’ CRS and returns nothing or nonsense. Set it explicitly per layer, as above.

  • Aggregating one day at a time is simple and slow. For a decade of history, aggregate the whole time axis in one pass with a chunked read instead — the per-call overhead dominates at 3,650 iterations.

  • The mean is the right statistic for rainfall and the wrong one for several other variables. Rainfall over a field is an area-weighted average by definition. Wind gusts, frost risk and heat stress are threshold questions where the field’s minimum or maximum cell matters more than its mean — a single frost-prone hollow decides whether a crop was damaged. Compute the statistic the decision needs rather than defaulting to the mean for every variable in the product.

This guide is part of Weather & Agronomic Data Integration for Field Models — see there for the water balance and the covariate table this feeds.

Parent section Weather & Agronomic Data Integration for Field Models Related Computing Growing Degree Days from Gridded Weather Data