Building a Band-Stacked GeoTIFF with Rasterio

Somebody opens a five-band reflectance raster from last season and needs the red edge band. Band three looks plausible. It is in fact near-infrared, because the sensor that produced it orders its output by wavelength and the one before it ordered by an internal channel number, and nothing in the file says which convention applied.

A band stack without descriptions is a file whose meaning lives in an email. This page covers writing one that carries its own meaning: band names, wavelengths, units, correct NoData for float reflectance, and compression that does not alter the values. It is the output stage of band alignment and stacking for multispectral sets.

Four decisions that make or break the file

Band identification. GeoTIFF supports a description per band, and almost nothing sets it. Writing the band name there, plus the centre wavelength as a tag, turns positional guessing into a lookup.

Data type and NoData. Reflectance is a float in roughly [0, 1]. Storing it as an integer with a scale factor saves space and reintroduces quantisation the calibration worked to avoid. NoData must be NaN rather than a sentinel, because every plausible sentinel — 0, −1, 255 — is either a valid reflectance or an integer the float band cannot represent exactly.

Compression. Reflectance rasters compress well with DEFLATE or ZSTD and a floating-point predictor. A lossy codec is never acceptable: an index is a difference of two bands, and a codec that alters each band by a fraction of a percent alters their difference by considerably more.

Tiling and overviews. A stacked raster is read band by band and window by window, which is what internal tiling is for. Without it, reading a small region of one band reads entire scanlines of all of them.

A band stack that carries its own meaning A five-band raster drawn as stacked planes. Each plane carries a description, a centre wavelength tag and a units tag: blue at four hundred and seventy-five nanometres, green at five hundred and sixty, red at six hundred and sixty-eight, red edge at seven hundred and seventeen, and near-infrared at eight hundred and forty. A panel to the side lists the file-level properties: float thirty-two data, NaN NoData, DEFLATE compression with predictor three, internal tiling and a band order tag. A note contrasts this with an unlabelled stack whose meaning lives in an email. blue · 475 nm · reflectance green · 560 nm · reflectance red · 668 nm · reflectance red edge · 717 nm · reflectance near-infrared · 840 nm · reflectance file-level properties dtype float32 nodata NaN compress DEFLATE, predictor 3 tiled 512 × 512 BAND_ORDER tag overviews, AVERAGE resampling Everything a reader needs is in the file. The alternative is a file whose meaning lives in an email nobody kept.

Figure 1 — What a self-describing band stack carries.

Minimal reproducible solution

import numpy as np
import rasterio


BAND_WAVELENGTHS_NM = {"blue": 475.0, "green": 560.0, "red": 668.0,
                       "red_edge": 717.0, "nir": 840.0}


def write_stack(bands: dict[str, np.ndarray], out_path: str, *,
                transform, crs, order: list[str] | None = None,
                wavelengths: dict[str, float] | None = None) -> dict:
    """Write a self-describing multi-band reflectance raster.

    Band order is explicit and recorded as a file tag as well as per-band
    descriptions, so a consumer that indexes positionally has something to
    check against rather than an assumption.
    """
    wavelengths = wavelengths or BAND_WAVELENGTHS_NM
    order = order or sorted(bands, key=lambda b: wavelengths[b])
    missing = [b for b in order if b not in bands]
    if missing:
        raise KeyError(f"bands {missing} were requested but not supplied")

    first = bands[order[0]]
    stack = np.stack([bands[b].astype("float32") for b in order])

    profile = {
        "driver": "GTiff", "height": first.shape[0], "width": first.shape[1],
        "count": len(order), "dtype": "float32", "crs": crs,
        "transform": transform, "nodata": np.nan,
        "compress": "deflate", "predictor": 3, "zlevel": 6,
        "tiled": True, "blockxsize": 512, "blockysize": 512,
        "BIGTIFF": "IF_SAFER",
    }

    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(stack)
        for i, band in enumerate(order, start=1):
            dst.set_band_description(i, band)
            dst.update_tags(i, WAVELENGTH_NM=f"{wavelengths[band]:.1f}",
                            UNITS="reflectance", VALID_RANGE="0,1")
        dst.update_tags(BAND_ORDER=",".join(order),
                        PRODUCT="calibrated surface reflectance")
        dst.build_overviews([2, 4, 8, 16], rasterio.enums.Resampling.average)

    return {"bands": order, "path": out_path}

predictor=3 is the floating-point predictor, and it is the difference between a stack that compresses to a third of its size and one that barely compresses at all. It is lossless — the predictor is a reversible transform applied before the entropy coder — so there is no trade-off to weigh.

Reading it back safely

import numpy as np
import rasterio


def read_bands(path: str, names: list[str]) -> dict[str, np.ndarray]:
    """Read named bands, failing loudly if the file does not identify them.

    Never index positionally. A stack written by a different sensor, or by an
    older version of the pipeline, may order its bands differently, and a
    positional read produces an index computed from the wrong pair with no
    error at all.
    """
    with rasterio.open(path) as src:
        descriptions = [d or "" for d in src.descriptions]
        if not any(descriptions):
            raise ValueError(f"{path} has no band descriptions — it cannot be "
                             "read safely by name")
        out = {}
        for name in names:
            if name not in descriptions:
                raise KeyError(f"band {name!r} not among {descriptions}")
            idx = descriptions.index(name) + 1
            out[name] = src.read(idx, masked=True).filled(np.nan)
        return out


def band_metadata(path: str) -> list[dict]:
    """Per-band descriptions and tags, for a manifest or a sanity check."""
    with rasterio.open(path) as src:
        return [{"index": i, "description": src.descriptions[i - 1],
                 **src.tags(i)} for i in range(1, src.count + 1)]
