Handling nodata and Scale Factors in Band Math

To compute an index correctly, mask the fill values first, convert digital numbers to physical units second, and only then do the arithmetic:

import rioxarray
import numpy as np

# masked=True promotes the declared nodata to NaN before anything else happens
red = rioxarray.open_rasterio("B04.tif", masked=True).squeeze()
nir = rioxarray.open_rasterio("B08.tif", masked=True).squeeze()

red = (red.astype("float32") * 0.0001).where(red.notnull())
nir = (nir.astype("float32") * 0.0001).where(nir.notnull())

ndvi = (nir - red) / (nir + red)     # NaN propagates through, as it should

This is the correctness layer under everything in Band Math Operations with xarray.

Why This Arises in Remote Sensing Workflows

Satellite reflectance products almost never store reflectance. They store scaled integers, because integers compress well and halve the storage of a float archive, and they mark absent pixels with a sentinel value chosen from the same integer range. Both facts are metadata, and both are routinely lost between the archive and the analysis.

The two failures compound. An unmasked fill value of 0 or −9999 enters the arithmetic as a real measurement; an unapplied scale factor leaves values in the thousands where the equations expect fractions. Neither raises an exception, and a normalised difference index is particularly good at hiding them, because dividing two large numbers still produces a number between −1 and 1.

The order matters as much as the operations. Masking is defined in the stored integer space, so it has to happen before any transformation of the values; scaling is defined on valid pixels only, so it has to happen after. Reversing them produces a mask that no longer matches its sentinel, and a fill value that has been quietly multiplied into a plausible-looking reflectance.

Three orderings, three different answers Starting from a stored value of minus 9999 and a valid value of 3140, masking then scaling yields NaN and 0.314. Scaling then masking leaves the sentinel scaled to minus 0.9999 and no longer matching the mask test, so it survives into the index. Skipping the mask entirely produces an index computed from the sentinel, which looks like a valid measurement. stored: nir = 3140, red = −9999 (fill) mask → scale → compute red = NaN nir = 0.3140 NDVI = NaN the pixel is honestly absent scale → mask → compute red = −0.9999 mask test == −9999 fails NDVI = −1.90 out of range, but finite no mask at all red = −9999 nir = 3140 NDVI = −1.91 and it will be averaged in Only the first column produces a value a downstream statistic can use. The other two are detectable by a range assertion — which is why that assertion belongs in the pipeline.

Environment & Setup

Package Version Why
rioxarray ≥0.15 open_rasterio(masked=True), CRS-aware writing
xarray ≥2023.6 Labelled arithmetic and where
rasterio ≥1.3.0 Underlying I/O, scales and offsets
numpy ≥1.23 dtype control and NaN handling 
pip install "rioxarray>=0.15" "xarray>=2023.6" "rasterio>=1.3.0"

Complete Working Example

This function reads a band, resolves its scaling from the file’s own metadata where possible, and returns a masked float array in physical units — the only form the arithmetic should ever see.

Where the scale factor hides, per product family Scaling can live in band attributes that rasterio exposes directly, in free-text tags whose names differ per provider, or nowhere at all — documented only in a product specification. Each case needs a different resolution strategy, and the third needs a project-level constant. Three places a scale factor can live location how to read it risk band scales/offsets src.scales, src.offsets low — travels with the band dataset or band tags src.tags(band) name differs per provider product spec only a constant in your code silently wrong after a baseline change Resolve it once at the read boundary and record what you applied in the output. 
import numpy as np
import rasterio
import rioxarray
import xarray as xr


def read_physical(
    path: str,
    *,
    default_scale: float = 1.0,
    default_offset: float = 0.0,
    band: int = 1,
) -> xr.DataArray:
    """Open a band as masked float32 in physical units.

    Scaling is taken from the file when it declares one, falling back to the
    caller's defaults, and the resolved values are recorded on the output.
    """
    with rasterio.open(path) as src:
        # rasterio exposes per-band scales/offsets when the file declares them
        scale = src.scales[band - 1] if src.scales else default_scale
        offset = src.offsets[band - 1] if src.offsets else default_offset
        tags = src.tags(band)
        # Some providers write the factor as a tag instead of a band attribute
        scale = float(tags.get("scale_factor", scale)) or default_scale
        offset = float(tags.get("add_offset", offset))

    da = rioxarray.open_rasterio(path, masked=True).squeeze(drop=True)  # fill → NaN
    physical = da.astype("float32") * np.float32(scale) + np.float32(offset)

    physical.attrs.update(
        scale_factor_applied=scale,
        add_offset_applied=offset,
        units="reflectance" if scale != 1.0 else da.attrs.get("units", "unknown"),
    )
    physical.rio.write_nodata(np.nan, inplace=True)
    return physical


def normalised_difference(a: xr.DataArray, b: xr.DataArray, *, name: str) -> xr.DataArray:
    """(a − b) / (a + b) with a guard on the zero denominator."""
    denom = a + b
    out = xr.where(denom != 0, (a - b) / denom, np.nan)
    out.name = name
    out.attrs.update(long_name=name, valid_min=-1.0, valid_max=1.0)
    return out.astype("float32")


