Packing Float Rasters into int16 with Scale and Offset

Map the product’s physical range onto the int16 range, keep −32768 for nodata, and record the transformation:

import numpy as np

lo, hi = -1.0, 1.0                                 # physical range of NDVI
scale = (hi - lo) / 65534                          # -32767..32767 spans the range
offset = (hi + lo) / 2

packed = np.where(np.isfinite(ndvi), np.round((ndvi - offset) / scale), -32768).astype("int16")
# physical = packed * scale + offset

Half the storage of float32, precision far finer than the measurement, and a file that tells every reader how to unpack it. This page belongs to raster dtypes, scaling and numerical precision in Core Raster Fundamentals & STAC Mapping.

How the Mapping Works

Physical range onto integer codes The physical range from minus one to plus one is mapped linearly onto the integers from minus 32767 to plus 32767. The integer minus 32768 is reserved for nodata and never produced by valid data. Each integer step corresponds to a physical step equal to the scale factor, about 0.00003 in this case. A linear map, with one code held back physical: -1.0 to +1.0 int16: -32767 to +32767 -32768 nodata offset maps to 0 One integer step = one scale unit ≈ 0.00003 in physical terms.

Using −32767 to +32767 rather than the full −32768 to +32767 costs one code and keeps the mapping symmetric around the offset, which makes the arithmetic easier to reason about and guarantees that no valid value ever lands on the nodata code.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Rounding, range checks and casting
rasterio >=1.3.0 Writing int16 with scales and offsets metadata
rioxarray >=0.15 Reading back with automatic unpacking 
pip install "numpy>=1.23" "rasterio>=1.3.0" "rioxarray>=0.15"

Complete Working Example

from dataclasses import dataclass

import numpy as np
import rasterio

NODATA = -32768
QMAX = 32767


@dataclass(frozen=True)
class Packing:
    scale: float
    offset: float

    @classmethod
    def for_range(cls, lo: float, hi: float) -> "Packing":
        if not hi > lo:
            raise ValueError("physical range must have hi > lo")
        return cls(scale=(hi - lo) / (2 * QMAX), offset=(hi + lo) / 2)

    def pack(self, arr: np.ndarray) -> np.ndarray:
        q = np.round((arr - self.offset) / self.scale)
        finite = np.isfinite(q)
        if finite.any() and (np.abs(q[finite]).max() > QMAX):
            raise OverflowError("values fall outside the declared physical range")
        return np.where(finite, q, NODATA).astype("int16")

    def unpack(self, packed: np.ndarray) -> np.ndarray:
        out = packed.astype("float32") * np.float32(self.scale) + np.float32(self.offset)
        out[packed == NODATA] = np.nan
        return out


def write_packed(arr: np.ndarray, profile: dict, path: str, packing: Packing,
                 *, description: str) -> None:
    profile = profile | {"dtype": "int16", "nodata": NODATA, "count": 1,
                         "compress": "zstd", "predictor": 2, "tiled": True}
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(packing.pack(arr), 1)
        dst.scales = (packing.scale,)
        dst.offsets = (packing.offset,)
        dst.set_band_description(1, description)


NDVI_PACKING = Packing.for_range(-1.0, 1.0)          # one packing for the whole product

Declaring one Packing per product, rather than computing one per scene, is the most important design decision here. Every file in the archive then shares a scale and offset, so any code that forgets to unpack still gets values on a consistent — if wrong — scale, and any two files can be compared as integers directly. predictor=2 is chosen because it suits integer data; for floats the predictor would be 3.

Choosing the Range, and What It Costs

Range decides precision NDVI over minus one to one packs to a step of about 0.00003, far finer than needed. Reflectance over minus 0.2 to 1.6 packs to about 0.00003 as well. Elevation over minus 500 to 9000 metres packs to about 0.15 metres, adequate for most DEMs. Biomass with no natural upper bound should stay float32. Step size for common products at int16 product physical range step verdict NDVI -1 to 1 0.00003 pack it reflectance -0.2 to 1.6 0.00003 pack it elevation -500 to 9000 m 0.15 m usually fine biomass unbounded — keep float32 Compare the step with the product's own uncertainty; pack when the step is ten times finer.

