⚙️ How-To Satellite imagery Raster & rasterio Satellite Imagery · step 8 Notebook By · Updated

How to Calculate NDVI from Sentinel-2 in Python

On this page

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.726 and water at −0.077. Done in the obvious order, the same pixels give 1.362.

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.935, median 0.710, 57.3% 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 has its own wrong answer1read the bandsuint162cast to float32or 1 − 1000 = 64,5373scale, checking the offsetor NDVI 1.3764mask cloud and shadowor shadow reads 0.6225divide, guarding zeroor infNone of the four failures raises an exception spatialworkflow.io
Each step out of order produces a specific wrong answer, and none of them raises an exception.

Cast to float, scale to reflectance, mask first

1. Cast to float before you touch the arrays

red = rasterio.open("B04.tif").read(1)
print(red.dtype, red.min())
print(red.min() - 1000)
uint16 1
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.362. See Radiometric levels explained.

3. Mask before dividing, not after

Cloud shadow over vegetation returns a median NDVI of 0.623 — 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 42.7% of pixels. Dilating the mask by two pixels to catch cloud edges leaves 39.7%.

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.726, cloud shadow 0.623, bare soil 0.362, cloud 0.154, water minus 0.077. measured NDVI medians, by scene-classification classvegetation0.726cloud shadow0.623bare soil0.362cloud0.15401Water sits at −0.077, off the bottom of this scale spatialworkflow.io
Measured over a real scene. Cloud shadow at 0.623 is why the mask is not optional.

NDVI with every check, and writing it out

The data used on this page

One Sentinel-2 L2A scene, S2A_30UVD_20250922_0_L2A (22 September 2025, processing baseline 05.11), from the Element 84 Earth Search sentinel-2-l2a collection, cut to a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa): x 422,880–432,880, y 5,875,440–5,885,440 in EPSG:32630. The setup writes the blue, green, red and near-infrared bands at 10 m (B02.tif, B03.tif, B04.tif, B08.tif) and the scene classification at 20 m (SCL.tif), each exactly as stored in the archive.

import os

import pystac_client
import rasterio
from rasterio.windows import from_bounds

BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
item = next(pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
            .search(ids=["S2A_30UVD_20250922_0_L2A"], collections=["sentinel-2-l2a"]).items())
for asset, name in (("blue", "B02"), ("green", "B03"), ("red", "B04"), ("nir", "B08"), ("scl", "SCL")):
    if os.path.exists(f"{name}.tif"):
        continue
    with rasterio.open(item.assets[asset].href) as src:
        window = from_bounds(*BOUNDS, transform=src.transform)
        profile = src.profile | {"width": int(window.width), "height": int(window.height),
                                 "transform": src.window_transform(window), "compress": "deflate",
                                 "driver": "GTiff", "tiled": False}
        with rasterio.open(f"{name}.tif", "w", **profile) as dst:
            dst.write(src.read(1, window=window), 1)
print({n: rasterio.open(f"{n}.tif").shape for n in ("B04", "B08", "SCL")})
{'B04': (1000, 1000), 'B08': (1000, 1000), 'SCL': (500, 500)}

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


ndvi2, profile2 = ndvi_from_scene("B04.tif", "B08.tif", "SCL.tif")
  NDVI median 0.716, 60.3% masked

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

The Snowdon window of 22 September 2025 by scene-classification class, with each class share and median NDVI. what the 10 km window containsSCL classsharemedian NDVIvegetation (4)36.8%0.726cloud (8, 9)31.2%0.154cloud shadow (3)23.2%0.623bare soil (5)2.6%0.362water (6)1.2%−0.077usable after masking: 42.7%; after two pixels of dilation: 39.7%Shadow is a quarter of the window and looks like healthy vegetation spatialworkflow.io
The mask is what separates the 0.726 from the 0.623; the NDVI value cannot.

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']}")


write_index("ndvi.tif", ndvi2, profile2, source="S2A_30UVD_20250922_0_L2A")
write_index("ndvi_int16.tif", ndvi2, profile2, source="S2A_30UVD_20250922_0_L2A", as_int16=True)
  wrote ndvi.tif as float32
  wrote ndvi_int16.tif as int16

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.

Example 4 — the numbers quoted on this page

