Computing NDWI and NBR for Water and Burn Mapping

Compute NDWI, MNDWI and NBR from Sentinel-2 or Landsat with rasterio and numpy: band selection, scale factors, safe division, and dNBR burn-severity differencing.

Both indices are normalised differences: NDWI = (green - nir) / (green + nir), NBR = (nir - swir2) / (nir + swir2). The arithmetic is trivial and the failure modes are not — offsets that do not cancel, division by zero over nodata, band numbers that differ between Sentinel-2 and Landsat, and burn-severity thresholds borrowed from a paper about a different biome.

This page covers water and burn indices specifically. For the shared machinery — reading bands, masking clouds, writing index rasters — see Raster Band Math and Index Calculation.

Why these particular band pairs separate water and burn Reflectance curves are drawn for water, healthy vegetation and burnt ground across the green, near infrared and shortwave infrared bands. Water is bright in green and dark in near infrared, which NDWI contrasts. Healthy vegetation is bright in near infrared and dark in shortwave infrared while burnt ground is the reverse, which NBR contrasts. Each index picks the two bands where its target class flips sign refl. greenNIRSWIR2 water vegetation burnt NDWI NBR

Why This Fails in Geospatial ML Pipelines

The first trap is the additive offset. A normalised difference is scale-invariant — multiplying both bands by 10,000 leaves the ratio unchanged — so teams learn that scale factors do not matter and generalise that to offsets. They are not equivalent. Sentinel-2 products from processing baseline 04.00 onwards carry a BOA_ADD_OFFSET of −1000 in digital-number space. Skip it and a reflectance of 0.05 is read as 0.15, which shifts NDWI by a large fraction of its useful range and makes a model trained on older scenes fail on newer ones with no code change anywhere.

The second trap is division. Over nodata, cloud shadow, or deep shade both bands can be zero, so the denominator is zero and NumPy emits inf or nan with a warning most pipelines have silenced. Those values then propagate into scaling and into the model. np.errstate plus an explicit where= mask is two lines and removes the whole class of problem.

Third, band numbering is not portable. Sentinel-2 SWIR2 is band 12; Landsat 8 SWIR2 is band 7; on a subset written by an earlier step the bands may be renumbered 1…n. Hard-coding indices produces an index that is arithmetically valid and physically meaningless — the classic case being NBR computed from SWIR1 instead of SWIR2, which still correlates with burn severity just weakly enough to look like a modelling problem rather than a bug.

Core Principles

  • Remove offsets before dividing; scale factors may cancel, offsets never do.
  • Resolve bands by name, not by number. Keep a per-collection mapping in configuration.
  • Divide safely. np.divide(..., where=denominator != 0) with a pre-filled nan output.
  • Mask clouds first. An index over cloud is a confident number about the atmosphere.
  • Prefer continuous indices over thresholded classes as model inputs — thresholds discard information the model can use.
  • Difference the same geometry. dNBR requires pre- and post-fire scenes on an identical grid.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
import rasterio

logger = logging.getLogger(__name__)

# Band name -> 1-based index, per collection. Never hard-code these inline.
BAND_MAP = {
    "sentinel-2-l2a": {"green": 3, "red": 4, "nir": 8, "swir1": 11, "swir2": 12},
    "landsat-c2-l2":  {"green": 3, "red": 4, "nir": 5, "swir1": 6, "swir2": 7},
}
# Additive offset in the product's own units (Sentinel-2 baseline >= 04.00).
BAND_OFFSET = {"sentinel-2-l2a": -1000.0, "landsat-c2-l2": 0.0}
BAND_SCALE = {"sentinel-2-l2a": 1e-4, "landsat-c2-l2": 2.75e-5}


def read_reflectance(path: str, collection: str, names: list[str]) -> dict[str, np.ndarray]:
    """Read named bands as physical reflectance, with nodata as NaN.

    Applies the collection's additive offset BEFORE the scale factor, which is the
    order the product specifications define.
    """
    if collection not in BAND_MAP:
        raise ValueError(f"unknown collection {collection!r}; add it to BAND_MAP")
    idx = BAND_MAP[collection]
    missing = [n for n in names if n not in idx]
    if missing:
        raise ValueError(f"{collection} has no band(s) {missing}")

    out = {}
    with rasterio.open(path) as src:
        for name in names:
            band = src.read(idx[name], masked=True).filled(np.nan).astype("float32")
            band = (band + BAND_OFFSET[collection]) * BAND_SCALE[collection]
            out[name] = band
    logger.info("read %s from %s as reflectance", names, collection)
    return out


