Vegetation Index Selection for Crop Stages

The wrong vegetation index at the wrong growth stage produces a map that looks authoritative and is quietly false. NDVI over a bare-soil V3 corn field reports soil brightness, not vigour. The same NDVI over a closed grain-fill canopy saturates and hides the nitrogen gradient a red-edge index would show plainly. This page is a decision reference: it tells you which index to trust for a given crop, phenological stage, canopy fraction, and sensor, and it ends with a single runnable dispatcher that computes the index you chose while refusing to fabricate bands the platform never captured.

The output is a defensible index choice per stage plus the exact reflectance formula and a validated raster, ready to feed threshold mapping for crop health detection or a variable-rate prescription. This guide is part of the Drone Imagery Processing & Vegetation Index Workflows section — see that overview for the full pipeline from ingest through temporal compositing. For the mechanics of the ratios themselves, the band math and raster algebra guide is the companion to this decision layer.

Prerequisites

Python packages (exact versions tested):

• rasterio>=1.3.9
• numpy>=1.24

Install with:

pip install "rasterio>=1.3.9" "numpy>=1.24"

Input data requirements:

• A multi-band reflectance raster (surface reflectance, not raw digital numbers) in a projected CRS — e.g. EPSG:32615 for UTM zone 15N. A geographic CRS such as EPSG:4326 is acceptable for the pure band math here, but any per-area statistics downstream require metres, so reproject early.
• Known band-to-index mapping. A MicaSense RedEdge-MX delivers Blue 475, Green 560, Red 668, RedEdge 717, and NIR 840 nm; a DJI P4 Multispectral is spectrally similar; a Parrot Sequoia gives Green, Red, RedEdge, NIR; an RGB Phantom 4 gives Blue, Green, Red only. Confirm the ordering with parsing Sentinel-2 vs drone multispectral bands in Python.
• Nodata set in the profile so that background and masked pixels do not enter the ratio.
• Radiometric calibration complete. Comparing an index across dates on uncalibrated data conflates illumination change with crop change; the ingesting multispectral drone imagery guide covers panel calibration.

This page assumes you already know how to open and window a raster. If not, work through the band math and raster algebra guide first.

1. Concept & Algorithm

A vegetation index is a ratio of reflectance in two or more bands, engineered so that its numerator moves with the biophysical property you care about and its denominator normalises out illumination and view geometry. Every index is a compromise: the band pair that is most sensitive to early biomass over bare soil is not the pair that stays responsive under a closed canopy at grain fill. Selection is therefore a function of three variables — growth stage (which sets canopy fraction and the dominant physiological signal), soil exposure (which sets how much background contaminates the pixel), and sensor bands (which sets what is even computable).

NDVI — canopy vigour, saturates at high LAI. The Normalized Difference Vegetation Index, (NIR − Red) / (NIR + Red), is the default vigour proxy. Red reflectance is driven down by chlorophyll absorption; NIR is driven up by leaf mesophyll scattering. NDVI is excellent from early canopy development until closure, then saturates: once leaf area index exceeds roughly 3–4 the red floor is near zero and the ratio flatlines near 0.85–0.9. Past that point NDVI cannot distinguish an adequately fed canopy from a nitrogen-luxury one. NDVI is also biased by soil brightness when cover is incomplete, which is why it is a poor early-season choice.

NDRE — chlorophyll and nitrogen status, mid-to-late season. The Normalized Difference Red Edge index, (NIR − RedEdge) / (NIR + RedEdge), swaps the saturated red band for a red-edge band near 717 nm. Red-edge reflectance keeps responding to leaf chlorophyll well after red has bottomed out, so NDRE holds its dynamic range through grain fill and reveals within-field nitrogen gradients that NDVI has already smeared to a uniform maximum. It requires a red-edge band — MicaSense RedEdge-MX, DJI P4 Multispectral, Parrot Sequoia, and Sentinel-2 have it; RGB drones do not.

