Computing NDVI Zonal Means per Parcel

One NDVI number per land parcel looks like a two-line calculation and almost never is: the division can overflow, cloud-masked pixels can arrive downstream as a convincing 0.0, and a parcel whose imagery is nine-tenths missing will hand you a mean that looks exactly like a real one. This guide is for anyone building an agronomy, subsidy-audit or crop-monitoring table and needing every row to be either trustworthy or visibly marked. It sits under Zonal Statistics & Raster Sampling in Spatial Analysis & Advanced Query Techniques, and assumes you can already open a scene with Rasterio.

Why This Approach / What Goes Wrong

NDVI is (NIR − red) / (NIR + red), and each of those three operations has a trap. The sum is where integer rasters break first: Sentinel-2 L2A bands are int16 and NAIP orthos are uint8, so nir + red is evaluated in the source dtype and can wrap silently — two 8-bit values of 200 sum to 400, which stores as 144, and the index that comes out is arithmetically valid and physically nonsense. Promoting both bands to float32 before any arithmetic removes the wrap and gives the division a sane denominator. The division itself then needs a guard, because in deep shadow or over calibrated black targets nir + red approaches zero and the quotient explodes to ±10⁵ or to inf.

The more insidious failure is the nodata path. A scene's fill value — 0 for Sentinel-2, or whatever the cloud mask wrote — is a perfectly ordinary integer as far as NumPy is concerned. If it flows into the ratio, the result is not an error but a number, usually something near 0.0, which downstream reads as bare soil rather than "no observation". By the time it reaches a parcel mean it is indistinguishable from data. The fix is to establish a boolean validity mask once, before the division, and to carry it all the way through to the output raster's own nodata value — never letting an unobserved pixel acquire a plausible number at any stage.

The last problem is geometric. Parcel boundaries do not follow the 10 m grid, so most edge pixels are partly inside and partly outside. Counting each touched pixel wholly in or wholly out biases the mean toward whatever the boundary happens to clip — badly for narrow strips and headlands. Weighting each pixel by the fraction of its cell that the polygon actually covers removes that bias, which is what exactextract computes exactly rather than by sampling; the trade-offs against the pure-Python alternative are laid out in zonal statistics with rasterstats vs exactextract.

Prerequisites

• rasterio>=1.3 — band reads with masked=True, the affine transform, and the GeoTIFF writer
• numpy>=1.24 — np.divide(..., where=...) and np.select for the flag ladder
• geopandas>=1.0 — vector I/O, to_crs, make_valid, and to_parquet with covering-bbox support
• exactextract>=0.2 — coverage-weighted statistics, accepting a GeoDataFrame directly
• pyarrow>=14 — the Parquet writer behind GeoPandas

conda install -c conda-forge "rasterio=1.3.*" "numpy=1.26.*" \
  "geopandas=1.0.*" "exactextract=0.2.*" "pyarrow=14.*"

Install everything from conda-forge in one solve: rasterio, geopandas and exactextract all bind GDAL and GEOS, and mixing wheels across channels produces the ABI mismatches that surface as import-time segfaults rather than clean errors.

Step-by-Step Implementation

1. Read both bands as masked arrays and prove they share a grid.

Band math is per-pixel, so the two rasters must be co-registered — same CRS, same affine transform, same shape. Overlapping is not enough. Sentinel-2 delivers B04 and B08 at 10 m on the identical tile grid, so the assertions pass; anything resampled or reprojected upstream must be aligned first with reproject_match on a raster cube.

import numpy as np
import rasterio

RED_PATH = "s2_b04_red_10m.tif"   # Sentinel-2 L2A red
NIR_PATH = "s2_b08_nir_10m.tif"   # Sentinel-2 L2A near-infrared

with rasterio.open(RED_PATH) as red_src, rasterio.open(NIR_PATH) as nir_src:
    assert red_src.crs == nir_src.crs, "Bands are in different CRSs"
    assert red_src.transform == nir_src.transform, "Bands sit on different grids"
    assert red_src.shape == nir_src.shape, "Bands differ in size"

    # masked=True honours the file's nodata value; astype keeps the mask
    red = red_src.read(1, masked=True).astype("float32")
    nir = nir_src.read(1, masked=True).astype("float32")

    raster_crs = red_src.crs
    profile = red_src.profile
    pixel_area_m2 = abs(red_src.transform.a * red_src.transform.e)  # 100.0 at 10 m

print(f"masked red pixels: {np.ma.getmaskarray(red).sum():,}")

2. Compute NDVI behind a validity gate, then clip.

filled(np.nan) turns the mask into NaN, which np.isfinite then folds into a single boolean array alongside the near-zero denominator test. The output buffer is pre-filled with the sentinel, and np.divide(..., where=valid) writes only where the gate passed — invalid cells keep -9999 because nothing ever touches them. Clipping afterwards is applied only to the valid subset, so the sentinel is not dragged up to -1.0.

NDVI_NODATA = np.float32(-9999.0)

