Debugging Raster and Band Math Errors

A vegetation-index pipeline rarely fails with a stack trace. It fails quietly: an NDVI raster that looks reasonable in a thumbnail but carries inf streaks along field edges, a mean index that has drifted because 40,000 nodata zeros were folded into the average, or a NDRE layer that is subtly wrong because src.read(4) returned the red edge band instead of NIR. This guide is a troubleshooting reference for the specific runtime faults that break band math on drone and satellite imagery in Python — divide-by-zero producing inf/NaN, NaN and nodata propagation through numpy arithmetic, integer dtype overflow, band-index off-by-one between rasterio and numpy, shape and broadcast mismatches, and unset nodata polluting index statistics. Each section pairs a diagnosis with a guarded fix, and the whole thing ends in a validation gate you can drop into any band math and raster algebra workflow.

This page is part of the Drone Imagery Processing & Vegetation Index Workflows guide. See that page for the full pipeline from ingest through temporal aggregation; this reference covers what to do when a stage in that pipeline throws or, worse, does not throw but returns corrupted numbers.

NaN and nodata propagation through band math, with guard interception points A single nodata pixel valued zero flows left to right through read, integer subtraction, division, and mean. The unguarded path produces overflow, then inf, then a corrupted mean. Three guard boxes below intercept the pixel at cast, at the errstate-plus-where division, and at the validation assert. Read pixel DN = 0 (nodata) uint16 subtract 0 − 3200 wraps divide by 0 → inf / NaN np.mean() → NaN result Corrupted index stats Guards that intercept the pixel: Guard A astype float32 nodata → NaN Guard B errstate + where denom==0 → NaN Guard C nanmean + range assert Each guard removes the fault before it reaches the next stage — fix upstream, not at the mean.

Prerequisites

Python packages (exact versions tested):

  • rasterio==1.3.10
  • numpy==1.26.4
  • scipy==1.13.0 (only for the optional morphological nodata infill)

Install with:

BASH 
pip install rasterio==1.3.10 numpy==1.26.4 scipy==1.13.0

Input data requirements:

  • A multi-band drone or satellite GeoTIFF in a projected CRS (for example EPSG:32615 for UTM zone 15N). A MicaSense RedEdge-MX stack (Blue 475, Green 560, Red 668, RedEdge 717, NIR 840 nm) or a DJI P4 Multispectral orthomosaic is assumed for the examples.
  • You must know two things about the file before doing arithmetic: the storage dtype (uint8, uint16, or float32) and the nodata value (commonly 0, 65535, or -9999). Both are read from src.profile in Step 1 — do not assume them.
  • The band-to-wavelength mapping. Band order is not standardised across export tools, so the 1-based rasterio index of red, red edge, and NIR must be confirmed, not guessed.

1. Concept: where band math actually breaks

Every vegetation index is a normalized difference of two reflectance bands — NDVI is (NIR - Red) / (NIR + Red), NDRE swaps red for the red edge band, SAVI adds a soil-adjustment constant. Mathematically trivial. The failures come from the fact that a raster is not a clean matrix of reflectance floats: it is an integer array with sentinel nodata values, read through a 1-based band API, that numpy will happily broadcast, overflow, and divide by zero without raising anything that stops your script.

Five distinct fault classes account for almost every corrupted index in practice.

Zero-denominator division. Over water, deep shadow, or nodata fill, NIR + Red is zero. 0 / 0 yields NaN; a division where the denominator is a tiny positive float yields a value that trends to inf. numpy emits a RuntimeWarning and keeps going, so the poison values land in your output raster and only surface later as a NaN mean or an inf-blown color ramp.

Integer dtype overflow. Raw drone bands are stored as uint8 (0–255) or uint16 (0–65535). Subtracting two uint16 bands, NIR - Red, when red exceeds NIR (bare soil, senescent canopy) wraps to a huge positive number near 65535 because unsigned integers cannot represent negatives. Addition can also overflow the type ceiling. The arithmetic must happen in float32.

