Python Script for SAVI Calculation on Drone Tiles

TL;DR: Read the NIR and Red bands from a multispectral drone GeoTIFF, apply SAVI = ((NIR - Red) / (NIR + Red + L)) * (1 + L) with L = 0.5, mask nodata pixels, and write a georeferenced float32 output — the complete runnable script is in Step 4 below.

Why SAVI Matters in Drone-Based Crop Monitoring

NDVI is the default vegetation index for most workflows, but it breaks down in the conditions that make drone imagery most valuable: early-season row crops with large gaps between plant rows, post-harvest residue fields, and newly transplanted orchards. In all three cases, the sensor is looking at a mixture of canopy and bare soil, and NDVI systematically overestimates green biomass because it treats bright soil in the NIR channel as vegetation signal.

The Soil Adjusted Vegetation Index introduces a correction factor L that shifts the NIR–Red denominator away from zero for bare soil, effectively flattening the soil line in NIR-Red feature space. Without this adjustment, a prescription map generated from an early-June flight of a corn field at V3–V4 growth stage will overstate canopy cover by 15–30%, leading to under-application of nitrogen in the zones that need it most. Running SAVI through the band math & raster algebra in Python pipeline before threshold mapping keeps the soil signal from corrupting management-zone boundaries.

The diagram below illustrates the role SAVI plays in the drone-to-prescription pipeline:

Prerequisites

This page assumes the environment described in the parent Band Math & Raster Algebra in Python cluster. The only additional requirement is that your tiles are already in surface reflectance (not raw digital numbers). If you are still working with DN outputs from a DJI P4 Multispectral or MicaSense RedEdge-MX, apply the sensor’s calibration coefficients before running the code below.

• rasterio >= 1.3
• numpy >= 1.24
• Input: single-tile GeoTIFF with at least NIR and Red bands, embedded CRS (EPSG code), nodata attribute set
• Band indices are 1-based (rasterio convention)

Step-by-Step Implementation

Step 1 — Identify Your Band Order

The band layout varies by sensor and export settings. Common configurations:

Always cross-check against the sensor’s band manifest or the GeoTIFF metadata before running any band math. Swapped NIR and Red bands will produce a SAVI raster with inverted values — vegetation appears as soil and vice versa, which is one of the most common and least obvious errors in production pipelines. The same issue applies when calculating NDVI and NDRE with rasterio, so establish a consistent band-ordering convention across your entire workflow before scaling.

Step 2 — Choose Your L Factor

The default L = 0.5 is appropriate for most mid-season flights. If you have ground-truth spectrometer readings or calibration panel measurements, use them to validate the L choice against measured greenness before locking in a value for a campaign.

Step 3 — Understand the Formula

SAVI = ((NIR - Red) / (NIR + Red + L)) * (1 + L)

The (1 + L) multiplier restores the dynamic range that the L offset would otherwise compress. When L = 0, the formula collapses to NDVI. The denominator (NIR + Red + L) cannot be zero for any L > 0 on physically valid reflectance data, but floating-point rounding after aggressive radiometric calibration can produce near-zero denominators — hence the guard in the script below.

Step 4 — Runnable Script

import rasterio
import numpy as np
from pathlib import Path
from typing import Optional


