Stacking Spectral Indices as Model Features

Declare each index as a named formula over band names, compute with a guarded division, and write the result with band descriptions:

import numpy as np

def nd(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    """Normalised difference with a guarded denominator."""
    den = a + b
    return np.divide(a - b, den, out=np.full(a.shape, np.nan, "float32"), where=den != 0)

features = {"NDVI": nd(B08, B04), "NDWI": nd(B03, B08), "NDMI": nd(B08, B11)}

The guard is the part that gets skipped and later produces a model that refuses to fit. This page belongs to feature engineering for pixel-based models in Raster Machine Learning & Model Inference.

Why Indices Help a Model That Already Has the Bands

A tree model can in principle learn any ratio from the raw bands, so it is reasonable to ask why an index helps at all. The answer is geometry: a decision tree splits on axis-aligned thresholds, and a ratio is a diagonal boundary in band space. Approximating a diagonal with axis-aligned steps takes many splits and generalises poorly beyond the training range.

A ratio is a diagonal a tree cannot draw In raw band space vegetated and non-vegetated pixels separate along a diagonal through the origin. A decision tree can only cut horizontally or vertically, so it approximates the diagonal with a staircase that mis-classifies pixels near every step. Supplying the ratio as a feature turns the same boundary into a single vertical cut. raw bands: staircase with NDVI: one cut red reflectance near infrared NDVI The staircase also fails outside the training range; the single cut does not.

The same argument explains why adding many indices stops helping. Once a diagonal has been supplied, further ratios of the same bands are largely determined by the ones already present, and the model gains columns without gaining separability.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Reading bands by description and writing the stack
numpy >=1.23 Guarded index arithmetic
xarray >=2023.12 Optional named-dimension handling for larger stacks 
pip install "rasterio>=1.3.0" "numpy>=1.23" "xarray>=2023.12"

Complete Working Example

from dataclasses import dataclass
from typing import Callable

import numpy as np
import rasterio


@dataclass(frozen=True)
class Index:
    name: str
    bands: tuple[str, ...]
    fn: Callable[..., np.ndarray]
    note: str


def _nd(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    den = a + b
    return np.divide(a - b, den, out=np.full(a.shape, np.nan, "float32"), where=den != 0)


INDICES = [
    Index("NDVI", ("B08", "B04"), _nd, "vegetation vigour"),
    Index("NDWI", ("B03", "B08"), _nd, "open water"),
    Index("NDMI", ("B08", "B11"), _nd, "canopy and soil moisture"),
    Index("NDBI", ("B11", "B08"), _nd, "built-up surfaces"),
]


def build_index_stack(src_path: str, dst_path: str,
                      indices: list[Index] = INDICES,
                      scale: float = 10_000.0) -> list[str]:
    with rasterio.open(src_path) as src:
        names = list(src.descriptions)
        missing = {b for ix in indices for b in ix.bands} - set(names)
        if missing:
            raise ValueError(f"source lacks bands {sorted(missing)}; has {names}")

        bands = {n: src.read(names.index(n) + 1).astype("float32") / scale
                 for n in {b for ix in indices for b in ix.bands}}

        profile = src.profile | {
            "count": len(indices), "dtype": "float32", "nodata": np.nan,
            "compress": "zstd", "zstd_level": 3, "tiled": True,
            "blockxsize": 512, "blockysize": 512,
        }
        profile.pop("photometric", None)

        with rasterio.open(dst_path, "w", **profile) as dst:
            for i, ix in enumerate(indices, start=1):
                dst.write(ix.fn(*(bands[b] for b in ix.bands)), i)
                dst.set_band_description(i, ix.name)
                dst.update_tags(i, formula="/".join(ix.bands), note=ix.note)

    return [ix.name for ix in indices]

Reading bands by description rather than by position is the single most valuable habit in this file. A stack built by someone else, or by an earlier version of your own pipeline, may order bands differently, and a positional read produces a beautifully coherent and entirely wrong feature raster.

Pruning Redundant Indices

Correlation between candidate indices A six by six correlation matrix over a sample of pixels shows NDBI and NDMI perfectly anti-correlated because one is the negation of the other, and SAVI strongly correlated with NDVI. Both pairs are candidates for pruning; the remaining indices are weakly correlated and each carries distinct information. Correlation over 50,000 sampled pixels NDVI NDWI NDMI NDBI SAVI NDVI NDWI NDMI NDBI SAVI r below -0.95 — drop one r above 0.9 — probably drop moderate — keep both weak — definitely keep NDBI is the exact negation of NDMI, so one of the pair is pure cost.

The check takes ten lines and is worth running whenever the index list changes.

import numpy as np
import rasterio

with rasterio.open("indices.tif") as src:
    names = list(src.descriptions)
    sample = src.read(window=((0, 1024), (0, 1024))).reshape(src.count, -1)

ok = np.isfinite(sample).all(axis=0)
corr = np.corrcoef(sample[:, ok])
for i in range(len(names)):
    for j in range(i + 1, len(names)):
        if abs(corr[i, j]) > 0.9:
            print(f"{names[i]} and {names[j]}: r = {corr[i, j]:+.3f}")

Two indices with an absolute correlation above 0.95 are the same feature twice. Above 0.9 they are close enough that dropping one usually costs nothing measurable and saves a band of read and compute on every pixel of every scene — a saving that compounds across a whole archive.

The general index definitions, and which ones suit which target, are covered in spectral index calculation pipelines; the point here is only which of them earn a place in a feature matrix.

Verifying the Stack before Training

Three properties are worth asserting every time the stack is rebuilt, because each corresponds to a failure that is invisible once the model has been fitted.

The shape of a healthy index histogram A healthy NDVI histogram is bimodal between minus zero point two and zero point nine, with a soil mode low and a vegetation mode high. A tall spike hard against minus one is the signature of unmasked nodata entering the ratio, and it is obvious in the histogram long before it is obvious in the model. NDVI distribution over one scene -1.0 0.1 0.6 0.95 nodata leaking into the ratio Soil mode near 0.1, vegetation mode near 0.6 — the shape a healthy scene should have. 
import numpy as np
import rasterio

with rasterio.open("indices.tif") as src:
    names = list(src.descriptions)
    arr = src.read()

assert all(names), "every feature band needs a description"
for name, band in zip(names, arr):
    finite = np.isfinite(band)
    assert finite.mean() > 0.3, f"{name} is mostly NaN"
    lo, hi = np.nanpercentile(band, [0.5, 99.5])
    assert -1.05 <= lo and hi <= 1.05, f"{name} outside the valid range: {lo:.2f}-{hi:.2f}"

The percentile bounds rather than the min and max are deliberate: a handful of extreme pixels is normal, and testing on the absolute extremes turns a healthy scene into a failed build. The finite fraction catches a band that failed to resample; the range check catches unscaled inputs.

Common Errors

RuntimeWarning: invalid value encountered in divide

The guarded form was not used, so a zero denominator produced a NaN with a warning rather than a controlled result. Use np.divide with out and where, and treat the warning as a defect rather than noise.

Index values outside the range −1 to 1

The inputs were not scaled to reflectance, so raw integer counts entered the ratio. Divide by the scale factor before computing, as covered in handling nodata and scale factors in band math.

The stack loads with empty band descriptions

set_band_description was called before the band was written, or the file was later rewritten by a tool that drops descriptions. Assert that every description is non-empty as part of the build.

Frequently Asked Questions

Q: How many indices should a feature stack contain? Three to six for most land cover problems. Each index is a ratio of bands the model already has, so beyond a handful they add correlated columns rather than new information, and every one is paid for on every pixel at inference time.

Q: Why does my NDVI contain infinities? Because red and near infrared were both zero, which happens over nodata or deep shadow, and the denominator vanished. Compute with an explicit where clause that returns NaN for a zero denominator rather than relying on the warning to tell you.

Q: Should indices be computed before or after resampling bands to a common grid? After. A ratio of two bands on different grids is meaningless, so every input must be on the reference grid first. Resampling the index instead of the bands also smooths a quantity that is already a ratio, which distorts its extremes.