Computing Vegetation Indices with Rasterio

An index is two bands and a division, which makes it look like the easiest stage in a multispectral pipeline. It is in fact the stage where the most values go silently wrong, because a division is exactly the operation that turns small problems in its inputs into large problems in its output.

A near-zero denominator produces an unbounded value. A NoData sentinel that was not masked produces a number computed from −9999. A shadowed pixel produces a ratio of two small noisy numbers. None of these raises an error, all of them land in the raster, and any statistic computed afterwards is dominated by whichever went furthest wrong.

This page covers computing indices at survey scale in a way that handles those cases explicitly, along with the masking and summarising that turn an index raster into a number an agronomist can act on. It follows the calibration and alignment stages in multispectral, thermal and index mapping pipelines.

Audience and prerequisites. Python 3.10+, rasterio, and a calibrated band-stacked reflectance raster with band descriptions. An index computed from uncalibrated digital numbers is not comparable between flights, whatever else this page does for it.

Prerequisites

Library / tool Minimum version Install command Role
rasterio ≥ 1.3 pip install "rasterio>=1.3" Windowed reads and writes
numpy ≥ 1.24 pip install numpy The arithmetic
shapely ≥ 2.0 pip install shapely Plot polygons
rasterstats ≥ 0.19 pip install rasterstats Zonal summaries (or hand-rolled)
scipy ≥ 1.10 pip install scipy Morphology for masking

Conceptual architecture

A normalised difference index is (A−B)/(A+B)(A-B)/(A+B) for two bands AA and BB. Its virtue is that multiplicative effects common to both bands — illumination, mostly — cancel, which is why indices work at all on imperfectly calibrated data. Its vice is the denominator.

Where both bands are large the denominator is well conditioned and small errors in the inputs produce small errors in the output. Where both are small — deep shadow, water, the NoData border of a warped frame — the denominator approaches zero and the index becomes an amplifier: a one percent error in a band becomes a large error in the ratio, and at the limit the value is unbounded.

The correct response is a floor on the denominator with masking beyond it, rather than clipping the output. Clipping produces a plausible value where no measurement exists; masking says so.

Index error amplification as the denominator approaches zero A curve of index error against the sum of the two bands, for a fixed one percent error in each band. At a band sum of one the index error is about zero point zero one. At zero point one it is about zero point one. At zero point zero two it is about zero point five, and below that it rises without bound. A shaded region below a band sum of zero point zero two is marked as the masked region, and annotations name the pixels that fall there: deep shadow, water and warped-frame borders. masked 0.005 0.02 0.1 0.4 1.0 sum of the two band values index error floor at 0.02 — error still ±0.5 shadow water warp borders A one percent band error becomes a fifty percent index error where the bands are small.

Figure 1 — Why the denominator floor is a physical decision rather than a numerical one.

Step 1: Compute windowed, so survey scale is not a problem

A five-band 30,000 × 30,000 reflectance raster is 18 GB as float32. Rasterio’s windowed reads make computing an index over it a matter of memory discipline rather than machine size.

import numpy as np
import rasterio
from rasterio.windows import Window


def compute_index(stack_path: str, out_path: str, *, band_a: str, band_b: str,
                  min_sum: float = 0.02, block: int = 2048) -> dict:
    """Normalised difference of two named bands, computed window by window.

    Reading by band description rather than index is what makes this safe
    across sensors, and the windowed loop is what makes it work on a raster
    larger than memory.
    """
    with rasterio.open(stack_path) as src:
        names = [d or "" for d in src.descriptions]
        for band in (band_a, band_b):
            if band not in names:
                raise KeyError(f"band {band!r} not in {names}")
        ia, ib = names.index(band_a) + 1, names.index(band_b) + 1

        profile = src.profile
        profile.update(count=1, dtype="float32", nodata=np.nan,
                       compress="deflate", predictor=3, tiled=True,
                       blockxsize=512, blockysize=512)

        masked_px = 0
        total_px = 0
        with rasterio.open(out_path, "w", **profile) as dst:
            for row in range(0, src.height, block):
                for col in range(0, src.width, block):
                    win = Window(col, row,
                                 min(block, src.width - col),
                                 min(block, src.height - row))
                    a = src.read(ia, window=win, masked=True).filled(np.nan)
                    b = src.read(ib, window=win, masked=True).filled(np.nan)

                    total = a + b
                    valid = np.isfinite(total) & (np.abs(total) > min_sum)
                    out = np.full(a.shape, np.nan, dtype="float32")
                    np.divide(a - b, total, out=out, where=valid)

                    masked_px += int(np.count_nonzero(~valid & np.isfinite(total)))
                    total_px += int(np.isfinite(total).sum())
                    dst.write(out, 1, window=win)

            dst.set_band_description(1, f"nd_{band_a}_{band_b}")
            dst.update_tags(1, FORMULA=f"({band_a} - {band_b}) / ({band_a} + {band_b})",
                            MIN_SUM=str(min_sum), UNITS="dimensionless")

    return {"masked_fraction": masked_px / max(total_px, 1),
            "output": out_path}