if __name__ == "__main__":
    red = read_physical("B04.tif", default_scale=0.0001)
    nir = read_physical("B08.tif", default_scale=0.0001)

    ndvi = normalised_difference(nir, red, name="ndvi")

    finite = ndvi.where(np.isfinite(ndvi))
    assert float(finite.min()) >= -1.0001 and float(finite.max()) <= 1.0001, "scaling mismatch"
    print("valid pixels:", int(finite.count()), "of", ndvi.size)

The assertion at the end is the cheapest available detector for every failure described above: a scaling applied twice, applied to one band only, or not applied at all all push the result outside ±1.

Variant Patterns

1. Baselines that add an offset

Some processing baselines introduce an additive offset alongside the multiplier, and a pipeline that handles only the multiplier produces values that are plausibly ranged and consistently wrong.

# Baseline with an offset: physical = (DN + offset) * scale
red = (rioxarray.open_rasterio("B04.tif", masked=True).squeeze(drop=True)
       .astype("float32")
       .pipe(lambda da: (da - 1000.0) * 0.0001))     # offset then scale, per the product spec

Because a normalised difference is invariant to a shared multiplicative scale but not to an additive offset, mixing baselines inside one time series produces a step change in the index that is easily mistaken for a real event. Record the baseline per scene and check it before stacking — the drift-detection argument from Extracting and Parsing Raster Metadata.

2. Combining a fill mask with a quality mask

Fill values and cloud flags are different kinds of invalidity, and both belong in the same mask before arithmetic.

import numpy as np
import rioxarray

scl = rioxarray.open_rasterio("SCL_20m.tif").squeeze(drop=True)
scl10 = scl.rio.reproject_match(red, resampling=rasterio.enums.Resampling.nearest)

invalid_classes = [0, 1, 3, 8, 9, 10, 11]
quality_ok = ~scl10.isin(invalid_classes)

red_clean = red.where(quality_ok)          # NaN where cloudy, shadowed or absent
nir_clean = nir.where(quality_ok)
ndvi = normalised_difference(nir_clean, red_clean, name="ndvi")

Nearest-neighbour is mandatory for the class layer, and applying the mask before the ratio rather than after is what keeps cloud edges out of the result — the argument developed in Masking Clouds with the Sentinel-2 SCL Band.

3. Writing a float result back into an integer file

Choosing the output encoding for an index float32 with NaN fill is the simplest and largest. int16 scaled by ten thousand halves the size and resolves the index to four decimal places, which exceeds the sensor's real precision. uint8 scaled to 0–250 quarters the size but resolves only to about 0.008, which is too coarse for change detection. Encoding an index that lives in [−1, 1] encoding size resolution fill value float32, NaN fill 440 MB / scene exact NaN — unambiguous int16 × 10000 220 MB / scene 0.0001 −32768, outside the range uint8, (x+1)×125 110 MB / scene 0.008 — coarse 255, at the range edge int16 × 10000 is the usual compromise: it out-resolves the sensor and halves the archive. Whatever you choose, write scale_factor and nodata into the file — the next reader has only the file. 
import numpy as np

SCALE, FILL = 10000, -32768

scaled = (ndvi * SCALE).round()
packed = scaled.fillna(FILL).astype("int16")
packed.rio.write_nodata(FILL, inplace=True)
packed.attrs.update(scale_factor=1 / SCALE, long_name="ndvi")
packed.rio.to_raster("ndvi_int16.tif", driver="COG", compress="DEFLATE", blocksize=512)

fillna before astype is essential: casting NaN to an integer dtype is undefined and produces whatever the platform’s conversion happens to yield, which is rarely the fill value you intended.

Common Errors

The index is uniformly close to −1 or +1

One band was scaled and the other was not, so the difference is dominated by the magnitude gap rather than by the spectral contrast. Assert that both bands share a scale factor before the arithmetic.

RuntimeWarning: invalid value encountered in divide

The denominator is zero somewhere — usually a pixel where both bands are zero after masking. Guard with xr.where(denom != 0, …, np.nan) rather than suppressing the warning.

Output opens with values in the tens of thousands

A packed integer result was written without scale_factor metadata, so the reader has no way to know it is scaled. Write the tag, or write float32.

Frequently Asked Questions

Q: Should I mask before or after scaling? Before. The fill value is expressed in the stored integer space, so comparing against it after scaling means comparing against a transformed sentinel — which usually no longer matches exactly, and silently leaves fill pixels in the data.

Q: Why do my NDVI values exceed 1 after applying a scale factor? A normalised difference is scale-invariant when both bands share the same factor, so a value outside ±1 means the two bands were scaled differently, or an additive offset was applied to one and not the other.

Q: How do I write a float result to an integer output? Pick a fill value that cannot occur in the data, replace NaN with it, scale the floats to the integer range, and record both the scale factor and the fill value in the output metadata.