Band Math & Raster Algebra in Python

Band math and raster algebra turn aligned multispectral bands — raw arrays of digital numbers — into agronomically interpretable metrics: NDVI canopy maps, NDRE early-stress detectors, SAVI soil-adjusted indices, and custom linear combinations for specific crop and soil conditions. This guide covers the complete production workflow, from pre-flight radiometric prerequisites through windowed pixel-wise computation to validated GeoTIFF output, producing a single-band reflectance index raster ready for threshold-based crop health mapping or input to a temporal vegetation index stack.

Prerequisites

Install the core stack:

pip install "rasterio>=1.3.9" "numpy>=1.24" "xarray>=2023.6" rioxarray "scikit-image>=0.21"

On macOS ARM64 or Windows, GDAL binary conflicts are common — use conda install -c conda-forge rasterio instead of pip for rasterio.

1. Concept & Algorithm

Vegetation indices are ratios or normalized differences of spectral bands chosen because plant tissues respond predictably at specific wavelengths. Healthy green vegetation absorbs red light (≈ 670 nm) for photosynthesis and strongly reflects near-infrared (≈ 840 nm) through cell wall scattering. Stressed or senescent tissue shows reduced NIR reflectance and increased red reflectance; bare soil shows both bands at roughly equal reflectance.

NDVI (Normalized Difference Vegetation Index) quantifies this contrast:

NDVI = (NIR − Red) / (NIR + Red)

Values range from −1 to +1. Actively growing canopy typically reads 0.6 – 0.85 on summer cereals; stressed canopy 0.3 – 0.55; bare soil 0.05 – 0.25.

NDRE (Normalized Difference Red Edge) substitutes the red-edge band (≈ 730 nm) for red. Red-edge reflectance is more sensitive to chlorophyll concentration at moderate-to-high LAI values — where NDVI saturates — making NDRE more discriminating in dense maize or soybean canopies mid-season. The MicaSense RedEdge-MX provides a dedicated red-edge band (band 4, ≈ 717 nm) specifically for this purpose.

SAVI (Soil-Adjusted Vegetation Index) adds an empirical soil factor L to reduce background soil reflectance influence in sparse canopy conditions, which are common early in the season or in semi-arid crops:

SAVI = ((NIR − Red) / (NIR + Red + L)) × (1 + L)

where L = 0.5 is the standard default. For a production-ready implementation of SAVI on drone tiles, see Python Script for SAVI Calculation on Drone Tiles.

All three indices share the same raster-algebra pattern: load co-registered bands, mask invalid pixels, apply the formula with safe division, and export with the original CRS and affine transform preserved. The critical agronomic constraint is that reflectance, not raw DN, must be the input — DN values are scene-specific and cannot be meaningfully compared across flight dates or camera units.

Before raw DN can be used here, they must be converted to surface reflectance. The ingesting multispectral drone imagery guide covers band layout and radiometric calibration for MicaSense, DJI, and Parrot sensors; the understanding CRS in precision agriculture guide covers how to confirm your raster is in a projected coordinate system (UTM) rather than geographic degrees before any area-weighted computation.

2. Step-by-Step Implementation

Step 1 — Ingest and validate band alignment

Open source rasters and assert that spatial metadata matches exactly. Any mismatch here will produce geometric artefacts in the computed index — pixels from two different locations get combined.

import rasterio
import numpy as np
from rasterio.crs import CRS

def open_and_validate(nir_path: str, red_path: str):
    """Open two single-band GeoTIFFs and verify they can be used together."""
    nir_src = rasterio.open(nir_path)
    red_src = rasterio.open(red_path)

    assert nir_src.crs == red_src.crs, (
        f"CRS mismatch: {nir_src.crs} vs {red_src.crs}"
    )
    assert nir_src.shape == red_src.shape, (
        f"Shape mismatch: {nir_src.shape} vs {red_src.shape}"
    )
    assert nir_src.transform == red_src.transform, (
        "Affine transform mismatch — bands are not pixel-aligned"
    )
    # Require a projected CRS so downstream area calculations are valid
    assert nir_src.crs.is_projected, (
        f"Expected a projected CRS (e.g. EPSG:32632), got {nir_src.crs}"
    )
    return nir_src, red_src

