Handling Partial Pixel Coverage in Zonal Statistics

Fix boundary bias in zonal statistics: compare all_touched, centre-in-polygon and exact area weighting, and choose the rule that matches your polygon size and pixel resolution.

Pixels do not respect polygon boundaries. Every zonal statistic therefore rests on a rule for what to do with the pixels that straddle the edge: include them whole (all_touched=True), include them only when their centre falls inside (the default), or weight them by the fraction of their area that lies inside (exact extraction). The three give measurably different answers, and which one is correct depends on how large your polygons are relative to your pixels.

This page compares the rules and shows how to pick. For the full zonal workflow, see Zonal Statistics for Polygon Aggregation.

Three answers to “does this edge pixel count?” The same polygon is overlaid on the same pixel grid three times. Under the centre-in-polygon rule only fully interior pixels are shaded. Under the all-touched rule every intersected pixel is shaded fully. Under exact weighting each edge pixel carries a fractional weight between zero and one that matches how much of it lies inside the polygon. Same polygon, same raster, three different means centre-in-polygon 4 pixels · mean 0.61 under-counts the edge all_touched 16 pixels · mean 0.54 pulls in outside values exact area weights 0.1 0.4 0.4 0.1 0.5 1.0 1.0 0.5 0.5 1.0 1.0 0.5 0.1 0.4 0.4 0.1 weighted · mean 0.58 matches the true area

Why This Fails in Geospatial ML Pipelines

The bias is systematic, not random, and that is what makes it dangerous. Under all_touched, every polygon acquires a ring of pixels from outside itself, so every zonal mean is pulled towards the surrounding landscape. Small polygons have proportionally more ring, so the contamination scales inversely with polygon size — which means the error correlates with a variable (parcel size) that is often itself predictive. The model then learns a relationship between size and target that is an artefact of the extraction rule.

Under the default centre rule the bias runs the other way: edge pixels are dropped, so a polygon’s statistic reflects only its core. For compact fields that is nearly harmless. For elongated features — a riparian strip, a hedgerow, a road verge — it can discard most of the polygon, and for anything narrower than a pixel it returns None, which most pipelines then impute with a mean and never mention again.

Third, the rule is rarely recorded. Two teams computing “mean NDVI per field” with different all_touched settings produce feature tables that differ by a few percent per row, with the difference concentrated in small fields. When one team’s model is evaluated on the other’s features, performance drops for reasons nobody can locate — the same class of silent contract mismatch as re-fitting a scaler, discussed in scaling features consistently between training and inference.

Core Principles

  • Measure the exposure first. Perimeter-to-area in pixel units tells you whether the rule matters at all.
  • Prefer exact area weighting when it does. exactextract is both more accurate and faster at scale.
  • Never let all_touched be the silent default for area classes. It inflates and contaminates.
  • Emit the pixel count as a feature. A mean over 6 pixels is not a mean over 6,000.
  • Record the rule in the output metadata. It is part of the feature’s definition.
  • Handle sub-pixel polygons explicitly rather than imputing their None away.

Production-Ready Code

from __future__ import annotations

import logging
import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio

logger = logging.getLogger(__name__)


def boundary_exposure(gdf: gpd.GeoDataFrame, cell_size: float) -> pd.Series:
    """Fraction of a polygon's area that sits within one pixel of its boundary.

    A useful proxy for how much the coverage rule can change the answer:
    below ~2% any rule is fine; above ~20% use exact weighting.
    """
    area = gdf.geometry.area
    inner = gdf.geometry.buffer(-cell_size).area.clip(lower=0)
    return ((area - inner) / area.replace(0, np.nan)).rename("boundary_frac")


def zonal_exact(gdf: gpd.GeoDataFrame, raster_path: str,
                stats: tuple = ("mean", "stdev", "count")) -> pd.DataFrame:
    """Area-weighted zonal statistics via exactextract.

    Each pixel contributes in proportion to the fraction of it inside the polygon,
    so boundary pixels are neither dropped nor counted whole.
    """
    from exactextract import exact_extract

    with rasterio.open(raster_path) as src:
        if gdf.crs != src.crs:
            logger.info("reprojecting polygons %s -> %s",
                        gdf.crs.to_string(), src.crs.to_string())
            gdf = gdf.to_crs(src.crs)

    result = exact_extract(raster_path, gdf, list(stats), output="pandas")
    result.index = gdf.index
    return result


