Zonal Statistics with rasterstats vs exactextract

Two libraries answer the same question — what is the mean raster value inside this polygon? — and on most zones they return the same number, which is exactly why the cases where they disagree catch people out. This guide is for anyone who has to defend a per-zone statistic to someone else: it sits under Zonal Statistics & Raster Sampling in Spatial Analysis & Advanced Query Techniques, and works from a six-by-six grid small enough to check by hand before scaling the same comparison to a real scene.

Why This Approach / What Goes Wrong

The disagreement is architectural, not a bug in either library. rasterstats builds a boolean mask: it burns the polygon onto the raster's own grid with the same rasterization GDAL uses, producing a True/False array, and then hands the selected values to NumPy. Every selected pixel carries weight 1 and every rejected pixel carries weight 0. exactextract never rasterizes. Its C++ core clips the polygon against each cell rectangle analytically and returns the fraction of that cell's area lying inside — a float in [0, 1] — then performs a weighted reduction over those fractions.

For a zone spanning hundreds of cells the distinction is cosmetic. A zone of N cells has roughly √N boundary cells, so the share of the sample that is ambiguous falls away as zones grow; at a thousand cells the two means typically differ in the fourth decimal. The trouble starts when the boundary is most of the zone: a 30 m riparian buffer on a 10 m grid, a residential parcel on Landsat, a road corridor, an administrative sliver left over from a polygon overlay. There, half the sample sits on the edge.

Worse, the resulting error is not noise that averages out. It is biased whenever pixel values correlate with position relative to the boundary — and in practice they nearly always do, because zone boundaries follow rivers, roads, field edges and coastlines, all of which are exactly where the raster changes. Include a boundary pixel whole and you import whatever lies outside the zone; drop it and you throw away the part that was genuinely inside. Under the centroid rule that choice is made by a coin-flip of a few metres.

Prerequisites

• rasterstats>=0.19 — pure Python over rasterio and NumPy; zonal_stats, point_query, categorical=True, add_stats
• exactextract>=0.2 — compiled C++ core with Python bindings; exact_extract returning a DataFrame
• rasterio>=1.3,<2 — the reader both libraries sit on, covered in Raster Data Handling with Rasterio
• geopandas>=1.0 and shapely>=2.0 — the zone layer and its geometries
• pandas>=2.0, numpy>=1.26 — the result tables and the hand-check arithmetic

conda install -c conda-forge "rasterstats=0.19.*" "exactextract=0.2.*" \
  "rasterio=1.3.*" "geopandas=1.0.*" "shapely=2.0.*" "pandas=2.*" "numpy=1.26.*"

Take both from the same channel. exactextract ships binary wheels for mainstream CPython versions, but if pip starts invoking CMake you have fallen off the wheel matrix and are about to build GEOS by hand — switch to conda-forge rather than fighting it.

Step-by-Step Implementation

1. Build a grid you can verify with a pencil.

A six-by-six float32 raster at 10 m in UTM zone 32N, with two low-NDVI columns on the left standing in for water and bare bank, and vegetation on the right. Every number printed later can be recomputed from this array.

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import from_origin
from shapely.geometry import box

ndvi = np.array([
    [0.12, 0.14, 0.61, 0.66, 0.68, 0.70],
    [0.13, 0.15, 0.62, 0.67, 0.69, 0.71],
    [0.11, 0.16, 0.60, 0.65, 0.70, 0.72],
    [0.10, 0.18, 0.59, 0.64, 0.71, 0.73],
    [0.12, 0.20, 0.58, 0.63, 0.72, 0.74],
    [0.13, 0.22, 0.57, 0.62, 0.73, 0.75],
], dtype="float32")

transform = from_origin(400000.0, 5650000.0, 10.0, 10.0)   # 10 m cells, UTM 32N
profile = dict(driver="GTiff", height=6, width=6, count=1, dtype="float32",
               crs="EPSG:32632", transform=transform, nodata=-9999.0)

with rasterio.open("ndvi_demo.tif", "w", **profile) as dst:
    dst.write(ndvi, 1)

