Raster Dtypes, Scaling and Numerical Precision

The most insidious errors in raster processing are numerical, because they do not raise. Two uint16 bands added together wrap around silently when the sum exceeds 65,535. A scale factor applied twice produces reflectance of 0.00004 that looks like a very dark pixel. A nodata value of zero survives into a ratio as a genuine zero reflectance. Each produces a raster that loads, plots and passes a schema check, and is wrong. This topic — part of Core Raster Fundamentals & STAC Mapping — is about making those errors impossible rather than merely unlikely.

The organising rule is simple: store compact, compute in float, and record the conversion between them in the file. Almost every numerical defect in a raster pipeline is a violation of one of those three clauses.

---

Prerequisites

pip install "rasterio>=1.3.0" "numpy>=1.23" "xarray>=2024.1" "rioxarray>=0.15"

The metadata inspection this builds on is covered in extracting nodata and dtype from a GeoTIFF, and the band-math context in handling nodata and scale factors in band math.

---

The Dtype Landscape

The pattern that falls out of that table is the one most mature archives use: pixels stored as scaled 16-bit integers, unpacked to float32 for any arithmetic, and repacked for storage. The rest of this topic is about doing each conversion without losing information or inventing it.

---

Step-by-Step Workflow

Step 1 — Read the storage contract before the values

import rasterio


def storage_contract(path: str) -> list[dict]:
    """Dtype, nodata, scale and offset for every band."""
    with rasterio.open(path) as src:
        return [
            {
                "band": i,
                "name": src.descriptions[i - 1],
                "dtype": src.dtypes[i - 1],
                "nodata": src.nodatavals[i - 1],
                "scale": src.scales[i - 1],
                "offset": src.offsets[i - 1],
            }
            for i in range(1, src.count + 1)
        ]


for band in storage_contract("S2A_36NYF_20260614_B04.tif"):
    print(band)

Many products do not set scales and offsets in the file and instead document them in a catalog record or a product specification — Sentinel-2 L2A before its processing-baseline change, most Landsat Collection 2 products, and a great many derived archives. src.scales then reports 1.0, which is technically correct and practically misleading. Where the file is silent, take the values from the catalog item and record where they came from.

Step 2 — Unpack to float with nodata masked first

import numpy as np
import rasterio


def read_physical(path: str, band: int = 1, *, scale: float | None = None,
                  offset: float | None = None) -> np.ndarray:
    """Read a band as float32 physical values with nodata as NaN."""
    with rasterio.open(path) as src:
        raw = src.read(band)
        nodata = src.nodata
        scale = scale if scale is not None else src.scales[band - 1]
        offset = offset if offset is not None else src.offsets[band - 1]

    out = raw.astype("float32")                  # convert BEFORE any arithmetic
    if nodata is not None:
        out[raw == nodata] = np.nan               # mask using the RAW values
    return out * np.float32(scale) + np.float32(offset)

Two ordering rules carry the whole function. Convert to float before applying the scale, or integer arithmetic truncates the result. And compare against nodata using the raw integer array, not the scaled one — after scaling, the nodata value is no longer the number you are looking for, and comparing floats for equality is fragile anyway.

Step 3 — Compute in float32

Once every input is float32 with NaN for missing, band arithmetic is safe by construction: sums cannot overflow, divisions by zero produce inf or nan that can be handled explicitly, and missing data propagates rather than being silently treated as zero.

import numpy as np

red = read_physical("B04.tif", scale=0.0001, offset=-0.1)
nir = read_physical("B08.tif", scale=0.0001, offset=-0.1)

den = nir + red
ndvi = np.divide(nir - red, den, out=np.full_like(den, np.nan), where=den != 0)

The offset=-0.1 in that example is not illustrative decoration: Sentinel-2 L2A products from processing baseline 04.00 onward apply an offset of −1000 digital numbers to reflectance, and ignoring it shifts every index. The version-specific details are covered in converting int16 reflectance to float safely.

Step 4 — Choose a storage dtype from the result’s range and precision

The middle branch is the one worth defaulting to for indices and reflectance. NDVI stored as int16 with a scale of 0.0001 represents every value from −1 to 1 at four decimal places — more precision than any sensor delivers — at half the size of float32. The packing arithmetic and its edge cases are the subject of packing float rasters into int16 with scale and offset.

Step 5 — Pack and record the conversion

import numpy as np
import rasterio


def write_packed(arr: np.ndarray, profile: dict, path: str, *,
                 scale: float = 0.0001, offset: float = 0.0,
                 nodata: int = -32768, description: str = "") -> None:
    packed = np.round((arr - offset) / scale)
    packed = np.where(np.isfinite(arr), packed, nodata)
    if np.nanmax(np.abs(packed[packed != nodata])) > 32767:
        raise OverflowError("values do not fit int16 at this scale")

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

Setting dst.scales and dst.offsets writes them into the GeoTIFF’s metadata, where rasterio, GDAL and rioxarray’s mask_and_scale will all find and apply them. A file that carries its own scale is self-describing; a file whose scale lives only in a README is a trap for the next person.

---

Nodata Is Part of the Dtype Decision

The nodata value has to be representable in the dtype, outside the range of valid data, and preserved by every operation between read and write. Those three constraints interact in ways that catch people repeatedly.

The conversion at the boundary is the whole discipline. Read a sentinel-coded integer band, turn the sentinel into NaN as part of the float conversion, compute with NaN throughout, and turn NaN back into the sentinel only when packing for storage. Any code path that does arithmetic on the raw sentinel-coded integers is a bug waiting for a scene with enough nodata to matter. The detailed choice of sentinel is covered in choosing a nodata value that survives band math.

---

Letting the Libraries Do the Unpacking

