Converting int16 Reflectance to Float Safely

Mask on the raw integers, cast, then scale and offset:

import numpy as np

raw = src.read(1)                                   # int16 / uint16 digital numbers
mask = raw == src.nodata                            # compare RAW values, not scaled ones
refl = raw.astype("float32") * np.float32(1e-4) + np.float32(offset)
refl[mask] = np.nan

That ordering — mask, cast, scale — is the whole technique. Every common error in this conversion is a violation of it. This page belongs to raster dtypes, scaling and numerical precision in Core Raster Fundamentals & STAC Mapping.

The Three Things That Go Wrong

Double scaling, missed offset, scaled nodata Applying the scale twice produces reflectance four orders of magnitude too small. Missing the processing-baseline offset shifts every value by 0.1, which moves indices in a way that looks like real change between old and new scenes. Scaling the nodata value along with the data turns minus 32768 into minus 3.2768, a plausible-looking number that slips past masks. Each error produces a raster that looks almost right double scaling reader applies scale, then code applies it again 0.00002 instead of 0.2 ratios survive; absolute thresholds all fail missed offset newer scenes carry -0.1, older ones carry none every value +0.1 too high shows up as false change at the baseline switch scaled nodata sentinel scaled with data then compared to itself -3.2768 slips past masks drags means and percentiles downward None raises an error. All three pass a dtype check and plot as a plausible image.

The missed offset is the most expensive of the three in practice, because it hides inside time series. An archive spanning the processing-baseline change shows an apparent jump in every index at the switch date — a jump that looks exactly like a regional environmental event and has sent more than one analysis in the wrong direction.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Raw reads, nodata and embedded scales
numpy >=1.23 Explicit casting and masking
pystac >=1.10 Reading scale and offset from the catalog item when the file lacks them 
pip install "rasterio>=1.3.0" "numpy>=1.23" "pystac>=1.10"

Complete Working Example

from dataclasses import dataclass

import numpy as np
import rasterio


@dataclass(frozen=True)
class Radiometry:
    scale: float
    offset: float
    nodata: int | None
    source: str                      # where the numbers came from, for the log


def radiometry_for(path: str, item=None, asset_key: str | None = None) -> Radiometry:
    """Prefer the file's own metadata; fall back to the catalog item's raster extension."""
    with rasterio.open(path) as src:
        scale, offset, nodata = src.scales[0], src.offsets[0], src.nodata
    if (scale, offset) != (1.0, 0.0):
        return Radiometry(scale, offset, nodata, "file")

    if item is not None and asset_key:
        bands = item.assets[asset_key].extra_fields.get("raster:bands", [{}])
        b = bands[0]
        if "scale" in b or "offset" in b:
            return Radiometry(b.get("scale", 1.0), b.get("offset", 0.0),
                              b.get("nodata", nodata), "catalog")

    raise ValueError(f"no scale/offset declared for {path}; refusing to guess")


def to_reflectance(path: str, radiometry: Radiometry, band: int = 1) -> np.ndarray:
    with rasterio.open(path) as src:
        raw = src.read(band)

    missing = np.zeros(raw.shape, dtype=bool)
    if radiometry.nodata is not None:
        missing = raw == radiometry.nodata               # 1. mask on RAW values

    out = raw.astype("float32")                          # 2. cast before arithmetic
    out *= np.float32(radiometry.scale)                  # 3. scale
    out += np.float32(radiometry.offset)                 #    then offset
    out[missing] = np.nan
    return out

Refusing to guess is the most important line in the file. A function that silently falls back to a scale of 1.0 when metadata is missing produces reflectance of 2,000 that downstream code may clip, threshold or average without complaint. Failing loudly forces the question to be answered once, where the answer can be recorded — and recording the source of the numbers makes a later audit possible.

Scenes Either Side of a Baseline Change

