Spectral Indices Explained: NDVI, NDWI and What They Assume

Problem statement

A spectral index turns several bands into one number, and that number gets treated as a measurement of something real: vegetation, water, built-up area, snow. It is not. It is a ratio, and a ratio is a contrast, not a quantity.

Two facts make this worth being careful about. First, indices that claim to measure the same thing disagree — measured over the same 466,954 vegetation pixels, NDVI and EVI correlate at only 0.612. Second, indices are easily fooled: cloud shadow over vegetation in the scene below has a median NDVI of 0.622, comfortably inside the range people call "healthy vegetation".

Quick answer

Every normalised-difference index is the same arithmetic with different bands:

import numpy as np


def normalised_difference(a, b):
    """(a - b) / (a + b), with the zero denominator handled."""
    a = np.asarray(a, dtype="float32")
    b = np.asarray(b, dtype="float32")
    total = a + b
    return np.where(np.abs(total) < 1e-6, np.nan, (a - b) / total)


ndvi = normalised_difference(nir, red)       # vegetation vs bare
ndwi = normalised_difference(green, nir)     # open water
mndwi = normalised_difference(green, swir)   # water against built-up
ndbi = normalised_difference(swir, nir)      # built-up
print(f"NDVI median {np.nanmedian(ndvi):.3f}, NDWI {np.nanmedian(ndwi):.3f}, "
      f"MNDWI {np.nanmedian(mndwi):.3f}, NDBI {np.nanmedian(ndbi):.3f} over every cell")

NDVI median 0.617, NDWI -0.581, MNDWI -0.439, NDBI -0.156 over every cell

Measured medians per surface class over a Sentinel-2 scene in Snowdonia:

What the division is for, what each index contrasts

1. Understand what the division is for

The numerator does the work: nir - red is large for vegetation because leaves absorb red and scatter near-infrared.

The denominator normalises. A slope facing away from the sun receives less light in every band, so both nir and red fall together, and their difference falls with them. Dividing by the sum cancels most of that, which is why an index is more comparable across terrain and across dates than a raw band.

That cancellation is also the limit. Anything that scales all bands equally is removed — and anything that does not, is not.

2. Know what each index actually contrasts

• NDVI (nir − red)/(nir + red) — the red-absorption/NIR-scattering gap. Not "greenness"; a stressed but structurally intact canopy can keep a high NDVI.
• NDWI (McFeeters) (green − nir)/(green + nir) — water absorbs NIR almost completely. Best index measured here for separating water from cloud shadow: 3.79 standard deviations against NDVI's 3.12.
• MNDWI (green − swir)/(green + swir) — better where buildings would otherwise be classed as water.
• NDBI (swir − nir)/(swir + nir) — built-up surfaces. Note it reads 0.287 over water in this scene, higher than over any land class. NDBI alone is not a built-up map.
• EVI — a soil- and atmosphere-corrected vegetation index that does not saturate as early as NDVI.
• SAVI — NDVI with a soil brightness term L, for sparse canopies.

3. Notice which indices are the same formula

On Sentinel-2, MNDWI and NDSI are both (green − swir16)/(green + swir16). They are the identical calculation with two different names, and in the table above they produce identical numbers to three decimal places.

They differ only in what you do with the result: high values are water in a coastal scene and snow in an alpine one. The arithmetic cannot tell you which — you need context, or a second index.

4. Check what the index does to your actual surfaces

The cloud shadow row is the whole argument for masking. Shadowed vegetation keeps a median NDVI of 0.622, because shadow suppresses red and near-infrared roughly proportionally and the ratio survives. NDVI cannot see the shadow; the scene classification band can.

NDVI over cloud shadow : 0.622
NDVI over vegetation   : 0.730

A composite that averages unmasked shadow into a vegetation time series drags it down by roughly 0.1 without ever producing an obviously wrong value.

5. Expect indices that "measure the same thing" to disagree

Over the 466,954 vegetation pixels:

       p5     p50     p95
NDVI  0.530  0.730  0.855
EVI   0.318  0.536  0.723
SAVI  0.300  0.493  0.619
correlation NDVI vs EVI: 0.612

A correlation of 0.612 means these indices rank pixels differently. If you threshold NDVI at 0.6 and someone else thresholds EVI at 0.4, you will disagree about a large fraction of the map and both be defensible.

Pick an index, state it, and keep it for the life of the analysis.

An index function that will not lie

The data used on this page

