Calculating NDVI from Sentinel-2 with rasterio

Turning Sentinel-2 Level-2A surface reflectance into a trustworthy Normalized Difference Vegetation Index (NDVI) raster is a narrow, repeatable task with three specific failure modes: integer truncation, a forgotten or double-applied radiometric offset, and an unmasked zero denominator. This page is the exact rasterio recipe for getting all three right, so the canopy-health values you hand to a forester or a conservation model are real signal rather than silent artifacts. It is one concrete workflow within Vegetation Index Calculation in Python, which sits inside the wider Ecological GIS Data Foundations in Python framework. The hard parts here are not the band ratio itself but the scaling, casting, and masking discipline that decides whether subtle reflectance gradients in low-biomass understories survive into the output.

The diagram below traces a single pixel through the transformation chain. Each stage is a place where a silent error can creep in — an integer that never becomes a float, an offset that is forgotten or double-applied, or a zero denominator that escapes the mask.

When to use this approach

Reach for this direct rasterio + NumPy recipe when you control the band files and want a transparent, auditable NDVI raster with explicit reflectance handling. Several alternatives exist, and the right one depends on how much of the pipeline you need to own:

If your bands are not yet on a common grid — different CRS, extent, or resolution — that is a prerequisite, not part of this recipe. Align them first via Coordinate Reference Systems for Forestry; NDVI arithmetic assumes the Red and NIR arrays already overlay pixel-for-pixel.

Data scaling and the radiometric offset

Sentinel-2 L2A surface reflectance products deliver the Near-Infrared (Band 8) and Red (Band 4) bands as 16-bit unsigned integers (uint16), scaled by a factor of 10,000 to preserve precision during atmospheric correction. Direct integer division truncates fractional differences, collapsing the NDVI range into binary artifacts that obscure subtle reflectance gradients. Arrays must therefore be cast to float32 before applying the standard formula (NIR - Red) / (NIR + Red). Refer to the official Sentinel-2 Level-2A processing documentation for authoritative reflectance scaling specifications.

The offset is the part most pipelines get wrong. Products generated under processing baseline 04.00 or later (acquisitions from 25 January 2022 onward) carry a radiometric BOA_ADD_OFFSET of -1000, so the correct conversion is reflectance = (DN - 1000) / 10000. Older SAFE-format products use the same scale with no offset. Always verify BOA_ADD_OFFSET and the scale in the product metadata (MTD_MSIL2A.xml) before computing indices — applying the offset to a pre-baseline product, or omitting it on a newer one, biases every reflectance value and drifts the whole NDVI surface.

Minimal reproducible example

The function below ingests separate Red (B04, 10 m) and NIR (B08, 10 m) GeoTIFFs, validates that they share a grid, applies the offset and scale, masks zero denominators, and writes a float32 NDVI raster. Sentinel-2 imagery is typically delivered as individual band GeoTIFFs within the granule’s IMG_DATA directory, not as a single multi-band file, which is why each band is opened separately.

Step 1 — Open both bands and validate the grid

Context-managed opening guarantees file descriptors close on exception, and the explicit CRS/shape checks catch the most common cause of garbage NDVI: two bands that look aligned but are not.

import rasterio
import numpy as np


def calculate_sentinel2_ndvi(red_path: str, nir_path: str, output_path: str) -> None:
    """
    Compute NDVI from separate Sentinel-2 L2A Red (B04) and NIR (B08) GeoTIFFs.
    Both bands must share identical CRS, extent, and resolution (10 m).
    """
    with rasterio.open(red_path) as red_src:
        red = red_src.read(1).astype("float32")
        out_meta = red_src.meta.copy()
        if red_src.crs is None:
            raise ValueError(f"Missing CRS in {red_path}")

    with rasterio.open(nir_path) as nir_src:
        nir = nir_src.read(1).astype("float32")
        if nir_src.crs != out_meta["crs"]:
            raise ValueError("Red and NIR bands have mismatched CRS.")
        if nir_src.shape != (out_meta["height"], out_meta["width"]):
            raise ValueError("Red and NIR bands have mismatched dimensions.")

    # Apply Sentinel-2 L2A radiometric offset (baseline >= 04.00) and scaling.
    # For pre-04.00 products set BOA_ADD_OFFSET = 0.0.
    BOA_ADD_OFFSET = -1000.0
    red = (red + BOA_ADD_OFFSET) / 10000.0
    nir = (nir + BOA_ADD_OFFSET) / 10000.0

    # Guard against zero denominators (clouds, water, no-data).
    denominator = nir + red
    mask = denominator == 0.0
    denominator[mask] = np.nan

    # Compute NDVI.
    ndvi = (nir - red) / denominator

    out_meta.update({
        "count": 1,
        "dtype": "float32",
        "nodata": np.nan,
    })

    with rasterio.open(output_path, "w", **out_meta) as dst:
        dst.write(ndvi, 1)

Step 2 — Mask cloud and shadow with the Scene Classification Layer

The zero-denominator guard only catches arithmetic blow-ups, not contaminated-but-finite pixels. For operational deployments, apply the Sentinel-2 Scene Classification Layer (SCL band, 20 m, resampled to 10 m) to drop cloud (SCL = 8, 9), cloud shadow (SCL = 3), and saturated/defective pixels (SCL = 1) before trusting the NDVI:

