NDVI stops responding in a dense canopy

Problem statement

The NDVI map of a field in mid-June is one shade of green. Fields that differ visibly in biomass have the same value, the within-field variation has vanished, and a time series through the peak is a flat line for three weeks.

This is saturation, and it is not a data problem. Chlorophyll absorbs essentially all incident red light once the leaf area index passes about three, so additional leaves change red reflectance almost not at all while near-infrared keeps rising slowly. The normalised difference asymptotes; the canopy does not.

Measured on a nearly cloud-free Sentinel-2 scene over Flevoland farmland on 12 June 2025: 498,177 of 982,656 vegetation pixels — 50.7% — were above NDVI 0.8, and across those pixels NDVI varied over a range of 0.199 with a coefficient of variation of 0.033.

Quick answer

Use a red-edge index over the saturated period:

import numpy as np

ndvi = (nir - red) / (nir + red)
ndre = (nir - rededge1) / (nir + rededge1)          # Sentinel-2 B8 and B5
evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)
vegetation = scl == 4                                # scene classification: vegetation

dense = vegetation & (ndvi > 0.8)
for name, arr in (("NDVI", ndvi), ("NDRE", ndre), ("EVI", evi)):
    v = arr[dense]
    print(f"{name:5} range {v.max()-v.min():.3f}  sd {v.std():.4f}  "
          f"CV {v.std()/abs(v.mean()):.3f}")

NDVI  range 0.199  sd 0.0282  CV 0.033
NDRE  range 0.521  sd 0.0608  CV 0.093
EVI   range 0.771  sd 0.0849  CV 0.113

NDRE gives 2.8 times the relative variation on identical ground. Across 379 arable parcels whose median NDVI exceeded 0.8, NDVI spanned 0.123 between fields and NDRE spanned 0.231.

Confirm saturation, measure its season, switch index

1. Confirm it is saturation and not a masking problem

Saturation shows as a high, narrow distribution: most vegetation pixels above 0.8 with a small spread. A masking problem shows as a bimodal distribution with a cloud or shadow population.

2. Measure how much of the season is affected

The field-median NDVI on the reference barley series stayed between 0.923 and 0.959 on ten dates from 26 April to 13 June — seven weeks with a total range of 0.036, while NDRE ranged over 0.157. That whole period contributes almost nothing to any analysis based on NDVI.

3. Switch to a red-edge index for the saturated window

Sentinel-2 band 5 at 705 nm sits on the steep part of the vegetation reflectance curve, where absorption is partial and further chlorophyll still changes reflectance.

4. Use the right red-edge band

Band 7 at 783 nm is already on the NIR plateau, so the normalised difference against NIR collapses towards zero. On the same scene, the median NDRE over the dense pixels was +0.645 with band 5 and +0.020 with band 7.

5. Consider EVI if you have a reliable blue band

EVI performed comparably to NDRE on the dense pixels — a range of 0.771 and a CV of 0.113 — and it is available on sensors with no red edge. Its weakness is that it depends on the blue band, whose atmospheric correction is the least reliable.

6. Fix the downstream consequences

Saturation does not only flatten the map. It flattens the peak of a phenology curve, so the peak date becomes arbitrary; it removes the boundaries a field delineation depends on; and it removes the feature that separates two crops at full canopy.

7. Keep NDVI for the parts of the season where it works

Emergence, early growth and senescence are all well below saturation, and NDVI's cross-sensor consistency there is worth keeping. The fix is to use both, not to abandon one.

Quantifying saturation on a real scene

The data used on this page

One real field: the largest spring-barley parcel (Gerst, zomer-) declared in 2025 in the Dutch parcel register BRP Gewaspercelen within 5.55–5.75° E, 52.48–52.56° N, east of Dronten in Flevoland. Its series comes from every Sentinel-2 L2A scene in Earth Search over the parcel from 1 March to 30 September 2025. For each date the setup records the share of the parcel that the scene classification calls vegetation or bare soil (clear_frac), and the median NDVI (B8, B4) and NDRE (B8, B5) over those pixels. Where two scenes share a date, the clearer one is kept. Dates with no clear pixel are dropped. Reading about a hundred scenes takes several minutes the first time; after that season.csv is reused.

import os

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.features import geometry_mask
from rasterio.windows import from_bounds

os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR",
                  CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif")

if not os.path.exists("field.gpkg"):
    url = "https://api.pdok.nl/rvo/gewaspercelen/ogc/v1/collections/brpgewas/items"
    feats, params = [], {"f": "json", "bbox": "5.55,52.48,5.75,52.56", "limit": 1000}
    while url:
        page = requests.get(url, params=params, timeout=120).json()
        feats += page["features"]
        url, params = next((l["href"] for l in page["links"] if l["rel"] == "next"), None), None
    brp = gpd.GeoDataFrame.from_features(feats, crs=4326)
    brp["ha"] = brp.to_crs(28992).area / 1e4
    brp.to_file("brp.gpkg")
    brp[brp.gewas.str.startswith("Gerst, zomer")].nlargest(1, "ha").to_file("field.gpkg")