b = {n: rasterio.open(f"{n}.tif").read(1).astype("float32") for n in ("B02", "B03", "B04", "B08")}
scl20 = rasterio.open("SCL.tif")
scl10 = scl20.read(1, out_shape=b["B04"].shape, resampling=Resampling.nearest)

water = scl10 == 6
print(f"water median reflectance: red  DN*1e-4 {np.median(b['B04'][water]) * 1e-4:.3f}, "
      f"DN*1e-4-0.1 {np.median(b['B04'][water]) * 1e-4 - 0.1:.3f}")


def nd(nir, red):
    with np.errstate(all="ignore"):
        return (nir - red) / (nir + red)


raw = nd(b["B08"] * 1e-4, b["B04"] * 1e-4)
offset = nd(b["B08"] * 1e-4 - 0.1, b["B04"] * 1e-4 - 0.1)
for code, name in [(4, "vegetation"), (5, "bare soil"), (6, "water"), (3, "cloud shadow"), ((8, 9), "cloud")]:
    m = np.isin(scl10, code)
    print(f"{name:13} {m.mean():6.1%} of the window, median NDVI {np.nanmedian(raw[m]):+.3f}")
print(f"vegetation with the offset applied again: {np.nanmedian(offset[scl10 == 4]):.3f}")
veg = scl10 == 4
print(f"vegetation reflectance: red {np.median(b['B04'][veg]) * 1e-4:.3f}, nir {np.median(b['B08'][veg]) * 1e-4:.3f}")
blue, red, nir = (b[k][veg] * 1e-4 for k in ("B02", "B04", "B08"))
with np.errstate(all="ignore"):
    evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)
ndvi_veg = raw[veg]
print("vegetation, 5th-95th percentile: NDVI", np.nanpercentile(ndvi_veg, [5, 95]).round(3),
      "| EVI", np.nanpercentile(evi, [5, 95]).round(3))

bad = np.isin(scl10, UNUSABLE)
print(f"usable: {1 - bad.mean():.1%}; after two pixels of dilation {1 - binary_dilation(bad, iterations=2).mean():.1%}")
water median reflectance: red  DN*1e-4 0.008, DN*1e-4-0.1 -0.092
vegetation     36.8% of the window, median NDVI +0.726
bare soil       2.6% of the window, median NDVI +0.362
water           1.2% of the window, median NDVI -0.077
cloud shadow   23.2% of the window, median NDVI +0.623
cloud          31.2% of the window, median NDVI +0.154
vegetation with the offset applied again: 1.362
vegetation reflectance: red 0.057, nir 0.362
vegetation, 5th-95th percentile: NDVI [0.522 0.852] | EVI [0.321 0.727]
usable: 42.7%; after two pixels of dilation 39.7%

What this scene says about the method

Over Snowdonia, vegetation reflected 0.057 in red and 0.362 in near-infrared, a factor of six — the contrast NDVI is built on, explained in vegetation indices explained. Dividing by the sum cancels anything that scales both bands together, which is what makes NDVI comparable across a hilly scene, and it is also why the mask matters more than the formula: shadow darkens both bands alike, so a shadowed canopy still read 0.623 against 0.726 in sun. Once an unmasked NDVI is averaged into a field statistic, the fact that part of the field was in shadow is gone.

The same ratio flattens over closed canopy. Here NDVI's 5th-to-95th percentile range over vegetation was 0.522 to 0.852, pressed against its ceiling, while EVI on the same pixels spanned 0.321 to 0.727 (Example 4) (NDVI saturates in dense canopy). And one range check catches almost every way to compute it wrongly — an offset applied twice, unsigned wraparound, mixed scaling — because each pushes values outside −1 to 1; telling them apart is NDVI values out of range.

Four NDVI failures and the single range assertion that catches three of them. two comparisons catch three of the four failuresfailurewhat the range check seesoffset applied twicevalues beyond ±1unsigned wraparoundvalues pinned near ±1mixed scaling between bandsa flat 0.999 everywherethe wrong band readnothing — still in rangeassert −1.001 ≤ ndvi ≤ 1.001 costs nothing and catches the common three spatialworkflow.io
One assertion, three of four failures. The fourth needs a class-median check.

Scales, offsets and nodata zero

  • 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.726. But cloud shadow reached 0.623, 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.