NDVI from Sentinel-2: A Step-by-Step Satellite Imagery Workflow

NDVI is probably the first index most people calculate from satellite imagery. The formula is short, and that makes it easy to get a number out and easy to get a wrong number out. Clouds, scaling and mismatched pixel grids all change the result without raising an error.

This topic walks through one complete workflow: a Sentinel-2 scene in, a cloud-masked NDVI raster out, with the checks in between.

What is NDVI?

NDVI (Normalized Difference Vegetation Index) compares how much near-infrared (NIR) and red light a surface reflects. Healthy green vegetation absorbs most red light and reflects a lot of NIR, so the two bands differ strongly over plants.

NDVI = (NIR − Red) / (NIR + Red)

Why this formula works

The formula follows from how different surfaces treat red and near-infrared light.

• Red light is absorbed by chlorophyll, the green pigment in leaves, which uses it for photosynthesis. A healthy leaf reflects very little red.
• Near-infrared light is not used by the plant. The internal cell structure of a healthy leaf scatters it, so a leaf reflects a lot of NIR.

A green plant is therefore dark in red and bright in NIR, and the gap between the two is large. Other surfaces have a smaller gap, or the gap runs the other way. Subtracting red from NIR measures that gap, and dividing by their sum scales it so that brighter or darker lighting changes the result less.

The reflectance values below are rounded, illustrative figures to show the pattern. Real values vary with the plant, the soil, the season and the sensor.

Two things follow from the table. First, a value near zero is ambiguous: bare soil and cloud can both give it, which is why you mask clouds before you interpret the index. Second, a low NDVI over a field tells you the plants are absent, sparse or stressed, but not which of those it is.

The result always falls between −1 and 1. As a rough guide, water and clouds give values near or below zero, bare soil and built-up areas give low positive values, and dense green vegetation gives high values. The exact thresholds depend on the sensor, the season and the surface, so treat any cut-off as something to check against your own area.

For Sentinel-2, the red band is B04 and the NIR band is B08. Both have a 10 m pixel size.

Why it matters

NDVI turns two bands into one layer that is easy to compare across dates. It is used to follow crop growth, spot vegetation stress, and measure how much green cover changed between two seasons or years.

It is also easy to misuse. A single cloudy pixel gives a meaningless value, and NDVI saturates over very dense vegetation, so it separates "moderately green" from "very green" less well than it separates "bare" from "green".

How it works

Satellite data goes through several stages before you calculate anything, and it helps to keep them separate.

• Level-1C is top-of-atmosphere reflectance. The light still includes the effect of the atmosphere.
• Level-2A is surface reflectance, after atmospheric correction. This is the usual starting point for NDVI.
• Orthorectification removes the distortion caused by terrain and viewing angle, so the image sits correctly on the map. Sentinel-2 Level-1C and Level-2A products are orthorectified, but you should not assume this for every satellite product. Check the documentation of the data you use.
• Cloud masking removes pixels you cannot trust. Sentinel-2 Level-2A includes a Scene Classification Layer (SCL) that labels each pixel as vegetation, water, cloud, cloud shadow and so on.

One detail trips people up: Sentinel-2 bands come at different resolutions. B04 and B08 are 10 m, while the SCL is 20 m. The arrays have different shapes, so you cannot combine them until you resample one to match the other.

Other indices and why their formulas work

NDVI is one of a family. Each index picks two or more bands where the surface you care about behaves differently from everything else, and combines them so that difference stands out. The reasoning is the same every time: find a band the target absorbs and a band it reflects.

These Sentinel-2 bands are used below.

Bands at 20 m must be resampled to the 10 m grid before you combine them with a 10 m band, in the same way as the SCL in the code later in this topic.

NDWI: open water

Formula: (Green − NIR) / (Green + NIR), which is (B03 − B08) / (B03 + B08).

Why it works: water absorbs NIR strongly and reflects a little green, while vegetation and soil reflect more NIR than green. So water gives positive values and land gives negative ones. This is the version by McFeeters; another index with the same name uses NIR and SWIR to measure water in leaves, so always state the formula you used.

Watch out for: built-up areas can also come out positive, because some surfaces reflect more green than NIR.

MNDWI: water in built-up areas

Formula: (Green − SWIR1) / (Green + SWIR1), which is (B03 − B11) / (B03 + B11).

Why it works: water absorbs SWIR almost completely, much more than it absorbs green. Buildings and bare soil reflect more SWIR than green, so they go negative. Replacing NIR with SWIR separates water from built-up land better than NDWI does.

NDBI: built-up land

Formula: (SWIR1 − NIR) / (SWIR1 + NIR), which is (B11 − B08) / (B11 + B08).

Why it works: built-up surfaces reflect somewhat more SWIR than NIR, while vegetation does the opposite, because leaf water and structure reflect NIR strongly. Positive values lean towards built-up.

