Calculating NDVI and NDRE with Rasterio Step by Step

TL;DR: Open a 5-band multispectral GeoTIFF with rasterio, cast the Red, NIR, and Red Edge bands to float32, apply np.divide(..., where=denominator > 0) to avoid zero-division warnings, clip to [-1, 1], replace NaN with -9999.0, and write two single-band GeoTIFFs that inherit the source CRS and affine transform unchanged.

Context

NDVI (Normalized Difference Vegetation Index) and NDRE (Normalized Difference Red Edge index) are the two most actionable spectral indices for in-season crop monitoring at field scale. NDVI uses the contrast between Red and NIR reflectance to track chlorophyll density; NDRE substitutes Red Edge for Red, making it sensitive to early-stage nitrogen stress before visible symptoms appear — a critical advantage for variable-rate fertiliser decisions.

Without correct band extraction and safe division, the resulting rasters carry silent errors: integer overflow when subtracting uint16 bands, RuntimeWarning floods from zero-denominator pixels, and metadata-stripped outputs that QGIS or downstream variable-rate export to ISOXML tooling cannot georeference. A single mismatched band index shifts NDVI into the range of bare-soil values and makes the entire flight unusable.

This workflow assumes imagery is already orthorectified and radiometrically calibrated to surface reflectance. If you are still working from raw digital numbers (DN), complete the calibration and stitching steps in orthomosaic stitching workflows before proceeding.

Prerequisites

These differ from the parent cluster only in band-count requirement:

  • rasterio ≥ 1.3.9 (ships float32 windowed write fixes; install via conda install -c conda-forge rasterio on macOS ARM64 or Windows to avoid GDAL DLL conflicts)
  • numpy ≥ 1.24
  • Input GeoTIFF: at minimum 5 spectral bands, uint16 or float32 DN/reflectance, with an embedded CRS (EPSG:32632 or similar UTM zone — never assume geographic WGS84/EPSG:4326 for pixel-level operations)
  • Sensor: MicaSense RedEdge-MX, DJI P4 Multispectral, or Parrot Sequoia with known band ordering

Band Ordering by Sensor

Before writing a line of code, confirm which array index maps to which wavelength. The table below covers the three most common field sensors:

Sensor Band 1 Band 2 Band 3 (Red) Band 4 (NIR) Band 5 (Red Edge)
MicaSense RedEdge-MX Blue 475 nm Green 560 nm Red 668 nm NIR 842 nm Red Edge 717 nm
DJI P4 Multispectral Blue 450 nm Green 560 nm Red 650 nm Red Edge 730 nm NIR 840 nm
Parrot Sequoia Green 550 nm Red 660 nm Red Edge 735 nm NIR 790 nm — (4 bands only)

Note that the DJI P4 Multispectral puts Red Edge at index 3 and NIR at index 4 — the reverse of MicaSense. Hardcoding nir = stack[3] without checking the sensor will produce inverted NDRE over every DJI flight.

NDVI and NDRE Calculation — Band Math Diagram

NDVI and NDRE calculation data flow A multispectral GeoTIFF with five bands splits into two paths: Red and NIR feed the NDVI formula, while NIR and Red Edge feed the NDRE formula. Both outputs are written as georeferenced single-band GeoTIFFs. Multispectral GeoTIFF (5 bands, float32) Red · NIR · Red Edge Red, NIR NIR, RedEdge NDVI (NIR − Red) ———————— (NIR + Red) NDRE (NIR − RedEdge) —————————— (NIR + RedEdge) ndvi.tif float32 · EPSG preserved nodata = −9999.0 ndre.tif float32 · EPSG preserved nodata = −9999.0 np.divide(..., where=denominator > 0) prevents zero-division warnings

Step-by-Step Implementation

The four steps below form a single, directly runnable pipeline. Each function validates its output before returning.

Step 1 — Install dependencies

BASH 
conda install -c conda-forge rasterio=1.3.10 numpy=1.26

Step 2 — Load bands and verify sensor band order

PYTHON 
import rasterio
import numpy as np
from pathlib import Path

# Band indices (0-based NumPy) for MicaSense RedEdge-MX
# DJI P4 Multispectral: swap RED_IDX=2, NIR_IDX=4, REDGE_IDX=3
RED_IDX   = 2   # Band 3, ~668 nm
NIR_IDX   = 3   # Band 4, ~842 nm
REDGE_IDX = 4   # Band 5, ~717 nm

def load_bands(src_path: str) -> tuple[np.ndarray, dict]:
    with rasterio.open(src_path) as src:
        # Assert CRS is present and projected (not geographic)
        assert src.crs is not None, "Input raster has no CRS — check your orthomosaic export settings."
        assert src.crs.is_projected, (
            f"Expected a projected CRS (e.g. UTM), got {src.crs}. "
            "Reproject before running band math."
        )
        assert src.count >= 5, (
            f"Expected at least 5 bands, found {src.count}. "
            "Confirm the sensor outputs a band-stacked GeoTIFF."
        )

        stack = src.read().astype(np.float32)   # shape: (bands, height, width)
        meta  = src.meta.copy()

        # Quick sanity check: NIR mean should exceed Red mean over natural scenes
        nir_mean = np.nanmean(stack[NIR_IDX])
        red_mean = np.nanmean(stack[RED_IDX])
        if nir_mean <= red_mean:
            raise ValueError(
                f"NIR mean ({nir_mean:.1f}) ≤ Red mean ({red_mean:.1f}). "
                "Band indices are likely swapped — verify against your sensor's spectral table."
            )
        return stack, meta

Step 3 — Compute NDVI and NDRE with safe division