zones = gpd.GeoDataFrame(
    {"zone_id": ["field_block", "riparian"]},
    geometry=[
        box(400010, 5649960, 400050, 5650000),   # grid-aligned: 4 x 4 whole cells
        box(400012, 5649982, 400041, 5650000),   # 29 m wide, clips six cells partly
    ],
    crs="EPSG:32632",                            # metric, matches the raster exactly
)

The CRS choice is deliberate: UTM 32N gives square 10 m cells of constant ground area, so an unweighted mean is already area-correct. In a geographic CRS the cells would be trapezoids and in EPSG:3857 they would grow with latitude, which makes every "area-weighted" statistic a fiction.

2. Run rasterstats — an integer pixel count under the centroid rule.

import pandas as pd
from rasterstats import zonal_stats

rs = zonal_stats(
    zones,                       # consumed as GeoJSON: the CRS is stripped on the way in
    "ndvi_demo.tif",             # a path, so reads stay windowed per feature
    stats=["count", "mean"],
    nodata=-9999.0,              # explicit, never trust the file to carry a tag
    all_touched=False,           # the default: pixel centre must fall inside
)
rs = pd.DataFrame(rs, index=zones.index).add_prefix("rs_")
print(rs)
#    rs_count  rs_mean
# 0        16   0.5281
# 1         6   0.4750

3. Run exactextract — the same call, weighted by coverage.

from exactextract import exact_extract

ee = exact_extract(
    "ndvi_demo.tif",
    zones,                              # a GeoDataFrame; a path or file object also works
    ["count", "mean"],
    include_cols=["zone_id"],           # carried through so the join key survives
    output="pandas",
)
print(ee)
#      zone_id  count    mean
# 0  field_block  16.00  0.5281
# 1     riparian   5.22  0.5044

Two things differ before the numbers do. exact_extract returns a DataFrame rather than a list of dicts, and include_cols carries the identifier through, so there is no positional join to get wrong. And count comes back as 16.00, not 16 — it is the sum of coverage fractions, so a float is the honest type. Print ee.columns once on any new source: a single-band raster gives plain operation names, while multi-band or explicitly named sources prefix them with the band name.

4. Put the two answers side by side.

comparison = zones[["zone_id"]].join(rs).merge(ee, on="zone_id", validate="one_to_one")
comparison["delta"] = comparison["mean"] - comparison["rs_mean"]
print(comparison[["zone_id", "rs_count", "rs_mean", "count", "mean", "delta"]])
#       zone_id  rs_count  rs_mean  count    mean   delta
# 0  field_block        16   0.5281  16.00  0.5281  0.0000
# 1     riparian         6   0.4750   5.22  0.5044  0.0294

On the grid-aligned block the two agree to the last printed digit, because every cell is either fully in or fully out and a coverage weight of exactly 1.0 reduces to a boolean mask. On the riparian strip they differ by 0.029 NDVI — about 6% — and the sign is not random: the centroid rule kept a low-value water-edge pixel at full weight while discarding a high-value vegetation pixel that was 10% inside.

5. When a number looks wrong, ask for the weights.

exactextract will hand back the raw per-cell arrays instead of a summary, which turns an argument about a suspicious mean into a table you can add up.

import numpy as np

cells = exact_extract(
    "ndvi_demo.tif",
    zones[zones["zone_id"] == "riparian"],
    ["coverage", "values"],          # array ops: one array per feature, not a scalar
    output="pandas",
)
row = cells.iloc[0]
print(list(zip(np.round(row["values"], 2), np.round(row["coverage"], 2))))
# [(0.14, 0.8), (0.61, 1.0), (0.66, 1.0), (0.68, 0.1),
#  (0.15, 0.64), (0.62, 0.8), (0.67, 0.8), (0.69, 0.08)]
print(float(np.sum(row["coverage"])))   # 5.22 — the count from step 3

The rasterstats equivalent is raster_out=True, which attaches the masked mini-array and its affine transform to every result dict — heavier, but it shows you the same pixels through the boolean lens.

6. Categorical rasters: a ragged dict versus fraction columns.