SAVI and OSAVI — sparse canopy, high soil exposure. The Soil-Adjusted Vegetation Index, ((NIR − Red) / (NIR + Red + L)) × (1 + L), introduces a soil-adjustment constant L that damps the soil-brightness contribution when bare ground is in the pixel. At L = 0.5 it is tuned for intermediate cover; at L = 0 it collapses to NDVI. OSAVI (Optimised SAVI) fixes L = 0.16, a value empirically optimal across a range of canopy densities, and removes the (1 + L) normalisation. Use SAVI/OSAVI from emergence through partial cover, then hand off to NDVI as rows close.

GNDVI — green-normalised vigour. The Green NDVI, (NIR − Green) / (NIR + Green), substitutes green for red. Green is less strongly absorbed by chlorophyll than red, so GNDVI saturates at a higher LAI than NDVI and correlates more tightly with canopy chlorophyll and nitrogen at moderate density. It is a useful mid-season alternative when a red-edge band is unavailable but a green band is.

CIre — early stress via red-edge chlorophyll. The red-edge Chlorophyll Index, (NIR / RedEdge) − 1, is a non-normalised ratio that is highly sensitive to small chlorophyll changes. Because it is unbounded it has more headroom than NDRE and can surface early stress before it is visible in a normalised index, making it valuable at tillering and for early nitrogen scouting. Its openness is also its weakness: it is noisier and harder to threshold consistently across fields.

Mapping these to the phenological timeline — emergence, tillering, canopy closure, grain fill, senescence — gives the selection matrix in Section 3. The governing rule is simple: soil-adjusted indices while soil shows, NDVI while the canopy is building and unsaturated, red-edge indices once it closes and nitrogen becomes the question.

2. Step-by-Step Implementation

Step 1 — Assign the growth stage and canopy fraction

Before touching pixels, pin down where the crop is. Canopy fraction is the single most predictive variable for index choice, and it tracks the phenological stage. The dictionary below encodes the default shortlist per stage; treat it as the machine-readable form of the selection matrix.

# Canonical stage → candidate index shortlist (primary first).
STAGE_INDEX_SHORTLIST = {
    "emergence":       ["savi", "osavi"],          # bare soil dominates the pixel
    "tillering":       ["savi", "ndvi", "cire"],    # canopy building, soil still visible
    "canopy_closure":  ["ndvi", "ndre"],            # NDVI peaks, about to saturate
    "grain_fill":      ["ndre", "gndvi"],           # NDVI saturated; nitrogen is the question
    "senescence":      ["ndre", "ndvi"],            # red-edge tracks chlorophyll breakdown
}

def shortlist_for_stage(stage: str) -> list:
    if stage not in STAGE_INDEX_SHORTLIST:
        raise KeyError(
            f"Unknown stage '{stage}'. "
            f"Expected one of {list(STAGE_INDEX_SHORTLIST)}"
        )
    return STAGE_INDEX_SHORTLIST[stage]

assert shortlist_for_stage("grain_fill")[0] == "ndre"
print("V-stage shortlist:", shortlist_for_stage("tillering"))

Step 2 — Resolve which indices the sensor can compute

The shortlist is aspirational; the sensor decides what is real. Encode each platform’s band set and intersect it with the bands each index needs, so a red-edge index is never requested from an RGB camera.

# Bands each platform delivers, mapped to logical names used by the indices.
SENSOR_BANDS = {
    "micasense_redege_mx": {"blue", "green", "red", "rededge", "nir"},
    "dji_p4_multispectral": {"blue", "green", "red", "rededge", "nir"},
    "parrot_sequoia":       {"green", "red", "rededge", "nir"},
    "sentinel2":            {"blue", "green", "red", "rededge", "nir"},
    "phantom4_rgb":         {"blue", "green", "red"},
}

