Fixing NoData and NaN Propagation in Band Math

TL;DR: Read bands with src.read(masked=True), convert the mask to NaN-float, guard the ratio with np.where inside np.errstate, set the output profile’s nodata to np.nan, and assert the output NaN count equals the input nodata count so no fill pixel ever reaches np.mean.

Why NoData Silently Corrupts Index Statistics

A drone orthomosaic is a rectangle, but the flight footprint inside it is an irregular polygon. Every pixel outside that polygon is fill — usually 0 for a uint16 product or -9999 for a float export. If the GeoTIFF profile never had its nodata value set, numpy has no way to know those pixels are not data. They are 0 reflectance, and 0 reflectance is a perfectly legal number.

The damage is entirely statistical and entirely silent. Nothing raises. The NDVI over the actual canopy is fine. But ndvi.mean() now averages tens of thousands of fill pixels — which compute to an NDVI of 0 or, over -9999 fill, to a wildly out-of-range value — into the field statistic. On a typical 40-hectare field where the footprint fills 60% of the raster bounding box, the unmasked mean can land 20–35% below the true canopy mean, and the histogram grows a false spike at index 0 that swamps the real vegetation distribution. Feed that biased mean into a management-zone threshold and a healthy field can be misclassified as stressed.

This is the most common of the faults catalogued in the parent guide, Debugging Raster and Band Math Errors, and unlike overflow or a broadcast error it never announces itself — you only notice when a zonal statistic looks wrong months later. The fix is to honor the nodata sentinel from the moment the band is read all the way to the written output.

Why an unmasked nodata sentinel poisons every statistic The upper path takes a raster whose nodata value is minus 9999 straight into band math; the sentinel survives the arithmetic and the resulting mean is minus 4123, a number with no physical meaning. The lower path masks nodata to NaN before the arithmetic and computes a nanmean of 0.62, the correct field mean. Reflectance raster nodata = −9999 18% of pixels outside the field band math, unmasked (nir − red) / (nir + red) sentinel survives index array holds −9999 and ±inf mean = −4123 no error raised anywhere mask first arr[arr == nodata] = np.nan NaN propagates safely errstate guards the zero denominator nanmean = 0.62 the field's real mean

Prerequisites

Only the packages used by the parent band math and raster algebra reference are needed:

TEXT 
rasterio==1.3.10
numpy==1.26.4

Install with:

BASH 
pip install rasterio==1.3.10 numpy==1.26.4

Input requirements: a two-or-more-band drone or satellite GeoTIFF in a projected CRS such as EPSG:32615. The nodata value should be recorded in the profile; if it is not, the script below accepts an explicit override so you can still mask correctly.

Step-by-Step

Five steps, and the last one is where most pipelines still lose Five steps for handling nodata safely: read the declared nodata value from the raster profile, cast to a float type so NaN can be represented, mask the sentinel to NaN before any arithmetic, guard the division against a zero denominator, and compute statistics with NaN-aware functions. Read the profile nodata is metadata Cast to float32 NaN needs a float type Mask to NaN before any arithmetic Guard division np.errstate Use nan-aware stats nanmean, not mean An integer raster cannot hold NaN at all, which is why the cast has to happen before the mask rather than after it.

Step 1 — Read bands as masked arrays

src.read(masked=True) returns a numpy.ma.MaskedArray whose mask is built from the dataset’s nodata value and internal mask band. This is the load-bearing call: it is the difference between numpy knowing which pixels are fill and numpy treating fill as reflectance. If the profile’s nodata is unset, pass an explicit value so the mask is still correct.

PYTHON 
import numpy as np
import rasterio

def read_masked(src, index: int, nodata_override=None):
    """Read a 1-based band as a masked array; honor an explicit nodata if the
    profile omits one. Returns a MaskedArray with fill pixels masked True."""
    band = src.read(index, masked=True)
    if src.nodata is None and nodata_override is not None:
        band = np.ma.masked_equal(band, nodata_override)
    return band

The dataset_mask() method gives the same footprint as a single uint8 array (0 = nodata, 255 = valid) if you prefer one shared mask across all bands rather than a per-band one — useful when only some bands carry the sentinel.

Step 2 — Convert to NaN-float and guard the division

Mixing a masked array with a plain array through np.where strips the mask and lets fill pixels back in, so convert the mask to NaN-float immediately with .filled(np.nan). Then run the normalized difference through the zero-denominator guard. np.where evaluates both branches, so the np.errstate context is required to suppress the divide-by-zero warning that fires before the NaN is substituted.

PYTHON 
def masked_to_nan(band) -> np.ndarray:
    """MaskedArray -> float32 array with masked pixels as NaN."""
    return band.astype("float32").filled(np.nan)


def safe_ndvi(nir: np.ndarray, red: np.ndarray) -> np.ndarray:
    """(NIR - Red) / (NIR + Red); NaN wherever a band is nodata or denom is 0."""
    denom = nir + red
    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = np.where(denom != 0, (nir - red) / denom, np.nan)
    return ndvi.astype("float32")

Because any arithmetic touching NaN yields NaN, a fill pixel masked in Step 1 stays NaN through the subtraction and division automatically — that is the propagation working for you instead of against you.

Step 3 — Propagate nodata into the output profile

A clean NaN index written to a profile that still says nodata=0 (or a uint16 dtype) casts NaN back to 0 on write, re-merging “undefined” with “index zero.” Update the profile to float32 with nodata=np.nan before writing.

Step 4 — Verify the NaN count and range

The decisive check: the number of NaN pixels in the output must equal the number of masked pixels on input (plus any genuine zero-denominator pixels). If they match, no fill pixel leaked into the valid data; if the output has fewer NaN, the mask was dropped somewhere.

The complete, directly runnable script:

PYTHON 
import numpy as np
import rasterio