red_f = red.filled(np.nan)      # plain ndarray; masked cells become NaN
nir_f = nir.filled(np.nan)

denom = nir_f + red_f
valid = np.isfinite(denom) & (np.abs(denom) > 1e-6)

ndvi = np.full(red_f.shape, NDVI_NODATA, dtype="float32")
np.divide(nir_f - red_f, denom, out=ndvi, where=valid)
ndvi[valid] = np.clip(ndvi[valid], -1.0, 1.0)

print(f"valid NDVI pixels: {valid.sum():,} of {valid.size:,}")

Because both bands were promoted to float32 first, the subtraction and the sum are floating point — no int16 wraparound, and no RuntimeWarning: invalid value encountered in divide, since the division never runs on the bad cells.

3. Write the NDVI surface as a tiled, compressed GeoTIFF.

predictor=3 is the floating-point predictor and typically halves the DEFLATE output on a smooth index surface. Declaring the sentinel as the file's nodata is the step that makes the mask legible to every downstream reader, including exactextract.

profile.update(
    driver="GTiff",
    dtype="float32",
    count=1,
    nodata=float(NDVI_NODATA),
    compress="deflate",
    predictor=3,
    tiled=True,
    blockxsize=512,
    blockysize=512,
    BIGTIFF="IF_SAFER",
)

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

Adding overviews here pays for itself if the same file also feeds a web map — the settings are covered in resampling and overviews when writing COGs.

4. Bring the parcels onto the raster's CRS.

exactextract compares polygon coordinates against the raster's grid directly and will not reproject for you; mismatched inputs return all-NaN rather than an error. Reproject the vector side, because the Sentinel-2 tile is already in a metric UTM zone and warping the raster would resample the values you are about to average. GeoPandas routes to_crs through pyproj with always_xy=True semantics, so you never hand-swap coordinates even though EPSG:4326's authority axis order is latitude-then-longitude — the details live in Coordinate Systems with PyProj.

import geopandas as gpd

parcels = gpd.read_file("parcels_2024.gpkg", layer="parcels")

if parcels.crs != raster_crs:
    parcels = parcels.to_crs(raster_crs)

assert parcels.crs.is_projected, "Coverage ratios need a projected, metric CRS"

parcels = parcels[parcels.geometry.notna() & ~parcels.geometry.is_empty].copy()
parcels["geometry"] = parcels.geometry.make_valid()   # bad rings break coverage math
parcels["parcel_area_m2"] = parcels.geometry.area

Self-intersecting rings make coverage fractions meaningless before they make anything crash, so repair them up front — see fixing self-intersecting polygons programmatically.

5. Run the coverage-weighted extraction.

exact_extract accepts a raster path and a GeoDataFrame and returns a DataFrame in feature order. mean is the coverage-weighted mean over cells that hold data, and count is the sum of coverage fractions across those same cells — an effective pixel count, so a parcel with 391.4 means it received the equivalent of 391.4 fully-covered valid cells. That second number is the one that makes the low-coverage check possible.

from exactextract import exact_extract

stats = exact_extract(
    "parcel_ndvi_surface.tif",
    parcels,
    ["mean", "count", "min", "max", "stdev"],
    include_cols=["parcel_id"],
    output="pandas",
)

print(stats.head(3).to_string(index=False))
# parcel_id     mean      count    min    max  stdev
#    101884 0.712043 391.402710 0.4120 0.8641 0.0714
#    101885 0.630918 118.994141 0.3805 0.7712 0.0663
#    101886 0.481277  27.331251 0.1002 0.7440 0.1918

For a single-band raster the columns are named after the operation; a multi-band input prefixes them band_1_, band_2_ and so on. Passing a rioxarray DataArray instead of a path skips the intermediate GeoTIFF entirely, which is worthwhile inside a raster cube workflow.

6. Turn coverage into an explicit status, and withhold the untrustworthy means.

Two thresholds do the work. A relative one — the share of the parcel area backed by valid pixels — catches cloud, shadow and scene-edge gaps. An absolute one catches slivers and access strips narrower than the grid, where a high coverage ratio can still rest on two or three cells. Neither threshold deletes anything: every parcel keeps its row, and ndvi_mean is populated only for the rows that earned it.

MIN_COVERAGE = 0.60          # share of parcel area backed by valid NDVI pixels
MIN_EFFECTIVE_PIXELS = 4.0   # guards slivers narrower than the 10 m grid

result = parcels.merge(stats, on="parcel_id", how="left", validate="one_to_one")

effective_cells = result["count"].fillna(0.0)
result["valid_coverage"] = effective_cells * pixel_area_m2 / result["parcel_area_m2"]

result["status"] = np.select(
    [
        effective_cells == 0,
        effective_cells < MIN_EFFECTIVE_PIXELS,
        result["valid_coverage"] < MIN_COVERAGE,
    ],
    ["no_data", "too_small", "low_coverage"],
    default="ok",
)
result["ndvi_mean"] = result["mean"].where(result["status"] == "ok")