INDEX_REQUIRES = {
    "ndvi":  {"red", "nir"},
    "ndre":  {"rededge", "nir"},
    "savi":  {"red", "nir"},
    "osavi": {"red", "nir"},
    "gndvi": {"green", "nir"},
    "cire":  {"rededge", "nir"},
}

def computable_indices(sensor: str, shortlist: list) -> list:
    have = SENSOR_BANDS[sensor]
    return [ix for ix in shortlist if INDEX_REQUIRES[ix] <= have]

# On an RGB drone at grain fill, the red-edge shortlist is empty — fall back.
print(computable_indices("phantom4_rgb", ["ndre", "gndvi"]))   # []
print(computable_indices("micasense_redege_mx", ["ndre", "gndvi"]))  # ['ndre', 'gndvi']

Step 3 — Compute the chosen index with the dispatcher

A single dispatcher keeps the band-availability guard, the L factor, and safe division in one place. It accepts named band arrays and returns a float32 index array with NaN where inputs are nodata. This is the “compute the chosen index” core referenced throughout the section.

import numpy as np

def compute_index(name: str, bands: dict, L: float = 0.5) -> np.ndarray:
    """
    Compute a vegetation index from a dict of reflectance arrays.

    bands : logical band name -> float32 reflectance array in [0, 1],
            e.g. {"red": ..., "nir": ..., "rededge": ..., "green": ...}
    L     : soil-adjustment factor, used by SAVI only.
    """
    name = name.lower()
    required = INDEX_REQUIRES[name]
    missing = required - bands.keys()
    if missing:
        raise ValueError(f"Index '{name}' needs bands {sorted(missing)} not supplied")

    # Promote to float32 once; safe division guards zero denominators.
    b = {k: bands[k].astype("float32") for k in required}

    with np.errstate(invalid="ignore", divide="ignore"):
        if name == "ndvi":
            num, den = b["nir"] - b["red"], b["nir"] + b["red"]
        elif name == "gndvi":
            num, den = b["nir"] - b["green"], b["nir"] + b["green"]
        elif name == "ndre":
            num, den = b["nir"] - b["rededge"], b["nir"] + b["rededge"]
        elif name == "savi":
            num = (b["nir"] - b["red"]) * (1.0 + L)
            den = b["nir"] + b["red"] + L
        elif name == "osavi":
            # OSAVI fixes L = 0.16 and drops the (1 + L) gain
            num, den = b["nir"] - b["red"], b["nir"] + b["red"] + 0.16
        elif name == "cire":
            # Non-normalised chlorophyll index: (NIR / RedEdge) - 1
            num, den = b["nir"], b["rededge"]
            out = np.where(den > 0, num / den, np.nan) - 1.0
            return out.astype("float32")
        else:
            raise ValueError(f"Unsupported index '{name}'")

        out = np.where(den != 0, num / den, np.nan)
    return out.astype("float32")

# Smoke test on synthetic reflectance
rng = np.random.default_rng(0)
demo = {b: rng.uniform(0.02, 0.6, (32, 32)).astype("float32")
        for b in ("blue", "green", "red", "rededge", "nir")}
ndvi = compute_index("ndvi", demo)
assert ndvi.shape == (32, 32) and ndvi.dtype == np.float32
assert np.nanmin(ndvi) >= -1.0 and np.nanmax(ndvi) <= 1.0
print("NDVI range:", round(float(np.nanmin(ndvi)), 3), "to", round(float(np.nanmax(ndvi)), 3))

Step 4 — Drive selection end-to-end from a raster

Wire the pieces together: read the stack, resolve the stage and sensor to a concrete index, compute it, and carry nodata through. This is the pattern you drop into a windowed loop for large orthomosaics, exactly as in the step-by-step NDVI and NDRE calculation with Rasterio guide.

import rasterio

# MicaSense RedEdge-MX 5-band stack, 1-based rasterio band indices
MICASENSE_BAND_INDEX = {"blue": 1, "green": 2, "red": 3, "rededge": 4, "nir": 5}