Recording the formula and the denominator floor as tags is what makes this index comparable with the next one. Two NDVI rasters computed with different floors are different products, and the difference shows up exactly where the crop is sparse.

Step 2: Mask before summarising

A plot polygon contains crop, soil, shadow and often a strip of track. Averaging all of it produces a number that moves with row spacing and sun angle as much as with the crop.

import numpy as np
from scipy import ndimage


def canopy_mask(index: np.ndarray, nir: np.ndarray, *,
                veg_threshold: float = 0.2,
                shadow_percentile: float = 10.0,
                erode_px: int = 2) -> np.ndarray:
    """Canopy pixels that are lit and away from the plot boundary.

    The shadow test uses near-infrared rather than brightness because
    vegetation is bright in NIR under any illumination; a low NIR value in a
    vegetated pixel therefore means shade rather than sparse cover.
    """
    vegetated = np.isfinite(index) & (index > veg_threshold)
    if not vegetated.any():
        return vegetated

    cutoff = np.percentile(nir[vegetated & np.isfinite(nir)], shadow_percentile)
    lit = vegetated & (nir > cutoff)
    if erode_px > 0:
        lit = ndimage.binary_erosion(lit, iterations=erode_px)
    return lit

Each of the three masks removes a different confounder, and the erosion is the one most often skipped — it removes the mixed pixels at plot edges where any residual band misalignment concentrates.

Step 3: Summarise robustly

import numpy as np


def plot_summary(index: np.ndarray, mask: np.ndarray,
                 *, min_pixels: int = 50) -> dict:
    """Robust summary statistics for one plot.

    Median and interquartile range rather than mean and standard deviation:
    index distributions within a plot are routinely skewed by a patch of poor
    establishment, and the median tracks the bulk of the crop while the
    spread describes the variability agronomists act on.
    """
    values = index[mask & np.isfinite(index)]
    if values.size < min_pixels:
        return {"pixels": int(values.size), "usable": False,
                "note": "too few valid canopy pixels to summarise"}

    q25, q50, q75 = np.percentile(values, [25, 50, 75])
    return {"pixels": int(values.size), "usable": True,
            "median": float(q50), "iqr": float(q75 - q25),
            "p10": float(np.percentile(values, 10)),
            "p90": float(np.percentile(values, 90)),
            "coverage": float(mask.sum() / max(np.isfinite(index).sum(), 1))}

Reporting the coverage — the share of the plot that survived masking — alongside the statistics is what keeps a comparison honest. A plot whose canopy coverage fell from 0.9 to 0.4 between flights has changed in a way the median index does not describe, and reporting both makes that visible.

What each mask removes from a plot, and how the summary changes A plot's pixel population broken down by what each mask removes. Of one hundred percent, bare soil between rows accounts for twenty-two percent, shadowed canopy for eleven percent, and the boundary buffer for six percent, leaving sixty-one percent as lit canopy. Beside it, the plot median index is shown for each successive masking stage, rising from zero point four one unmasked to zero point six eight after all three masks, with the interquartile range narrowing from zero point three one to zero point zero nine. bare soil 22 % shadow 11 % boundary 6 % lit canopy 61 % none −soil −shadow −edge 0.41 0.68 plot median index The unmasked median measures the plot's geometry as much as its crop.

Figure 2 — What masking removes, and what it does to the number the agronomist reads.

Step 4: Attach the plot geometry without rasterising twice

