Computing EVI and SAVI for Closed Canopy in Python

NDVI is the default vegetation index, but over dense, closed-canopy forest it saturates — once the canopy is nearly continuous, adding leaf area barely moves the value, so NDVI flattens near its ceiling and stops discriminating condition. The Enhanced Vegetation Index (EVI) was built precisely to stay sensitive in that high-biomass regime, while the Soil-Adjusted Vegetation Index (SAVI) corrects for bright soil background where the canopy is not yet closed. This guide computes both with rasterio and numpy and shows when each earns its place over NDVI. It sits within the Vegetation Index Calculation in Python workflow, part of the broader Ecological GIS Data Foundations in Python framework.

When to Prefer EVI or SAVI over NDVI

All three indices contrast near-infrared against red reflectance; they differ in what confounder they neutralise. NDVI normalises illumination but has no defence against two things that dominate forest scenes: saturation in dense canopy, and bright soil showing through a sparse one. EVI adds a blue-band aerosol correction and a canopy background term to resist saturation; SAVI adds a soil-brightness term $L$ to suppress background reflectance.

The short rule: saturation problem → EVI; soil-background problem → SAVI; neither → NDVI. A practical caveat is that EVI depends on a trustworthy blue band, which is the noisiest, most aerosol-sensitive band on most sensors — if the blue band is compromised, SAVI is the more robust choice even in denser canopy. For the canonical NDVI walkthrough and its scaling and masking details, see the sibling Calculating NDVI from Sentinel-2 with rasterio guide.

The Formulas

NDVI, for reference, is the plain normalised difference:

$$\text{NDVI}=\frac{{\rho }_{\text{NIR}}-{\rho }_{\text{Red}}}{{\rho }_{\text{NIR}}+{\rho }_{\text{Red}}}$$

SAVI inserts a soil-brightness constant $L$ into the denominator and rescales by $(1+L)$ so the output keeps a comparable range. With $L=0.5$ for intermediate canopy cover ($L\to 0$ recovers NDVI for dense canopy; $L\to 1$ for very sparse cover):

$$\text{SAVI}=\frac{{\rho }_{\text{NIR}}-{\rho }_{\text{Red}}}{{\rho }_{\text{NIR}}+{\rho }_{\text{Red}}+L}(1+L)$$

EVI adds a gain $G$, a blue-band aerosol correction with coefficients $C_1$ and $C_2$, and a canopy background term $L$. The standard Sentinel-2/MODIS coefficients are $G=2.5$, $C_1=6$, $C_2=7.5$, $L=1$:

$$\text{EVI}=G\cdot \frac{{\rho }_{\text{NIR}}-{\rho }_{\text{Red}}}{{\rho }_{\text{NIR}}+C_1{\rho }_{\text{Red}}-C_2{\rho }_{\text{Blue}}+L}$$

All three require surface reflectance, not raw digital numbers. EVI in particular is defined on reflectance in $[0,1]$; feeding it unscaled integers produces values far outside its nominal range. Convert digital numbers to reflectance before any of this arithmetic.

Minimal Reproducible Example

The example reads red, NIR, and blue bands from a reflectance raster, computes SAVI and EVI with guarded denominators, and writes each as a float32 GeoTIFF. The denominator guard matters most for EVI, whose denominator can pass through zero over dark or masked pixels where NDVI’s cannot.

import numpy as np
import rasterio


def savi(nir, red, L=0.5):
    """Soil-Adjusted Vegetation Index; L=0.5 for partial canopy cover."""
    den = nir + red + L
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(den == 0, np.nan, (nir - red) / den * (1 + L))


def evi(nir, red, blue, G=2.5, C1=6.0, C2=7.5, L=1.0):
    """Enhanced Vegetation Index; stays sensitive in high-biomass forest."""
    den = nir + C1 * red - C2 * blue + L
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(den == 0, np.nan, G * (nir - red) / den)


