Avoiding Integer Overflow in Index Calculations

Cast to float before any arithmetic, and the problem cannot occur:

import numpy as np

red = red_uint16.astype("float32")
nir = nir_uint16.astype("float32")
den = nir + red
ndvi = np.divide(nir - red, den, out=np.full_like(den, np.nan), where=den != 0)

Every integer overflow bug in band math is a version of doing the arithmetic before that astype. This page belongs to raster dtypes, scaling and numerical precision in Core Raster Fundamentals & STAC Mapping.

What Wrapping Looks Like

Unsigned subtraction wraps around Over vegetation red is around 400 and near infrared around 3000. Computing red minus near infrared in uint16 should give minus 2600, but unsigned integers cannot be negative, so the result wraps to 62936. The same wrap happens for sums exceeding 65535, where two bright bands add to a small number. red − nir in uint16 over vegetation 0 65535 red 400 62936 should be -2600; falls off the left end and reappears on the right No warning, no error — just a large positive number where a negative one belongs.

The resulting NDVI is not merely inaccurate; it is a finite number that passes every range check once the ratio is taken, because the wrapped numerator and a possibly wrapped denominator combine into something between −1 and 1. That is why overflow bugs survive into published products.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Dtype control and guarded division
rasterio >=1.3.0 Reading bands in their stored dtype
pytest >=8.0 Extreme-value tests for index functions 
pip install "numpy>=1.23" "rasterio>=1.3.0" "pytest>=8.0"

Complete Working Example

import numpy as np


def as_float(*bands: np.ndarray) -> list[np.ndarray]:
    """Cast every input to float32 before any arithmetic touches it."""
    return [b.astype("float32", copy=False) for b in bands]


