Mapping Burned Area with NBR Differencing

To map fire extent and severity, compute the Normalised Burn Ratio before and after the fire and difference them:

def nbr(nir, swir22):
    denom = nir + swir22
    return (nir - swir22) / denom

dnbr = nbr(nir_pre, swir_pre) - nbr(nir_post, swir_post)   # positive = burn severity

dNBR is the best-established application of the method described in Change Detection and Differencing Workflows, and the one case where published class breaks make a derived threshold optional.

Why This Arises in Remote Sensing Workflows

Fire changes vegetation in a way the spectrum reports unusually clearly. Combustion removes chlorophyll and water, so near-infrared reflectance falls sharply, while char and exposed soil raise short-wave infrared reflectance. NBR combines exactly those two bands, so it moves further and more consistently after fire than any general vegetation index — which is why burned-area mapping is one of the few change applications where a single date pair is genuinely adequate.

The operational demand is immediate: agencies need extent within days for response and severity within weeks for rehabilitation planning. That timeline rules out composite pairs and time-series fitting, and it means the pipeline has to be robust with whatever two clear scenes exist around the event.

The complication is look-alikes. Anything that removes vegetation and exposes soil looks like fire in NBR: clear-felling, ploughing, a drained reservoir, a landslide. Separating them is the part of the workflow that needs judgement rather than arithmetic.

Why NBR responds so strongly to fire Healthy vegetation is bright in the near-infrared and dark in the long short-wave infrared, giving a high NBR. After fire the near-infrared collapses and the short-wave infrared rises because char and soil are exposed, so NBR falls sharply. The difference between the two states is larger and more consistent than for a vegetation index alone. Reflectance before and after fire 0.4 0 B08 / B8A — near-infrared B12 — SWIR 2.2 µm 0.36 0.16 0.10 0.28 pre-fire post-fire NBR: 0.38 → −0.47 · dNBR = 0.85 Both bands move, and they move in opposite directions — which is what makes the ratio so sensitive.

Environment & Setup

Package Version Why
rioxarray ≥0.15 Masked reads, alignment, COG output
xarray ≥2023.6 Labelled arithmetic
rasterio ≥1.3.0 Resampling and I/O
scikit-image ≥0.22 Minimum mapping unit cleanup
numpy ≥1.23 Classification and statistics 
pip install "rioxarray>=0.15" "xarray>=2023.6" "rasterio>=1.3.0" "scikit-image>=0.22"

Complete Working Example

This function computes dNBR, classifies it with the widely used severity breaks, removes water and small objects, and returns both the continuous and classified layers.

dNBR severity classes and their usual interpretation The widely used dNBR breaks map to unburned, low, moderate-low, moderate-high and high severity. The class boundaries were calibrated on particular ecosystems, so they are a starting point rather than a universal truth, and the low boundary is where area estimates are decided. Standard dNBR class breaks dNBR range class typical appearance < 0.10 unburned indistinguishable from background 0.10 – 0.27 low severity surface scorch, canopy intact 0.27 – 0.44 moderate-low partial canopy loss 0.44 – 0.66 moderate-high most canopy consumed > 0.66 high severity complete consumption, bare soil Recalibrate the boundaries against local field severity data before publishing area by class. 
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from rasterio.enums import Resampling
from skimage.morphology import remove_small_objects

# Widely used dNBR severity breaks (USGS convention)
SEVERITY_BREAKS = [
    (-np.inf, 0.10, 0, "unburned"),
    (0.10, 0.27, 1, "low"),
    (0.27, 0.44, 2, "moderate-low"),
    (0.44, 0.66, 3, "moderate-high"),
    (0.66, np.inf, 4, "high"),
]


def nbr(nir: xr.DataArray, swir: xr.DataArray) -> xr.DataArray:
    denom = nir + swir
    return xr.where(denom != 0, (nir - swir) / denom, np.nan).astype("float32")


def map_burn(
    pre: dict, post: dict, *, mmu_pixels: int = 6, water_ndwi_threshold: float = 0.0
) -> tuple[xr.DataArray, xr.DataArray, dict]:
    """Return (dnbr, severity, stats) for one fire event."""
    # B8A and B12 share the 20 m grid — no resampling needed between them
    nir_pre = rioxarray.open_rasterio(pre["b8a"], masked=True).squeeze(drop=True).astype("float32")
    swir_pre = rioxarray.open_rasterio(pre["b12"], masked=True).squeeze(drop=True).astype("float32")
    nir_post = rioxarray.open_rasterio(post["b8a"], masked=True).squeeze(drop=True)
    swir_post = rioxarray.open_rasterio(post["b12"], masked=True).squeeze(drop=True)

    nir_post = nir_post.rio.reproject_match(nir_pre, resampling=Resampling.bilinear).astype("float32")
    swir_post = swir_post.rio.reproject_match(nir_pre, resampling=Resampling.bilinear).astype("float32")

    nbr_pre, nbr_post = nbr(nir_pre, swir_pre), nbr(nir_post, swir_post)
    dnbr = (nbr_pre - nbr_post).astype("float32")

    # Water look-alike removal: exclude pixels that are water on either date
    green_pre = rioxarray.open_rasterio(pre["b03"], masked=True).squeeze(drop=True).astype("float32")
    green_pre = green_pre.rio.reproject_match(nir_pre, resampling=Resampling.bilinear)
    ndwi_pre = xr.where((green_pre + nir_pre) != 0,
                        (green_pre - nir_pre) / (green_pre + nir_pre), np.nan)
    dnbr = dnbr.where(ndwi_pre < water_ndwi_threshold)

    severity = xr.zeros_like(dnbr, dtype="uint8")
    for lo, hi, code, _label in SEVERITY_BREAKS:
        severity = xr.where((dnbr >= lo) & (dnbr < hi), code, severity)
    severity = severity.where(dnbr.notnull())

    burned = (severity.fillna(0).values >= 1)
    burned = remove_small_objects(burned, min_size=mmu_pixels)
    severity = severity.where(xr.DataArray(burned | (severity.fillna(0).values == 0),
                                           coords=severity.coords, dims=severity.dims), 0)

    pixel_area_ha = abs(dnbr.rio.resolution()[0] * dnbr.rio.resolution()[1]) / 10_000
    stats = {
        "burned_pixels": int(burned.sum()),
        "burned_area_ha": round(float(burned.sum()) * pixel_area_ha, 2),
        "high_severity_pixels": int((severity == 4).sum()),
        "date_pre": pre["date"],
        "date_post": post["date"],
        "mmu_pixels": mmu_pixels,
    }
    return dnbr, severity, stats