Band-index off-by-one. rasterio.open(...).read(n) is 1-based: read(1) is the first band. A numpy array indexed stack[0] is the first band. Move a band read from rasterio into a numpy slice without adjusting and you read the neighbouring wavelength — a bug that produces a perfectly valid-looking but agronomically wrong index, because red edge and NIR are adjacent in most stacks.

Shape and broadcast mismatch. Reading two bands with different windows, or mixing a full-resolution NIR with a resampled 20 m SWIR, gives arrays of different shape. numpy either raises ValueError: operands could not be broadcast or, more dangerously, broadcasts a (1, N) row against an (N, N) grid and returns a wrong-shaped result that still writes to disk.

Unset nodata leaking into statistics. If the GeoTIFF profile has no nodata value, the 0 fill pixels outside the flight footprint are real data as far as numpy is concerned. They pull the field-mean NDVI toward zero and add a spurious spike at index 0 in every histogram. This one never raises — it just quietly biases every zonal statistic downstream.

Why it matters agronomically: an NDVI shifted by even 0.05 can move a pixel across a management-zone boundary in a variable-rate map. A single inf in the array blows the min/max used to scale a crop-health threshold, and a nodata-biased mean can flip a whole field’s fertility classification. Correctness here is not cosmetic.

2. Step-by-Step Diagnosis and Fix

One check splits every band-math bug into two families A decision diagram for a wrong-looking index. If the values fall outside the physically possible range, the cause is arithmetic — integer types, an unmasked nodata sentinel or a zero denominator. If the range is plausible but the map is wrong, the cause is band identity or missing calibration. Index output looks wrong is the range inside −1 to 1? yes Suspect band identity or calibration the arithmetic is sound, so check band order, then whether the input is reflectance at all no Suspect dtype or nodata values outside the physical range come from integer arithmetic, an unmasked sentinel or a zero denominator, in that order of likelihood

Step 1 — Inspect dtype, band count, and nodata first

Before writing a single arithmetic expression, read the three properties that determine which faults are possible. This is the cheapest debugging you will ever do.

PYTHON 
import rasterio
import numpy as np

INPUT_PATH = "rededge_mx_ortho_field7.tif"

# 1-based rasterio indices for a MicaSense RedEdge-MX stack.
# Confirm against YOUR export — band order is not standardised.
BANDS = {"blue": 1, "green": 2, "red": 3, "rededge": 4, "nir": 5}

with rasterio.open(INPUT_PATH) as src:
    print(f"CRS: EPSG:{src.crs.to_epsg()}")
    print(f"Band count: {src.count}")
    print(f"Dtype: {src.dtypes[0]}")
    print(f"Nodata: {src.nodata}")

    # Guard the off-by-one class: the highest index we will read must exist.
    assert src.count >= max(BANDS.values()), (
        f"Stack has {src.count} bands but band map needs index "
        f"{max(BANDS.values())} — check the band order."
    )
    # Flag the silent-leakage class early.
    if src.nodata is None:
        print("WARNING: nodata is unset — 0 pixels will pollute index stats")

If src.nodata is None and the dtype is an unsigned integer, you already know two of the five faults are live. If the dtype is uint8 or uint16, overflow is live the moment you subtract.

Step 2 — Cast to float32 and mask nodata to NaN

This single step neutralises both dtype overflow and nodata leakage. Casting to float32 makes negative differences and large sums representable; converting the sentinel value to NaN means it can no longer masquerade as reflectance zero.

PYTHON 
def read_band_as_float(src, index: int, nodata=None) -> np.ndarray:
    """Read a 1-based band, cast to float32, and set nodata pixels to NaN."""
    arr = src.read(index).astype("float32")     # Guard A: no more overflow
    nd = src.nodata if nodata is None else nodata
    if nd is not None:
        arr[arr == nd] = np.nan                  # sentinel can no longer pose as data
    return arr


with rasterio.open(INPUT_PATH) as src:
    red = read_band_as_float(src, BANDS["red"])
    nir = read_band_as_float(src, BANDS["nir"])

# Verify the cast actually happened and shapes agree.
assert red.dtype == np.float32 and nir.dtype == np.float32
assert red.shape == nir.shape, f"Shape mismatch: {red.shape} vs {nir.shape}"

