⚙️ How-To Satellite imagery Raster & rasterio Satellite Imagery · step 8 of 16 Updated

How to Calculate NDVI from Sentinel-2 in Python

Problem statement

NDVI is two bands and one line of arithmetic, which is why it is the first thing anyone computes from satellite imagery and the thing most often computed wrongly.

The formula is not the hard part. The hard parts are the four things that happen around it:

  • the bands are unsigned integers on a scale nobody told you about
  • half the pixels are cloud, cirrus or cloud shadow
  • a shadowed canopy still returns a healthy-looking NDVI of 0.62
  • the denominator can be zero, or negative, and NumPy will not stop you

Done properly on a real Sentinel-2 scene over Snowdonia, vegetation comes out at a median NDVI of 0.730 and water at −0.096. Done in the obvious order, the same pixels give 1.376.

Quick answer 

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

UNUSABLE = (0, 1, 3, 8, 9, 10)     # nodata, saturated, shadow, cloud x2, cirrus

with rasterio.open("B04.tif") as ds:                  # red, 10 m
    red = ds.read(1).astype("float32") * 1e-4
    shape, profile = ds.shape, ds.profile
with rasterio.open("B08.tif") as ds:                  # nir, 10 m
    nir = ds.read(1).astype("float32") * 1e-4
with rasterio.open("SCL.tif") as ds:                  # classification, 20 m
    scl = ds.read(1, out_shape=shape, resampling=Resampling.nearest)

usable = ~np.isin(scl, UNUSABLE)
red = np.where(usable, red, np.nan)
nir = np.where(usable, nir, np.nan)

with np.errstate(divide="ignore", invalid="ignore"):
    total = nir + red
    ndvi = np.where(np.abs(total) < 1e-6, np.nan, (nir - red) / total)

print(f"NDVI {np.nanmin(ndvi):.3f} to {np.nanmax(ndvi):.3f}, "
      f"median {np.nanmedian(ndvi):.3f}, {np.isnan(ndvi).mean():.1%} masked")
NDVI -1.000 to 0.942, median 0.715, 51.2% masked
The five steps of an NDVI calculation in order: read, cast to float, scale, mask, then divide — with the failure each wrong order produces.
Each step out of order produces a specific wrong answer, and none of them raises an exception.

Step-by-step solution

1. Cast to float before you touch the arrays 

red = ds.read(1)
print(red.dtype, red.min())      # uint16 1
print(red.min() - 1000)          # 64537

Unsigned integer arithmetic wraps. Any subtraction on a raw band — an offset, a difference, a dark-object correction — must happen after .astype("float32").

2. Scale to reflectance, and check the scale

NDVI is a ratio, so a pure multiplicative scale cancels out: (nir − red)/(nir + red) is identical whether both bands are in DN or reflectance. The scale only matters if there is an offset, and then it matters enormously.

Check against water before trusting a declared offset:

water median reflectance
  DN * 1e-4          0.008     plausible
  DN * 1e-4 - 0.1   -0.092     impossible

The archive measured here declares offset: -0.1 in STAC metadata but ships COGs with the offset already applied. Applying it again gives vegetation an NDVI of 1.376. See Radiometric levels explained.

3. Mask before dividing, not after

Cloud shadow over vegetation returns a median NDVI of 0.622 — comfortably inside the range people call healthy. A ratio cancels the proportional darkening that a shadow causes, so nothing about the NDVI value reveals the pixel is unusable.

usable = ~np.isin(scl, (0, 1, 3, 8, 9, 10))

Over this scene that leaves 48.9% of pixels. Dilating the mask by two pixels to catch cloud edges leaves 44.0%.

4. Handle the zero denominator explicitly

nir + red is zero wherever both bands are zero — no-data fill, or the deep shadow of a scene edge. In NumPy that gives nan with a RuntimeWarning, or inf if only one is zero.

with np.errstate(divide="ignore", invalid="ignore"):
    total = nir + red
    ndvi = np.where(np.abs(total) < 1e-6, np.nan, (nir - red) / total)

Returning NaN is right; returning 0 is wrong, because 0 is a legitimate NDVI meaning "as much red as near-infrared".

5. Check the output range before you use it 

assert np.nanmin(ndvi) >= -1.0001 and np.nanmax(ndvi) <= 1.0001

NDVI is bounded by −1 and 1 whenever both inputs are non-negative. Anything outside means the inputs were not non-negative, which means the scaling is wrong. This assertion costs nothing and catches the single most common error in the whole workflow.