Watch out for: dry bare soil also reflects more SWIR than NIR and looks similar. NDBI on its own does not cleanly separate built-up from bare soil.

NDMI: moisture in vegetation

Formula: (NIR − SWIR1) / (NIR + SWIR1), which is (B08 − B11) / (B08 + B11).

Why it works: water in leaves absorbs SWIR. A well-watered canopy reflects less SWIR and gives a higher value, and a dry or stressed one reflects more SWIR and gives a lower value. It can pick up water stress before the canopy turns less green, which is when NDVI starts to fall.

NBR: burned areas

Formula: (NIR − SWIR2) / (NIR + SWIR2), which is (B08 − B12) / (B08 + B12).

Why it works: fire removes leaves, so NIR reflectance falls, and it exposes dry soil and ash, so SWIR2 reflectance rises. A healthy canopy gives a high NBR and a recently burned surface gives a low or negative one. To measure the burn, calculate NBR before and after the fire and subtract: dNBR = NBR(before) − NBR(after).

SAVI: sparse vegetation over soil

Formula: (NIR − Red) / (NIR + Red + L) × (1 + L), where L is a soil adjustment factor. L = 0.5 is the usual default.

Why it works: where plants are sparse, the soil between them reflects light and shifts NDVI up or down depending on how bright the soil is. The L term reduces that effect. With L = 0 the formula becomes NDVI. Use it for early-season crops and arid areas.

EVI: dense vegetation

Formula: 2.5 × (NIR − Red) / (NIR + 6 × Red − 7.5 × Blue + 1).

Why it works: NDVI levels off over very dense canopies, so thick forest and very thick crops look alike. EVI keeps responding there. The blue band is used to correct for haze (atmospheric scattering), and the 1 in the denominator is a canopy background adjustment. EVI expects reflectance on a 0 to 1 scale, so apply the scale factor first.

NDRE: chlorophyll in dense crops

Formula: (NIR − Red edge) / (NIR + Red edge), which is (B08 − B05) / (B08 + B05).

Why it works: the red edge is the steep rise in reflectance between red and NIR. Its position and level respond to chlorophyll content, and red edge light penetrates deeper into a dense canopy than red light. NDRE is useful when NDVI has already levelled off, such as mid-season crops.

Comparing them on the same surfaces

The table uses one set of rounded, illustrative reflectances for each surface and calculates every index from them. The values only show the pattern. Real values depend on the sensor, the date, the soil and the vegetation type.

The reflectances behind it, from blue to SWIR2 (B02, B03, B04, B05, B08, B11, B12):

Read each column to see which surface the index singles out: NDWI and MNDWI peak on water, NBR drops on the burned area, NDBI is positive on built-up land but also on bare soil and the burned area, and NDVI, SAVI, EVI and NDRE peak on healthy vegetation. Each index is built for one target, and it can be fooled by surfaces that behave similarly in the bands it uses. NDMI and NBR are meant for vegetated land, so ignore their water values and mask water first.

Calculating any normalized difference index

NDVI, NDWI, MNDWI, NDBI, NDMI, NBR and NDRE all share one shape: (A − B) / (A + B). One small function covers all of them.

# One function for every normalized difference index.
# a and b are NumPy arrays of reflectance on the same grid.

import numpy as np

def normalized_difference(a, b):
    # Step 1: Add the two bands. Where the sum is zero, leave the result empty.
    total = a + b
    result = np.full(a.shape, np.nan, dtype="float32")

    # Step 2: Divide only where the sum is not zero
    np.divide(a - b, total, out=result, where=total != 0)
    return result

# Step 3: Examples (arrays already read and scaled to reflectance)
# ndvi = normalized_difference(nir, red)
# ndwi = normalized_difference(green, nir)
# mndwi = normalized_difference(green, swir1)
# ndbi = normalized_difference(swir1, nir)
# ndmi = normalized_difference(nir, swir1)
# nbr = normalized_difference(nir, swir2)
# ndre = normalized_difference(nir, red_edge)

SAVI and EVI have extra terms, so write each as its own line from the formulas above.

Practical GIS example

Take an agricultural area in Maharashtra and compare a cloud-free scene from the growing season with one from the dry season. The aim is a cloud-masked NDVI raster for each date, which you can then subtract or classify.

Getting the Sentinel-2 files

The guide How to download Sentinel-2 imagery shows two routes: the Copernicus Browser, where you click and download a zip, and a Python script that reads only your area. Choose Level-2A, preview the scene over your own area, and check the cloud before you download.

If you downloaded the zip, the bands you need are inside the unzipped .SAFE folder, in GRANULE, in a subfolder named after the tile:

• IMG_DATA/R10m/..._B04_10m.jp2 is the red band.
• IMG_DATA/R10m/..._B08_10m.jp2 is the NIR band.
• IMG_DATA/R20m/..._SCL_20m.jp2 is the scene classification layer.