The range should be the product’s physical range with a little headroom, not the observed range of one scene. Reflectance can exceed 1.0 over bright cloud and snow and dip below zero after atmospheric correction over water, so a range of −0.2 to 1.6 covers the real distribution without clipping. A range tuned to one scene’s minimum and maximum fails the first time a brighter scene arrives, and the overflow check exists precisely to make that failure loud.

Interaction with Compression

Packing and compression reinforce each other. A float32 array of smoothly varying values has noisy low-order mantissa bits that compress poorly; the same values rounded to integers have none of that noise, and horizontal differencing turns smooth fields into runs of small numbers that DEFLATE and ZSTD compress very well.

In practice a packed int16 NDVI raster with ZSTD and predictor=2 is often a quarter to a sixth of the size of the equivalent float32 file with default compression — the halving from the dtype, and then another factor of two to three from compressibility. Across an archive that is a large saving for a transformation that loses nothing measurable, and it compounds with the choices set out in choosing COG compression: ZSTD vs DEFLATE.

Verification

Round-trip error stays within half a step After packing and unpacking, the difference from the original values is spread evenly between minus half a scale step and plus half a scale step. Any error outside that band indicates a wrong scale, a wrong offset, or a value clipped by the range. Error after pack then unpack uniform within ±half step -scale/2 +scale/2 Anything outside the dashed lines means the packing is wrong, not just imprecise. 
import numpy as np
import rioxarray

original = ndvi                                        # float32 with NaN
back = rioxarray.open_rasterio("ndvi_packed.tif", mask_and_scale=True).squeeze().values

both = np.isfinite(original) & np.isfinite(back)
err = np.abs(back[both] - original[both])
assert err.max() <= NDVI_PACKING.scale / 2 + 1e-7, f"max error {err.max()}"
assert (np.isnan(back) == np.isnan(original)).all(), "nodata footprint changed"

Reading back through mask_and_scale rather than through the Packing.unpack method is deliberate: it confirms that a third-party reader, using only the metadata written into the file, recovers the same values. That is the property that makes the file safe to hand to anyone.

Common Errors

OverflowError on a new scene

The scene contains values outside the declared physical range. Widen the range for the product, not for the scene.

Readers see the raw integers

They are not applying band metadata. Use mask_and_scale=True, or document that the file must be unpacked.

Different files in one archive have different scales

Packing was computed per scene. Declare one packing per product and use it everywhere.

The packed file is not much smaller

predictor=2 and a strong compressor were omitted. The dtype halves the size; compression does the rest.

Values near the ends of the range lose precision

They do not — the step is uniform across the range. If the extremes look coarse, the range was declared too wide, spreading the codes over values that never occur.

Frequently Asked Questions

Q: How much precision does int16 packing lose? At most half a scale step per value. For NDVI packed with a scale of 0.0001 that is 0.00005 — far below any sensor’s radiometric uncertainty — so nothing that was ever real is lost.

Q: Should the scale come from this scene’s minimum and maximum? No. Derive it from the product’s physical range, the same for every scene. Per-scene scales make each file individually optimal and mutually inconsistent, and an operation that forgets to apply each file’s own scale produces a patchwork.

Q: Do readers apply the scale automatically? Readers that honour band metadata do — rioxarray with mask_and_scale, xarray for CF-encoded Zarr and NetCDF, and GDAL’s unscaling options. Plain rasterio reads return the packed integers, so code reading raw values must apply the scale itself.

Q: Is uint16 ever better than int16 for packing? When the product is strictly non-negative, yes — it gives twice the resolution for the same range, with 0 or 65535 reserved for nodata. For anything that can dip below zero, including atmospherically corrected reflectance, int16 is the safer default.