21 min read 4701 words

Zonal Statistics of GHI over Candidate Parcels with rasterstats

The scenario: a portfolio of 3,200 parcels is ranked by mean plane-of-array irradiance, and 400 of them come back with identical values to four decimal places. They are all smaller than one raster cell, so each took the value of the single cell its centroid fell in, and the ranking among them is an artefact of the grid rather than of the resource. This page computes zonal statistics that say so, and it extends solar irradiance raster processing.

Root-cause analysis

Three defaults produce a zonal statistic that looks precise and is not.

  1. Centroid-based sampling on parcels smaller than a cell. With all_touched=False and a parcel below the cell size, the statistic is one cell’s value. That is not wrong so much as unquantified — the parcel’s true mean could differ by whatever the local gradient is.
  2. Unweighted means over partially covered cells. A cell half inside the parcel contributes as much as one wholly inside it. On a compact parcel spanning many cells the bias is negligible; on a long, thin parcel it is not, and long thin parcels are common along transmission corridors.
  3. Nodata treated as zero. An undeclared fill value of −9999 dragged into a mean produces a number that is obviously wrong; a fill of 0 produces one that is plausibly wrong, which is worse.
Centroid, all-touched and coverage-weighted sampling Three copies of the same irregular parcel drawn over a raster grid. In the first, centroid sampling, only the single cell containing the parcel centroid is highlighted. In the second, all-touched, every cell the parcel intersects is fully highlighted including cells only clipped at a corner. In the third, coverage weighting, each intersecting cell is shaded in proportion to the fraction of it the parcel covers, from nearly white at the edges to solid in the interior. Each panel reports the resulting mean: 1,842, 1,829 and 1,836 kilowatt-hours per square metre per year. Same parcel, same grid, three sampling rules centroid 1 842 kWh/m²·yr one cell — no variation all-touched 1 829 kWh/m²·yr edge cells count fully coverage-weighted 1 836 kWh/m²·yr each cell in proportion The three differ by less than a percent here and by several percent on a long, thin parcel — which is exactly the shape a transmission-corridor site takes.

Pre-flight validation

The decisive number is how many cells each parcel actually covers. Below about ten, the statistic is dominated by the grid and should be labelled as such.

python
import geopandas as gpd
import rasterio


def parcel_cell_coverage(parcels: gpd.GeoDataFrame, raster_path: str) -> gpd.GeoDataFrame:
    """How many raster cells each parcel spans — the honesty check for a zonal mean."""
    with rasterio.open(raster_path) as src:
        cell_area = abs(src.transform.a * src.transform.e)
        crs = src.crs
    p = parcels.to_crs(crs)
    out = p.copy()
    out["cells_covered"] = p.geometry.area / cell_area
    out["statistic_quality"] = out["cells_covered"].map(
        lambda n: "grid-dominated" if n < 10 else ("usable" if n < 100 else "well-resolved")
    )
    return out[["cells_covered", "statistic_quality", "geometry"]]

Running this before the statistics turns “3,200 parcels ranked by irradiance” into “2,800 ranked meaningfully and 400 whose ranking is grid noise”, which is a different report.

Fix implementation

python
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.features import geometry_mask, rasterize


