Computing SAVI and MSAVI for Sparse Vegetation

Both indices need reflectance in the 0–1 range, not scaled integers:

import numpy as np

nir = nir_dn * 1e-4          # Sentinel-2 L2A scaling (apply offset for newer baselines)
red = red_dn * 1e-4
L = 0.5
savi = (1 + L) * (nir - red) / (nir + red + L)
msavi2 = (2 * nir + 1 - np.sqrt((2 * nir + 1) ** 2 - 8 * (nir - red))) / 2

Where vegetation is sparse, the soil between plants dominates the pixel, and NDVI responds to soil brightness as much as to plants. This page belongs to spectral index calculation pipelines in Satellite Processing Workflows & Index Pipelines.

Why NDVI Struggles over Bare Soil

The soil line and index isolines In red–NIR space, bare soils of different brightness lie along a soil line. NDVI isolines all pass through the origin, so dark and bright soils with the same plant cover fall on different NDVI values. SAVI shifts the convergence point away from the origin so its isolines run closer to parallel with the soil line, reducing the soil effect. Red–NIR space for sparse cover soil line NDVI isolines meet at origin SAVI isolines closer to parallel

Every NDVI isoline passes through the origin of red–NIR space, so a dark wet soil and a bright dry soil with identical plant cover sit on different isolines and get different NDVI values. The soil-adjustment constant L moves the convergence point to the negative side of the origin, rotating the isolines so they run closer to parallel with the soil line. The practical result is that soil brightness changes — after rain, after tillage, across soil types — move the index less.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Index arithmetic
xarray >=2023.1 Labelled arrays for time series
rasterio >=1.3.0 Reading bands and scale metadata 
pip install "numpy>=1.23" "xarray>=2023.1" "rasterio>=1.3.0"

Complete Working Example

import numpy as np
import xarray as xr


def to_reflectance(dn: xr.DataArray, scale: float = 1e-4, offset: float = 0.0) -> xr.DataArray:
    r = dn.astype("float32") * scale + offset
    return r.where((r > -0.1) & (r < 1.5))           # drop impossible values


def savi(nir, red, L: float = 0.5) -> xr.DataArray:
    return ((1 + L) * (nir - red) / (nir + red + L)).rename("savi")


def msavi2(nir, red) -> xr.DataArray:
    a = 2 * nir + 1
    disc = np.maximum(a ** 2 - 8 * (nir - red), 0)    # guard against tiny negatives
    return ((a - np.sqrt(disc)) / 2).rename("msavi2")


def ndvi(nir, red) -> xr.DataArray:
    return ((nir - red) / (nir + red)).rename("ndvi")


def soil_adjusted_indices(ds: xr.Dataset, scale=1e-4, offset=0.0) -> xr.Dataset:
    nir = to_reflectance(ds.nir, scale, offset)
    red = to_reflectance(ds.red, scale, offset)
    return xr.merge([ndvi(nir, red), savi(nir, red), msavi2(nir, red)])

The single most common bug with SAVI is feeding it integer digital numbers. The L constant is designed for reflectance between 0 and 1; added to values in the thousands it is negligible, and SAVI silently becomes NDVI multiplied by 1.5. Converting to reflectance first — with the product’s scale and offset, which changed for Sentinel-2 processing baseline 04.00 — is not optional. The index functions follow the same shape as those in writing reusable index functions with xarray apply_ufunc.

Choosing L, or Not Choosing It

L by cover, or MSAVI2 to skip the choice An L of 1 suits very sparse cover, 0.5 intermediate cover and 0.25 dense cover; with L of zero SAVI equals NDVI. Because cover varies across a scene and through a season, a single L is always a compromise. MSAVI2 computes an effective L per pixel from the data itself, so no constant is needed. The soil adjustment factor L = 1 very sparse L = 0.5 intermediate L = 0.25 dense MSAVI2 L chosen per pixel L = 0.5 is the conventional default; MSAVI2 is the safer choice when cover varies widely.