Sentinel-2 L2A scene S2A_30UVD_20250922_0_L2A from Earth Search, a partly cloudy late-September pass over Snowdonia, 4.16–4.00° W, 53.03–53.12° N. The setup reads the blue, green, red and near-infrared bands on their 10 m grid, brings the 20 m shortwave infrared (B11) onto it bilinearly and the scene classification band (scl) by nearest neighbour, and divides the digital numbers by 10,000 — Earth Search's copy needs no further offset.

import numpy as np
import rasterio
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import transform_bounds

AOI = (-4.16, 53.03, -4.00, 53.12)
item = next(Client.open("https://earth-search.aws.element84.com/v1").search(
    ids=["S2A_30UVD_20250922_0_L2A"], collections=["sentinel-2-l2a"]).items())
with rasterio.open(item.assets["red"].href) as ds:
    bounds = transform_bounds("EPSG:4326", ds.crs, *AOI)
    w = rasterio.windows.from_bounds(*bounds, transform=ds.transform).round_offsets().round_lengths()
    shape = (int(w.height), int(w.width))


def read(asset, resampling=Resampling.bilinear):
    with rasterio.open(item.assets[asset].href) as ds:
        win = rasterio.windows.from_bounds(*bounds, transform=ds.transform).round_offsets().round_lengths()
        return ds.read(1, window=win, out_shape=shape, resampling=resampling)


blue, green, red, nir, swir = (read(a).astype("float32") / 1e4 for a in ("blue", "green", "red", "nir", "swir16"))
scl = read("scl", Resampling.nearest)
print(f"{shape[0]} x {shape[1]} cells of 10 m; scene classes present: {sorted(np.unique(scl).tolist())}")

1017 x 1087 cells of 10 m; scene classes present: [2, 3, 4, 5, 6, 7, 8, 9, 10]

Example 1 — an index function that will not lie to you

import numpy as np


def spectral_index(bands, formula, valid_mask=None, clip=True):
    """Compute an index with masking applied first and the range checked after."""
    arrays = {k: np.asarray(v, dtype="float32") for k, v in bands.items()}

    if valid_mask is not None:
        for key in arrays:
            arrays[key] = np.where(valid_mask, arrays[key], np.nan)

    with np.errstate(divide="ignore", invalid="ignore"):
        out = formula(**arrays)

    finite = out[np.isfinite(out)]
    if finite.size:
        lo, hi = float(finite.min()), float(finite.max())
        print(f"  range {lo:.3f} to {hi:.3f}, "
              f"{np.isnan(out).sum():,} masked of {out.size:,}")
        if clip and (lo < -1.001 or hi > 1.001):
            raise ValueError(
                f"normalised index outside [-1, 1] ({lo:.3f} to {hi:.3f}) — "
                "check the band scaling before trusting this"
            )
    return out


ndvi_masked = spectral_index(
    {"nir": nir, "red": red},
    lambda nir, red: (nir - red) / (nir + red),
    valid_mask=~np.isin(scl, [0, 3, 8, 9, 10]),
)

  range -1.000 to 0.942, 565,605 masked of 1,105,479

Half the scene is gone — 51% of these pixels are cloud, cirrus or cloud shadow, which is what a single September date over Snowdonia looks like once you are honest about it.

Masking before the arithmetic rather than after matters: np.nan propagates through the division, but a cloud pixel that survives into the array will be averaged into every zonal statistic downstream.

Example 2 — comparing indices on the same pixels

import numpy as np
import pandas as pd


def compare_indices(indices, classes, class_names, min_pixels=500):
    """One row per index, one column per surface class: the table above."""
    rows = []
    for name, values in indices.items():
        row = {"index": name}
        for code, label in class_names.items():
            mask = (classes == code) & np.isfinite(values)
            row[label] = round(float(np.median(values[mask])), 3) \
                if mask.sum() >= min_pixels else None
        rows.append(row)

    frame = pd.DataFrame(rows).set_index("index")
    print(frame.to_string())

    # the separation that actually matters: can this index tell A from B?
    return frame


def separation(values, classes, code_a, code_b):
    """Distance between two class medians in pooled standard deviations."""
    a = values[(classes == code_a) & np.isfinite(values)]
    b = values[(classes == code_b) & np.isfinite(values)]
    pooled = np.sqrt(0.5 * (a.std() ** 2 + b.std() ** 2))
    return float(abs(np.median(a) - np.median(b)) / pooled)


evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)
savi = 1.5 * (nir - red) / (nir + red + 0.5)
ndsi = normalised_difference(green, swir)
indices = {"NDVI": ndvi, "NDWI": ndwi, "MNDWI": mndwi, "NDSI": ndsi, "NDBI": ndbi, "EVI": evi}
compare_indices(indices, scl, {4: "vegetation", 5: "bare soil", 6: "water", 9: "cloud", 3: "cloud shadow"})
print("water vs cloud shadow, separation in standard deviations")
for name in ("NDVI", "NDWI", "MNDWI"):
    print(f"  {name:6} {separation(indices[name], scl, 6, 3):.2f}")

veg = scl == 4
print(f"{veg.sum():,} vegetation pixels     p5     p50     p95")
for name, values in (("NDVI", ndvi), ("EVI", evi), ("SAVI", savi)):
    print(f"  {name:5} " + "  ".join(f"{v:.3f}" for v in np.percentile(values[veg], [5, 50, 95])))
print(f"  correlation NDVI vs EVI: {np.corrcoef(ndvi[veg], evi[veg])[0, 1]:.3f}")

       vegetation  bare soil  water  cloud  cloud shadow
index                                                   
NDVI        0.730      0.353 -0.096  0.083         0.622
NDWI       -0.689     -0.358  0.350 -0.081        -0.570
MNDWI      -0.539     -0.352  0.093 -0.001        -0.410
NDSI       -0.539     -0.352  0.093 -0.001        -0.410
NDBI       -0.225     -0.028  0.287 -0.082        -0.190
EVI         0.536      0.258 -0.003  0.341         0.225
water vs cloud shadow, separation in standard deviations
  NDVI   3.12
  NDWI   3.79
  MNDWI  2.60
466,954 vegetation pixels     p5     p50     p95
  NDVI  0.530  0.730  0.855
  EVI   0.318  0.536  0.723
  SAVI  0.300  0.493  0.619
  correlation NDVI vs EVI: 0.612

Class medians alone are misleading — two classes can have distant medians and still overlap heavily. Separation in pooled standard deviations is what tells you whether a threshold will work.

Example 3 — choosing a threshold from the data, not from a paper

import numpy as np


def otsu_threshold(values, bins=256):
    """The split that minimises within-class variance. No magic constant."""
    values = values[np.isfinite(values)]
    counts, edges = np.histogram(values, bins=bins)
    centres = (edges[:-1] + edges[1:]) / 2

    weight_low = np.cumsum(counts)
    weight_high = counts.sum() - weight_low
    valid = (weight_low > 0) & (weight_high > 0)

    mean_low = np.cumsum(counts * centres) / np.maximum(weight_low, 1)
    total = (counts * centres).sum()
    mean_high = (total - np.cumsum(counts * centres)) / np.maximum(weight_high, 1)

    between = weight_low * weight_high * (mean_low - mean_high) ** 2
    between[~valid] = -np.inf
    return float(centres[int(np.argmax(between))])


clear = ~np.isin(scl, [0, 1, 3, 8, 9, 10])
for name, values in (("NDWI", ndwi), ("NDVI", ndvi)):
    t = otsu_threshold(values[clear])
    print(f"Otsu threshold for {name} on clear pixels: {t:+.3f}")
water_by_otsu = (ndwi > otsu_threshold(ndwi[clear])) & clear
print(f"NDWI above it: {water_by_otsu.sum():,} cells; the classification band's water: {(scl == 6).sum():,}")

Otsu threshold for NDWI on clear pixels: -0.209
Otsu threshold for NDVI on clear pixels: +0.453
NDWI above it: 23,270 cells; the classification band's water: 17,834

Published thresholds — "NDWI > 0" for water, "NDVI > 0.3" for vegetation — are defaults from one sensor over one landscape. Deriving the threshold from your own histogram, and recording what it came out as, is both more accurate and more honest. On this scene Otsu put the NDWI split at −0.209 and the NDVI split at +0.453, well away from both defaults; NDWI above −0.209 picked out 23,270 clear cells against the classification band's 17,834 water cells, so even a derived threshold wants checking against a second source.

The normalised difference form and its limits

Why the normalised difference form is everywhere

(a − b)/(a + b) has three properties that make it hard to beat:

1. Bounded. Whenever a and b are non-negative, the result lies in [−1, 1]. That makes maps comparable and gives an immediate sanity check — a value outside that range means the inputs were not non-negative, which almost always means the scaling is wrong.
2. Scale invariant. Multiply both bands by any positive constant and the answer does not change. Illumination differences from slope, aspect and sun angle largely cancel.
3. Cheap. Two bands, three operations, no calibration, no training data.

The cost of that last property is that it has no idea what it is looking at.