def zonal_ghi(
    parcels: gpd.GeoDataFrame,
    raster_path: str,
    *,
    id_field: str = "parcel_id",
    supersample: int = 4,
) -> gpd.pd.DataFrame:
    """Area-weighted zonal mean with a coverage fraction, computed per parcel window."""
    rows = []
    with rasterio.open(raster_path) as src:
        p = parcels.to_crs(src.crs)
        for _, parcel in p.iterrows():
            window = rasterio.windows.from_bounds(*parcel.geometry.bounds, transform=src.transform)
            window = window.round_offsets().round_lengths()
            if window.width < 1 or window.height < 1:
                window = rasterio.windows.Window(
                    int(window.col_off), int(window.row_off), max(1, int(window.width)),
                    max(1, int(window.height)),
                )
            data = src.read(1, window=window, masked=True).astype("float32")
            transform = src.window_transform(window)

            # Supersampled coverage: what fraction of each cell the parcel actually covers.
            fine = rasterize(
                [(parcel.geometry, 1)],
                out_shape=(data.shape[0] * supersample, data.shape[1] * supersample),
                transform=transform * rasterio.Affine.scale(1 / supersample, 1 / supersample),
                fill=0, dtype="uint8", all_touched=True,
            )
            weights = fine.reshape(
                data.shape[0], supersample, data.shape[1], supersample
            ).mean(axis=(1, 3)).astype("float32")

            valid = weights * (~data.mask).astype("float32")
            total_weight = float(valid.sum())
            if total_weight <= 0:
                rows.append({id_field: parcel[id_field], "mean_ghi": np.nan,
                             "coverage_fraction": 0.0, "cells": 0})
                continue

            mean = float((data.filled(0) * valid).sum() / total_weight)
            rows.append({
                id_field: parcel[id_field],
                "mean_ghi": mean,
                "min_ghi": float(data.min()) if data.count() else np.nan,
                "max_ghi": float(data.max()) if data.count() else np.nan,
                "coverage_fraction": total_weight / float(weights.sum()) if weights.sum() else 0.0,
                "cells": int((weights > 0).sum()),
            })
    return gpd.pd.DataFrame(rows)

The supersampled weight array is what turns “all-touched or not” — a binary choice that is wrong in one direction or the other — into a continuous coverage fraction. At a supersample of four, a cell half inside the parcel contributes about half, which is the answer the question actually wants.

Sampling error against the number of cells a parcel spans A chart with the number of raster cells a parcel spans on a logarithmic horizontal axis from one to one thousand, and the difference between centroid and coverage-weighted means on the vertical. The curve falls steeply: 12 percent at one cell, 2.4 percent at ten, 0.6 percent at forty and under 0.1 percent at four hundred. A shaded region below ten cells is labelled grid-dominated, and a note records that 400 of 3,200 parcels in the worked portfolio fall in it. Below about ten cells, the grid decides the ranking grid-dominated 0% 5% 10% 1 10 100 1000 cells spanned → 12.0% 2.4% 0.25% 400 of 3 200 parcels in the worked portfolio span under ten cells. Their mutual ranking is an artefact of the grid, and the cell count is what says so.

Fallback routing and performance tuning

  • Window per parcel, never read the whole raster. A national GHI grid read in full for each of 3,200 parcels is the usual reason a zonal run takes hours; windowed reads make it seconds.
  • Sort parcels by tile before iterating. Reading windows in spatial order keeps the GDAL block cache warm and can halve wall-clock on a tiled source.
  • Use rasterstats for the simple case. Its zonal_stats with all_touched=True is fine when parcels span many cells; the supersampled weighting above earns its complexity on small or thin parcels.
  • Keep the supersample modest. Four is enough for a coverage fraction; sixteen costs sixteen times the rasterisation for a difference below the raster’s own uncertainty.
  • Batch by raster, not by parcel. For a multi-band stack, read the window once and reduce every band from it rather than reopening per band.

Downstream validation

python
import numpy as np


def assert_zonal_sane(stats, *, min_coverage: float = 0.9) -> None:
    """A zonal statistic must be inside the source range and adequately covered."""
    finite = stats.dropna(subset=["mean_ghi"])
    assert not finite.empty, "every parcel returned NaN — check the CRS and the extents overlap"
    assert (finite["mean_ghi"] >= finite["min_ghi"] - 1e-6).all(), "mean below the observed minimum"
    assert (finite["mean_ghi"] <= finite["max_ghi"] + 1e-6).all(), "mean above the observed maximum"
    poor = finite[finite["coverage_fraction"] < min_coverage]
    assert poor.empty or len(poor) / len(finite) < 0.05, (
        f"{len(poor)} parcels below {min_coverage:.0%} coverage — nodata or an extent mismatch"
    )
    tiny = finite[finite["cells"] < 10]
    if len(tiny):
        print(f"note: {len(tiny)} parcels span under 10 cells — their ranking is grid-dominated")

Reading a zonal result honestly

Three columns turn a zonal table from a number into a claim that can be checked.