def index_for_scene(path, stage, sensor="micasense_redege_mx",
                    band_index=None, L=0.5):
    band_index = band_index or MICASENSE_BAND_INDEX
    chosen = computable_indices(sensor, shortlist_for_stage(stage))
    if not chosen:
        raise RuntimeError(
            f"No index in the {stage} shortlist is computable on {sensor}. "
            "Use a red-edge capable sensor or fall back to NDVI."
        )
    name = chosen[0]
    with rasterio.open(path) as src:
        assert src.crs is not None, "Scene has no CRS — reproject before indexing"
        needed = INDEX_REQUIRES[name]
        bands = {b: src.read(band_index[b]).astype("float32") / 10000.0
                 for b in needed}   # /10000 scales DN → reflectance for many products
    arr = compute_index(name, bands, L=L)
    print(f"{stage} on {sensor}: chose {name.upper()} "
          f"(median {np.nanmedian(arr):.3f})")
    return name, arr

3. Key Parameters & Tuning

The selection matrix condenses Sections 1–2 into a lookup. Read down the stage column to the row of the sensor-available index you can compute.

The parameters that most change an index’s behaviour:

For the full soil-adjusted derivation and windowed implementation, see the Python script for SAVI calculation on drone tiles.

4. Edge Cases & Failure Modes

NDVI reported as “healthy” over bare soil. At emergence, a moist dark seedbed pushes NDVI to 0.2–0.35 with no crop present, and a downstream threshold map reads that as sparse vegetation. Any stage below ~15% cover must use SAVI or OSAVI; the soil-adjustment term is precisely the fix for this artefact.

NDVI plateau mistaken for uniform vigour. Once a corn or wheat canopy closes, NDVI saturates and a genuinely variable nitrogen field maps as a flat 0.87 sheet. Teams then conclude the field is uniform and skip variable-rate topdress. The tell is a compressed NDVI histogram (interquartile range < 0.03) on a full canopy — switch to NDRE, whose histogram will re-open the gradient.

Red-edge requested from a platform that lacks it. Copying an NDRE recipe onto a Phantom 4 RGB flight throws a KeyError at best and silently reuses the red band as red-edge at worst, producing a plausible-looking but meaningless raster. The computable_indices() guard in Step 2 refuses this case explicitly; never index a stack without asserting the band set first.

Cross-date drift from uncalibrated reflectance. An index computed from raw DN on two dates conflates a sunnier afternoon with a greener crop. This is a leading cause of false stress alerts in time series; calibrate against reflectance panels during ingest and, for satellite stacks, mask contaminated pixels using cloud masking for agricultural imagery before any index comparison.

SAVI L frozen across the season. Holding L = 0.5 from emergence to closure over-damps the signal once the canopy fills, flattening real variation. Either step L down with canopy fraction or adopt OSAVI’s fixed 0.16 for a defensible constant.

Mixed-pixel edges at field boundaries. Half-vegetation, half-road pixels along a headland skew any index. Clip to the field polygon and handle the border, as covered in handling edge effects in raster index generation, before deriving zone statistics.

5. Verification & Output Validation

An index raster is only trustworthy if its range, distribution, and stage-appropriateness check out. Run these guards after index_for_scene().

import numpy as np

def validate_index(arr: np.ndarray, name: str) -> dict:
    """Confirm an index array is physically plausible for its type."""
    name = name.lower()
    finite = arr[np.isfinite(arr)]
    assert finite.size > 0, "Index is all-NaN — check band scaling and nodata"

    # Normalised-difference indices live in [-1, 1]; CIre is unbounded but > -1
    if name in {"ndvi", "ndre", "savi", "osavi", "gndvi"}:
        assert finite.min() >= -1.001 and finite.max() <= 1.001, (
            f"{name.upper()} outside [-1, 1]: {finite.min():.3f}..{finite.max():.3f} "
            "— reflectance likely unscaled (raw DN?)"
        )
    elif name == "cire":
        assert finite.min() >= -1.001, "CIre below -1 is impossible; check NIR/RedEdge"

    report = {
        "index": name,
        "valid_fraction": float(finite.size / arr.size),
        "median": float(np.median(finite)),
        "iqr": float(np.subtract(*np.percentile(finite, [75, 25]))),
    }
    # A saturated NDVI on a closed canopy shows a collapsed IQR — warn to switch to NDRE
    if name == "ndvi" and report["median"] > 0.8 and report["iqr"] < 0.03:
        print("WARNING: NDVI appears saturated (median>0.8, IQR<0.03) — use NDRE for this stage")
    return report

