Computing Weighted Zonal Statistics for Partial Pixels

Weight each pixel by the fraction of it the polygon covers, rather than counting it wholly in or wholly out:

import numpy as np

# coverage: fraction 0..1 of each pixel inside the polygon; values: the raster window
w = np.where(np.isfinite(values), coverage, 0.0)
weighted_mean = float((values * w).sum() / w.sum())
weighted_sum = float(np.nansum(values * coverage))       # e.g. total biomass

For small or thin zones the difference from the pixel-centre rule is not a rounding error — it can be a tenth of the answer. This page belongs to zonal statistics and vector–raster integration in Core Raster Fundamentals & STAC Mapping.

Why the Centre Rule Biases Small Zones

Boundary pixels are most of a small zone A small field polygon overlaps a five by four pixel neighbourhood. Only four interior pixels are fully covered. Twelve boundary pixels are covered by between ten and eighty per cent. Under the centre rule, some boundary pixels count fully and others not at all, so the zone's area and total are misstated; weighting by coverage fraction counts each boundary pixel in proportion. Coverage fraction per pixel 0.55 0.90 0.60 0.85 1.00 0.95 0.45 0.80 0.50 Centre rule: 7 pixels counted. Weighted: 6.6 pixel-equivalents. The difference grows as zones shrink.

The error is not random: it depends on how the zone happens to sit on the grid, so two identical fields a few metres apart can report different areas. For national statistics the errors average out across millions of zones. For a per-field product — yield per hectare, carbon per plot — they do not, and each field inherits its own misstatement.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Windowed reads and supersampled rasterization
numpy >=1.23 Coverage aggregation and weighted reductions
geopandas >=0.14 Zones and reprojection
shapely >=2.0 Geometry bounds and validity 
pip install "rasterio>=1.3.0" "numpy>=1.23" "geopandas>=0.14" "shapely>=2.0"

Complete Working Example

import numpy as np
import rasterio
from rasterio.features import rasterize
from rasterio.transform import Affine
from rasterio.windows import from_bounds


def coverage_fraction(geom, window_transform: Affine, shape: tuple[int, int],
                      factor: int = 10) -> np.ndarray:
    """Fraction of each pixel covered by geom, via supersampled rasterization."""
    fine_transform = window_transform * Affine.scale(1 / factor)
    fine = rasterize([(geom, 1)], out_shape=(shape[0] * factor, shape[1] * factor),
                     transform=fine_transform, fill=0, dtype="uint8")
    return fine.reshape(shape[0], factor, shape[1], factor).mean(axis=(1, 3))


def weighted_zonal(raster_path: str, geom, *, band: int = 1) -> dict[str, float]:
    with rasterio.open(raster_path) as src:
        win = from_bounds(*geom.bounds, transform=src.transform).round_offsets().round_lengths()
        win = win.intersection(rasterio.windows.Window(0, 0, src.width, src.height))
        values = src.read(band, window=win, masked=True).astype("float64")
        wt = src.window_transform(win)
        px_area = abs(src.transform.a * src.transform.e)

    cov = coverage_fraction(geom, wt, values.shape)
    valid = ~np.ma.getmaskarray(values)
    w = np.where(valid, cov, 0.0)
    v = values.filled(0.0)

    covered = float(w.sum())
    return {
        "weighted_mean": float((v * w).sum() / covered) if covered else float("nan"),
        "weighted_sum": float((v * w).sum()),
        "covered_area_m2": covered * px_area,
        "centre_rule_pixels": int((cov >= 0.5).sum()),     # approximate centre rule
        "valid_fraction": float(w.sum() / max(cov.sum(), 1e-9)),
    }

valid_fraction is the diagnostic worth reporting alongside every weighted statistic: the share of the zone’s coverage that fell on valid pixels. A zone half-covered by cloud or nodata has a weighted mean that describes only its visible half, and that caveat belongs next to the number. The unweighted, library-based approach this refines is in computing zonal statistics with rasterstats.