Step 2 — Retrieve per-band nodata values

MicaSense and DJI P4 Multispectral exports set nodata = 0 by convention; Sentinel-2 Level-2A uses 0 or 65535 depending on the processing baseline. Always read nodata from the file metadata rather than hard-coding it.

def get_nodata(src: rasterio.DatasetReader, fallback: float = 0.0) -> float:
    nd = src.nodata
    return float(nd) if nd is not None else fallback

Step 3 — Build a block iterator for windowed I/O

UAV orthomosaics at 2 cm/px over a 50 ha field exceed 10 GB. Loading the entire raster into RAM triggers MemoryError on standard 16 GB instances. Rasterio’s internal block structure (typically 256×256 or 512×512 tiles for COGs) provides a natural iteration unit.

def iter_windows(src: rasterio.DatasetReader):
    """Yield (window, transform) pairs over the raster's internal block grid."""
    for _, window in src.block_windows(1):
        yield window

Step 4 — Execute pixel-wise algebra with safe division

The following function computes a normalized difference index (configurable to NDVI, NDRE, or any band pair) using windowed reads, NaN-safe division, and atomic writes.

def compute_normalized_difference(
    nir_path: str,
    secondary_path: str,
    output_path: str,
    nir_nodata: float = 0.0,
    sec_nodata: float = 0.0,
) -> str:
    """
    Compute (NIR - B) / (NIR + B) with windowed I/O and safe division.

    Parameters
    ----------
    nir_path       : path to the NIR single-band GeoTIFF
    secondary_path : path to Red, RedEdge, or any target band
    output_path    : destination GeoTIFF path
    nir_nodata     : nodata value for NIR band (default 0.0)
    sec_nodata     : nodata value for secondary band (default 0.0)

    Returns
    -------
    output_path on success
    """
    with rasterio.open(nir_path) as nir_src, \
         rasterio.open(secondary_path) as sec_src:

        # Validate before allocating any memory
        assert nir_src.crs == sec_src.crs, "CRS mismatch"
        assert nir_src.shape == sec_src.shape, "Shape mismatch"
        assert nir_src.transform == sec_src.transform, "Transform mismatch"

        profile = nir_src.profile.copy()
        profile.update(
            dtype="float32",
            count=1,
            compress="deflate",
            predictor=2,
            nodata=np.nan,
        )

        with rasterio.open(output_path, "w", **profile) as dst:
            for window in iter_windows(nir_src):
                nir = nir_src.read(1, window=window).astype("float32")
                sec = sec_src.read(1, window=window).astype("float32")

                # Mask nodata pixels
                valid = (nir != nir_nodata) & (sec != sec_nodata)

                result = np.full_like(nir, np.nan, dtype="float32")

                with np.errstate(divide="ignore", invalid="ignore"):
                    num = nir[valid] - sec[valid]
                    den = nir[valid] + sec[valid]
                    result[valid] = np.where(den != 0, num / den, np.nan)

                dst.write(result, 1, window=window)

    return output_path

Assertion: confirm a non-trivial valid pixel fraction after the first window

# Quick sanity check on window [0,0] before committing to the full raster
with rasterio.open(output_path) as check:
    w0 = next(iter_windows(check))
    sample = check.read(1, window=w0)
    valid_pct = np.sum(~np.isnan(sample)) / sample.size
    assert valid_pct > 0.01, (
        f"< 1 % valid pixels in first window — check nodata value or CRS"
    )

Step 5 — SAVI variant with L factor

SAVI requires a trivial modification to the denominator. Keep the windowed structure identical; only the algebra changes.