The read() call also fixes the band-index class if — and only if — the BANDS map is correct, which is why Step 1 asserts on it. For the deeper nodata story, including reading through src.read(masked=True) and the dataset mask, see fixing nodata and NaN propagation in band math.

Step 3 — Guard the denominator with errstate and where

Never divide two bands directly. Wrap the division in np.errstate so the RuntimeWarning does not spam your logs, and use np.where to replace the result with NaN wherever the denominator is zero. This is the fix for the inf/NaN class.

PYTHON 
def safe_normalized_difference(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    """(a - b) / (a + b) with zero-denominator pixels set to NaN, not inf."""
    denom = a + b
    with np.errstate(divide="ignore", invalid="ignore"):
        out = np.where(denom != 0, (a - b) / denom, np.nan)
    return out.astype("float32")


ndvi = safe_normalized_difference(nir, red)

# There must be no infinities left, though NaN is expected over nodata/water.
assert not np.isinf(ndvi).any(), "inf survived — denominator guard is not covering all pixels"
print(f"NaN pixels (nodata + zero-denom): {np.isnan(ndvi).sum()}")

Two subtleties. First, np.where still evaluates (a - b) / denom for every pixel including the zero-denominator ones — that is why the errstate context is required, to suppress the warning from the division that happens before where selects. Second, NaN is the correct sentinel for “no valid index here,” not 0: a zero NDVI is a real value (bare soil, water) and must stay distinguishable from “undefined.” This is exactly the pattern used in calculating NDVI and NDRE with rasterio step by step.

Step 4 — Align band indices and shapes before broadcasting

The off-by-one and broadcast classes both surface here. Read every band through the same named map and the same window, and assert equal shapes before any binary operation.

PYTHON 
from rasterio.windows import Window

win = Window(col_off=0, row_off=0, width=2048, height=2048)

with rasterio.open(INPUT_PATH) as src:
    # Same window for every band — mismatched windows are the usual cause
    # of "operands could not be broadcast together".
    red = src.read(BANDS["red"], window=win).astype("float32")
    nir = src.read(BANDS["nir"], window=win).astype("float32")
    rededge = src.read(BANDS["rededge"], window=win).astype("float32")

# Broadcast guard: fail loudly rather than let numpy silently stretch a row.
for name, band in [("nir", nir), ("rededge", rededge)]:
    assert band.shape == red.shape, (
        f"{name} shape {band.shape} != red {red.shape}; "
        "check window and resampling before band math"
    )

ndvi = safe_normalized_difference(nir, red)
ndre = safe_normalized_difference(nir, rededge)

If a band lives at a coarser native resolution (a resampled SWIR, for example), resample it to the target grid with an explicit out_shape and Resampling method during the read call rather than letting it reach the arithmetic at the wrong shape. A coarser band read without resampling is the single most common source of the broadcast ValueError in multi-sensor stacks discussed in the band math and raster algebra reference.

Step 5 — Validate the output before writing

Poison values that slipped through every guard get caught here. Assert on the physically valid range and on a sane NaN fraction, and use np.nanmean so the field statistic ignores undefined pixels.

PYTHON 
def validate_index(index: np.ndarray, name: str = "NDVI") -> dict:
    finite = np.isfinite(index)
    assert not np.isinf(index).any(), f"{name} contains inf"

    vals = index[finite]
    assert vals.min() >= -1.0001 and vals.max() <= 1.0001, (
        f"{name} out of [-1, 1]: min {vals.min():.3f}, max {vals.max():.3f} "
        "— likely dtype overflow or wrong band pairing"
    )

    nan_frac = float((~finite).sum() / index.size)
    if nan_frac > 0.5:
        print(f"WARNING: {nan_frac:.0%} of {name} is NaN — check nodata handling")

    return {"mean": float(np.nanmean(index)), "nan_fraction": nan_frac}


print(validate_index(ndvi, "NDVI"))
print(validate_index(ndre, "NDRE"))

An NDVI outside [-1, 1] is proof of an unresolved overflow or a red/NIR swap — the range check turns a silent agronomic error into a hard failure.

3. Key Parameters & Tuning

These are the values that determine whether band math is numerically safe. Set them explicitly; the defaults that libraries pick for you are usually wrong for agricultural imagery.

Parameter Type Default Agronomic Effect
nodata value int / float None (unset) Must be set. 0 for uint8/uint16 masks the flight-footprint fill; -9999 is the float convention. Left unset, footprint zeros bias every zonal mean toward 0 and add a false histogram spike at index 0.
working dtype numpy dtype source dtype Force float32. Doing math in uint8/uint16 overflows on NIR - Red over bare soil and wraps to ~65535, throwing NDVI far outside [-1, 1]. float64 doubles memory for no precision gain at reflectance scale.
errstate policy context numpy warns, continues Set divide="ignore", invalid="ignore" around the ratio only. Silencing globally hides real bugs; leaving it default floods logs and still writes inf to disk.
denominator guard np.where vs masked array none np.where(denom != 0, ratio, np.nan) keeps a plain array; numpy.ma masked arrays preserve the mask through later ops but many rasterio writers drop it. Pick one convention per pipeline.
output clip range tuple none Clip valid index to [-1, 1] only after validation, never before — clipping first hides an overflow that the range assert would otherwise catch.
band map indices dict (1-based) assumed order Record the 1-based rasterio index per wavelength and assert src.count. Guessing red-edge vs NIR order silently produces a wrong-but-plausible index.

4. Edge Cases & Failure Modes

A uint16 subtraction that wraps instead of going negative. Over senescing or bare-soil pixels red reflectance exceeds NIR, so NIR - Red is negative. In uint16 there is no negative — the result wraps to roughly 65500 and the NDVI numerator becomes enormous. The tell is an NDVI histogram with a spike near +1 exactly where the crop is dead. Casting to float32 in Step 2 is the fix; the Step 5 range assert is the detector.

Four symptoms and the fastest confirming test for each A table of four observations from a bad index raster — uniformly high values, an inverted map, a rim of apparent stress at the field edge, and scattered impossible values — with the most likely cause and the quickest test that confirms it. Observation Most likely cause Confirm by Index is uniformly high everywhere Input is digital numbers, not reflectance Check the value range of a band Index is a mirror image sign inverted Red and NIR swapped Median red versus median NIR A rim of stress at the edge field boundary Mixed pixels or zero padding Buffer inward and recompute Scattered impossible values isolated Zero denominator on nodata Count non-finite values

Masked-array versus np.where mismatch. src.read(masked=True) returns a numpy.ma.MaskedArray. If you then combine it with a plain array via np.where, the mask is stripped and nodata pixels re-enter as their raw fill value. Either stay in masked-array land end to end (np.ma.masked_where, .filled(np.nan)) or convert to NaN-float immediately on read. Mixing the two is a leading cause of nodata reappearing after you thought it was masked — the dedicated pattern is in fixing nodata and NaN propagation in band math.

NaN written to an integer output profile. You compute a clean float NDVI with NaN nodata, then write it to a profile still set to uint8 — NaN casts to 0, silently merging “undefined” with “index zero.” Update the output profile to dtype="float32", nodata=np.nan (or scale to a known integer nodata like 255) before writing.

Cloud-contaminated reflectance producing valid-but-wrong denominators. A cloud pixel has high, roughly equal reflectance across bands, so NIR + Red is large and the ratio is defined — no divide-by-zero, no NaN, just a wrong NDVI near zero over what should be canopy. The denominator guard cannot catch this; it must be removed upstream with cloud masking for agricultural imagery before band math runs.

Broadcast that succeeds silently. Reading a band that ended up shape (1, 2048) (a single-row window bug) against a (2048, 2048) grid broadcasts to (2048, 2048) without error — every row is identical. The Step 4 shape assert is the only thing standing between this and a striped index raster that passes every CRS check.

5. Verification & Output Validation

Correctness for band math is a small set of properties: no infinities, values in physical range, nodata accounted for, and a NaN fraction that matches the expected off-crop area. The guard functions above cover the compute path; this end-to-end check confirms the written file.

PYTHON 
import numpy as np
import rasterio

def audit_index_raster(path: str, valid_range=(-1.0, 1.0)) -> dict:
    """Post-write audit: catches overflow, inf leakage, and nodata bias."""
    with rasterio.open(path) as src:
        assert src.crs is not None, "index raster lost its CRS"
        data = src.read(1, masked=True).filled(np.nan).astype("float32")

    assert not np.isinf(data).any(), "inf present — denominator guard failed"

    finite = data[np.isfinite(data)]
    lo, hi = valid_range
    assert finite.min() >= lo - 1e-4 and finite.max() <= hi + 1e-4, (
        f"values outside {valid_range}: "
        f"[{finite.min():.3f}, {finite.max():.3f}] — overflow or band swap"
    )

    nan_frac = float(np.isnan(data).sum() / data.size)
    return {
        "mean": round(float(np.nanmean(data)), 4),
        "std": round(float(np.nanstd(data)), 4),
        "nan_fraction": round(nan_frac, 4),
        "min": round(float(finite.min()), 4),
        "max": round(float(finite.max()), 4),
    }


report = audit_index_raster("ndvi_field7.tif")
print(report)
assert report["nan_fraction"] < 0.9, "almost everything is NaN — masking is too aggressive"

A healthy corn-canopy NDVI at mid-season lands with a mean around 0.6–0.85 and a modest NaN fraction over headlands and water. A mean near zero with a large NaN fraction, or a max at exactly 1.0 with a spike, points straight back to the overflow or nodata classes above.

6. Integration with the Pipeline

These guards are not a separate stage — they wrap the arithmetic wherever it happens in your workflow.

Upstream. The band-index map and nodata value both come from ingest. If the stack was built by ingesting multispectral drone imagery with a consistent band order and a nodata value baked into the profile, Steps 1–2 collapse to a formality. Where cloud contributes valid-but-wrong denominators, run cloud masking for agricultural imagery first so those pixels arrive as NaN.

Core arithmetic. Every ratio in the band math and raster algebra reference and the step-by-step NDVI and NDRE guide should route its division through safe_normalized_difference and its output through validate_index.

The nodata deep-dive. The most common silent corruption — nodata pixels dragging index means and histograms — has its own focused walkthrough with the src.read(masked=True) and dataset-mask patterns in fixing nodata and NaN propagation in band math.

Downstream. A validated index raster with NaN nodata flows cleanly into temporal aggregation of vegetation indices, which relies on np.nanmean-style compositing ignoring undefined pixels rather than averaging in zeros.

Frequently Asked Questions

Why is my NDVI array full of inf and NaN values?

Both cases come from the denominator. Where NIR plus red equals zero the division is undefined and numpy returns NaN with a runtime warning; where the sum is a tiny non-zero float the result explodes toward inf. Cast both bands to float32, wrap the division in a numpy errstate context, and use np.where to substitute NaN wherever the denominator is zero.

Why do my index values wrap around to huge or negative numbers?

You are doing arithmetic in integer dtype. A uint8 or uint16 band overflows when you add or subtract bands because the result cannot exceed the type maximum and wraps modulo the range. Convert every band to float32 with astype before subtracting or adding, so NIR minus red can go negative and the sum can exceed 65535 without wrapping.

Why does rasterio read the wrong band from my stack?

rasterio band indices are 1-based, so src.read(1) returns the first band, whereas numpy array axes are 0-based. Mixing the two conventions reads red where you expected NIR and silently produces a plausible but wrong index. Keep a named band map that records the 1-based rasterio index for each wavelength and assert the band count before reading.

Parent section Drone Imagery Processing & Vegetation Index Workflows Related Band Math & Raster Algebra in Python Related Batch and Async Orthomosaic Processing Related Cloud Masking for Agricultural Imagery Related Radiometric Calibration & Reflectance Conversion Related Temporal Aggregation of Vegetation Indices Related Threshold Mapping for Crop Health: Field Zone Classification Related Vegetation Index Selection for Crop Stages