How Much It Matters, by Zone Size

Centre-rule error shrinks with zone size For zones of about ten pixels the centre rule's area error is commonly fifteen to twenty per cent. At a hundred pixels it falls to around five per cent, at a thousand to under two, and beyond ten thousand pixels it is negligible. The weighted method's error stays near zero across the range. Typical area error against zone size 20% 0% centre rule coverage-weighted 10 px 100 1,000 10,000 Zone size in pixels (log scale). Smallholder fields at 10 m sit squarely in the left half.

The practical threshold is a few hundred pixels. Above it the centre rule is fine for most purposes and much faster; below it weighting is worth the cost. Smallholder agriculture at 10 m resolution — fields of a hectare or less, which is 100 pixels or fewer — sits firmly in the range where weighting changes the answer.

Performance at Scale

Supersampled rasterization is simple and exact enough, but it costs factor² times the pixels per zone, which becomes expensive across hundreds of thousands of polygons. For that workload, dedicated libraries compute exact geometric coverage fractions directly from the polygon edges without supersampling, and they are an order of magnitude faster; the approach is covered in zonal statistics over millions of polygons with exactextract.

The implementation above remains useful as a reference: it is transparent, depends only on rasterio, and gives a way to check that a faster tool is producing the numbers you expect on a sample of zones before trusting it with the full set.

Verification

Coverage must add up to the polygon's area The sum of coverage fractions multiplied by pixel area should equal the polygon's own geometric area to within about one per cent at a supersampling factor of ten. A larger discrepancy means the polygon and raster are in different CRSs or the window missed part of the polygon. An area identity that should hold Σ coverage × pixel area from the raster side polygon.area from the vector side ≈ Within 1% at factor 10. A larger gap is a CRS or windowing bug. 
stats = weighted_zonal("biomass.tif", field_geom)
rel = abs(stats["covered_area_m2"] - field_geom.area) / field_geom.area
assert rel < 0.01, f"coverage area off by {rel:.1%} — CRS or window mismatch"

The area identity is a strong check because it involves no raster values at all: it tests the geometry pipeline — CRS, window and supersampling — in isolation. If it holds, the weights are right, and any remaining question is about the values being weighted.

Common Errors

Weighted area is far from the polygon area

The polygon is in a different CRS from the raster. Reproject the zones to src.crs before computing.

The weighted mean is NaN

Every covered pixel is nodata, so the total weight is zero. Report the zone as unobserved rather than zero.

Results differ slightly from a library tool

Supersampling approximates coverage to about one per cent at factor ten. Raise the factor, or treat the difference as within tolerance.

Computation is slow for many zones

Supersampling costs factor-squared per zone. Use an exact-coverage library for bulk work.

Zones near the raster edge report too little area

The window was clipped to the raster, so coverage outside it is missing. That is correct for the statistic but should be flagged; compare against the polygon’s full area and report the observed share.

Frequently Asked Questions

Q: When does partial-pixel weighting matter? When zones are small relative to pixels or elongated, so that boundary pixels are a large share of the total. For a field of a few dozen pixels the pixel-centre rule can misstate area by ten per cent or more; for a province of millions of pixels it makes no practical difference.

Q: Does weighting change the mean or only the sum? Both, but it matters far more for sums and totals — area, biomass, population — because every boundary pixel is either wholly counted or wholly ignored under the centre rule. Means shift less, and only when boundary pixels differ systematically from interior ones.

Q: How precise do coverage fractions need to be? A supersampling factor of ten per axis gives fractions to the nearest one per cent, which is well within the uncertainty of the vector boundary itself. Exact geometric coverage is available from dedicated libraries when many zones must be processed.

Q: Should the weights also account for pixel area varying with latitude? Only in a geographic CRS, where pixel area shrinks toward the poles. In a projected equal-area or UTM grid, pixel area is effectively constant within a zone and the coverage fraction alone is the right weight.