if __name__ == "__main__":
    dnbr, severity, stats = map_burn(
        {"b8a": "B8A_pre.tif", "b12": "B12_pre.tif", "b03": "B03_pre.tif", "date": "2023-08-04"},
        {"b8a": "B8A_post.tif", "b12": "B12_post.tif", "date": "2023-08-21"},
    )
    print(stats)
    severity.rio.to_raster("burn_severity.tif", driver="COG", compress="DEFLATE", dtype="uint8")

Using B8A rather than B08 is a small decision with a real payoff: it is already on the 20 m grid that B12 uses, so the two bands need no relative resampling and the index is computed from unmodified observations.

Variant Patterns

1. Relativised dNBR for mixed vegetation

Absolute dNBR depends on how much vegetation there was to burn, so the same fire severity produces a smaller dNBR in sparse shrubland than in closed forest. RdNBR normalises by pre-fire condition.

import numpy as np

denom = np.sqrt(np.abs(nbr_pre.where(np.abs(nbr_pre) > 0.001)))
rdnbr = (dnbr / denom).astype("float32")

Use RdNBR when one map spans several vegetation types and the severity classes must mean the same thing across them; use dNBR when the area is homogeneous and interpretability matters more.

2. Offsetting for phenology between the dates

Correcting dNBR for background change Unburned ground between two dates several weeks apart is not exactly unchanged: seasonal drying alone can shift dNBR by 0.03 to 0.05. Subtracting the median dNBR of unburned surroundings recentres the distribution so the published severity breaks apply as intended. dNBR inside and outside the perimeter unburned surroundings inside the perimeter offset ≈ +0.04 −0.1 0.3 1.0 dNBR Subtract the unburned median before classifying, or every class boundary sits 0.04 too low. The correction matters most for the low-severity boundary, which is where area estimates are decided. 
outside = dnbr.where(~perimeter_buffer)          # ring around the fire, known unburned
offset = float(outside.median(skipna=True))
dnbr_corrected = dnbr - offset

3. Separating harvest from fire

Clear-felling produces a dNBR signature close to moderate-severity fire. Three discriminators usually resolve it: shape, since harvest blocks have straight edges and fires do not; the short-wave infrared magnitude, which is higher over char than over slash; and independent context such as an active-fire detection during the window. Where the distinction matters, encode the shape test — a compactness or edge-straightness measure on each connected component — rather than leaving it to visual inspection, and record which components were reclassified.

Validating a Burn Map

Fire is one of the better-served applications for validation because independent evidence usually exists.

Compare the mapped perimeter against an independent product if one is published for the event, and report the agreement as an area ratio rather than as a subjective assessment. Differences at the low-severity boundary are expected; differences in the high-severity core are a warning.

Check that severity decreases outward. A real fire almost always shows a severity gradient from the core to the perimeter, so a map with high-severity pixels scattered outside a low-severity ring is reporting noise, not structure.

Confirm the unburned control. Sample the surroundings and confirm their dNBR is centred near zero after the offset correction — if it is not, the correction did not work and the class boundaries are misplaced, which is discussed in Thresholding Change Maps and Removing Noise.

Common Errors

Lakes and rivers appear as high severity

Water was not masked, or the water mask came from only one date. Compute NDWI on both dates and exclude pixels that are water on either — the index is covered in Computing EVI and NDWI from Sentinel-2 Bands.

The whole scene shows low severity

Phenological drying between the dates has shifted the background. Apply the unburned-surroundings offset before classifying.

Severity classes look wrong for this vegetation type

The published breaks were calibrated on different ecosystems. Either recalibrate against local field severity data, or switch to RdNBR and say so in the product metadata.

Frequently Asked Questions

Q: Which Sentinel-2 bands does NBR use? B08 near-infrared and B12 long short-wave infrared for the standard formulation. B8A is often preferred over B08 because it is narrower and shares the 20 m grid with B12, which removes a resampling step.

Q: Is dNBR or RdNBR better? dNBR measures absolute change and is easier to interpret; RdNBR normalises by pre-fire condition and compares better across vegetation types with different pre-fire biomass. Sparse pre-fire vegetation is where the two disagree most.

Q: Why does water show up as high severity? Water has strongly negative NBR, so a lake that was dry before and wet after produces a large positive dNBR. Mask water on both dates with an index such as NDWI before classifying.