Why NDVI saturates and EVI does not

NDVI's denominator grows with the near-infrared signal. Once a canopy is dense, extra leaf area adds near-infrared without removing much more red — red absorption is already near complete — so NDVI creeps towards its ceiling and stops responding.

That is visible in the percentiles above: NDVI's 5th-to-95th spread over vegetation is 0.53–0.855, compressed against the top of its range, while EVI spans 0.318–0.723 with room above.

EVI adds a blue term to correct for aerosols and a soil-adjustment constant in the denominator, which keeps it responsive at high biomass. It also needs a well-calibrated blue band, which is the band atmospheric correction gets wrong most often.

Why an index is not a classifier

An index is one number per pixel; a class is a decision. Turning the first into the second needs a threshold, and the table above shows why a fixed one is fragile: NDBI reads 0.287 over water, higher than over bare soil. Threshold NDBI at 0.2 to find buildings and you will map every lake in the scene.

The fix is not a better index. It is either a second index that separates the confused pair — NDWI would remove the water immediately — or a classifier that sees all the bands at once. See How to classify land cover from satellite imagery in Python.

Why masking has to come first

Cloud shadow at NDVI 0.622 is the measurement that makes this concrete. Indices are ratios, and a ratio survives the thing that made the pixel dark. The index cannot detect the shadow, so a downstream threshold cannot either.

The scene classification band, or a dedicated cloud mask, is the only part of the pipeline that knows the pixel is unusable. If it runs after the index, its information is already lost into an average.

MNDWI against NDSI, division by zero and casting

• MNDWI and NDSI are the same formula on Sentinel-2. Green and SWIR-1. Interpretation, not arithmetic, tells you whether high values are water or snow.
• Divide by zero happens. Where nir + red is 0 the index is undefined; return NaN rather than 0, which is a legitimate index value.
• Cast to float before the arithmetic. Integer division truncates; unsigned subtraction wraps.
• An index outside [−1, 1] is a scaling bug, not an interesting pixel. See NDVI values are above 1 or below −1.
• Resample coarse bands before mixing them. MNDWI needs a 10 m green and a 20 m SWIR on one grid.
• NDBI reads high over water (0.287 here). Always pair it with a water index.
• EVI needs the blue band, which atmospheric correction handles worst. Over hazy scenes EVI is noisier than NDVI even though it is theoretically better.
• Thresholds do not transfer between sensors, seasons or regions. Derive them, record them, and re-derive when the input changes.

Internal links

• Spectral bands explained: what satellite imagery actually measures — the reflectance shapes these indices contrast
• How to calculate NDVI from Sentinel-2 in Python — the implementation, with masking in the right order
• NDVI values are above 1 or below −1: how to fix it — what an out-of-range index means
• Cloud masking explained — why cloud shadow reaches NDVI 0.62
• How to extract a vegetation index time series for a polygon — indices over time rather than over space
• How to classify land cover from satellite imagery in Python — when one threshold stops being enough
• Radiometric levels explained — why an index changes if the level does
• How to detect change between two satellite images in Python — differencing indices across dates
• NDVI stops responding in a dense canopy — saturation, measured

FAQ

What is a spectral index?

A ratio between bands that emphasises one contrast — usually the normalised difference (a − b)/(a + b), which is bounded between −1 and 1 and largely cancels illumination differences.

What is a good NDVI value for healthy vegetation?

Over this scene, vegetation had a median NDVI of 0.73 with a 5th-to-95th percentile range of 0.53 to 0.855. But cloud shadow reached 0.622, so a threshold alone cannot separate them.

Why do NDVI and EVI disagree?

They contrast different band combinations and respond differently at high biomass. Measured on the same 466,954 vegetation pixels they correlate at 0.612 — related, but not interchangeable.

Which water index should I use?

NDWI separated water from cloud shadow best here (3.79 pooled standard deviations against MNDWI's 2.60). MNDWI is preferable where built-up surfaces would otherwise be classed as water.

Is NDBI a reliable built-up index?

Not on its own. It read 0.287 over water in this scene, higher than over any land class. Pair it with a water index before thresholding.

Should I mask clouds before or after computing an index?

Before. An index is a ratio, so it survives shadow — cloud shadow here has a median NDVI of 0.622. Once the index is computed the information that the pixel was unusable is gone.

Can I compare NDVI between two satellites?

Only roughly. Band widths and centre wavelengths differ, so the same ground gives slightly different NDVI. For a time series across sensors, either apply a published cross-calibration or treat each sensor as a separate series.