def zonal_with_rule(gdf: gpd.GeoDataFrame, raster_path: str,
                    all_touched: bool) -> pd.DataFrame:
    """rasterstats zonal means under an explicit coverage rule, for comparison."""
    from rasterstats import zonal_stats

    with rasterio.open(raster_path) as src:
        nodata = src.nodata
        if gdf.crs != src.crs:
            gdf = gdf.to_crs(src.crs)

    rows = zonal_stats(gdf, raster_path, stats=["mean", "count"],
                       nodata=nodata, all_touched=all_touched)
    out = pd.DataFrame(rows, index=gdf.index)
    return out.rename(columns={"mean": f"mean_at{int(all_touched)}",
                               "count": f"count_at{int(all_touched)}"})


def compare_rules(gdf: gpd.GeoDataFrame, raster_path: str,
                  cell_size: float) -> pd.DataFrame:
    """Side-by-side of the three rules plus the exposure metric."""
    exposure = boundary_exposure(gdf, cell_size)
    strict = zonal_with_rule(gdf, raster_path, all_touched=False)
    loose = zonal_with_rule(gdf, raster_path, all_touched=True)
    exact = zonal_exact(gdf, raster_path).rename(columns={"mean": "mean_exact"})

    out = pd.concat([exposure, strict, loose, exact[["mean_exact"]]], axis=1)
    out["spread"] = out[["mean_at0", "mean_at1", "mean_exact"]].max(axis=1) - \
                    out[["mean_at0", "mean_at1", "mean_exact"]].min(axis=1)
    return out

Step-by-Step Walkthrough

Step 1 — Measure how exposed your polygons are

import geopandas as gpd
import rasterio

fields = gpd.read_file("fields.gpkg")
with rasterio.open("ndvi.tif") as src:
    cell = abs(src.transform.a)
    fields = fields.to_crs(src.crs)

exposure = boundary_exposure(fields, cell)
print(exposure.describe())
print(f"{(exposure > 0.2).mean():.1%} of polygons are more than 20% boundary")

If the 75th percentile is under about 0.05, the coverage rule barely matters for this dataset and you can use whichever tool is convenient. If a meaningful share exceed 0.2, exact weighting is not an optimisation — it is a correctness requirement.

Step 2 — Compare the three rules on a sample

sample = fields.sample(min(2_000, len(fields)), random_state=0)
cmp = compare_rules(sample, "ndvi.tif", cell)

print(cmp[["boundary_frac", "mean_at0", "mean_at1", "mean_exact", "spread"]].describe())
print("\nlargest disagreements:")
print(cmp.nlargest(5, "spread")[["boundary_frac", "mean_at0", "mean_at1", "mean_exact"]])

The correlation between boundary_frac and spread is the diagnostic: if it is strong, the rule is injecting a size-dependent bias into the features.

Step 3 — Compute the production statistics with exact weights

stats = zonal_exact(fields, "ndvi.tif", stats=("mean", "stdev", "count", "coverage_fraction"))
fields = fields.join(stats.add_prefix("ndvi_"))

fields["ndvi_pixels"] = fields["ndvi_count"]
fields["ndvi_thin"] = fields["ndvi_count"] < 10        # thin-evidence flag
print(f"{fields['ndvi_thin'].mean():.1%} of polygons rest on fewer than 10 pixels")

Step 4 — Deal with sub-pixel polygons deliberately

tiny = fields.loc[fields["ndvi_count"].isna() | (fields["ndvi_count"] < 1)]
if len(tiny):
    logger.warning("%d polygons are smaller than one pixel", len(tiny))
    centroid_vals = point_sample("ndvi.tif", tiny.geometry.centroid)
    fields.loc[tiny.index, "ndvi_mean"] = centroid_vals
    fields.loc[tiny.index, "ndvi_thin"] = True