def normalised_difference(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    """(a - b) / (a + b), overflow-proof for any integer or float input."""
    a, b = as_float(a, b)
    den = a + b
    return np.divide(a - b, den, out=np.full(den.shape, np.nan, dtype="float32"),
                     where=den != 0)


def evi(nir: np.ndarray, red: np.ndarray, blue: np.ndarray, *,
        scale: float = 1e-4) -> np.ndarray:
    """EVI needs reflectance, not raw counts — convert, then compute."""
    nir, red, blue = (x * np.float32(scale) for x in as_float(nir, red, blue))
    den = nir + 6 * red - 7.5 * blue + 1
    return np.divide(2.5 * (nir - red), den,
                     out=np.full(den.shape, np.nan, dtype="float32"), where=den != 0)


if __name__ == "__main__":
    red = np.array([[400, 60000]], dtype="uint16")
    nir = np.array([[3000, 62000]], dtype="uint16")
    print("unsafe:", ((nir - red) / (nir + red)))          # wrapped
    print("safe:  ", normalised_difference(nir, red))

EVI makes the point that casting is necessary but not sufficient. Its formula contains an additive constant of 1 that only makes sense in reflectance units, so feeding it raw digital numbers — even as floats — produces a result dominated by the wrong term. Convert to the physical units the formula assumes, then compute; the scale conventions are set out in converting int16 reflectance to float safely.

Where Overflow Hides

Which operations wrap, and when Subtraction on uint16 wraps whenever the result would be negative, which is routine for index numerators. Addition wraps when the sum exceeds the maximum, which happens for bright pixels and for multi-band sums. Multiplication by a constant wraps quickly. Squaring for a standard deviation wraps almost immediately. Division alone is safe because it promotes to float in NumPy. Integer operations on reflectance-scale data operation when it wraps uint16 subtraction any time the result is negative — every vegetated pixel in red − nir squaring (variance) any value above 255 in uint16, above 181 in int16 addition of bands bright pixels, or sums over many bands division safe: NumPy true division promotes to float Division being safe is why the bug hides: the final step looks fine; the step before it was not.

Squaring deserves special mention because it appears inside standard deviations and texture measures, often buried in a helper function. uint16 values above 255 overflow when squared, which is essentially every reflectance value; a focal standard deviation computed on raw counts is meaningless everywhere. The focal statistics in feature engineering for pixel-based models are written against float input for exactly this reason.

Making It Impossible Rather Than Unlikely

The durable fix is to make integer arrays unable to reach arithmetic in the first place. Two habits achieve that.

The first is a single read function for every band in the pipeline that always returns float32 with nodata as NaN, so there is no code path that yields raw integers to a caller. The second is a type annotation discipline — npt.NDArray[np.float32] on every index function’s parameters — together with a runtime assertion at the top of each one. Neither is elaborate; together they mean a new contributor adding an index function receives floats whether or not they know about overflow.

For large pipelines, add a lint-style check that searches for arithmetic on arrays read with src.read() directly. It is crude, but it catches the one pattern responsible for nearly every overflow in practice: a quick script that reads two bands and subtracts them.

Overflow in Other Places

Band math is the most common site, but not the only one. Three others catch people regularly.

Accumulators in reductions inherit the input dtype unless told otherwise. np.sum over a large uint16 array promotes to a wider type on most platforms, but hand-written loops, Dask reductions with an explicit dtype, and some library helpers do not, and a sum over a whole scene overflows long before it finishes. Passing dtype="float64" for sums and dtype="float32" for means removes the question.

Histogram bins are a quieter case. Binning uint16 data into 256 bins with integer arithmetic on the edges truncates the bin width, so the top bin silently collects everything above the last full edge. Using float edges from np.linspace rather than integer division avoids it.

Resampling with integer output is the third. A bilinear resample of uint8 data can compute an intermediate slightly above 255 at a sharp edge; whether that clips or wraps depends on the library and its version. For continuous data, resample in float and cast afterwards; the resampling choices themselves are covered in choosing the right resampling method for Sentinel-2.

Verification

Test at the corners of the input space An index function is tested with inputs at zero, at mid-range and at the dtype maximum for each band, in every combination. The outputs must all lie within the index's valid range or be NaN, which catches every overflow path in a handful of cases. Nine cases catch every wrap 0 / 0 0 / mid 0 / max mid / 0 mid / mid mid / max max / 0 max / mid max / max max / max is the case unguarded addition fails first. 
import itertools

import numpy as np
import pytest


@pytest.mark.parametrize("dtype", ["uint16", "int16"])
def test_normalised_difference_never_leaves_range(dtype):
    info = np.iinfo(dtype)
    levels = [0, info.max // 2, info.max]
    for a, b in itertools.product(levels, repeat=2):
        out = normalised_difference(np.array([a], dtype=dtype), np.array([b], dtype=dtype))
        assert np.isnan(out[0]) or -1.0 <= out[0] <= 1.0, (a, b, out)

Nine input combinations at zero, mid-range and the dtype maximum exercise every overflow path a two-band index has, and the test runs in milliseconds. Adding it for every index function in a pipeline is the cheapest insurance available against this entire class of bug.

Common Errors

NDVI around 0.95 over bare soil

The numerator wrapped from negative to a large positive number. Cast before subtracting.

A standard deviation raster that looks like noise

Values were squared in their integer dtype. Cast to float32 before any focal or variance calculation.

Sums over many bands come out small

The accumulator overflowed. Accumulate in float32, or pass dtype="float32" to the reduction.

xarray arithmetic still overflows

The DataArray holds raw integers. Open with mask_and_scale=True or call .astype("float32") first.

Frequently Asked Questions

Q: Why doesn’t NumPy raise an error on overflow? Because integer array arithmetic in NumPy follows fixed-width machine semantics for speed: the result wraps modulo the dtype’s range. Python integers grow without limit, but array elements do not, and there is no check on every addition.

Q: Is subtraction on uint16 also dangerous? More so. Unsigned subtraction that would go negative wraps to a large positive number, so red minus near infrared over vegetation becomes something near 65,000. The numerator of NDVI is exactly that subtraction.

Q: Does xarray protect against overflow? Only when the data has been decoded to float, for instance by mask_and_scale. If the DataArray still holds raw uint16 values, xarray arithmetic wraps exactly as NumPy does.

Q: Is float16 a safe compromise for arithmetic? No. Its maximum is about 65,500 and it has only three significant digits, so band sums of raw counts overflow to infinity and small differences vanish. Use float32 for computation and reserve compact types for storage.