field = gpd.read_file("field.gpkg")


def window(href, bounds, shape=None):
    with rasterio.open(href) as src:
        win = from_bounds(*bounds, transform=src.transform)
        extra = {"out_shape": shape} if shape else {}
        return src.read(1, window=win, boundless=True, **extra).astype("float32"), src.window_transform(win)


if not os.path.exists("season.csv"):
    body = {"collections": ["sentinel-2-l2a"], "intersects": field.geometry.iloc[0].__geo_interface__,
            "datetime": "2025-03-01T00:00:00Z/2025-09-30T23:59:59Z", "limit": 100}
    items = []
    while body:
        page = requests.post("https://earth-search.aws.element84.com/v1/search", json=body, timeout=120).json()
        items += page["features"]
        body = next((l["body"] for l in page["links"] if l["rel"] == "next"), None)
    rows = []
    for day in sorted({i["properties"]["datetime"][:10] for i in items}):
        best = None
        for item in (i for i in items if i["properties"]["datetime"][:10] == day):
            a = item["assets"]
            utm = field.to_crs(item["properties"]["proj:epsg"])
            b = utm.total_bounds + np.array([-20, -20, 20, 20])
            red, transform = window(a["red"]["href"], b)
            inside = ~geometry_mask(utm.geometry, red.shape, transform)
            scl, _ = window(a["scl"]["href"], b, red.shape)
            clear = inside & np.isin(scl, [4, 5])          # vegetation, bare soil
            if best is None or clear.sum() > best[1].sum():
                best = (item, clear, inside, b, red)
        item, clear, inside, b, red = best
        if not clear.any():
            continue
        # DN / 10,000 is reflectance: no offset to remove, whatever the item's
        # boa_offset_applied flag says (water reads 1.5–2.6% on 7 and 10 March either way)
        nir, _ = window(item["assets"]["nir"]["href"], b, red.shape)
        re1, _ = window(item["assets"]["rededge1"]["href"], b, red.shape)
        r, n, e = red[clear], nir[clear], re1[clear]
        rows.append({"date": day, "scene": item["id"],
                     "clear_frac": round(float(clear.sum() / inside.sum()), 3),
                     "ndvi_med": round(float(np.median((n - r) / (n + r))), 3),
                     "ndre_med": round(float(np.median((n - e) / (n + e))), 3)})
    pd.DataFrame(rows).to_csv("season.csv", index=False)

season = pd.read_csv("season.csv", parse_dates=["date"])
print(f"{field.gewas.iloc[0]} parcel, {field.ha.iloc[0]:.1f} ha; "
      f"{len(season)} dates with a clear view of it, {season.date.min().date()} to {season.date.max().date()}")

Gerst, zomer- parcel, 11.4 ha; 47 dates with a clear view of it, 2025-03-07 to 2025-09-30

The scene-wide figures use one whole Sentinel-2 scene over the same 5.55–5.75° E, 52.48–52.56° N box: S2A_31UFU_20250612_0_L2A, 12 June 2025, a nearly cloud-free tile. The setup reads its blue, red, NIR, band 5 and band 7 and the scene classification, all resampled to the 10 m grid, and converts DN to reflectance by dividing by 10,000.

from affine import Affine

SCENE = "S2A_31UFU_20250612_0_L2A"
if not os.path.exists("scene_20250612.npz"):
    item = requests.get("https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/"
                        + SCENE, timeout=60).json()
    corners = gpd.GeoSeries.from_xy([5.55, 5.75], [52.48, 52.56], crs=4326).to_crs(item["properties"]["proj:epsg"])
    b = (corners.x.min(), corners.y.min(), corners.x.max(), corners.y.max())
    bands = {}
    bands["red"], transform = window(item["assets"]["red"]["href"], b)
    for name in ("blue", "nir", "rededge1", "rededge3", "scl"):
        bands[name], _ = window(item["assets"][name]["href"], b, bands["red"].shape)
    np.savez_compressed("scene_20250612.npz", transform=np.array(transform)[:6],
                        epsg=item["properties"]["proj:epsg"], **bands)

z = np.load("scene_20250612.npz")
blue, red, nir, rededge1, rededge3 = (z[k] / 10_000 for k in ("blue", "red", "nir", "rededge1", "rededge3"))
scl, transform, epsg = z["scl"], Affine(*z["transform"]), int(z["epsg"])
brp = gpd.read_file("brp.gpkg")
print(f"{SCENE}: {red.shape[1]} x {red.shape[0]} pixels, {(scl == 4).mean():.1%} vegetation; "
      f"{len(brp):,} BRP parcels in the box")