Sampling the centroid is an honest answer for a sub-pixel polygon — it is literally the only information the raster has — provided the ndvi_thin flag travels with it so the model can discount it.

Step 5 — Record the rule with the features

fields.attrs["zonal_rule"] = "exactextract area-weighted"
fields.attrs["raster"] = "ndvi.tif"
fields.attrs["cell_size_m"] = cell
fields.to_parquet("field_features.parquet")
The rule matters exactly where polygons are small A scatter plot places boundary exposure on the horizontal axis and the spread between the three coverage rules on the vertical axis. Points cluster near zero spread at low exposure and fan upwards steeply as exposure rises past twenty percent. Disagreement between rules vs boundary exposure fraction of area within one pixel of the edge → spread ~20% exposure any rule is fine here exact weighting required Because exposure tracks polygon size, the bias correlates with a real predictor.

Verification

import numpy as np
from shapely.geometry import box


def test_exact_weights_match_the_geometric_fraction():
    """A polygon covering exactly half a pixel must give that pixel weight 0.5."""
    # 10 m pixels; polygon spans one full pixel and half of the next.
    poly = box(0, 0, 15, 10)
    weights = pixel_coverage(poly, transform=from_origin(0, 10, 10, 10), shape=(1, 2))
    assert np.allclose(weights, [[1.0, 0.5]], atol=1e-6)


def test_all_touched_never_reduces_the_pixel_count():
    strict = zonal_with_rule(fields, "ndvi.tif", all_touched=False)
    loose = zonal_with_rule(fields, "ndvi.tif", all_touched=True)
    assert (loose["count_at1"] >= strict["count_at0"]).all()


def test_exact_mean_lies_between_the_two_rules():
    """Area weighting should not fall outside the bracket the two rules define."""
    cmp = compare_rules(fields.head(200), "ndvi.tif", cell_size=10.0)
    lo = cmp[["mean_at0", "mean_at1"]].min(axis=1) - 1e-6
    hi = cmp[["mean_at0", "mean_at1"]].max(axis=1) + 1e-6
    ok = cmp["mean_exact"].between(lo, hi)
    assert ok.mean() > 0.95, "exact means fall outside the bracket too often"


def test_subpixel_polygon_is_flagged_not_silently_dropped():
    tiny = fields_with_a_2m_polygon()
    out = zonal_exact(tiny, "ndvi_10m.tif")
    assert out["count"].iloc[0] < 1

A useful field check: compute the same feature under both rules, fit the model twice, and compare held-out scores under spatial cross-validation. If the scores differ materially, the extraction rule is a hyperparameter and deserves to be recorded as one.

Choosing the coverage rule from the geometry A decision flow begins with the boundary exposure measurement. Low exposure leads to the default centre-in-polygon rule. Moderate to high exposure leads to exact area weighting. Polygons narrower than one pixel lead to a recommendation to use a finer raster or a different analysis unit. measure boundary exposure area within one pixel of the edge below ~5% → default rule is fine 5% to 40% → exact area weighting narrower than a pixel → finer raster or new unit Whichever branch you take, emit the pixel count with the statistic. It is the only way a model can tell a solid estimate from a thin one.

FAQ

When does partial pixel coverage actually matter?

When polygons are small relative to pixels. A 100 ha field on a 10 m raster has a boundary effect well under 1%; a 0.4 ha allotment is about 40 pixels, half of which touch the edge. The boundary_exposure metric above turns that intuition into a number you can threshold.

Is exactextract always better than rasterstats?

It is more accurate for small or elongated polygons, and faster on large jobs. For large compact polygons the two agree to a fraction of a percent, so the simpler dependency is perfectly reasonable — see computing zonal statistics with rasterstats.

How should I handle a polygon smaller than one pixel?

Accept that the statistic is one pixel’s value, record the pixel count so the model knows the evidence is thin, and ask whether a finer raster or a different analysis unit is the real fix.

Part of: Zonal Statistics for Polygon Aggregation Part of: Spatial Feature Engineering for Machine Learning