def compute_indices(src_path, savi_path, evi_path):
    """Read reflectance bands (1=Red, 2=NIR, 3=Blue) and write SAVI and EVI."""
    with rasterio.open(src_path) as src:
        red = src.read(1).astype("float32")
        nir = src.read(2).astype("float32")
        blue = src.read(3).astype("float32")
        profile = src.profile.copy()

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

    savi_arr = savi(nir, red).astype("float32")
    evi_arr = evi(nir, red, blue).astype("float32")

    with rasterio.open(savi_path, "w", **profile) as dst:
        dst.write(savi_arr, 1)
    with rasterio.open(evi_path, "w", **profile) as dst:
        dst.write(evi_arr, 1)


compute_indices("reflectance_rgbnir.tif", "savi.tif", "evi.tif")

If your source is raw Sentinel-2 L2A digital numbers rather than reflectance, scale first — reflectance = (DN - 1000) / 10000 for baseline 04.00+ — then clip to [0, 1] before calling savi or evi, exactly as the parent Vegetation Index Calculation in Python pipeline does on ingest.

Parameter Reference

The coefficients and their ecological effect. Change the soil and aerosol terms deliberately; the gain and blue-correction constants are conventionally fixed.

SAVI’s $L$ is the one parameter you should tune to the scene: $0.25$ for near-closed canopy, $0.5$ for intermediate cover, and up to $1.0$ where soil dominates. EVI’s constants are a fixed convention — changing them breaks comparability with published EVI products.

Expected Output and Verification

SAVI, like NDVI, is bounded roughly to $[-1,1]$ (the $(1+L)$ rescaling keeps it in range). EVI is nominally in $[-1,1]$ for valid vegetated reflectance but can exceed it slightly over bright or noisy pixels, so verification asserts a slightly wider physical band and flags an excessive out-of-range fraction rather than a hard fail.

import numpy as np
import rasterio


def verify_index(path, lo=-1.0, hi=1.0, max_out_frac=0.02):
    with rasterio.open(path) as src:
        arr = src.read(1)
    valid = arr[np.isfinite(arr)]
    assert valid.size > 0, "index is entirely NaN — check scaling and band order"
    out_of_range = np.mean((valid < lo) | (valid > hi))
    assert out_of_range <= max_out_frac, (
        f"{out_of_range:.1%} of pixels outside [{lo}, {hi}] — "
        "likely unscaled DN or a bad blue band"
    )
    return {"min": float(valid.min()), "max": float(valid.max()),
            "mean": float(valid.mean()), "out_of_range_frac": float(out_of_range)}


print(verify_index("savi.tif"))
print(verify_index("evi.tif"))

A useful cross-check over dense forest: EVI should retain visible contrast where NDVI is flat. Extract both indices at plots spanning a canopy-density gradient and confirm EVI’s range across those plots is wider than NDVI’s — that widening is the saturation resistance you switched indices to get. Over young or post-fire stands, SAVI should sit below NDVI because it has discounted the soil signal NDVI counted as vegetation; that gap feeds naturally into recovery tracking alongside Burn Severity Mapping with Python.

Common Pitfalls

• Running EVI on unscaled digital numbers. EVI’s coefficients assume reflectance in $[0,1]$; on raw uint16 DN the denominator explodes and the output is meaningless. Convert to reflectance and clip before computing.
• A noisy blue band poisoning EVI. The blue band is the most aerosol- and haze-sensitive; a bad blue value drives EVI’s denominator to extremes. If the blue band is unreliable, use SAVI instead of forcing EVI.
• Leaving SAVI’s L at 0.5 everywhere. L = 0.5 suits intermediate cover; over near-closed canopy it under-corrects and over very sparse cover it over-corrects. Tune L to the scene’s canopy fraction.
• Integer truncation. Computing any of these on integer arrays floors the result. Cast every band to float32 before the arithmetic — the single most common silent corruption in index work.

Frequently Asked Questions

Related

• Vegetation Index Calculation in Python — the parent workflow covering NDVI, EVI, SAVI, and NDWI end to end.
• Calculating NDVI from Sentinel-2 with rasterio — the scaling, casting, and masking details for the canonical index.
• Burn Severity Mapping with Python — index-based change detection for post-fire recovery.
• Ecological GIS Data Foundations in Python — the reproducible-pipeline framework this guide belongs to.