Writing a Safe NDVI with NoData Handling

The NDVI raster looks right in the viewer and its mean comes back as nan. Or it comes back as 0.43 when every visible pixel is clearly around 0.7. Or the histogram has a spike at exactly 1.0 and another at −1.0 that correspond to nothing on the ground.

All three come from the same source: an NDVI written as (nir - red) / (nir + red) on arrays that contain values the formula was not designed for. The formula is correct; the guards around it are missing.

This page covers the five guards a production NDVI needs, why each is necessary, and how to test that each works. It is the defensive detail behind computing vegetation indices with rasterio.

The five ways a naive NDVI goes wrong

Sentinel NoData treated as data. A band with nodata = -9999 read without masking produces an NDVI of roughly −1 or +1 over every NoData pixel, and those pixels are then indistinguishable from genuine extremes.

Division by a near-zero denominator. Where both bands are small, the ratio amplifies band noise without bound. In floating point the result is a large finite number rather than an error.

Division by exactly zero. Produces inf or nan depending on the numerator, and NumPy emits a warning that is routinely suppressed at the top of the file.

Integer arithmetic. (nir - red) on unsigned integer arrays wraps around when red exceeds nir, producing a huge positive number where a small negative one belongs. This one is spectacular and mercifully obvious.

NaN propagating into statistics. A single NaN makes np.mean return NaN for the whole array, which is at least loud — unlike the cases above.

Five input conditions and what a naive NDVI returns for each A table of five problem inputs against the naive result and the guarded result. A sentinel NoData of minus nine thousand nine hundred and ninety-nine yields a naive value near plus one and a guarded value of NaN. A near-zero denominator yields an unbounded naive value and a guarded NaN. An exactly zero denominator yields infinity or NaN with a suppressed warning, and a guarded NaN. Unsigned integer subtraction wraps to a large positive value naively and is guarded by casting to float. A NaN input propagates into every statistic naively and is excluded by nan-aware reductions. naive result guarded result sentinel −9999 denominator ≈ 0 denominator = 0 uint subtraction NaN in the array ≈ +1.00, looks real NaN, masked unbounded, finite NaN, below the floor inf or NaN, warning hidden NaN, never evaluated 65 000-ish, wrapped correct, cast to float every statistic is NaN excluded by nan-aware reductions Only the last one announces itself. The other four produce plausible numbers.

Figure 1 — Four silent failures and one loud one. The loud one is the least dangerous.

Minimal reproducible solution

import numpy as np


def ndvi(nir: np.ndarray, red: np.ndarray, *,
         nir_nodata: float | None = None, red_nodata: float | None = None,
         min_sum: float = 0.02, valid_range: tuple = (-0.05, 1.5)) -> np.ndarray:
    """Normalised difference vegetation index with every guard in place.

    Five defences, in the order they must be applied:
      1. cast to float, so integer subtraction cannot wrap
      2. replace sentinels with NaN before any arithmetic
      3. reject physically implausible reflectance
      4. require the denominator to clear a floor
      5. evaluate the division only where all of the above hold
    """
    a = np.asarray(nir, dtype=np.float64)
    b = np.asarray(red, dtype=np.float64)

    if nir_nodata is not None:
        a = np.where(a == nir_nodata, np.nan, a)
    if red_nodata is not None:
        b = np.where(b == red_nodata, np.nan, b)

    lo, hi = valid_range
    a = np.where((a < lo) | (a > hi), np.nan, a)
    b = np.where((b < lo) | (b > hi), np.nan, b)

    total = a + b
    valid = np.isfinite(total) & (total > min_sum)

    out = np.full(a.shape, np.nan, dtype=np.float32)
    np.divide(a - b, total, out=out, where=valid)
    return out

np.divide with a where argument is the key construct: it evaluates the division only at the positions the mask selects, so the invalid positions are never computed and no warning is raised. Writing np.where(valid, (a - b) / total, np.nan) looks equivalent and is not — the division is evaluated everywhere first, warnings and all, and only then discarded.

Requiring total > min_sum rather than abs(total) > min_sum is deliberate for reflectance. A negative sum means at least one band is negative, which is a calibration fault rather than a dark pixel, and NDVI is not defined there in any useful sense.

Testing each guard

A function with five guards deserves five tests, and they are short enough that there is no excuse.

import numpy as np


def test_guards():
    """Each guard, exercised in isolation."""
    # 1. integer inputs must not wrap
    nir = np.array([[100]], dtype=np.uint16)
    red = np.array([[200]], dtype=np.uint16)
    assert ndvi(nir, red, valid_range=(-1, 1000))[0, 0] < 0

    # 2. sentinels become NaN
    out = ndvi(np.array([[-9999.0]]), np.array([[0.4]]), nir_nodata=-9999)
    assert np.isnan(out[0, 0])

    # 3. implausible reflectance is rejected
    assert np.isnan(ndvi(np.array([[12.0]]), np.array([[0.2]]))[0, 0])

    # 4. a tiny denominator is masked rather than amplified
    assert np.isnan(ndvi(np.array([[0.004]]), np.array([[0.004]]))[0, 0])

    # 5. a valid pair computes correctly
    got = ndvi(np.array([[0.6]]), np.array([[0.2]]))[0, 0]
    assert abs(got - 0.5) < 1e-6

    # and no warnings are emitted anywhere
    with np.errstate(all="raise"):
        ndvi(np.array([[0.0, 0.5]]), np.array([[0.0, 0.1]]))