NDVI medians by surface class: vegetation 0.730, cloud shadow 0.622, bare soil 0.353, cloud 0.083, water minus 0.096.
Measured over a real scene. Cloud shadow at 0.622 is why the mask is not optional.

Code examples

Example 1 — NDVI with every check built in 

import numpy as np
import rasterio
from rasterio.enums import Resampling
from scipy.ndimage import binary_dilation

UNUSABLE = (0, 1, 3, 8, 9, 10)


def ndvi_from_scene(red_path, nir_path, scl_path=None, scale=1e-4,
                    offset=0.0, dilate_px=2):
    """NDVI as float32 with NaN where the pixel is not usable ground."""
    with rasterio.open(red_path) as ds:
        red = ds.read(1)
        shape, profile = ds.shape, ds.profile
    with rasterio.open(nir_path) as ds:
        nir = ds.read(1)

    if nir.shape != red.shape:
        raise ValueError(
            f"red is {red.shape} and nir is {nir.shape} — "
            "resample onto one grid before computing an index"
        )

    # cast first: these are unsigned and any offset would wrap
    red = red.astype("float32") * scale + offset
    nir = nir.astype("float32") * scale + offset

    if scl_path:
        with rasterio.open(scl_path) as ds:
            scl = ds.read(1, out_shape=shape, resampling=Resampling.nearest)
        bad = np.isin(scl, list(UNUSABLE))
        if dilate_px:
            bad = binary_dilation(bad, iterations=dilate_px)
        red = np.where(bad, np.nan, red)
        nir = np.where(bad, np.nan, nir)

    with np.errstate(divide="ignore", invalid="ignore"):
        total = nir + red
        ndvi = np.where(np.abs(total) < 1e-6, np.nan,
                        (nir - red) / total).astype("float32")

    finite = ndvi[np.isfinite(ndvi)]
    if finite.size and (finite.min() < -1.001 or finite.max() > 1.001):
        raise ValueError(
            f"NDVI outside [-1, 1] ({finite.min():.3f} to {finite.max():.3f}) — "
            f"check scale={scale} offset={offset} against a known dark target"
        )

    profile.update(dtype="float32", count=1, nodata=np.nan, compress="deflate")
    print(f"  NDVI median {np.nanmedian(ndvi):.3f}, "
          f"{np.isnan(ndvi).mean():.1%} masked")
    return ndvi, profile
  NDVI median 0.720, 54.0% masked

Dilating the mask by two pixels drops another 2.8% of the scene and lifts the median from 0.715 to 0.720 — the cloud-edge pixels it removed were dragging the answer down. Both numbers are the mask doing its job.

Example 2 — writing NDVI out so it stays readable 

import numpy as np
import rasterio


def write_index(path, values, profile, name="NDVI", source=None,
                as_int16=False):
    """Write an index raster with the metadata needed to interpret it later."""
    profile = dict(profile)
    if as_int16:
        # 10,000x scaling keeps four decimals in half the bytes
        data = np.where(np.isnan(values), -32768,
                        np.round(values * 10000)).astype("int16")
        profile.update(dtype="int16", nodata=-32768)
    else:
        data = values.astype("float32")
        profile.update(dtype="float32", nodata=np.nan)

    profile.update(count=1, compress="deflate", tiled=True,
                   blockxsize=512, blockysize=512)

    with rasterio.open(path, "w", **profile) as dst:
        dst.write(data, 1)
        dst.update_tags(index=name, scale=0.0001 if as_int16 else 1.0,
                        source=source or "", valid_range="-1 to 1")
        dst.build_overviews([2, 4, 8, 16])
    print(f"  wrote {path} as {profile['dtype']}")

Store the scale in the tags. An int16 NDVI raster with no recorded scale factor is a raster of numbers between −10,000 and 10,000 that somebody will plot as-is.

Example 3 — NDVI over a stack, keeping only real observations 

import numpy as np


def ndvi_stack_median(scenes):
    """Per-pixel median NDVI across dates, ignoring masked observations."""
    stack = np.stack([ndvi for ndvi, _ in scenes])          # (time, y, x)
    n_obs = np.isfinite(stack).sum(axis=0)

    with np.errstate(all="ignore"):
        median = np.nanmedian(stack, axis=0)

    # a pixel seen once is not a median; say so rather than pretending
    median = np.where(n_obs >= 3, median, np.nan)

    print(f"  {len(scenes)} scenes, observations per pixel: "
          f"min {n_obs.min()}, median {np.median(n_obs):.0f}, max {n_obs.max()}")
    print(f"  pixels with fewer than 3 clear looks: "
          f"{(n_obs < 3).mean():.1%}")
    return median, n_obs