A false step at the baseline switch Mean red reflectance over a stable area is plotted across two years. Without applying the per-scene offset, values jump upward by about 0.1 on the date the processing baseline changed and stay high afterwards. Applying each scene's own offset removes the step entirely and leaves only the seasonal cycle. Stable target, two treatments baseline change offset ignored: false +0.1 step offset per scene: no step older baseline, no offset newer baseline, offset -0.1 A change detector run on the red series would flag the whole region at the switch date.

The safe pattern is per-scene radiometry: read each scene’s scale and offset from its own metadata, never from a constant defined once for the archive. Catalogs that expose the raster extension make this easy, since every asset carries its own values. The catalog mechanics are covered in querying STAC catalogs programmatically, and the sensor-level details of the correction in comparing L1C and L2A Sentinel-2 products.

One Place to Apply the Scale

Most double-scaling bugs come from two well-meaning layers each applying the conversion. A reader configured with mask_and_scale=True returns physical values, and then a helper written for raw data multiplies by the scale again. Neither piece of code is wrong on its own.

The durable fix is organisational rather than technical: decide that reflectance is converted in exactly one function at the edge of the pipeline, and that everything past that function works in physical units. Name variables to match — raw_b04 versus refl_b04 — so a reviewer can see at a glance which side of the boundary a value is on. When a library is doing the conversion, disable any hand-written scaling rather than trying to detect whether it has already happened; detection by inspecting value ranges is fragile, especially for dark scenes.

Verification

Where converted reflectance should land For a typical land scene the first percentile of converted reflectance sits just below zero and the ninety-ninth below about 0.8. Values in the thousands mean no scaling, values near 0.00001 mean double scaling, and a first percentile far below minus 0.1 means a nodata sentinel leaked into the data. Plausibility bands for the 1st and 99th percentiles -0.05 to 0.8: healthy below -0.1: sentinel leaked above 1.5: unscaled Everything around 1e-5: scaled twice. 
import numpy as np

refl = to_reflectance("B04.tif", radiometry_for("B04.tif", item, "red"))
valid = refl[np.isfinite(refl)]
p1, p99 = np.percentile(valid, [1, 99])
print(f"p1 {p1:.4f}  p99 {p99:.4f}")
assert p1 > -0.1, "sentinel or offset problem at the low end"
assert p99 < 1.5, "values look unscaled"
assert p99 > 0.01, "values look scaled twice"

Three thresholds catch all three failure modes, and they are loose enough that a legitimately unusual scene — mostly water, or mostly snow — still passes. Run them per scene in the ingest step and the class of bugs described on this page cannot reach anything downstream.

Common Errors

Values around 0.00002

The scale was applied twice. Remove the manual scaling or disable mask_and_scale, not both.

An abrupt shift in every index on one date

The processing baseline changed and the offset was not applied per scene. Read the offset from each scene’s metadata.

Very negative outliers in statistics

The nodata sentinel was scaled with the data. Mask on the raw integers before scaling.

Integer results from a conversion

An integer offset was added before casting. Cast to float32 first, then scale, then offset.

Frequently Asked Questions

Q: What scale and offset does Sentinel-2 L2A use? A scale of 1/10000 throughout, and from processing baseline 04.00 onward an additive offset of minus 1000 digital numbers, which becomes minus 0.1 in reflectance. Older scenes have no offset, so the offset must be read per scene rather than assumed for the whole archive.

Q: Why cast to float before multiplying by the scale? Because NumPy follows the input dtype for some operations, and integer arithmetic truncates or overflows. Multiplying an int16 array by 0.0001 does promote to float, but adding an integer offset first, or dividing with an integer, can silently stay integer. Casting first removes the question.

Q: Can negative reflectance be real? Small negative values appear over dark targets such as deep water after atmospheric correction, and the offset in newer products exists precisely so they can be stored. They are noise around zero rather than physics, but clipping them to zero biases dark-target statistics, so keep them unless a downstream step requires non-negative input.

Q: What about Landsat Collection 2? The same pattern with different numbers: a multiplicative scale and an additive offset for surface reflectance, documented per product and exposed in its catalog items. The function above handles it unchanged, because it reads the values rather than assuming them.