def compute_savi(
    nir_path: str,
    red_path: str,
    output_path: str,
    L: float = 0.5,
    nir_nodata: float = 0.0,
    red_nodata: float = 0.0,
) -> str:
    """
    SAVI = ((NIR - Red) / (NIR + Red + L)) * (1 + L)
    L=0.5 is the standard default for moderate vegetation density.
    """
    with rasterio.open(nir_path) as nir_src, \
         rasterio.open(red_path) as red_src:

        assert nir_src.crs == red_src.crs
        assert nir_src.shape == red_src.shape

        profile = nir_src.profile.copy()
        profile.update(dtype="float32", count=1,
                       compress="deflate", predictor=2, nodata=np.nan)

        with rasterio.open(output_path, "w", **profile) as dst:
            for window in iter_windows(nir_src):
                nir = nir_src.read(1, window=window).astype("float32")
                red = red_src.read(1, window=window).astype("float32")

                valid = (nir != nir_nodata) & (red != red_nodata)
                result = np.full_like(nir, np.nan, dtype="float32")

                with np.errstate(divide="ignore", invalid="ignore"):
                    den = nir[valid] + red[valid] + L
                    result[valid] = np.where(
                        den != 0,
                        ((nir[valid] - red[valid]) / den) * (1 + L),
                        np.nan,
                    )

                dst.write(result, 1, window=window)

    return output_path

3. Key Parameters & Tuning

4. Handling Edge Cases & Failure Modes

Zero denominator in normalized differences. Any pixel where NIR + Red = 0 produces a division by zero. This legitimately occurs over deep water, shadows, or sensor-saturated pixels. Using np.errstate(divide='ignore', invalid='ignore') suppresses the warning; the np.where(den != 0, ...) guard ensures those pixels become NaN rather than Inf, which downstream readers handle gracefully.

Nodata value leaking into algebra. If the nodata value (0.0) is within the valid reflectance range — which it is not for reflectance-calibrated data, but can be for raw uint16 DN — a pixel with a legitimate zero reading gets masked incorrectly. Always check the sensor’s actual nodata convention. Some MicaSense RedEdge-MX exports set nodata to 65535 (the uint16 maximum); if your profile says dtype=uint16, nodata=65535 and you cast to float32 before masking, 65535 persists as a very large float and corrupts the algebra. Mask before casting, or mask after casting by comparing to float(src.nodata).

UTM zone boundary crossings. A flight that straddles a UTM zone boundary (e.g. EPSG:32632 / EPSG:32633) will have mismatched transforms if the two bands were exported separately in different zones. Reproject both to a common UTM zone before any alignment check. The understanding CRS in precision agriculture guide covers reprojection with explicit EPSG codes.

Tile seam lines in windowed output. Reading tiles without overlap can cause 1-pixel seam artefacts if any smoothing kernel was applied before export. For morphological operations applied per-tile, add a overlap = 8 border using rasterio.windows.Window with boundless=True and trim the overlap region before writing.

Radiometric inconsistency across flight lines. NDVI values along adjacent flight-line boundaries can differ by 0.05–0.15 if vignetting or BRDF effects were not corrected during orthomosaic stitching. This is not a band-math bug — it is a pre-processing gap in the orthomosaic stitching workflow. Verify that the orthomosaic was built with a radiometric balancing step before running index computation.

Cloud shadows misidentified as crop stress. Cloud shadows reduce both NIR and Red reflectance, often producing NDVI values in the 0.2–0.4 range that mimic moderate crop stress. Apply cloud masking for agricultural imagery before band math, not after. Post-hoc shadow masking based on index thresholds is unreliable because stressed crops occupy the same index range.

5. Verification & Output Validation

A computed index raster is not finished until the output distribution has been confirmed against expected agronomic ranges.

Histogram inspection

import rasterio
import numpy as np

def check_index_histogram(index_path: str, expected_min: float = -0.2,
                           expected_max: float = 0.95):
    """
    Confirm the index distribution falls within agronomically plausible bounds.
    Raises AssertionError if the 2nd or 98th percentile lies outside the
    expected range — a common indicator of nodata leakage or DN input.
    """
    values = []
    with rasterio.open(index_path) as src:
        for window in iter_windows(src):
            tile = src.read(1, window=window).ravel()
            values.append(tile[~np.isnan(tile)])

    all_vals = np.concatenate(values)
    p2, p98 = np.percentile(all_vals, [2, 98])

    assert p2 >= expected_min, (
        f"2nd percentile {p2:.3f} below expected minimum {expected_min} "
        "— possible DN input or nodata leakage"
    )
    assert p98 <= expected_max, (
        f"98th percentile {p98:.3f} above expected maximum {expected_max} "
        "— check for saturated pixels or reflectance > 1.0"
    )
    print(f"Index histogram OK: p2={p2:.3f}, p98={p98:.3f}, n={len(all_vals):,}")
    return p2, p98