Returning n_obs alongside the composite is what makes it interpretable. A median of three observations and a median of thirty look identical in the output raster and are not the same measurement.

Explanation

Why NDVI works at all

Chlorophyll absorbs red light strongly for photosynthesis. The internal structure of leaf mesophyll scatters near-infrared light strongly, because there is nothing in a leaf that absorbs at 0.84 µm.

The result is the largest contrast in the whole visible-to-near-infrared range for any natural surface. Measured over this scene, vegetation reflectance was 0.056 in red and 0.361 in near-infrared — a factor of six. Bare soil rises gently across both, giving a much smaller difference, and water absorbs both.

The normalised form divides by the sum, which cancels anything that scales both bands equally: sun angle, slope, aspect, thin haze. That is what makes NDVI comparable across a hilly scene rather than just across a flat one.

Why the mask matters more than the formula

That same cancellation is why cloud shadow gets through. A shadow reduces the illumination on both bands by roughly the same factor, and the ratio does not change much. Measured: 0.622 in shadow against 0.730 in full sun — a real drop, but nowhere near enough to identify shadow from the NDVI value alone.

So the mask is the only part of the pipeline that knows. Once NDVI has been computed on an unmasked array and averaged into a field statistic, the information that a third of the field was in shadow is gone.

Why the range check catches nearly everything

Every common failure in this workflow produces an out-of-range NDVI:

  • an offset applied twice makes the denominator negative → values beyond ±1
  • unsigned wraparound makes one band enormous → values approaching ±1 everywhere
  • mixing a scaled band with an unscaled band → arbitrary values
  • reading the wrong band → usually still in range, and the only one this misses

Two comparisons catch three of the four. That is a good return on one assertion.

Why NDVI saturates, and when to use something else

NDVI's response flattens once a canopy is closed: red absorption is already nearly complete, so more leaf area adds near-infrared but cannot remove more red. Measured over vegetation here, NDVI's 5th-to-95th percentile range was 0.53 to 0.855, compressed against its ceiling, while EVI on the same pixels spanned 0.318 to 0.723.

If your question is about dense forest biomass, NDVI is the wrong instrument. If it is about where vegetation is and roughly how much, NDVI is hard to beat for two bands and one division.

Four NDVI failure modes against a single range assertion, which catches three of them.
One assertion, three of four failures. The fourth needs a class-median check.

Edge cases or notes

  • The scale cancels; the offset does not. A pure multiplicative scale leaves NDVI unchanged. Only add an offset if you have verified it is needed.
  • nodata is 0, which is a legal reflectance value. Mask on the classification band, not on zero.
  • np.nanmedian over an all-NaN column warns and returns NaN. Wrap stack statistics in np.errstate and filter by observation count.
  • NDVI of exactly −1 or 1 means one band was zero. Those are edge or fill pixels, not extreme vegetation.
  • Do not store NDVI as uint8. Scaling −1..1 into 0..255 loses three decimal places and hides small changes.
  • Bilinear resampling of NDVI is not the same as NDVI of bilinear bands. Compute the index at native resolution, then resample the result if you must.
  • Snow gives NDVI near 0 and can be mistaken for bare soil. Class 11 in the scene classification separates them.
  • A single-date NDVI over a cloudy region is mostly holes. Over Snowdonia the median scene was 0.2% usable — use a composite.

FAQ

How do I calculate NDVI in Python?

Read the red and near-infrared bands, cast to float, put them on the same grid, mask cloud and shadow using the scene classification band, then compute (nir − red) / (nir + red) with the zero denominator handled.

Do I need to convert to reflectance first?

Not for NDVI alone — a multiplicative scale cancels in the ratio. You do need to get any offset right, and the range check will tell you if you have not.

Why is my NDVI above 1?

The denominator went negative, almost always from applying a radiometric offset that was already applied, or from unsigned integer wraparound. See NDVI values are above 1 or below −1.

What NDVI value means healthy vegetation?

Over this scene, vegetation had a median of 0.730. But cloud shadow reached 0.622, so a threshold on its own cannot separate vegetation from shadowed vegetation.

Should I mask clouds before or after computing NDVI?

Before. NDVI is a ratio and largely cancels the darkening a shadow causes, so it cannot detect what the mask missed.

How do I save NDVI to a GeoTIFF?

As float32 with NaN nodata, or as int16 scaled by 10,000 with the scale recorded in the file tags. Never as uint8.

Why does my NDVI raster have holes?

Because those pixels were cloud, cirrus, shadow or no data. On a single Sentinel-2 date over Snowdonia that was 51% of the scene. Build a composite from several dates.