Metadata a band-stacked file must carry to stay usable Four rows. Band descriptions naming each band and its centre wavelength prevent the most common downstream error, which is computing an index from the wrong two bands. A nodata value per band keeps masked areas out of every statistic that follows. The data type and any scale factor record whether values are reflectance, scaled integers or raw digital numbers, which nothing else in the file states. The calibration provenance — panel, irradiance sensor or none — is what makes a value comparable with another flight rather than only with itself. band descriptions name and centre wavelength — prevents an index from the wrong two bands per-band nodata keeps masked areas out of every downstream statistic data type and scale reflectance, scaled integers or raw DN — nothing else states it calibration provenance panel, irradiance sensor or none — what makes it comparable A stack without these is a valid GeoTIFF and an unusable deliverable.

Figure 3 — Four fields that decide whether the file is usable by anyone else.

Edge-case matrix

Situation Consequence Handling
No band descriptions Positional guessing downstream Refuse to read; re-write with descriptions
Integer dtype with a scale Quantisation reintroduced Write float32
NoData of 0 Valid reflectance treated as missing NaN NoData
Lossy compression Index errors larger than band errors DEFLATE or ZSTD only
No predictor File three times larger predictor=3 for float
Not tiled Windowed reads pull whole scanlines Internal tiling
Bands of different shapes Stack fails or silently pads Assert shapes before stacking
Band order differs between flights Positional reads break Read by name; record BAND_ORDER

Verification snippet

import numpy as np
import rasterio


def verify_stack(path: str, expected_bands: list[str]) -> dict:
    """Structural and value checks a delivered reflectance stack must pass."""
    problems = []
    with rasterio.open(path) as src:
        descriptions = [d or "" for d in src.descriptions]
        if descriptions != expected_bands:
            problems.append(f"band order is {descriptions}, expected {expected_bands}")
        if src.dtypes[0] != "float32":
            problems.append(f"dtype is {src.dtypes[0]}, expected float32")
        if src.nodata is not None and not np.isnan(src.nodata):
            problems.append(f"nodata is {src.nodata}, expected NaN")
        if not src.profile.get("tiled"):
            problems.append("raster is not internally tiled")
        if src.overviews(1) == []:
            problems.append("no overviews built")

        sample = src.read(1, masked=True).filled(np.nan)
        finite = sample[np.isfinite(sample)]
        if finite.size:
            out_of_range = float(np.count_nonzero((finite < -0.05) | (finite > 1.1))
                                 / finite.size)
            if out_of_range > 0.01:
                problems.append(f"{out_of_range:.1%} of band 1 is outside [0, 1] — "
                                "check the calibration")
    return {"bands": descriptions, "problems": problems, "ok": not problems}

The value-range check is the one that catches calibration errors at the point where they are still cheap to fix. Reflectance outside [0, 1] beyond a fraction of a percent means something upstream is wrong, and finding it in the writer is much better than finding it in a plot statistic.

File size and index error by compression choice Four compression options compared for a five-band float reflectance raster. Uncompressed is one hundred percent of size with no error. DEFLATE without a predictor is eighty-eight percent with no error. DEFLATE with predictor three is thirty-one percent with no error. A lossy JPEG-style codec is nine percent but introduces a band error of half a percent, which becomes a two point four percent error in a normalised difference index. A note states that the lossy option is never acceptable for a measurement product. no compression 100 % DEFLATE 88 % DEFLATE with predictor 3 31 % · lossless lossy codec 9 % · but 2.4 % index error The predictor is free. The lossy option costs the measurement. A half-percent band error becomes several times that in a normalised difference.

Figure 2 — Why one of these four is the answer and one of them is disqualified.

Stacking indices alongside the bands, or not

A recurring question is whether computed indices belong in the same file as the reflectance bands. There is a good case each way and the decision is worth making deliberately rather than by habit.

Keeping them separate is the cleaner default. Reflectance is a measurement; an index is a derived product with its own definition, version and masking policy, and those change more often than the bands do. A separate index raster can be recomputed and reissued without touching the reflectance product, and its own metadata can record the formula and the mask thresholds that produced it.

Keeping them together suits delivery to a client who will open one file in a desktop GIS and expects everything in it. Where that is the audience, append the indices as additional bands with descriptions such as ndvi and ndre, and tag each with the formula used — the same discipline that makes the reflectance bands readable applies to them.

What does not work is a file containing indices and no record of how they were computed. An NDVI with an unspecified denominator floor, computed from an unspecified band pair, is a raster of numbers rather than a measurement, and it will be compared against next season’s NDVI computed slightly differently.

def index_band_tags(name: str, formula: str, *, mask: dict,
                    source_bands: list[str]) -> dict:
    """Tags that make a derived index band self-describing."""
    return {"PRODUCT": "vegetation index", "INDEX_NAME": name,
            "FORMULA": formula, "SOURCE_BANDS": ",".join(source_bands),
            "MASK": ";".join(f"{k}={v}" for k, v in mask.items()),
            "UNITS": "dimensionless", "VALID_RANGE": "-1,1"}

When to escalate

  • The stack must interoperate with software that cannot read float NaN. Some older packages cannot. Supply a companion integer product with a documented scale factor rather than compromising the primary one.
  • Band descriptions are lost by a downstream tool. Assert them on every read and re-set them on every write; a tool that drops them is a fact to work around rather than to argue with.
  • The file is too large for the delivery channel. Overviews and tiling make range-based access practical; consider serving it rather than shipping it, as in serving COG range requests without a tile server.

← Band Alignment and Stacking for Multispectral Sets