PYTHON 
def calculate_ndvi_ndre(
    stack: np.ndarray,
    red_idx: int   = RED_IDX,
    nir_idx: int   = NIR_IDX,
    redge_idx: int = REDGE_IDX,
) -> tuple[np.ndarray, np.ndarray]:
    red      = stack[red_idx]
    nir      = stack[nir_idx]
    red_edge = stack[redge_idx]

    # Pre-allocate with NaN so uncomputed pixels are never zero
    ndvi = np.full(nir.shape, np.nan, dtype=np.float32)
    ndre = np.full(nir.shape, np.nan, dtype=np.float32)

    denom_ndvi = nir + red
    denom_ndre = nir + red_edge

    # np.divide with `where` avoids RuntimeWarning and leaves masked pixels as NaN
    np.divide(nir - red,      denom_ndvi, out=ndvi, where=denom_ndvi > 0)
    np.divide(nir - red_edge, denom_ndre, out=ndre, where=denom_ndre > 0)

    # Clip to physically valid range (values outside indicate calibration drift)
    np.clip(ndvi, -1.0, 1.0, out=ndvi)
    np.clip(ndre, -1.0, 1.0, out=ndre)

    # Inline verification: valid pixel fraction must exceed 50%
    valid_ndvi = np.sum(~np.isnan(ndvi))
    total      = ndvi.size
    assert valid_ndvi / total > 0.50, (
        f"Only {valid_ndvi/total:.1%} of NDVI pixels are valid — "
        "check for nodata extent or incorrect band masking."
    )
    return ndvi, ndre

Step 4 — Export georeferenced GeoTIFFs

PYTHON 
def export_indices(
    ndvi: np.ndarray,
    ndre: np.ndarray,
    meta: dict,
    out_dir: str,
    nodata: float = -9999.0,
) -> None:
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)

    # Update only the fields that change; CRS and transform are inherited
    meta.update({
        "count":    1,
        "dtype":    "float32",
        "compress": "deflate",
        "predictor": 2,        # floating-point delta predictor
        "tiled":    True,
        "blockxsize": 256,
        "blockysize": 256,
        "nodata":   nodata,
    })

    for filename, data in {"ndvi.tif": ndvi, "ndre.tif": ndre}.items():
        # Replace NaN with the explicit nodata sentinel before writing
        export_data = np.where(np.isnan(data), nodata, data).astype(np.float32)
        out_path = out_dir / filename
        with rasterio.open(out_path, "w", **meta) as dst:
            dst.write(export_data, 1)
        print(f"Written: {out_path}  ({export_data[export_data != nodata].mean():.4f} mean)")

# --- Run the pipeline ---
if __name__ == "__main__":
    stack, meta = load_bands("farm_north_reflectance.tif")
    ndvi, ndre  = calculate_ndvi_ndre(stack)
    export_indices(ndvi, ndre, meta, out_dir="indices/")

Inline verification: after the export, confirm outputs are correctly georeferenced:

PYTHON 
with rasterio.open("indices/ndvi.tif") as chk:
    print(chk.crs)        # should match source EPSG
    print(chk.nodata)     # should be -9999.0
    arr = chk.read(1)
    valid = arr[arr != -9999.0]
    print(f"NDVI range: {valid.min():.3f} – {valid.max():.3f}")
    assert -1.0 <= valid.min() and valid.max() <= 1.0, "NDVI out of valid range."

Gotchas and Edge Cases

  • Band index is 1-based in rasterio metadata, 0-based in the returned NumPy array. src.read(4) (1-based) returns the same data as src.read()[3] (0-based array). Mixing the two conventions in the same function is the most common source of band-swap bugs.
  • DJI P4 Multispectral puts Red Edge at Band 4 and NIR at Band 5 — the reverse of MicaSense RedEdge-MX. Always override NIR_IDX and REDGE_IDX constants at the top of your script rather than hardcoding them inside functions.
  • np.nan as nodata is not reliably recognised by GDAL-backed tools. QGIS, ArcGIS Pro, and most field-equipment platforms (John Deere Operations Center, Climate FieldView) expect a numeric nodata tag in the GeoTIFF header. Write NaN during intermediate calculations, then substitute -9999.0 before the final rasterio.open write call.
  • uint16 source data subtracted before float cast causes integer underflow. If Red > NIR (which can happen in water or deep shadow), uint16(Red) - uint16(NIR) wraps around to ~65000. Always cast with .astype(np.float32) before the subtraction step.

Frequently Asked Questions

Why does my NDVI output contain only zeros or NaN across the entire raster?

This almost always means the wrong bands were assigned. Confirm by printing stack[NIR_IDX].mean() — for reflectance-calibrated MicaSense imagery it should be roughly 0.4–0.6 over a green crop canopy. If it reads near zero, NIR and Red are swapped.

Should I use np.nan or -9999.0 as the nodata sentinel?

Use np.nan during in-memory computation so masked pixels never pollute arithmetic. Replace with -9999.0 before writing to disk. Many GDAL-backed tools do not recognise IEEE NaN as a GeoTIFF nodata value unless TIFFTAG_GDAL_NODATA is explicitly set, and numeric sentinels are universally portable across field platforms.

My NDRE values look inverted — high over bare soil, low over crops.

NIR and Red Edge are swapped. Assert np.nanmean(stack[NIR_IDX]) > np.nanmean(stack[REDGE_IDX]) at load time; if it fails, swap the index constants. This check is already included in the load_bands function above (comparing NIR vs Red, which provides the same directional guard for NDVI).

This guide is part of Band Math & Raster Algebra in Python — see there for the full pipeline context including radiometric calibration, tiled/windowed processing for large orthomosaics, and chaining multiple indices in a single I/O pass.

Parent section Band Math & Raster Algebra in Python Related Python Script for SAVI Calculation on Drone Tiles