Copy these three into your working folder. They are JPEG 2000 (.jp2) files, which Rasterio opens like any other raster, as long as your GDAL build includes the JPEG 2000 driver (the standard Rasterio installs do).

Workflow

1. Download a Level-2A scene with the guide above, and copy B04, B08 and SCL into your working folder.
2. Check that the files share a CRS. Bands from one Sentinel-2 tile are in the same UTM zone. If you combine tiles, confirm the CRS of each (see the blog post on CRS errors that can ruin a GIS analysis).
3. Read the bands and convert to reflectance. The files hold integer values, and you scale them to reflectance.
4. Resample the SCL to the 10 m grid and build a mask (the guide How to mask clouds in Sentinel-2 explains the classes and what each choice costs).
5. Calculate NDVI and set masked pixels to no-data.
6. Check the result visually and with a few numbers before using it.

# Cloud-masked NDVI from a Sentinel-2 Level-2A scene.
# Needs the three files you copied from the downloaded .SAFE folder.

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

# Step 0: Put your own file names here. The real names are long and include
# the tile and date, for example T43QCV_20240115T050109_B04_10m.jp2.
RED_FILE = "B04_10m.jp2"
NIR_FILE = "B08_10m.jp2"
SCL_FILE = "SCL_20m.jp2"

# Step 1: Scaling values. Level-2A stores reflectance x 10000.
# Scenes processed with baseline 04.00 or later (from January 2022) also
# carry an offset of -1000. Check your scene's metadata and set this to 0
# for older scenes.
SCALE = 10000.0
OFFSET = -1000.0

with rasterio.open(RED_FILE) as red_src, \
     rasterio.open(NIR_FILE) as nir_src, \
     rasterio.open(SCL_FILE) as scl_src:

    # Step 2: Read the red and NIR bands as decimal numbers
    red = (red_src.read(1).astype("float32") + OFFSET) / SCALE
    nir = (nir_src.read(1).astype("float32") + OFFSET) / SCALE

    # Step 3: Read the 20 m SCL onto the 10 m grid.
    # Nearest neighbour keeps class numbers intact (no averaging of classes).
    scl = scl_src.read(
        1,
        out_shape=(red_src.height, red_src.width),
        resampling=Resampling.nearest,
    )

    profile = red_src.profile

# Step 4: Mark the pixels we do not trust.
# SCL classes: 1 saturated, 3 cloud shadow, 8 cloud (medium),
# 9 cloud (high), 10 thin cirrus
bad = np.isin(scl, [1, 3, 8, 9, 10])

# Step 5: Calculate NDVI. Where red + NIR is zero, leave the result empty.
total = nir + red
ndvi = np.full(red.shape, np.nan, dtype="float32")
np.divide(nir - red, total, out=ndvi, where=total != 0)

# Step 6: Remove masked pixels
ndvi[bad] = np.nan

# Step 7: Save the result as a single-band float raster
profile.update(dtype="float32", count=1, nodata=np.nan)
with rasterio.open("ndvi.tif", "w", **profile) as dst:
    dst.write(ndvi, 1)

# Step 8: Quick check of the values
print("Min:", np.nanmin(ndvi), "Max:", np.nanmax(ndvi))
print("Masked share:", round(bad.mean() * 100, 1), "%")

The script does the arithmetic on whole arrays at once, which is why it stays short. The where= argument stops a divide-by-zero warning on empty pixels.

Common mistakes

• Skipping the cloud mask. Cloud pixels have NDVI near zero or below, which looks like bare ground or water in the result.
• Forgetting the scale factor or offset. The ratio is not the same with and without the offset, so NDVI values drift away from other people's results.
• Mixing resolutions. Combining a 10 m band and a 20 m layer without resampling either fails with a shape error, or worse, gets silently resized by a tool in a way you did not choose.
• Resampling class layers with bilinear or cubic methods. Averaging class numbers produces class numbers that do not exist. Use nearest neighbour for categories.
• Comparing dates without checking the season. An NDVI drop between two dates can be a harvest, a dry spell or a different stage of the crop calendar. The index does not tell you which.
• Reading NDVI as a measurement of yield or health. It is an index of greenness. Link it to ground observations before drawing conclusions.

Key takeaways

• NDVI = (NIR − Red) / (NIR + Red); for Sentinel-2 that is (B08 − B04) / (B08 + B04).
• Use Level-2A surface reflectance, and apply the scale factor and any offset from the metadata.
• Mask clouds and shadows with the SCL before calculating, and resample the SCL with nearest neighbour.
• Keep preprocessing steps separate in your notes: atmospheric correction, orthorectification and cloud masking are different operations.
• Check the output visually and against a known area before comparing dates.