Four columns every zonal table should carry A worked zonal table for four parcels. Parcel A: mean 1,836 kilowatt-hours per square metre per year, coverage 1.00, 412 cells, range 1,801 to 1,874 — well resolved. Parcel B: mean 1,842, coverage 0.62, 88 cells, range 1,812 to 1,871 — usable but partly nodata. Parcel C: mean 1,829, coverage 1.00, 6 cells, range 1,826 to 1,833 — grid-dominated. Parcel D: mean not available, coverage 0.00, zero cells — outside the raster extent. Each row carries a quality label derived from the coverage and cell count rather than from the mean. The mean alone does not say how much to trust it parcel mean coverage cells range quality A 1 836 1.00 412 1 801–1 874 well resolved B 1 842 0.62 88 1 812–1 871 partly nodata C 1 829 1.00 6 1 826–1 833 grid-dominated D — 0.00 0 — outside the extent Parcel C has the tightest range and the least information: six cells cannot resolve within-parcel variation, so its narrow range is the grid speaking rather than the resource.

Coverage fraction says how much of the parcel had valid data behind it. A parcel at 0.62 coverage has a mean over the 62 percent that was not nodata, which may be perfectly representative or may be systematically biased if the missing part is a lake or a cloud-masked region.

Cell count says whether the statistic is resolving anything. Ten cells is a coarse average; several hundred is a meaningful distribution, and only then do the minimum and maximum carry information about within-parcel variation.

Range — the minimum and maximum alongside the mean — is what shows whether the parcel is uniform. A 40-hectare parcel whose GHI ranges by 3 percent is a different siting proposition from one that ranges by 0.2 percent, and the mean alone hides that entirely.

Publishing the three alongside the mean costs nothing, because the same read produced them, and it answers the question a reviewer asks first: how much should I trust the fourth decimal place.

Frequently asked questions

Should all_touched be True or False?

Neither, for parcels near the cell size — that is the point of the coverage weighting. Where the simpler API is being used, all_touched=True is the safer default because it never returns an empty result for a small parcel, at the cost of including cells that barely overlap. all_touched=False silently returns nothing for a parcel that contains no cell centre.

Does the raster need to be in the same CRS as the parcels?

They need to be reconciled, and reprojecting the parcels is almost always the cheaper direction — thousands of geometries against billions of pixels. The exception is when several rasters are being combined, where a common grid matters more and the vectors follow it.

How should cloud-masked or seasonal nodata be handled?

As a coverage question rather than a value question. Compute the statistic over valid cells and report the fraction; substituting a fill value or an interpolated estimate hides the gap and biases the mean toward whatever the substitute was. A parcel with 40 percent coverage in one month deserves a flag, not an invented value.

Is a mean the right statistic for siting?

For a first-pass ranking, yes. For anything downstream, the percentiles matter more: a parcel whose tenth percentile is high is a better site than one with the same mean and a long low tail, because the layout will not use the whole parcel. Computing a small set of percentiles from the same window read costs nothing extra.

How do I make the run reproducible?

Record the raster path and its checksum, the supersample factor, the parcel layer vintage and the CRS the statistic was computed in. Two zonal tables computed with different supersampling or different all-touched settings are not comparable, and nothing in the numbers themselves says which was used.

Can the same code summarise a whole hourly stack?

Yes, and it is the efficient shape: read the parcel window once across every band, compute the weights once, and reduce each band against them. That turns 8,760 separate zonal runs into one window read per parcel, which is the difference between a minute and most of a day.

Can zonal statistics be computed on a cloud-hosted raster without downloading it?

Yes, and it is the normal case for national products. A Cloud-Optimised GeoTIFF served over HTTP supports windowed reads, so each parcel fetches only the tiles it overlaps — typically kilobytes. What breaks it is a striped, uncompressed GeoTIFF, where every window read pulls whole rows and the transfer dwarfs the computation. Checking that the source is tiled before running a portfolio is a one-line check that saves hours.

What happens when a parcel spans two raster tiles?

Nothing special, provided the read is done through the dataset rather than per file: rasterio resolves the window across internal tiles transparently. It matters when the source is a set of separate files rather than one mosaic, in which case the parcel needs a VRT or a merged source, or the statistic is silently computed over whichever tile the code happened to open.

Related articles

Back to Solar Irradiance Raster Processing