Everything above can be done by hand, and it is worth understanding at that level. In day-to-day work, though, the safer route is to let a library apply the storage contract consistently, because hand-rolled unpacking is exactly where scale factors get applied twice.

import rioxarray

# mask_and_scale applies nodata -> NaN, then scale and offset, returning float
red = rioxarray.open_rasterio("B04.tif", mask_and_scale=True).squeeze("band")
nir = rioxarray.open_rasterio("B08.tif", mask_and_scale=True).squeeze("band")
ndvi = (nir - red) / (nir + red)

mask_and_scale=True reads the file’s own nodata, scale and offset and applies them in the right order, returning float values with NaN for missing data. It does exactly what read_physical above does, with one important caveat: it can only apply what the file declares. For products that document their scale outside the file, it will silently apply 1.0 — correctly according to the metadata and wrongly according to the product.

The robust pattern is therefore to decide, per product family, where the contract lives. If the file carries it, use mask_and_scale and nothing else. If it does not, write the scale and offset into your own copy of the file once, at ingest, so that everything downstream can rely on the file again. Mixing the two — some code applying the scale by hand, some reading it from the file — is how a pipeline ends up with a scale applied twice in one branch and not at all in another.

The same principle applies to xarray’s CF conventions when reading Zarr or NetCDF: the scale_factor, add_offset and _FillValue attributes are applied automatically on read and reversed on write. That is what makes the scaled-int16 encoding in writing an xarray datacube to Zarr invisible to the code that uses the cube.

---

Precision Budgets in Practice

It helps to be concrete about how much precision each quantity actually needs, because that decides whether a compact dtype loses anything.

Surface reflectance from a well-calibrated sensor has radiometric uncertainty of a few thousandths. Storing it at a resolution of 0.0001 — which scaled int16 does — keeps a decimal place more than the measurement supports. NDVI inherits that uncertainty and is similarly well served. Elevation from a global DEM has vertical errors of metres; storing it to a centimetre in int16 with a scale of 0.01 covers −327 to +327 m, which is too narrow for most terrain, so elevation typically uses int16 in whole metres or float32. Derived continuous products such as biomass or canopy height have uncertainties of tens of per cent and no natural bound, and float32 is the honest choice there.

The general test is to compare the quantisation step with the measurement’s own uncertainty. When the step is ten times finer than the uncertainty, the compact dtype loses nothing that was ever real. When it approaches the uncertainty, use a wider type. The storage-side payoff of getting this right compounds with compression, as the size comparisons in choosing COG compression: ZSTD vs DEFLATE show.

---

Parameter Reference

---

Verification & Testing

import numpy as np

phys = read_physical("B04.tif", scale=0.0001, offset=-0.1)
valid = phys[np.isfinite(phys)]

assert valid.size > 0, "every pixel is nodata — check the nodata value"
assert -0.2 < np.percentile(valid, 1) and np.percentile(valid, 99) < 1.5, \
    "reflectance outside plausible range — scale or offset wrong"
assert phys.dtype == np.float32
print(f"median reflectance {np.median(valid):.4f}")

The plausibility range is the check that catches most scaling bugs. Surface reflectance lives between roughly 0 and 1, with a little overshoot from cloud and snow; values in the thousands mean the scale was never applied, values around 0.00001 mean it was applied twice, and a median shifted by exactly 0.1 means an offset was missed.

A second check worth running once per product type: read the same pixel through two independent paths — plain rasterio with manual scaling, and rioxarray with mask_and_scale=True — and confirm they agree. Disagreement means one path is reading metadata the other ignores.

---

Troubleshooting

Index values like 3.2 or −40 appear in an NDVI raster

Integer arithmetic overflowed or truncated. Convert to float32 before any band arithmetic; the full explanation is in avoiding integer overflow in index calculations.

Reflectance values are around 0.00001

The scale factor was applied twice — once by a reader with mask_and_scale and once by hand. Pick one place to apply it.

Indices are systematically shifted between two dates

One scene has an offset in its processing baseline and the other does not. Apply the offset per scene from its own metadata.

Masked pixels reappear after writing

NaN was cast straight to an integer dtype, which produces an arbitrary value. Replace NaN with the sentinel before casting.

Statistics include values near −3276.8

The sentinel was scaled along with valid data. Mask the raw sentinel before applying the scale.

---

Frequently Asked Questions

Q: Why do satellite products store reflectance as integers? Because reflectance needs about four significant digits and an int16 or uint16 holds that exactly with a scale factor, at half the size of float32. Across a petabyte archive that halving is worth the extra step of applying the scale on read.

Q: What is the safest dtype for doing arithmetic? float32 for almost everything. It represents scaled reflectance and indices with ample precision, never overflows on band sums, and propagates NaN for missing data. float64 doubles memory for precision no raster product needs.

Q: Why is my index full of strange large numbers? Integer overflow. Adding two uint16 bands whose sum exceeds 65,535 wraps around to a small number, and the resulting ratio is wrong without any error. Convert to float32 before any arithmetic.

Q: Is float16 ever a good idea? For storage of display renditions or model inputs where three significant digits suffice, occasionally. It saves half again over int16-with-scale’s equivalent float, but support across GIS tools is patchy and it cannot hold values much above 65,000, so scaled int16 is the safer compact choice.

---

Related

• Converting int16 Reflectance to Float Safely — scales, offsets and processing baselines.
• Choosing a Nodata Value That Survives Band Math — sentinels, NaN and the conversion boundary.
• Packing Float Rasters into int16 with Scale and Offset — compact storage without precision loss.
• Avoiding Integer Overflow in Index Calculations — the silent wrap-around and how to prevent it.
• Band Math Operations with xarray — the arithmetic these conventions protect.