Known-value assertion on calibration targets. MicaSense and DJI P4M flights should include a calibration panel of known reflectance (e.g. 0.50 in NIR, 0.10 in Red → expected NDVI ≈ 0.667). Extract the panel pixels from the output index raster and assert the computed value is within ±0.03 of the theoretical value.

CRS and transform round-trip check

with rasterio.open(output_path) as out:
    assert out.crs.to_epsg() is not None, "CRS must be expressible as EPSG code"
    assert out.nodata is not None or np.isnan(out.nodata), "nodata must be set"
    print(f"CRS: EPSG:{out.crs.to_epsg()}, shape: {out.shape}, nodata: {out.nodata}")

Visual spot-check. Load the output in QGIS or rasterio.plot.show() and overlay a field boundary from the field boundary extraction workflow. Visually confirm that: (1) the high-NDVI region spatially aligns with known crop areas, (2) field edges match the boundary polygon, and (3) no rectangular tile artefacts are visible.

6. Integration with the Broader Pipeline

Band math sits at the centre of the drone imagery processing pipeline. Its upstream dependencies and downstream consumers are:

Upstream (must be complete before running band math):

• Multispectral bands co-registered and exported as reflectance-calibrated GeoTIFFs — see ingesting multispectral drone imagery
• Orthomosaic stitched and radiometrically balanced — see orthomosaic stitching workflows
• Cloud and shadow pixels masked out — see cloud masking for agricultural imagery
• Coordinate reference system validated as a projected UTM zone — see understanding CRS in precision agriculture

Downstream (consumes this page’s output):

• Threshold mapping reads the single-band index raster and applies agronomic cut-offs to classify field zones — see threshold mapping for crop health
• Temporal aggregation stacks index rasters from multiple flight dates using xarray to track seasonal trajectories — see temporal aggregation of vegetation indices
• Prescription generation clips the index raster to field boundaries and converts zone polygons to variable-rate ISOXML files — see variable-rate export to ISOXML

The output format matters here. Export index rasters as float32 GeoTIFF with compress=deflate, predictor=2, tiled=True, blockxsize=512, blockysize=512. This produces a Cloud-Optimized GeoTIFF (COG) structure that downstream rasterio reads can access with windowed I/O without full-file downloads — critical when rasters are stored in object storage. For the full pipeline architecture, see Drone Imagery Processing & Vegetation Index Workflows.

Frequently Asked Questions

Why does my NDVI output contain unexpected NaN or Inf values? NaN or Inf in NDVI typically means the denominator (NIR + Red) hit zero — most commonly over water bodies, deep shadows, or pixels where the nodata value was not masked before division. Wrap the division in np.errstate(divide='ignore', invalid='ignore') and pre-mask any pixels equal to the source nodata value.

Do I need to convert DN to reflectance before computing vegetation indices? Yes. Vegetation indices computed directly on raw digital numbers are scene-dependent and cannot be compared across flight dates, times of day, or sensor units. Convert to at-sensor radiance (using the sensor’s gain/offset coefficients) and then to surface reflectance using a panel-based empirical line or atmospheric correction before applying any normalized difference formula.

Why do seam lines appear between tiles in my windowed output raster? Seam lines in windowed output occur when adjacent tiles are read without overlap. Read each window with a 1–2 pixel boundary buffer using boundless=True and discard the overlap before writing. For large orthomosaics, a 16-pixel overlap buffer eliminates seam artefacts caused by edge effects in resampling or smoothing operations applied before export.

Related

• Calculating NDVI and NDRE with Rasterio Step by Step — full walkthrough of NDVI and NDRE computation with band-selection logic for MicaSense and DJI P4M sensors
• Python Script for SAVI Calculation on Drone Tiles — complete tiled SAVI implementation with L-factor tuning guidance
• Cloud Masking for Agricultural Imagery — upstream masking step that must precede band math on satellite imagery
• Threshold Mapping for Crop Health — consumes the index rasters produced here to classify management zones
• Temporal Aggregation of Vegetation Indices — stacks single-date index rasters from this workflow into multi-flight time series

This guide is part of Drone Imagery Processing & Vegetation Index Workflows — see there for the full pipeline context, from orthomosaic ingestion through prescription export.