S2A_31UFU_20250612_0_L2A: 1324 x 939 pixels, 79.0% vegetation; 6,464 BRP parcels in the box

Example 1 — quantify the saturation on your own scene

import numpy as np

def saturation_report(ndvi, ndre, evi, vegetation_mask, threshold=0.8):
    veg = vegetation_mask & np.isfinite(ndvi) & np.isfinite(ndre)
    dense = veg & (ndvi > threshold)
    out = {"vegetation_pixels": int(veg.sum()),
           "dense_pixels": int(dense.sum()),
           "dense_share": round(float(dense.sum() / veg.sum()), 3)}
    for name, arr in (("ndvi", ndvi), ("ndre", ndre), ("evi", evi)):
        v = arr[dense]
        v = v[np.isfinite(v)]
        out[f"{name}_range"] = round(float(v.max() - v.min()), 3)
        out[f"{name}_sd"] = round(float(v.std()), 4)
        out[f"{name}_cv"] = round(float(v.std() / abs(v.mean())), 3)
    out["ndre_advantage"] = round(out["ndre_cv"] / out["ndvi_cv"], 2)
    return out

print(saturation_report(ndvi, ndre, evi, scl == 4))

dense = (scl == 4) & (ndvi > 0.8)
ndre7 = (nir - rededge3) / (nir + rededge3)          # band 7, 783 nm, against B8
print(f"median NDRE on dense pixels: band 5 {np.median(ndre[dense]):+.3f}, "
      f"band 7 {np.median(ndre7[dense]):+.3f}")

{'vegetation_pixels': 982656, 'dense_pixels': 498177, 'dense_share': 0.507, 'ndvi_range': 0.199, 'ndvi_sd': 0.0282, 'ndvi_cv': 0.033, 'ndre_range': 0.521, 'ndre_sd': 0.0608, 'ndre_cv': 0.093, 'evi_range': 0.771, 'evi_sd': 0.0849, 'evi_cv': 0.113, 'ndre_advantage': 2.82}
median NDRE on dense pixels: band 5 +0.645, band 7 +0.020

Example 2 — separation between fields, which is what usually matters

import numpy as np, pandas as pd

def field_separation(field_table, threshold=0.8):
    hi = field_table[field_table["ndvi"] > threshold]
    print(f"{len(hi)} of {len(field_table)} fields above NDVI {threshold}")
    for col in ("ndvi", "ndre"):
        print(f"  {col.upper():5} spread between fields "
              f"{hi[col].max() - hi[col].min():.3f}, sd {hi[col].std():.4f}")
    print(f"  within-field CV: NDVI median {hi['ndvi_cv'].median():.3f}, "
          f"NDRE median {hi['ndre_cv'].median():.3f}")
    return hi


from rasterio.features import rasterize

parcels = brp[(brp.category == "Bouwland") & (brp.ha >= 5)].to_crs(epsg).reset_index(drop=True)
ids = rasterize(((g, i + 1) for i, g in enumerate(parcels.geometry)),
                out_shape=scl.shape, transform=transform, fill=0, dtype="int32")
rows = []
for i in range(1, len(parcels) + 1):
    inside = ids == i
    veg = inside & (scl == 4)
    if not inside.any() or veg.sum() < 0.9 * inside.sum():      # skip bare or cloudy parcels
        continue
    a, b = ndvi[veg], ndre[veg]
    rows.append({"ndvi": np.median(a), "ndre": np.median(b),
                 "ndvi_cv": a.std() / a.mean(), "ndre_cv": b.std() / b.mean()})
field_table = pd.DataFrame(rows)
print(f"{len(parcels)} arable parcels of 5 ha or more; {len(field_table)} at least 90% vegetation")
hi = field_separation(field_table)

754 arable parcels of 5 ha or more; 544 at least 90% vegetation
379 of 544 fields above NDVI 0.8
  NDVI  spread between fields 0.123, sd 0.0250
  NDRE  spread between fields 0.231, sd 0.0576
  within-field CV: NDVI median 0.057, NDRE median 0.069

Between fields, NDRE spreads 1.9 times wider (0.231 against 0.123). Within fields the two are similar — median CVs of 0.057 and 0.069 — which is worth knowing: saturation costs most when comparing fields, and less when mapping variation inside one.

Example 3 — a season stack that switches index by growth stage

import numpy as np, pandas as pd

def blended_index(dates, ndvi, ndre, saturate_at=0.75):
    """Use NDVI while it is informative and NDRE once it saturates."""
    out, source = [], []
    for d, a, b in zip(dates, ndvi, ndre):
        if np.nanmedian(a) > saturate_at:
            out.append(b)
            source.append("ndre")
        else:
            out.append(a)
            source.append("ndvi")
    return np.array(out), pd.Series(source, index=pd.DatetimeIndex(dates))