INPUT_PATH = "rededge_mx_ortho_field7.tif"
OUTPUT_PATH = "ndvi_field7.tif"

# 1-based rasterio band indices for a MicaSense RedEdge-MX stack.
RED, NIR = 3, 5
# If the source profile omits nodata, state it explicitly (0 for uint16 fill).
NODATA_OVERRIDE = 0


def read_masked(src, index, nodata_override=None):
    band = src.read(index, masked=True)
    if src.nodata is None and nodata_override is not None:
        band = np.ma.masked_equal(band, nodata_override)
    return band


def safe_ndvi(nir, red):
    denom = nir + red
    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = np.where(denom != 0, (nir - red) / denom, np.nan)
    return ndvi.astype("float32")


with rasterio.open(INPUT_PATH) as src:
    assert src.crs is not None and src.crs.is_projected, "reproject to a metric CRS first"

    red_m = read_masked(src, RED, NODATA_OVERRIDE)
    nir_m = read_masked(src, NIR, NODATA_OVERRIDE)

    # Count fill pixels on input (a pixel is invalid if either band is masked).
    input_nodata = int((np.ma.getmaskarray(red_m) | np.ma.getmaskarray(nir_m)).sum())

    red = red_m.astype("float32").filled(np.nan)
    nir = nir_m.astype("float32").filled(np.nan)

    ndvi = safe_ndvi(nir, red)

    profile = src.profile.copy()
    profile.update(count=1, dtype="float32", nodata=np.nan)

    with rasterio.open(OUTPUT_PATH, "w", **profile) as dst:
        dst.write(ndvi, 1)
        dst.set_band_description(1, "NDVI (NaN = nodata)")

# ── Inline verification ───────────────────────────────────────────────────
with rasterio.open(OUTPUT_PATH) as out:
    result = out.read(1)

output_nan = int(np.isnan(result).sum())
valid = result[np.isfinite(result)]

# No fill pixel leaked into valid data: output NaN >= input nodata,
# the surplus being genuine zero-denominator pixels (water, deep shadow).
assert output_nan >= input_nodata, (
    f"nodata leaked: {output_nan} output NaN < {input_nodata} input nodata"
)
assert valid.min() >= -1.0001 and valid.max() <= 1.0001, (
    f"NDVI out of [-1, 1]: [{valid.min():.3f}, {valid.max():.3f}]"
)
print(f"input nodata: {input_nodata:,}  output NaN: {output_nan:,}")
print(f"masked field mean NDVI: {np.nanmean(result):.4f}")

The two printed numbers are your proof. If np.nanmean(result) is meaningfully higher than a naive result.mean() on the same data with nodata unmasked, you have quantified exactly how much the fill pixels were biasing the field statistic.

Gotchas & Edge Cases

  • np.where still divides by zero. It evaluates the true branch for every pixel before selecting, so without the np.errstate wrapper you get a RuntimeWarning flood and, on older numpy, inf written before the NaN substitution. Keep the divide="ignore", invalid="ignore" context around the ratio.
  • Mixing a masked array with a plain array drops the mask. The instant you combine src.read(masked=True) output with an unmasked array via np.where or +, the mask is gone and fill pixels return as raw values. Convert to NaN-float with .filled(np.nan) at the read boundary and stay in float land thereafter.
  • -9999 fill without a set nodata is worse than 0 fill. A -9999 reflectance produces a huge-magnitude NDVI that blows the array min/max and the color ramp, whereas 0 fill quietly biases the mean. Both are fixed by masking on read, but only the range assert in Step 4 catches the -9999 case.
  • Writing NaN to an integer profile. float NaN cast to uint8/uint16 becomes 0. Always set dtype="float32", nodata=np.nan, or scale to an explicit integer sentinel like 255 and record it as the profile nodata.
Nodata is a per-dtype decision, and mixing them is the bug A table of four data types used in an imagery pipeline with the nodata convention each supports and what breaks when conventions are mixed: unsigned integers need an explicit sentinel, floats should use NaN, class rasters usually reserve zero, and a pipeline mixing them forces every consumer to remember both. Data type Nodata convention What breaks if mixed uint16 reflectance no NaN available An explicit sentinel, often 0 The sentinel enters the arithmetic float32 index NaN is available NaN, and nothing else A sentinel needs masking everywhere uint8 zone classes 0 usually means none 0 as 'unclassified' Averaging includes class zero Mixed in one pipeline the usual case Convert at the boundary Every consumer must remember both

Frequently Asked Questions

Why is my field mean NDVI lower than the values I see on the map?

The nodata fill pixels outside the flight footprint are being averaged in as real zeros. If the raster profile has no nodata value set, numpy treats those 0 or minus 9999 pixels as data and drags the mean down. Read the bands as masked arrays, convert the mask to NaN, and compute the field statistic with np.nanmean so only valid canopy contributes.

Should I use zero or NaN as the nodata value for a computed index?

Use NaN for the in-memory float index, because zero is a real NDVI value for bare soil and water and must stay distinguishable from undefined. When writing the output, set the profile nodata to NaN for a float32 raster, or scale to an integer sentinel such as 255 if the downstream tool needs an integer type.

Does np.where evaluate the divide-by-zero before it selects?

Yes. np.where computes both the true and false branches for every pixel and only then chooses, so the division runs on the zero-denominator pixels and emits a runtime warning before the NaN is substituted. Wrap the expression in a numpy errstate context with divide and invalid set to ignore so the warning is suppressed while the guard still does its job.

Parent Guide

This guide is part of Debugging Raster and Band Math Errors — see there for the full set of band-math failure classes, from dtype overflow and band-index off-by-one to broadcast mismatches and the shared diagnosis-to-fix workflow.

Parent section Debugging Raster and Band Math Errors