Compute NDVI from Sentinel-2 imagery with Python

The Normalized Difference Vegetation Index, or NDVI, is probably the most widely used number in remote sensing. It compares how much near-infrared and red light a surface reflects. Healthy plants absorb most red light for photosynthesis and bounce back a great deal of near-infrared, so they score high. Bare soil, water and buildings score low.

In this tutorial you will compute NDVI from a Sentinel-2 image, remove clouds and shadows, save the result as a GeoTIFF and draw it as a map. The whole script is under 50 lines. You will also learn how to read the result sensibly, because a number you cannot interpret is not much use.

What you need

• Python 3.9 or newer.
• Three single-band GeoTIFF files from the same Sentinel-2 Level-2A scene: B04 (red, 10 m), B08 (near-infrared, 10 m) and SCL (scene classification, 20 m).
• The packages rasterio, numpy and matplotlib.

pip install rasterio numpy matplotlib

Sentinel-2 images are free. You can download them from the Copernicus Data Space browser or fetch them from a public STAC catalogue such as Earth Search. Choose Level-2A products, which are already corrected for the atmosphere, and pick a scene with low cloud cover over your area. If you download whole scene archives, each band is a separate JP2 or TIFF file inside the folder named for the resolution.

Step 1: read the bands

We read red and near-infrared as 32-bit floating point numbers. This matters because the files store whole numbers, and dividing whole numbers can silently round your result. The SCL layer has half the resolution of the other bands, so we ask rasterio to resample it to match using nearest-neighbour, which keeps the class codes intact.

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

with rasterio.open("B04.tif") as red_src, rasterio.open("B08.tif") as nir_src:
    red = red_src.read(1).astype("float32")
    nir = nir_src.read(1).astype("float32")
    profile = red_src.profile

with rasterio.open("SCL.tif") as scl_src:
    # SCL is 20 m; resample it onto the 10 m grid without blending class codes
    scl = scl_src.read(1, out_shape=red.shape, resampling=Resampling.nearest)

Step 2: calculate NDVI

The formula is one line. The tiny extra value in the denominator avoids dividing by zero over pixels where both bands are empty, such as no-data edges.

ndvi = (nir - red) / (nir + red + 1e-6)

Because it is a ratio, NDVI always falls between −1 and 1, and it is less affected by overall brightness than raw band values. That is why it works well for comparing different dates and places.

Step 3: mask clouds and shadows

Clouds are bright in both bands, so they produce NDVI values near zero that look like bare ground. If you leave them in, your averages and maps will be wrong. The SCL layer labels every pixel with a class, which lets you remove the bad ones.

• 0 no data, 1 saturated or defective pixels
• 3 cloud shadow
• 8 and 9 cloud (medium and high probability), 10 thin cirrus
• 11 snow or ice
• Useful classes to keep: 4 vegetation, 5 not vegetated, 6 water, 7 unclassified

bad = np.isin(scl, [0, 1, 3, 8, 9, 10, 11])
ndvi = np.where(bad, np.nan, ndvi).astype("float32")

print(f"Masked {bad.mean():.1%} of pixels")

Step 4: save the result

We copy the metadata of the red band, so the output keeps the same location and projection, and write a single band of 32-bit floats. Setting the no-data value to NaN tells GIS software to treat masked pixels as empty.

profile.update(dtype="float32", count=1, nodata=np.nan)

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

You can now open ndvi.tif in QGIS or any other GIS program, overlay it on a basemap and style it with a green colour ramp.

Step 5: draw a quick map

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 8))
im = ax.imshow(ndvi, cmap="RdYlGn", vmin=-0.2, vmax=0.9)
ax.set_title("NDVI")
ax.axis("off")
fig.colorbar(im, ax=ax, shrink=0.7, label="NDVI")
plt.savefig("ndvi.png", dpi=150, bbox_inches="tight")

How to read the result

• Below zero: water, snow or clouds that slipped through the mask.
• About 0 to 0.2: bare soil, rock, sand or built-up land.
• About 0.2 to 0.5: sparse vegetation, grassland or young crops.
• Above about 0.5: dense, healthy vegetation such as mature crops and forest.

Treat these bands as a guide, not a rule. The same crop looks different at each growth stage, and a pale, dry grassland in winter can score lower than a healthy crop in spring while being perfectly normal for its season. The most reliable way to use NDVI is comparison: this field against its neighbours, or this month against the same month last year.

Common problems

• Everything is near zero or looks washed out: check you used red and near-infrared (B04 and B08), not two other bands.
• Values look odd or banded: confirm you converted to float before dividing.
• The SCL array has a different shape: you skipped the resampling step.
• The map is mostly empty: your scene is too cloudy. Pick another date.

Where to go next

• Compute the average NDVI inside a field boundary, using a polygon from a shapefile.
• Repeat the process for every clear image across a season to draw a growth curve.
• Try the water index NDWI, which uses green and near-infrared, with the same pattern.
• Feed NDVI as an extra feature into a land cover classifier, which usually improves accuracy.

Going further: an NDVI time series

A single date tells you how a field looks today. A series of dates tells you how it is behaving, which is far more useful for judging crop health. The loop below computes the average NDVI of a scene for every date in a folder. It assumes files named like 20260612_B04.tif and 20260612_B08.tif.

import glob
import numpy as np
import rasterio

means = {}
for red_path in sorted(glob.glob("scenes/*_B04.tif")):
    date = red_path.split("/")[-1][:8]
    nir_path = red_path.replace("_B04", "_B08")
    with rasterio.open(red_path) as r, rasterio.open(nir_path) as n:
        red = r.read(1).astype("float32")
        nir = n.read(1).astype("float32")
    ndvi = (nir - red) / (nir + red + 1e-6)
    means[date] = float(np.nanmean(ndvi))

for date, value in means.items():
    print(date, round(value, 3))

Plot these values against the dates and you have a growth curve. In a healthy crop, NDVI rises through the early season, peaks and then falls at ripening and harvest. A field whose curve lags behind its neighbours is worth a visit. Remember to apply the cloud mask from step 3 to each date, or one cloudy scene will pull a dip into the line.

Where NDVI falls short

• Saturation: in dense, leafy canopies, NDVI levels off and stops distinguishing between dense and very dense vegetation. Indices such as EVI are designed to keep their sensitivity there.
• Soil background: in sparse crops, bright or dark soil shifts the values. SAVI includes a soil adjustment for this.
• Chlorophyll sensitivity: red-edge indices, such as NDRE, respond earlier to nitrogen stress in many crops.
• Season and species: the same value means different things for rice, maize and forest, so compare like with like.

Frequently asked questions

Why is my NDVI always between 0 and 0.3?

Check that you used red and near-infrared bands, B04 and B08, and not two visible bands. A mostly urban or bare scene will also score low.

Do I need Level-2A data?

It is strongly recommended. Level-2A is corrected for the atmosphere, so values are comparable between dates. Level-1C values vary with haze and sun angle.

Can I compute NDVI in QGIS instead?

Yes. The raster calculator can apply the same formula, and it is a good way to check your Python result.

Why does the SCL band have a different size?

Sentinel-2 delivers different bands at 10, 20 and 60 metres. The scene classification is at 20 metres, so it must be resampled to 10 metres, using nearest neighbour to keep the class codes intact.

You now have a repeatable recipe for one of the most useful measurements in remote sensing. The same pattern of read, calculate, mask and save applies to almost every index you will meet.