For land cover the two libraries diverge in output shape, not just in weighting. rasterstats returns one dictionary per zone whose keys are only the classes actually present, so the result is ragged and has to be normalised before it becomes a table.

classes = np.array([
    [80, 80, 40, 40, 10, 10],
    [80, 80, 40, 40, 10, 10],
    [80, 40, 40, 10, 10, 10],
    [80, 40, 40, 10, 10, 10],
    [40, 40, 40, 10, 10, 10],
    [40, 40, 10, 10, 10, 10],
], dtype="uint8")

with rasterio.open("landcover_demo.tif", "w", **{**profile, "dtype": "uint8", "nodata": 0}) as dst:
    dst.write(classes, 1)

tally = zonal_stats(
    zones, "landcover_demo.tif",
    categorical=True,
    category_map={10: "tree_cover", 40: "cropland", 80: "water"},
    nodata=0,
)
print(tally)
# [{'water': 2, 'cropland': 8, 'tree_cover': 6}, {'water': 2, 'cropland': 4}]

shares = pd.DataFrame(tally, index=zones.index).fillna(0)   # absent class = 0, not NaN

Look at the second zone: tree_cover is missing from the dict entirely. It is genuinely present in the riparian strip — two cells are 10% and 8% inside — but no tree_cover pixel centre falls in the zone, so the class vanishes from the tally rather than appearing with a small share. exactextract reports it, because a fraction of a cell is still a fraction:

frac = exact_extract(
    "landcover_demo.tif", zones,
    ["frac", "majority", "variety"],
    include_cols=["zone_id"],
    output="pandas",
)
print(frac)
#       zone_id  frac_10  frac_40  frac_80  majority  variety
# 0  field_block   0.3750   0.5000   0.1250        40        3
# 1     riparian   0.0345   0.6897   0.2759        40        3

frac expands to one frac_<value> column per class encountered and the row sums to 1; majority is the coverage-weighted modal class and variety the class count, which is rasterstats' unique stat under another name. Neither library should ever be asked for a mean on this raster — the average of class codes 10 and 80 is 45, which is not a class.

7. Know what changes at ten thousand zones.

rasterstats executes a Python loop: per feature it computes a window, reads it, rasterizes the geometry and runs a NumPy reduction. That is a handful of Python-level calls per zone plus one interpreter round-trip, and it is why cost scales cleanly but with a large constant. exactextract pushes the whole loop into C++ and returns once, so the per-zone constant is far smaller.

Two levers matter once the job is long enough to notice. exactextract accepts strategy="raster-sequential", which walks the raster once and updates every overlapping zone's accumulators, instead of the default "feature-sequential" that revisits the file per feature — a large win when zones are dense and the raster is chunked, at the cost of holding more state; max_cells_in_memory bounds that state. On the rasterstats side there is no such switch, so the fix is to shard the zone layer across processes with joblib and concatenate, giving each worker its own dataset handle because rasterio datasets are not fork-safe.

8. Apply the decision rule.

The rule that survives contact with real work is short. Reach for exactextract when the zones are small relative to the cell, when there are tens of thousands of them, or when the number has to be defensible to a regulator or a client — coverage weighting is the only choice you never have to justify. Reach for rasterstats when the statistic is unusual enough to need add_stats, when the raster has no nodata tag and you must supply one in the call, when you also need point_query, or when the deployment target cannot take a compiled wheel. For a few hundred large zones, either is correct and the one already in your environment wins.

Verification

The comparison is only meaningful if the agreement case really agrees. Assert both halves — identical on the grid-aligned zone, materially different on the sliver — so a future library upgrade that changes a default breaks the test rather than the report.

merged = comparison.set_index("zone_id")

# Whole cells only: coverage weight 1.0 collapses to a boolean mask, so the means match
block = merged.loc["field_block"]
assert abs(block["mean"] - block["rs_mean"]) < 1e-5, "aligned zone must agree exactly"
assert block["count"] == float(block["rs_count"])          # 16.00 == 16