import numpy as np
import rasterio
from rasterio.enums import Resampling


def apply_scl_mask(ndvi: np.ndarray, scl_path: str, shape: tuple[int, int]) -> np.ndarray:
    """Set NDVI to NaN where the SCL flags cloud, shadow, or saturation."""
    invalid_classes = {1, 3, 8, 9, 10, 11}  # sat/defect, shadow, cloud med/high, cirrus, snow
    with rasterio.open(scl_path) as scl_src:
        scl = scl_src.read(
            1, out_shape=shape, resampling=Resampling.nearest
        )
    contaminated = np.isin(scl, list(invalid_classes))
    ndvi = ndvi.copy()
    ndvi[contaminated] = np.nan
    return ndvi

The SCL is categorical, so it is the one place nearest-neighbour resampling is correct — never bilinear, which would invent fractional class codes. Continuous bands are the opposite case (see the pitfalls below).

Step 3 — Process full scenes with windowed reads

Full-scene Sentinel-2 tiles (100 km × 100 km) at 10 m resolution produce arrays of roughly 10,000 × 10,000 pixels, about 400 MB per band as float32. Iterating block_windows() keeps peak memory bounded by one block while preserving spatial continuity:

import numpy as np
import rasterio


def calculate_ndvi_windowed(red_path: str, nir_path: str, output_path: str) -> None:
    """Memory-efficient NDVI via block-wise processing."""
    with rasterio.open(red_path) as red_src, rasterio.open(nir_path) as nir_src:
        meta = red_src.meta.copy()
        meta.update(count=1, dtype="float32", nodata=np.nan)

        with rasterio.open(output_path, "w", **meta) as dst:
            for _, window in red_src.block_windows(1):
                # Offset (-1000) applies to baseline >= 04.00 products; use 0 for older ones.
                red = (red_src.read(1, window=window).astype("float32") - 1000.0) / 10000.0
                nir = (nir_src.read(1, window=window).astype("float32") - 1000.0) / 10000.0
                denom = nir + red
                with np.errstate(divide="ignore", invalid="ignore"):
                    ndvi = np.where(denom == 0, np.nan, (nir - red) / denom)
                dst.write(ndvi, 1, window=window)

To clip the output to an administrative boundary or watershed polygon — reducing storage and accelerating downstream Raster-Vector Overlay Techniques — pass the geometry through rasterio.mask.mask after writing.

Parameter reference

Expected output and verification

A correct run produces a single-band float32 GeoTIFF on the same grid as the inputs, with NDVI values in [-1, 1] over land and NaN only where pixels were masked. Vegetated canopy typically lands in 0.3–0.9; bare soil and senescent vegetation near 0.1–0.25; water and built surfaces at or below zero. Validate before trusting it downstream:

import numpy as np
import rasterio

with rasterio.open("ndvi.tif") as src:
    ndvi = src.read(1)
    assert src.count == 1 and src.dtypes[0] == "float32"
    assert src.crs is not None, "NDVI is missing a CRS — check the input bands"

    finite = ndvi[np.isfinite(ndvi)]
    assert finite.min() >= -1.0 and finite.max() <= 1.0, "values outside [-1, 1]"

    # Quantify masking so a silently over-aggressive mask is visible.
    nan_fraction = np.isnan(ndvi).mean()
    print(f"masked fraction: {nan_fraction:.1%}")
    print(f"median vegetated NDVI: {np.nanmedian(finite[finite > 0]):.3f}")

A masked fraction far above the scene’s true cloud cover usually means the offset was wrong (shifting healthy pixels out of range) or the SCL mask was too broad. The np.isnan(ndvi).sum() count is also the right provenance figure to record before these outputs feed Spatial Plot Sampling Design workflows, where data loss per plot must be auditable.

Common pitfalls

• Skipping the float32 cast. uint16 arithmetic truncates (NIR - Red) / (NIR + Red) to integers, flattening NDVI to 0, 1, or -1 and erasing every gradient. Cast both bands before any division.
• Mishandling BOA_ADD_OFFSET. Applying -1000 to a pre-baseline-04.00 product (or omitting it on a newer one) biases reflectance and drifts NDVI by a consistent amount that is easy to miss. Read the baseline from MTD_MSIL2A.xml rather than assuming.
• Leaving the denominator unmasked. Cloud, water, and no-data pixels sum to zero after scaling, raising RuntimeWarning: invalid value encountered in divide and writing ±inf that ecological models cannot ingest. Force those pixels to NaN.
• Bilinear-resampling the SCL or nearest-resampling the NDVI. The categorical SCL must use nearest-neighbour; the continuous NDVI must use bilinear when reprojected. Swapping them either invents fake class codes or stamps step artifacts into the canopy gradient.

Frequently Asked Questions

Related

• Vegetation Index Calculation in Python — the parent workflow this recipe belongs to
• Ecological GIS Data Foundations in Python — the full data-foundations framework
• Coordinate Reference Systems for Forestry — align Red and NIR bands before computing the index
• Raster-Vector Overlay Techniques — clip NDVI to watersheds and management boundaries
• Spatial Plot Sampling Design — intersect validated NDVI with field plots