The np.errstate(all="raise") block at the end is the test that catches a regression back to an unguarded division. Any future edit that evaluates the ratio at an invalid position will raise there, in a test, rather than emitting a suppressed warning in production.

Four conditions a naive NDVI expression does not survive Four rows. A zero denominator, where red and near-infrared sum to zero, produces a division by zero that NumPy reports as infinity rather than raising. Nodata pixels participating in the arithmetic produce a plausible finite value from meaningless inputs, which is worse than an error because nothing flags it. Integer inputs cause the division to truncate before the ratio is formed, producing a raster of zeros and ones. A scaled reflectance product, where values are stored as integers multiplied by ten thousand, yields a correct-looking ratio only by coincidence and an incorrect one whenever the two bands carry different scales. a zero denominator NumPy returns infinity rather than raising — nothing stops the pipeline nodata in the arithmetic a plausible finite value from meaningless inputs, unflagged integer inputs the division truncates before the ratio forms — a raster of zeros and ones differently scaled bands a correct-looking ratio by coincidence, wrong whenever scales differ Each produces output. None produces an error. That is what makes them expensive.

Figure 3 — Four conditions that all return a raster.

Edge-case matrix

Input Naive Guarded
nir=-9999, red=0.3 ≈ −1.0 NaN
nir=0.004, red=0.004 0.0 or unbounded NaN
nir=0.0, red=0.0 NaN plus a warning NaN, no warning
uint16(100) - uint16(200) 65 436 correct negative
nir=12.0 (uncalibrated DN) plausible ratio NaN
nir=0.6, red=-0.05 > 1.0 NaN (band negative)
Whole array valid correct correct
One NaN in a big array mean is NaN excluded by nan-aware reductions

The fifth row deserves emphasis. Running an NDVI on uncalibrated digital numbers produces values that look entirely reasonable — the ratio form cancels the scale — which is why the mistake survives. The valid_range guard catches it, and the choice to define that range in reflectance units is what makes the function refuse data it should not be given.

Verification snippet

import numpy as np


def ndvi_report(index: np.ndarray) -> dict:
    """Distribution checks that catch a guard that did not fire."""
    finite = index[np.isfinite(index)]
    problems = []
    if finite.size == 0:
        return {"problems": ["no finite values"]}

    if np.any(np.abs(finite) > 1.0001):
        problems.append("values outside [-1, 1] — a band was negative")

    # A spike exactly at the extremes usually means unmasked NoData.
    at_extreme = float(np.count_nonzero(np.abs(np.abs(finite) - 1.0) < 1e-4)
                       / finite.size)
    if at_extreme > 0.001:
        problems.append(f"{at_extreme:.2%} of pixels sit exactly at ±1 — "
                        "NoData is probably being treated as data")

    masked = 1 - finite.size / index.size
    return {"masked_fraction": masked, "at_extreme_fraction": at_extreme,
            "median": float(np.median(finite)),
            "problems": problems}

The spike test is the most useful diagnostic here. Genuine NDVI values are distributed continuously; a sharp concentration at exactly ±1 means a population of pixels where one band was zero or a sentinel survived into the arithmetic, and it is visible in a histogram at a glance.

NDVI histogram with and without the guards Two overlaid histograms of NDVI values. The guarded distribution is a smooth bimodal curve with a soil peak near zero point one five and a canopy peak near zero point seven. The unguarded distribution adds two narrow spikes at exactly minus one and plus one, containing about four percent of pixels between them, which come from NoData and zero-denominator pixels. A note states that the spikes are invisible in a map view because they are scattered, and unmistakable in a histogram. −1.0 −0.2 0.3 0.7 +1.0 NDVI value soil canopy NoData spike zero-denominator spike Scattered across the map and therefore invisible there; unmistakable in a histogram.

Figure 2 — The histogram is the cheapest guard-failure detector there is.

Generalising beyond NDVI

Every normalised difference index shares this structure, so the guards belong in one function rather than being re-implemented per index.

def normalised_difference(a: np.ndarray, b: np.ndarray, **kwargs) -> np.ndarray:
    """The guarded form, for any band pair."""
    return ndvi(a, b, **kwargs)


INDEX_DEFINITIONS = {
    "ndvi": ("nir", "red"),
    "ndre": ("nir", "red_edge"),
    "gndvi": ("nir", "green"),
    "ndwi": ("green", "nir"),
}


def compute_named_index(bands: dict, name: str, **kwargs):
    """Compute any defined index from a band dictionary, by name."""
    if name not in INDEX_DEFINITIONS:
        raise KeyError(f"unknown index {name!r}; known: {sorted(INDEX_DEFINITIONS)}")
    a_name, b_name = INDEX_DEFINITIONS[name]
    missing = [n for n in (a_name, b_name) if n not in bands]
    if missing:
        raise KeyError(f"{name} needs bands {missing}, which are not present")
    return normalised_difference(bands[a_name], bands[b_name], **kwargs)

Defining the index by band names rather than positions keeps the failure mode explicit: asking for NDRE from a sensor without a red edge band raises rather than silently computing something else.

When to escalate

  • A large fraction of the raster is masked. The guards are working and the data is the problem. Look at the band statistics: usually the calibration produced negatives, or the mosaic has large NoData regions.
  • Values are outside [−1, 1] after guarding. A band is negative and the range guard was set too permissively. Fix the calibration rather than widening the guard.
  • The index is being computed on digital numbers deliberately. That is defensible for within-flight ranking and nothing else; state it in the deliverable, and set the valid range to match the sensor’s dynamic range rather than to reflectance.

← Computing Vegetation Indices with Rasterio