# Partial cells: the answers must differ, and exactextract must count fewer cells
strip = merged.loc["riparian"]
assert strip["count"] < strip["rs_count"]                  # 5.22 < 6
assert abs(strip["mean"] - strip["rs_mean"]) > 0.01, "no divergence — check the geometry"

# The coverage fractions are the mean's denominator: reproduce it by hand
weights = np.array([0.80, 1.00, 1.00, 0.10, 0.64, 0.80, 0.80, 0.08])
values = np.array([0.14, 0.61, 0.66, 0.68, 0.15, 0.62, 0.67, 0.69])
assert abs(np.average(values, weights=weights) - strip["mean"]) < 1e-4

print(f"aligned {block['mean']:.4f} | strip rasterstats {strip['rs_mean']:.4f} "
      f"vs exactextract {strip['mean']:.4f}")
# aligned 0.5281 | strip rasterstats 0.4750 vs exactextract 0.5044

The third assertion is the one worth keeping: it proves the reported mean is np.average(values, weights=coverage) and nothing more mysterious, which is usually the end of the argument about which number to publish.

Edge Cases & Debugging

• rasterstats returns {'count': 0, 'mean': None} for narrow zones. No pixel centre fell inside — the zone is thinner than a cell. Switch to exactextract, which still returns a coverage-weighted mean, or accept all_touched=True and its outward bias.
• exactextract means look contaminated but rasterstats looks fine. The file carries no nodata tag, so there is nothing for exactextract to honour while rasterstats took your nodata= argument. Write the tag with gdal_edit.py -a_nodata -9999 scene.tif, or pass a masked rioxarray DataArray (da.where(da != -9999)) as the raster source instead of the path.
• KeyError from include_cols. The named column is not on the GeoDataFrame — usually renamed by an upstream merge. Check zones.columns before the call; exact_extract will not invent the key.
• Columns come back as band_1_mean instead of mean. The source is multi-band or explicitly named, so operations are band-prefixed. Print ee.columns and select by suffix rather than hard-coding names.
• A MergeError when joining results back. validate="one_to_one" caught a duplicated zone_id, typically from an earlier explode() on multipart zones. Dissolve back to one row per id first, as in Computing NDVI Zonal Means per Parcel.

Frequently Asked Questions

Should every project just default to exactextract? Not automatically, but it is the better default for polygon zonal means. The cases where rasterstats is genuinely required are specific: a statistic no operation covers (add_stats takes any callable over the masked array), a raster whose nodata you must override at call time, point sampling through point_query, or an environment that cannot install compiled wheels. Everything else — smaller error, honest counts, several times the throughput, a DataFrame straight out — favours exactextract.

Can I make rasterstats reproduce exactextract's number? Not exactly, only asymptotically. The boolean mask has no representation for "48% of this pixel", so the two converge only as the zone grows relative to the cell. You can bracket the exact answer by running all_touched=False and all_touched=True — the coverage-weighted mean lies between them, as it does in the figure above — and you can shrink the gap by upsampling the raster before masking, which trades memory for a finer approximation of the same integral exactextract solves analytically. If the bracket is tight enough for your purposes, the question is moot; if it is not, that is the signal to switch.

What happens to a zone smaller than a single pixel? exactextract returns the containing cell's value with a count below 1 — a mean resting on one measurement, which is at least visible as such in the count column. rasterstats returns count: 0 and None for every statistic unless the centroid happens to fall inside, and all_touched=True promotes the same one or two cells to full weight. Either way the honest reading is that the raster does not resolve the zone; filter on count and report those zones as unmeasured rather than shipping a number with no support behind it.

Do the two libraries agree on sum the way they agree on mean? No, and the difference is larger than for the mean. rasterstats sums the selected pixel values; exactextract sums value × coverage fraction, which is the fraction-weighted total. For an extensive quantity — biomass, population, rainfall volume — the exactextract sum multiplied by the cell ground area is the quantity actually inside the polygon, while the pixel sum over-counts every boundary cell it accepted. Compare sum only after deciding whether you want "the total over the pixels I selected" or "the total inside the zone"; they are different questions.