Choosing a Nodata Value That Survives Band Math

Store with a sentinel the data can never take, compute with NaN, and convert only at the boundaries:

import numpy as np

SENTINEL = -32768                                  # int16 storage

x = raw.astype("float32")
x[raw == SENTINEL] = np.nan                        # read boundary: sentinel -> NaN
result = (a - b) / (a + b)                         # NaN propagates on its own
packed = np.where(np.isfinite(result), np.round(result / 1e-4), SENTINEL).astype("int16")

Nothing in the middle of the pipeline ever sees the sentinel, so nothing in the middle can mistake it for data. This page belongs to raster dtypes, scaling and numerical precision in Core Raster Fundamentals & STAC Mapping.

How a Sentinel Poisons a Computation

One unmasked pixel, two wrong answers A window of nine reflectance pixels contains one nodata sentinel of minus 32768. Averaged unmasked, the window's mean becomes about minus 3640 instead of about 0.2. In a band ratio the sentinel produces a finite value near one that looks like a perfectly ordinary index and is thresholded like one. A single sentinel in a 3x3 window -32768 eight pixels near 0.2 mean, unmasked: about -3640 obvious if anyone looks; invisible inside a larger aggregate ratio, unmasked: a finite value near 1 looks like dense vegetation and passes every range check With NaN instead of the sentinel, both results are NaN for that pixel and correct for the rest.

The ratio case is the dangerous one. When both bands carry the same sentinel, (a − b) / (a + b) becomes 0 / −65536 = 0; when only one does, the result is close to ±1. Either way the output is a finite number inside the valid index range, and no plausibility check can distinguish it from a real pixel.

Environment & Setup

Package Version pin Used for
numpy >=1.23 NaN handling and nan-aware reductions
rasterio >=1.3.0 Reading and declaring nodata per band
xarray >=2024.1 where, fillna and CF _FillValue handling 
pip install "numpy>=1.23" "rasterio>=1.3.0" "xarray>=2024.1"

Complete Working Example

import numpy as np
import rasterio

SENTINELS = {
    "uint8": 0,          # acceptable only when 0 is impossible as data (class maps)
    "uint16": 65535,
    "int16": -32768,
    "int32": -2147483648,
    "float32": np.nan,
}


def read_for_compute(path: str, band: int = 1) -> np.ndarray:
    """float32 array with every declared or conventional sentinel turned into NaN."""
    with rasterio.open(path) as src:
        raw = src.read(band)
        declared = src.nodata
        dtype = src.dtypes[band - 1]

    out = raw.astype("float32")
    sentinel = declared if declared is not None else SENTINELS.get(dtype)
    if sentinel is not None and not (isinstance(sentinel, float) and np.isnan(sentinel)):
        out[raw == sentinel] = np.nan
    return out


def write_from_compute(arr: np.ndarray, profile: dict, path: str, *,
                       dtype: str = "int16", scale: float = 1e-4) -> None:
    """Pack a float array with NaN into an integer dtype with a sentinel."""
    sentinel = SENTINELS[dtype]
    packed = np.round(arr / scale)
    packed = np.where(np.isfinite(packed), packed, sentinel)
    info = np.iinfo(dtype)
    valid = packed != sentinel
    if valid.any() and (packed[valid].min() <= info.min or packed[valid].max() > info.max):
        raise OverflowError("valid values collide with the sentinel or overflow the dtype")

    with rasterio.open(path, "w", **(profile | {"dtype": dtype, "nodata": sentinel,
                                                "count": 1})) as dst:
        dst.write(packed.astype(dtype), 1)
        dst.scales = (scale,)

The overflow check does double duty: it catches values too large for the dtype, and it catches a valid value landing exactly on the sentinel, which would silently turn a real pixel into missing data. The <= on the minimum is deliberate for that reason — with int16 and a sentinel of −32768, the lowest valid packed value must be −32767 or higher.

NaN-Aware Reductions