Per-plot statistics need the index values inside each polygon. The naive implementation rasterises the polygon over the whole raster extent once per plot, which on a trial with four hundred plots is four hundred full-extent rasterisations of a grid that is mostly empty.

The efficient version reads only the window each polygon covers and rasterises within it.

import numpy as np
import rasterio
from rasterio.features import geometry_mask
from rasterio.windows import from_bounds


def plot_windows(index_path: str, plots: list[dict]) -> list[dict]:
    """Per-plot index statistics, reading only each plot's own window.

    `plots` carry an id and a shapely geometry in the raster's CRS. Reading a
    window per plot turns an operation that scales with plots times raster
    area into one that scales with the total plot area.
    """
    results = []
    with rasterio.open(index_path) as src:
        for plot in plots:
            geom = plot["geometry"]
            try:
                win = from_bounds(*geom.bounds, transform=src.transform)
            except ValueError:
                results.append({"id": plot["id"], "usable": False,
                                "note": "plot lies outside the raster"})
                continue

            data = src.read(1, window=win, masked=True, boundless=True,
                            fill_value=np.nan).filled(np.nan)
            win_transform = src.window_transform(win)
            inside = ~geometry_mask([geom], out_shape=data.shape,
                                    transform=win_transform, invert=False)

            values = data[inside & np.isfinite(data)]
            results.append({"id": plot["id"], "pixels": int(values.size),
                            "usable": values.size >= 50,
                            "median": float(np.median(values)) if values.size else None})
    return results

boundless=True handles plots that extend past the raster edge, which happens on trial sites where the flight boundary was drawn tight. Without it, those plots raise rather than returning the partial coverage they actually have — and the partial coverage, correctly reported, is more useful than an exception.

Step 5: Compare flights, not rasters

The output a programme actually consumes is a time series per plot, and assembling it correctly needs two things the per-flight statistics do not carry.

A consistent plot geometry. Re-digitising plot boundaries between flights introduces a change of a few percent in every plot mean, indistinguishable from a real effect. Store the polygons once with the trial and reuse them.

A comparability check per pair of flights. Two plot medians are comparable only if the index was computed the same way, the masks used the same thresholds and the calibration was consistent. Those facts live in the tags written earlier; checking them is mechanical.

import rasterio


COMPARABILITY_TAGS = ("FORMULA", "MIN_SUM")


def flights_comparable(path_a: str, path_b: str,
                       masks_a: dict, masks_b: dict) -> dict:
    """Are two index rasters and their masking policies like-for-like?"""
    def tags(path):
        with rasterio.open(path) as src:
            return src.tags(1)

    ta, tb = tags(path_a), tags(path_b)
    problems = []
    for key in COMPARABILITY_TAGS:
        if ta.get(key) != tb.get(key):
            problems.append(f"{key} differs: {ta.get(key)!r} vs {tb.get(key)!r}")
    for key in set(masks_a) | set(masks_b):
        if masks_a.get(key) != masks_b.get(key):
            problems.append(f"mask {key} differs: {masks_a.get(key)} vs {masks_b.get(key)}")

    return {"comparable": not problems, "problems": problems}

A series assembled without that check will, sooner or later, contain a step caused by somebody adjusting a threshold. The check costs nothing and makes the step impossible to introduce silently.

Step 6: Decide what the index is being asked to do

Three quite different uses get called “computing an index”, and they have different tolerances.

Relative ranking within one flight — which plot is doing best today — is the least demanding. Calibration errors common to the whole flight cancel in the ranking, so even an uncalibrated index ranks plots correctly provided the illumination was stable during the flight. This is where most agronomic scouting sits, and it is why uncalibrated indices remain popular despite everything on this page.

Change within a plot between flights is more demanding, because the flights differ. Here calibration matters, masking consistency matters, and the plot geometry must not move. This is the usual monitoring case and the one the pipeline above is built for.

Absolute values compared against a published threshold is the most demanding and the least often achievable. A statement such as “NDVI above 0.75 indicates canopy closure” comes from a study with its own sensor, its own band centres and its own calibration, and a value from a different rig is not directly comparable to it — band centres alone shift an index by several hundredths.