blend, src = blended_index(season.date, season.ndvi_med, season.ndre_med)
print(src.value_counts().to_string())

plateau = season[(season.clear_frac >= 0.6) & (season.ndvi_med > 0.9)]
print(f"NDVI above 0.9 on {len(plateau)} dates, {plateau.date.min():%d %b} to {plateau.date.max():%d %b}: "
      f"NDVI range {plateau.ndvi_med.max() - plateau.ndvi_med.min():.3f} "
      f"({plateau.ndvi_med.min():.3f}–{plateau.ndvi_med.max():.3f}), "
      f"NDRE range {plateau.ndre_med.max() - plateau.ndre_med.min():.3f}")

ndvi    30
ndre    17
NDVI above 0.9 on 10 dates, 26 Apr to 13 Jun: NDVI range 0.036 (0.923–0.959), NDRE range 0.157

On the barley series, 17 of the 47 dates had a field median above 0.75 and switched to NDRE. Blending is a pragmatic choice rather than a principled one — the two indices are on different scales — so it suits a segmentation or a classifier feature far better than it suits a published time series. For a published series, report both.

Absorption is the limit; the red edge keeps responding

Why absorption is the limit

Red light is absorbed by chlorophyll. A canopy with a leaf area index of one already absorbs most of the red incident on it; by three, essentially all of it. Further leaves therefore cannot lower red reflectance any further, so the numerator nir − red grows only as fast as NIR, and the denominator grows at the same rate. The ratio flattens.

Why the red edge keeps responding

Between roughly 700 and 730 nm, chlorophyll absorption is falling steeply and is only partial even in a dense canopy. Additional chlorophyll still reduces reflectance measurably there, so a normalised difference against NIR keeps varying. That is also why red-edge indices track leaf nitrogen: the relationship is with chlorophyll concentration rather than merely with cover.

Why band 7 does not work

By 783 nm the leaf's absorption has essentially stopped and reflectance is on the NIR plateau, close to band 8's value. The normalised difference between two similar values is near zero and carries little information — hence the +0.020 median against +0.645 for band 5.

Why the downstream effects are larger than the map

A flat peak makes the peak date arbitrary, which propagates into every phenology metric. It removes the between-field contrast that a delineation needs, causing adjacent fields to merge. And it removes the biomass differences that separate two cereals at full canopy, which is where a crop classifier makes most of its errors. The flat map is the least of it.

Thresholds, shadow and band resolutions

• The threshold is not universal. 0.8 suits temperate arable; check your own distribution.
• Shadow looks like dense canopy to every index. Mask it.
• Band resolutions differ. Red and NIR are 10 m; the red edges are 20 m.
• Not every sensor has a red edge. Landsat does not; use EVI there.
• EVI depends on the blue band, whose atmospheric correction is least reliable.
• Saturation is worse at nadir and varies with sun angle.
• NDVI is still right for emergence and senescence. Use both.
• State which index each figure used. They are not on the same scale.

Internal links

• Vegetation indices explained: NDVI, EVI, NDRE and when each fails — the whole family
• Crop phenology and growing seasons explained — the flat peak
• How to extract phenology metrics from an NDVI time series — the arbitrary peak date
• Delineated field boundaries merge fields or swallow roads — the missing boundaries
• A crop classifier confuses two crops every year — the missing separation
• Spectral indices explained — the wider context
• How to calculate NDVI in Python — the basic computation
• NDVI values are out of range — a different NDVI problem
• Urban cooling and green space explained — why trees and grass cool differently

FAQ

Why is my NDVI map all one colour in summer?

Saturation. On a real June scene, 50.7% of vegetation pixels were above 0.8 and varied over a range of only 0.199.

What should I use instead?

A red-edge index — NDRE from Sentinel-2 bands 8 and 5 — or EVI. On the same pixels NDRE varied with 2.8 times NDVI's coefficient of variation, and EVI with 3.4 times.

Which red-edge band?

Band 5 at 705 nm. Band 7 at 783 nm is already on the NIR plateau and gave a median NDRE of +0.020 against +0.645 for band 5.

Does saturation matter within a field or between fields?

Mostly between. Across 379 dense arable parcels, NDRE spread 1.9 times wider than NDVI; within fields the two were similar.

What if my sensor has no red-edge band?

Use EVI, which was comparable to NDRE on the dense pixels — at the cost of depending on the blue band's atmospheric correction.

Should I stop using NDVI?

No. It works well at emergence and senescence and has unmatched cross-sensor consistency. Use it where it responds and something else where it does not.