mean versus nanmean Over a window where two of nine pixels are NaN, numpy mean returns NaN for the whole window, discarding seven good values. nanmean returns the average of the seven valid pixels. The nan-aware family, including nanmedian and nanpercentile, is what makes NaN a practical nodata representation. Two missing pixels out of nine np.mean → nan seven good values discarded np.nanmean → 0.213 average of the seven valid pixels nanmean · nanmedian · nanpercentile · nanstd · nansum · nanmax xarray's reductions skip NaN by default; plain NumPy does not Which one you want depends on whether a partly missing window should produce a value at all.

Choosing between them is a real decision rather than a technicality. np.mean returning NaN for any window with a missing pixel is conservative — it refuses to summarise incomplete data — and is sometimes exactly right. np.nanmean summarises whatever is present, which is usually what a composite or a zonal statistic wants. For those, record how many valid pixels went into each result, because a mean of one pixel and a mean of nine are not equally trustworthy; the same reasoning drives the observation-count feature in building temporal feature stacks from image time series.

Carrying Masks Through Resampling

Resampling is where nodata handling most often breaks after it has been done correctly everywhere else. A bilinear kernel over a region containing a sentinel averages the sentinel into its neighbours, producing a ring of corrupted pixels along every nodata edge; the same kernel over NaN produces NaN in the whole kernel footprint, which erodes the valid region by a pixel.

The robust approach is to pass the nodata value explicitly to the warp — src_nodata and dst_nodata in rasterio’s reproject — so GDAL excludes missing pixels from each kernel and fills uncovered output cells cleanly. That is covered in depth in reprojecting a raster from UTM to WGS84 with rasterio, and it applies to every resampling step in a pipeline, not only reprojection.

Verification

The footprint must survive the round trip The set of missing pixels in the input is compared with the set in the output. For an operation that does not change geometry, they should be identical; any pixel that was missing and is now valid means data was invented, and any pixel that was valid and is now missing means data was lost. Compare footprints, not values input: missing corner output: identical missing corner = A pixel-wise operation must neither create nor destroy missing pixels. 
import numpy as np

a = read_for_compute("B04.tif")
b = read_for_compute("B08.tif")
idx = (b - a) / (b + a)

missing_in = ~np.isfinite(a) | ~np.isfinite(b)
missing_out = ~np.isfinite(idx)
assert (missing_in <= missing_out).all(), "a missing input produced a valid output"
extra = missing_out & ~missing_in
print(f"{extra.sum()} pixels became missing (zero denominators): acceptable")

The one-directional comparison is intentional. A band ratio can legitimately create new missing pixels where the denominator is zero, but it must never turn a missing input into a valid output. That asymmetry is exactly what an unmasked sentinel violates.

Common Errors

A valid dark pixel disappears

Its packed value equals the sentinel. Choose a sentinel outside the packed valid range and check for collisions before writing.

NaN turns into a random integer on write

NaN was cast directly to an integer dtype. Replace it with the sentinel before astype.

The nodata region grows after resampling

NaN inside a kernel spreads to the whole kernel. Pass src_nodata and dst_nodata to the warp instead.

A file reports no nodata, but has obvious fill

The value was never declared. Declare it on a copy at ingest so every reader masks it.

Frequently Asked Questions

Q: Why is zero a bad nodata value? Because zero is a plausible measurement — dark water, deep shadow, and many index values near zero — so masking it removes real data, and failing to mask it treats missing pixels as real dark ones. A value outside the physically possible range removes the ambiguity.

Q: Should NaN be used as nodata on disk? For float32 files it is the cleanest choice, since NaN cannot be confused with a measurement and every nan-aware function skips it. For integer files it is impossible, so a sentinel is needed there, converted to NaN when the data is read for computation.

Q: What happens to nodata in a ratio of two bands? With NaN, the result is NaN wherever either input is missing, which is exactly right. With an unmasked sentinel, the ratio is a finite and entirely meaningless number that will be averaged and thresholded like real data.

Q: What should a class raster use as nodata? Zero, reserved and never assigned to a real class. Class codes start at one, which leaves zero free and makes the convention easy to remember and to check. It is the one place zero is the right choice.