The original SAVI paper recommends L = 0.5 as a general-purpose value, and most software uses it. That works well for rangeland and cropland mid-season. But a scene containing both bare fields and closed canopy has no single right L, and a time series spanning bare soil to full canopy changes the ideal value every few weeks. MSAVI2 solves this by computing, in closed form, the L that each pixel’s own red and NIR imply. It has no parameter to tune, which also makes it easier to compare across studies.

When the Adjustment Is Worth It

Soil adjustment matters most where fractional cover is below roughly 40%: drylands, rangeland, orchards and vineyards with bare inter-rows, and cropland early in the season. There, NDVI of the same crop can differ by 0.1 between a light sandy field and a dark clay field, and SAVI or MSAVI narrow that gap substantially. In dense forest and mature crops, the soil contribution is small and all three indices tell the same story; SAVI then offers little except a slightly lower saturation point. For bare-soil and water mapping, none of them is the right tool — use indices designed for those surfaces, such as those in computing NDWI and MNDWI for water mapping.

Adding Them to a Multi-Index Pipeline

SAVI and MSAVI2 slot into a configuration-driven pipeline the same way any other index does: declare the bands, the formula and any parameter, and let the pipeline handle loading, scaling and writing. Keeping L as a configuration value rather than a literal in the code means it can be set per region, and recording it in the output’s metadata means a SAVI map can always be traced back to the constant that produced it. The pattern is described in building a YAML-driven multi-index pipeline. Write the result as float32 or, if storage matters, as int16 with an explicit scale factor of 1e-4, which preserves four decimal places — far more precision than the index’s physical meaning supports.

Verification

Dark and bright soils should read alike For two bare soils of different brightness, NDVI differs noticeably between them, SAVI less so, and MSAVI2 least. A pipeline whose SAVI shows the same gap as NDVI is probably fed unscaled integers. Index gap between dark and bright bare soil NDVI SAVI MSAVI2 Smaller is better: less sensitivity to soil brightness. 
import numpy as np

idx = soil_adjusted_indices(ds)
dark, bright = soil_masks["dark"], soil_masks["bright"]      # bare-soil polygons rasterised
for name in ("ndvi", "savi", "msavi2"):
    gap = abs(float(idx[name].where(bright).mean()) - float(idx[name].where(dark).mean()))
    print(f"{name:7s} soil gap {gap:.3f}")
assert float(idx.savi.max()) <= 1.5, "SAVI above 1.5 means inputs were not reflectance"

Common Errors

SAVI is exactly 1.5 × NDVI

Inputs are digital numbers, so L is negligible. Convert to reflectance first.

MSAVI2 returns NaN in a few pixels

Floating-point error makes the discriminant slightly negative. Clip it at zero before the square root.

Values drift after early 2022

Sentinel-2 added a −1000 offset at processing baseline 04.00. Apply the offset from the metadata.

SAVI looks noisier than NDVI

Negative or near-zero reflectance in water and shadow. Mask those surfaces first.

Frequently Asked Questions

Q: What L value should I use for SAVI? 0.5 is the standard default and works for intermediate cover. Use values near 1 for very sparse cover and near 0.25 for dense cover, or use MSAVI2, which chooses the adjustment per pixel.

Q: What is the difference between MSAVI and MSAVI2? MSAVI originally required an iterative computation of L; MSAVI2 is the closed-form solution of that iteration and is what almost all software implements under either name.

Q: Why is my SAVI identical to scaled NDVI? Because the bands were integers rather than reflectance, making the L term negligible. Apply the scale factor and offset first.

Q: Is SAVI better than NDVI everywhere? No. It helps where soil is visible between plants. In dense vegetation all indices agree, and NDVI is more widely used for comparison.

Q: Can I use SAVI with Landsat? Yes. Use the surface reflectance scale and offset for Collection 2, which differ from Sentinel-2’s, then apply the same formula.