7. Write the result as GeoParquet.

Parquet keeps float64 statistics at full precision and does not truncate valid_coverage to a ten-character field name, which is the practical reason this table is not a shapefile — the full comparison is in GeoParquet vs Shapefile for storage. The covering-bbox column lets readers push spatial filters down without decoding geometry.

out = result[[
    "parcel_id", "parcel_area_m2", "ndvi_mean", "mean", "min", "max",
    "stdev", "count", "valid_coverage", "status", "geometry",
]].rename(columns={"mean": "ndvi_mean_raw", "count": "valid_cell_equivalents"})

out.to_parquet(
    "parcel_ndvi_2024.parquet",
    index=False,
    compression="zstd",
    geometry_encoding="WKB",
    schema_version="1.1.0",
    write_covering_bbox=True,
)

Verification

Check the raster side and the table side separately: the sentinel must be untouched by the clip, and the parcel count must be identical before and after the join.

import geopandas as gpd

# Raster: valid cells are inside the physical range, invalid cells are exactly the sentinel
assert ndvi[valid].min() >= -1.0 and ndvi[valid].max() <= 1.0
assert (ndvi[~valid] == NDVI_NODATA).all(), "Sentinel was overwritten by the clip"

check = gpd.read_parquet("parcel_ndvi_2024.parquet")

assert len(check) == len(parcels), "Row count changed — the merge fanned out"
assert check.crs == parcels.crs, "CRS lost on the Parquet round-trip"

usable = check[check["status"] == "ok"]
assert usable["ndvi_mean"].between(-1.0, 1.0).all(), "Nodata leaked into a mean"
assert check.loc[check["status"] != "ok", "ndvi_mean"].isna().all()

print(check["status"].value_counts().to_string())
# ok              2283
# low_coverage      96
# no_data           38
print(f"mean NDVI over usable parcels: {usable['ndvi_mean'].mean():.3f}")
# mean NDVI over usable parcels: 0.634

If low_coverage dominates on a clear-sky scene, the thresholds are not the problem — the parcels are probably straddling the tile edge, or one of the two bands was masked more aggressively than the other.

Edge Cases & Debugging

• Every mean is NaN. The parcels are not in the raster's CRS; exactextract never reprojects, so nothing intersects. Re-run step 4 and confirm parcels.crs == raster_crs.
• NDVI is biased low across the whole scene. Sentinel-2 processing baseline 04.00 added BOA_ADD_OFFSET = -1000; the offset does not cancel in a ratio. Apply red_f -= 1000.0 and nir_f -= 1000.0 before computing denom.
• NaN used as the raster nodata value. GDAL stores it, but nan != nan, so masked-value comparison silently fails in some readers. Use a finite sentinel such as -9999.0, as above.
• A parcel smaller than one pixel returns a mean anyway. Coverage weighting happily averages the two or three cells it touches. MIN_EFFECTIVE_PIXELS catches these; for genuinely point-like features, sample raster values at point locations instead.
• Extraction crawls on a national scene. The default strategy="feature-sequential" re-reads raster blocks per feature. Pass strategy="raster-sequential" when parcels are dense relative to the raster, or restrict the read window first as in windowed reads from a Cloud-Optimized GeoTIFF.

Frequently Asked Questions

Should I average NDVI per parcel, or compute NDVI from averaged bands? Average the NDVI. The index is a nonlinear ratio, so the mean of the ratio and the ratio of the means are different numbers, and the second is pulled toward bright, high-reflectance pixels — roads, roofs and dry headlands inside the parcel boundary. Averaging per-pixel NDVI is the convention every agronomic reference assumes, so mixing the two makes your values incomparable with everyone else's.

Why weight by coverage instead of just masking the raster? A binary mask forces a yes-or-no decision on every boundary pixel, and both defaults are wrong: centroid-based inclusion under-counts narrow fields, while all_touched over-counts by pulling in cells barely clipped by the edge. Coverage weighting assigns each cell the exact fraction the polygon covers, so an edge cell half inside contributes half as much. The bias only disappears for large, compact parcels — for strips and small fields it is the dominant error, as rasterstats vs exactextract quantifies.

What coverage threshold should I actually use? There is no universal number; it depends on parcel size relative to pixel size and on what the mean feeds. At 10 m over roughly one-hectare parcels, 0.60 is a reasonable default because a parcel that size holds around 100 cells, so 60 % still rests on a real sample. The important part is not the threshold but shipping valid_coverage as a column, so a downstream consumer can raise the bar without re-running the extraction.

Can I run this over a time series without writing a GeoTIFF per date? Yes — exact_extract accepts an in-memory rioxarray DataArray, so a stacked cube can loop over the time dimension and extract each slice without touching disk, appending a date column per pass. That path is described in xarray & rioxarray raster cubes; the intermediate GeoTIFF is only worth writing when the NDVI surface itself is a deliverable.