Zonal Statistics over Millions of Polygons with exactextract

For very large polygon sets, compute exact coverage-weighted statistics directly and stream the results:

from exactextract import exact_extract

df = exact_extract(
    "ndvi_2026.tif", "fields.gpkg",
    ["mean", "count", "stdev", "frac"],        # coverage-weighted by default
    include_cols=["field_id"],
    output="pandas",
)
df.to_parquet("field_ndvi.parquet")

Coverage-weighted, exact, and fast enough for national parcel datasets on a single machine. This page belongs to zonal statistics and vector–raster integration in Core Raster Fundamentals & STAC Mapping.

Where the Time Goes in Zonal Statistics

At national scale the statistics are almost free and the overhead is everything. Five million parcels at a millisecond of overhead each is well over an hour before any value is summed. Removing per-polygon mask allocation is what makes the difference, and exact coverage — the weighting described in computing weighted zonal statistics for partial pixels — comes along for free.

Environment & Setup

pip install "exactextract>=0.2" "geopandas>=0.14" "pyarrow>=15.0" "rasterio>=1.3.0"

Complete Working Example

from concurrent.futures import ProcessPoolExecutor
from pathlib import Path

import geopandas as gpd
import pandas as pd
import rasterio
from exactextract import exact_extract

OPS = ["mean", "stdev", "count", "frac", "median"]


def partition(polygons_path: str, raster_path: str, tile_px: int = 4096) -> dict[tuple, gpd.GeoDataFrame]:
    """Group polygons by the raster tile containing their centroid."""
    gdf = gpd.read_file(polygons_path)
    with rasterio.open(raster_path) as src:
        if gdf.crs != src.crs:
            gdf = gdf.to_crs(src.crs)
        inv = ~src.transform
    cols, rows = zip(*[inv * (p.x, p.y) for p in gdf.geometry.centroid])
    gdf["_tile"] = [(int(r) // tile_px, int(c) // tile_px) for r, c in zip(rows, cols)]
    return {k: g.drop(columns="_tile") for k, g in gdf.groupby("_tile")}


def run_partition(args) -> str:
    key, part, raster_path, out_dir = args
    df = exact_extract(raster_path, part, OPS, include_cols=["field_id"], output="pandas")
    out = Path(out_dir) / f"part_{key[0]:03d}_{key[1]:03d}.parquet"
    df.to_parquet(out, index=False)
    return str(out)


def zonal_at_scale(polygons_path: str, raster_path: str, out_dir: str,
                   *, workers: int = 8) -> pd.DataFrame:
    Path(out_dir).mkdir(parents=True, exist_ok=True)
    parts = partition(polygons_path, raster_path)
    jobs = [(k, g, raster_path, out_dir) for k, g in parts.items()]
    with ProcessPoolExecutor(max_workers=workers) as pool:
        files = list(pool.map(run_partition, jobs))
    return pd.concat((pd.read_parquet(f) for f in files), ignore_index=True)

Partitioning by the raster’s tiles is what keeps the reads efficient. Each worker reads one compact region — a few internal tiles of a COG — rather than every worker seeking across the whole file, which on object storage turns scattered requests into cached, near-sequential ones. The parallelism model follows the scene-per-task pattern in scaling raster processing with Dask, expressed here with a plain process pool.

Choosing Operations Carefully

For categorical rasters — land cover per parcel, for instance — the useful operations are per-class fractions rather than means. Asking for the coverage-weighted fraction of each class gives a composition per polygon that sums to one, which is almost always what a downstream analysis wants and is far more informative than a single mode.

Memory and Output at Scale

Two things break naive runs at this size. Holding every result in memory until the end fails once the result table itself runs to gigabytes; writing one Parquet file per partition and concatenating lazily avoids that entirely, and leaves a restartable record if the run is interrupted. And loading every polygon into every worker multiplies memory by the worker count; passing only each partition’s polygons to its worker, as above, keeps each process’s footprint proportional to its own share.

A third, quieter issue is polygon validity. National parcel datasets routinely contain self-intersections and slivers, and a single invalid geometry can raise inside a worker and lose a whole partition. Running make_valid once before partitioning, and dropping empty geometries, is cheap insurance. The general validity handling is covered in rasterizing vector polygons onto a raster grid.

Verification

result = zonal_at_scale("fields.gpkg", "ndvi_2026.tif", "out/zonal")
n_in = len(gpd.read_file("fields.gpkg", columns=["field_id"]))
assert len(result) == n_in, f"{n_in - len(result)} polygons missing from the output"
assert result["field_id"].is_unique, "a polygon was counted in two partitions"

Row-count reconciliation is the check that catches a failed partition, which in a parallel run is easy to miss because the other partitions succeed and the concatenation looks complete. Follow it with a spot check of twenty polygons against the transparent reference method so the values themselves are confirmed, not only their number.

Common Errors

Some polygons are missing from the output

A partition failed, often on an invalid geometry. Validate geometries before partitioning and check row counts after.

Workers run out of memory

Each received the whole polygon set. Pass only the partition to each worker, never the full layer.

Every value is NaN for a region

The polygons are in a different CRS from the raster, or that region of the raster is entirely nodata. Reproject before partitioning, and check the count column to tell the two cases apart.

Median makes the run very slow

It collects and sorts every polygon’s pixels. Drop it unless a robust statistic is genuinely needed.

Results differ from rasterstats for small polygons

That is expected: rasterstats uses the pixel-centre rule, exactextract weights by coverage. The difference shrinks as polygons grow, and for small parcels the weighted value is the more faithful one.

Frequently Asked Questions

Q: Why is exactextract faster than rasterize-based approaches? Because it computes each pixel’s coverage fraction directly from the polygon edges rather than rasterizing a mask per polygon, and it walks the raster once per polygon’s bounding box. For many small polygons that removes most of the per-polygon overhead that dominates simpler tools.

Q: How should the polygons be partitioned? By the raster’s own tiling — assign each polygon to the tile containing its centroid. Each job then reads one compact region of the raster from disk or object storage, which keeps reads sequential and caches effective.

Q: What about polygons that cross tile boundaries? Assigning by centroid keeps each polygon in exactly one partition, and the job reads whatever part of the raster its polygons need, including across the tile edge. The partition only groups work; it does not clip geometry.

Q: Can this run against a remote COG? Yes, and partitioning by tile matters even more there, since it turns each worker’s reads into a compact set of range requests. Set the GDAL options described in diagnosing slow COG reads with GDAL VSI logging in each worker’s environment.

Related

• Zonal Statistics and Vector–Raster Integration — the parent topic.
• Computing Zonal Statistics with rasterstats — the simpler tool for smaller jobs.
• Computing Weighted Zonal Statistics for Partial Pixels — the weighting exactextract performs.
• Processing COGs on AWS Batch with Docker — running partitions as batch jobs.