def normalised_difference(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    """(a - b) / (a + b), NaN where the denominator vanishes or inputs are invalid."""
    denom = a + b
    out = np.full(a.shape, np.nan, dtype="float32")
    valid = np.isfinite(a) & np.isfinite(b) & (np.abs(denom) > 1e-6)
    with np.errstate(invalid="ignore", divide="ignore"):
        np.divide(a - b, denom, out=out, where=valid)
    return out


def ndwi(bands: dict) -> np.ndarray:
    """McFeeters NDWI — open water against vegetation."""
    return normalised_difference(bands["green"], bands["nir"])


def mndwi(bands: dict) -> np.ndarray:
    """Modified NDWI — suppresses built-up surfaces that NDWI reads as water."""
    return normalised_difference(bands["green"], bands["swir1"])


def nbr(bands: dict) -> np.ndarray:
    """Normalised Burn Ratio — healthy vegetation positive, burnt ground negative."""
    return normalised_difference(bands["nir"], bands["swir2"])


def dnbr(pre_path: str, post_path: str, collection: str) -> np.ndarray:
    """Burn severity as pre-fire NBR minus post-fire NBR, on a verified shared grid."""
    with rasterio.open(pre_path) as a, rasterio.open(post_path) as b:
        if a.crs != b.crs or (a.width, a.height) != (b.width, b.height) \
                or not a.transform.almost_equals(b.transform, precision=1e-6):
            raise ValueError(
                "pre- and post-fire scenes are not on the same grid; reproject or "
                "re-window one onto the other before differencing")

    pre = nbr(read_reflectance(pre_path, collection, ["nir", "swir2"]))
    post = nbr(read_reflectance(post_path, collection, ["nir", "swir2"]))
    return (pre - post).astype("float32")

Step-by-Step Walkthrough

Step 1 — Read bands as reflectance and check the range

bands = read_reflectance("S2_20260714.tif", "sentinel-2-l2a",
                         ["green", "nir", "swir1", "swir2"])
for name, arr in bands.items():
    print(f"{name:>6}: {np.nanmin(arr):.3f} to {np.nanmax(arr):.3f}")

Surface reflectance should sit roughly in [0, 1], with a small tail above 1 over bright cloud and snow. Values clustered around 1,000–8,000 mean the scale factor was not applied; values shifted by about 0.1 mean the offset was not.

Step 2 — Mask clouds before computing anything

from rasterio import open as rio_open

with rio_open("S2_20260714_SCL.tif") as scl_src:
    scl = scl_src.read(1)

# Sentinel-2 scene classification: 3 shadow, 8-10 cloud, 11 snow.
cloudy = np.isin(scl, [3, 8, 9, 10, 11])
for name in bands:
    bands[name] = np.where(cloudy, np.nan, bands[name])
print(f"masked {cloudy.mean():.1%} of pixels")

The masking machinery is covered in full in masking clouds before computing NDVI; the same mask applies to every index.

Step 3 — Compute the indices and inspect their distributions

water = ndwi(bands)
water_urban = mndwi(bands)
burn = nbr(bands)

for name, arr in [("NDWI", water), ("MNDWI", water_urban), ("NBR", burn)]:
    finite = arr[np.isfinite(arr)]
    print(f"{name:>6}: median {np.median(finite):+.3f}, "
          f"p5 {np.percentile(finite, 5):+.3f}, p95 {np.percentile(finite, 95):+.3f}")
    assert np.nanmin(arr) >= -1.0001 and np.nanmax(arr) <= 1.0001, f"{name} out of range"

Every normalised difference is bounded to [-1, 1] by construction. A value outside that range is proof of a bug — usually a negative reflectance from an unapplied offset.

Step 4 — Difference two dates for burn severity

severity = dnbr("S2_prefire.tif", "S2_postfire.tif", "sentinel-2-l2a")
print(f"dNBR p99 = {np.nanpercentile(severity, 99):.3f}")

Keep severity continuous when it feeds a model. Thresholding into severity classes throws away resolution the model would otherwise exploit, and imported thresholds rarely transfer between biomes.

Step 5 — Write a multi-band index stack

with rasterio.open("S2_20260714.tif") as src:
    profile = src.profile.copy()
profile.update(count=4, dtype="float32", nodata=np.nan, compress="deflate", tiled=True)

with rasterio.open("indices.tif", "w", **profile) as dst:
    for i, (arr, name) in enumerate(
            zip([water, water_urban, burn, severity], ["ndwi", "mndwi", "nbr", "dnbr"]), 1):
        dst.write(arr, i)
        dst.set_band_description(i, name)
The scale factor cancels; the offset does not Two conversion paths are compared for the same pair of digital numbers. In the first the scale factor is applied without removing the additive offset, giving inflated reflectance values and an index that is biased away from the true value. In the second the offset is removed first and the index matches the true value. Sentinel-2 baseline 04.00: BOA_ADD_OFFSET = −1000 offset ignored DN green 1500, DN nir 1200 × 1e−4 → 0.150, 0.120 NDWI = +0.111 offset removed first DN green 1500, DN nir 1200 (DN − 1000) × 1e−4 → 0.050, 0.020 NDWI = +0.429 The same pixel: water in one pipeline, not water in the other. A model trained on pre-baseline scenes silently mispredicts on newer ones.

Verification

import numpy as np


def test_indices_are_bounded():
    rng = np.random.default_rng(0)
    a, b = rng.uniform(0, 1, 10_000), rng.uniform(0, 1, 10_000)
    nd = normalised_difference(a, b)
    assert np.nanmin(nd) >= -1.0 and np.nanmax(nd) <= 1.0


def test_zero_denominator_yields_nan_not_inf():
    a = np.array([0.0, 0.4], dtype="float32")
    b = np.array([0.0, 0.2], dtype="float32")
    nd = normalised_difference(a, b)
    assert np.isnan(nd[0]) and np.isfinite(nd[1])


def test_scale_factor_cancels_but_offset_does_not():
    g, n = 0.15, 0.12
    assert np.isclose(normalised_difference(np.array([g * 10]), np.array([n * 10]))[0],
                      normalised_difference(np.array([g]), np.array([n]))[0])
    assert not np.isclose(normalised_difference(np.array([g + 0.1]), np.array([n + 0.1]))[0],
                          normalised_difference(np.array([g]), np.array([n]))[0])


def test_known_water_pixel_is_positive_ndwi():
    """Open water: bright green, near-zero NIR."""
    bands = {"green": np.array([[0.09]], "float32"), "nir": np.array([[0.02]], "float32")}
    assert ndwi(bands)[0, 0] > 0.3


def test_burnt_pixel_is_negative_nbr():
    bands = {"nir": np.array([[0.12]], "float32"), "swir2": np.array([[0.28]], "float32")}
    assert nbr(bands)[0, 0] < -0.3

The last two tests encode domain knowledge rather than arithmetic, and they are the ones that catch a swapped band mapping — the failure that arithmetic tests cannot see.

Why a continuous index beats a thresholded class An NDWI histogram is bimodal, with a large land mode below zero and a smaller water mode above zero. The region between the modes is shaded to show shallow water, wet soil and mixed pixels, where any fixed threshold makes errors that a model given the continuous value could avoid. NDWI distribution over a mixed scene −1.00.0+1.0 land open water wet soil, shallow water, mixed pixels Any single cut through the shaded band throws away information the model could use.

FAQ

What is the difference between NDWI and MNDWI?

NDWI contrasts green with near-infrared and separates open water from vegetation. MNDWI swaps near-infrared for shortwave infrared, which suppresses built-up surfaces that NDWI frequently reads as water. In urban catchments MNDWI is usually the better predictor; in rural ones they behave similarly.

Do I need to apply a scale factor before computing an index?

The multiplicative scale cancels in a normalised difference, so it is optional for the index alone. The additive offset does not cancel — Sentinel-2 baseline 04.00 and later carry one, and ignoring it biases every value.

How should dNBR thresholds be chosen?

Locally, from field-validated plots, or not at all. Published breakpoints are tied to a sensor, a biome and a time offset. Feeding continuous dNBR to the model usually outperforms any thresholding, and avoids importing someone else’s calibration.

Part of: Raster Band Math and Index Calculation Part of: Spatial Feature Engineering for Machine Learning