demo_ndvi = compute_index("ndvi", demo)
print(validate_index(demo_ndvi, "ndvi"))

Beyond the numeric guards, cross-check the choice against ground truth: overlay the index on a true-colour composite in QGIS and confirm high-index zones correspond to visibly denser canopy, and that a scouted stress patch appears in the red-edge index at grain fill even when NDVI shows nothing there. If a stage-appropriate index and NDVI disagree on a closed canopy, trust the red-edge index — that disagreement is the saturation you selected around.

6. Integration with the Pipeline

Index selection sits between calibrated imagery and every decision product downstream.

Upstream. Selection assumes calibrated, cloud-clean reflectance. Ingest and calibration come from ingesting multispectral drone imagery; for satellite stacks, pixels must first pass through cloud masking for agricultural imagery so masked values never enter the ratio.

Computation. The dispatcher here is the decision wrapper around the raw formulas in band math and raster algebra in Python; for the windowed NDVI/NDRE and SAVI implementations, use calculating NDVI and NDRE with Rasterio step by step and the SAVI calculation on drone tiles script.

Downstream. The chosen index feeds threshold mapping for crop health detection to delineate management zones, and its multi-date form feeds temporal aggregation of vegetation indices for season-long composites. For a worked crop-specific application of everything here, see NDVI vs NDRE vs SAVI for corn growth stages, which walks the same decisions through concrete V- and R-stages.

Frequently Asked Questions

Why does NDVI stop responding after canopy closure?

NDVI saturates once leaf area index rises above roughly 3 to 4 because the red band reflectance floor is already near zero and the NIR term dominates the ratio. Additional biomass no longer moves the value, so a full high-vigour canopy and a moderately dense one both read near 0.85. Switch to a red-edge index such as NDRE or CIre after closure, since red-edge reflectance keeps changing with chlorophyll content well past the point where NDVI flatlines.

Which index should I use when the soil is still visible between rows?

Use SAVI or its optimised variant OSAVI while bare soil is in the pixel footprint. Plain NDVI is biased by soil brightness at low canopy cover, reading artificially high over dark moist soil and low over bright dry soil. The L soil-adjustment factor damps that background contribution; set L near 0.5 for partial cover and drop toward 0.16 with OSAVI as the canopy closes.

Can I compute NDRE without a red-edge camera?

No. NDRE requires a discrete red-edge band around 717 nanometres, which RGB payloads and standard four-band NIR cameras do not capture. Sensors such as the MicaSense RedEdge-MX, DJI P4 Multispectral, Parrot Sequoia, and Sentinel-2 provide it, but a DJI Phantom 4 RGB does not. Without red-edge, fall back to NDVI or GNDVI for vigour and treat late-season nitrogen status as unobservable from that platform.

Related

• NDVI vs NDRE vs SAVI for Corn Growth Stages — a concrete V-stage and R-stage walkthrough applying this decision matrix to maize
• Band Math & Raster Algebra in Python — the raw index formulas and windowed raster patterns behind the dispatcher
• Calculating NDVI and NDRE with Rasterio Step by Step — production windowed implementations of the two most-used indices
• Python Script for SAVI Calculation on Drone Tiles — the soil-adjusted derivation and tile-based SAVI code
• Threshold Mapping for Crop Health Detection — turning the selected index raster into management zones