def index_use_requirements(use: str) -> dict:
    """What a given use of an index actually requires of the pipeline."""
    return {
        "ranking": {"calibration": "optional", "masking": "consistent within flight",
                    "geometry": "consistent within flight",
                    "note": "illumination must be stable during the flight"},
        "change": {"calibration": "required", "masking": "identical across flights",
                   "geometry": "fixed across flights",
                   "note": "the usual monitoring case"},
        "absolute": {"calibration": "required", "masking": "identical",
                     "geometry": "fixed",
                     "note": "band centres must match the reference study; "
                             "usually they do not"},
    }[use]

Being explicit about which of the three a deliverable serves prevents the most common disappointment in multispectral work: a carefully calibrated programme whose results are then compared against a literature threshold measured on a different sensor.

Parameter deep-dive

Parameter Type Default Valid range Effect
min_sum float 0.02 0.005–0.10 Denominator floor; below it the index is an amplifier
veg_threshold float 0.2 0.1–0.35 Separates canopy from soil; crop and stage dependent
shadow_percentile float 10.0 0–25 Share of the darkest canopy treated as shaded
erode_px int 2 0–5 Removes mixed pixels at the plot boundary
min_pixels int 50 20–500 Below this the summary is not stable
block int 2048 512–8192 Windowed read size; memory against overhead
Band pair names red/nir — Must be read by name, never by index
Statistic enum median median/mean Median is robust to a poor patch

Verification and output inspection

import numpy as np
import rasterio


def verify_index_raster(path: str, *, max_masked: float = 0.25) -> dict:
    """Checks an index raster should pass before anybody summarises it."""
    with rasterio.open(path) as src:
        arr = src.read(1, masked=True).filled(np.nan)
        tags = src.tags(1)

    finite = arr[np.isfinite(arr)]
    problems = []
    if finite.size == 0:
        problems.append("no finite values at all")
        return {"problems": problems}

    masked_fraction = 1 - finite.size / arr.size
    if masked_fraction > max_masked:
        problems.append(f"{masked_fraction:.0%} of the raster is masked")
    if np.any(np.abs(finite) > 1.0001):
        problems.append("values outside [-1, 1] — the denominator floor is too low")
    if "FORMULA" not in tags:
        problems.append("no formula recorded — this raster is not comparable")
    if "MIN_SUM" not in tags:
        problems.append("no denominator floor recorded")

    return {"masked_fraction": masked_fraction,
            "median": float(np.median(finite)),
            "p01": float(np.percentile(finite, 1)),
            "p99": float(np.percentile(finite, 99)),
            "problems": problems}

The out-of-range check is the definitive one for a normalised difference: the formula cannot produce a value outside [−1, 1] unless a band was negative, which means the calibration produced negative reflectance somewhere and that is a calibration problem rather than an index problem.

What an index value can and cannot support Three rows. Within one calibrated flight, an index supports relative comparison between areas, ranking of plots and locating anomalies, which is what most agronomic use actually needs. Across flights of the same sensor with consistent calibration, it supports trend over time, provided the calibration method did not change between them. Across sensors it supports almost nothing without a cross-calibration, because band centres and widths differ and two sensors' bands of the same name measure different things. within one flight relative comparison, ranking, anomaly location — most agronomic use across flights, same sensor trend over time, provided the calibration method did not change across sensors almost nothing without cross-calibration — band centres and widths differ Most disappointed index users were making a comparison the third row does not support.

Figure 3 — Three scopes, and only the first is free.

Troubleshooting

The index raster is full of NaN. Either the bands were read as a sentinel that was not masked, or the denominator floor is far too high. Check the band statistics before the index.

Values outside −1 to 1. A band is negative. That comes from the dark-current subtraction in radiometric calibration in Python; clip at zero there rather than here.

Plot medians move between flights with no agronomic cause. Compare canopy coverage as well as the median. A change in what survived masking looks like a change in the crop.

Extreme values ring every field boundary. Band misregistration, not index arithmetic. See fixing band misregistration artifacts in index rasters.

Two indices from the same data disagree about which plot is best. Expected, and informative: different indices saturate at different canopy densities. See choosing between NDVI, NDRE and GNDVI.

Computation is slow on a large raster. Check that the source is internally tiled; a windowed read of a striped raster reads whole scanlines, so a small window costs almost as much as the full row.

A plot returns no usable pixels. Either it lies outside the flight extent, or masking removed everything — which on a plot of bare soil is the correct answer rather than a fault.

← Multispectral, Thermal & Index Mapping Pipelines