def calculate_savi_tile(
    input_path: str,
    output_path: str,
    nir_band_idx: int,
    red_band_idx: int,
    l_factor: float = 0.5,
    nodata_out: float = -9999.0,
    nodata_override: Optional[float] = None,
) -> None:
    """
    Compute SAVI on a single multispectral drone tile and write a
    georeferenced float32 GeoTIFF.

    Band indices follow rasterio's 1-based convention.
    Input bands must be in surface reflectance units.
    """
    src_path = Path(input_path)
    if not src_path.exists():
        raise FileNotFoundError(f"Input tile not found: {src_path}")

    with rasterio.open(src_path) as src:
        if src.crs is None:
            raise ValueError(
                f"Tile has no CRS. Reproject to a defined coordinate system "
                f"(e.g. EPSG:32632) before running band math."
            )
        if not (1 <= nir_band_idx <= src.count):
            raise ValueError(
                f"nir_band_idx={nir_band_idx} out of range (tile has {src.count} bands)"
            )
        if not (1 <= red_band_idx <= src.count):
            raise ValueError(
                f"red_band_idx={red_band_idx} out of range (tile has {src.count} bands)"
            )

        # Read only the required bands to minimise RAM footprint
        nir = src.read(nir_band_idx).astype(np.float32)
        red = src.read(red_band_idx).astype(np.float32)

        # Build validity mask: explicit nodata + physically impossible reflectance
        src_nodata = nodata_override if nodata_override is not None else src.nodata
        invalid = np.zeros(nir.shape, dtype=bool)
        if src_nodata is not None:
            invalid |= (nir == src_nodata) | (red == src_nodata)
        invalid |= (nir < 0) | (red < 0) | (nir > 1.5) | (red > 1.5)

        # SAVI with zero-denominator guard
        denom = nir + red + l_factor
        with np.errstate(divide="ignore", invalid="ignore"):
            savi = np.where(
                denom != 0,
                ((nir - red) / denom) * (1.0 + l_factor),
                nodata_out,
            ).astype(np.float32)

        savi[invalid] = nodata_out

        out_meta = src.meta.copy()
        out_meta.update(
            dtype="float32",
            count=1,
            nodata=nodata_out,
            compress="deflate",
            predictor=2,   # floating-point delta encoding
            tiled=True,
            blockxsize=256,
            blockysize=256,
        )

        Path(output_path).parent.mkdir(parents=True, exist_ok=True)
        with rasterio.open(output_path, "w", **out_meta) as dst:
            dst.write(savi, 1)

    print(f"SAVI written → {output_path}  (nodata={nodata_out}, L={l_factor})")


# Quick usage example (DJI P4 Multispectral, 5-band stack, UTM projection):
# calculate_savi_tile(
#     input_path="flight_2026-06-15/tile_042.tif",
#     output_path="indices/tile_042_savi.tif",
#     nir_band_idx=5,
#     red_band_idx=3,
#     l_factor=0.5,
# )

Inline verification — run this immediately after the function returns:

with rasterio.open("indices/tile_042_savi.tif") as check:
    data = check.read(1)
    valid = data[data != -9999.0]
    assert valid.min() >= -1.0 and valid.max() <= 1.0, (
        f"SAVI out of range: min={valid.min():.4f}, max={valid.max():.4f} — "
        "check band order and reflectance calibration"
    )
    print(f"SAVI range: {valid.min():.3f} – {valid.max():.3f}, "
          f"mean={valid.mean():.3f}, nodata pixels: {(data == -9999.0).sum()}")

Gotchas and Edge Cases

• Band index is 1-based in rasterio.open but 0-based in the returned NumPy array. src.read(4) returns band four, but that data lands at array[0] (the only array dimension) not array[3]. Never mix the two conventions in the same expression.

• Reflectance values above 1.0 are physically valid but rare. Over bright artificial surfaces or calibration panels, reflectance can reach 1.2–1.4. The script’s upper guard (nir > 1.5) avoids masking those pixels. If your calibration produces values well above 1.5 in vegetated areas, the calibration coefficients are likely wrong — cross-check against the manufacturer’s panel reflectance spec.

• CRS must be set on the tile before band math. The function raises immediately if src.crs is None. If you are processing tiles exported without a projection, reproject them to the appropriate UTM zone first using the approach in understanding CRS in precision agriculture before running SAVI. Never embed a reprojection inside the index calculation loop — resampling artifacts will distort spectral values.

• nodata_out = -9999.0 must match what downstream tools expect. QGIS and most Python raster libraries read the embedded nodata value from the file profile, but some ISOXML export tools expect NaN. Pass nodata_out=float("nan") if your prescription pipeline requires it, and add np.isnan to the verification assertion.

This guide is part of Band Math & Raster Algebra in Python — see there for the full pipeline context, including windowed I/O patterns, nodata handling conventions, and integration with cloud-masking and threshold-mapping stages.

Related

• Calculating NDVI and NDRE with Rasterio Step by Step — apply normalized difference formulas for canopy stress and red-edge reflectance on the same tile stack
• Band Math & Raster Algebra in Python — parent guide covering the full raster algebra pipeline from band alignment to windowed export
• Threshold Mapping for Crop Health — use your SAVI output to delineate management zones and drive variable-rate prescriptions
• Handling Edge Effects in Raster Index Generation — prevent seam-line artefacts when stitching SAVI tiles into a full-field mosaic