Thresholding dNBR to Standard Burn Severity Classes in Python

A dNBR raster is a continuous number that nobody manages by. Severity becomes actionable only once it is a class — unburned, low, moderate, high — and the break values decide how many hectares fall into each. Published thresholds are widely used and were calibrated in a specific forest type on a specific sensor, so applying them unexamined to another system quietly moves hectares between classes. This guide covers doing it defensibly, inside Burn Severity Mapping with Python, part of Fire Risk & Fuel Assessment.

When to use each index

Use dNBR in a homogeneous forest where every stand had roughly the same pre-fire biomass. Use RdNBR or RBR when the fire crossed a mosaic — closed forest, open woodland, shrub — because absolute change conflates severity with how much vegetation existed to begin with. That distinction and its arithmetic are set out in Burn Severity Mapping with Python, which owns the index definitions this page thresholds.

Offset correction comes before any threshold

Two scenes a year apart differ in phenology, illumination and atmosphere even where nothing burned. That difference — the offset — shifts the whole dNBR distribution and therefore every class boundary. Correcting it is a single subtraction and it is not optional.

$${\mathrm{dNBR}}_{\text{corrected}}=\left({\mathrm{NBR}}_{\text{pre}}-{\mathrm{NBR}}_{\text{post}}\right)\times 1000-\mathrm{offset}$$

where the offset is the median scaled dNBR of unburned pixels just outside the perimeter, in the same vegetation types that burned.

Minimal reproducible example

import numpy as np
import geopandas as gpd
import rasterio
from rasterio.mask import mask

# Widely used scaled-dNBR breaks. Treat them as a starting point, not a standard.
BREAKS = {
    "unburned": (-np.inf, 99),
    "low": (100, 269),
    "moderate": (270, 659),
    "high": (660, np.inf),
}


def compute_offset(dnbr_path: str, perimeter_path: str,
                   buffer_m: float = 2000.0, inner_m: float = 200.0) -> float:
    """Median scaled dNBR in an unburned ring around the fire perimeter."""
    perim = gpd.read_file(perimeter_path)
    with rasterio.open(dnbr_path) as src:
        geom = perim.to_crs(src.crs).geometry.union_all()
        ring = geom.buffer(buffer_m).difference(geom.buffer(inner_m))
        arr, _ = mask(src, [ring], crop=True, filled=True, nodata=np.nan)
    vals = arr[0][np.isfinite(arr[0])]
    return float(np.median(vals))


def classify_severity(dnbr: np.ndarray, offset: float = 0.0) -> np.ndarray:
    """Ordinal severity classes from offset-corrected scaled dNBR."""
    d = dnbr - offset
    cls = np.zeros(d.shape, dtype="uint8")      # 0 = nodata
    valid = np.isfinite(d)
    cls[valid] = 1                              # 1 = unburned / very low
    cls[valid & (d >= 100)] = 2                 # 2 = low
    cls[valid & (d >= 270)] = 3                 # 3 = moderate
    cls[valid & (d >= 660)] = 4                 # 4 = high
    return cls

The inner buffer matters. Taking the offset from pixels immediately outside the perimeter contaminates it with edge effects — scorch, partly burned ground, and perimeter mapping error — all of which bias the offset upward and shift every class boundary in the same direction.

Calibrate the breaks where you can

Published breaks were derived by regressing dNBR against field Composite Burn Index measurements in particular forests. Where you have CBI plots, fit local breaks; the fitted values often differ from the published ones by 50–100 units, which is several percent of the scar.

import numpy as np


def calibrate_breaks(dnbr_at_plots: np.ndarray, cbi: np.ndarray,
                     cbi_breaks=(0.1, 1.25, 2.25)) -> dict:
    """Map CBI class boundaries onto dNBR by a monotone fit through the plot data."""
    ok = np.isfinite(dnbr_at_plots) & np.isfinite(cbi)
    d, c = dnbr_at_plots[ok], cbi[ok]
    order = np.argsort(c)
    d_sorted, c_sorted = d[order], c[order]
    fitted = {}
    for b in cbi_breaks:
        # Interpolate the dNBR value at each CBI boundary.
        fitted[b] = float(np.interp(b, c_sorted, d_sorted))
    return {
        "n_plots": int(ok.sum()),
        "dnbr_at_cbi": {str(k): round(v) for k, v in fitted.items()},
        "spearman": float(np.corrcoef(np.argsort(np.argsort(d)),
                                      np.argsort(np.argsort(c)))[0, 1]),
    }

A Spearman correlation below about 0.7 between dNBR and CBI means the index is not tracking severity well at this site — usually because pre-fire vegetation was heterogeneous, which is the signal to switch to RdNBR rather than to adjust the breaks.

Parameter reference

Expected output and verification

The checks that catch real errors are about distribution shape and about the relationship to independent evidence:

import numpy as np


def severity_report(cls: np.ndarray, dnbr: np.ndarray, offset: float) -> dict:
    """Class areas plus the diagnostics that reveal a bad offset or bad breaks."""
    valid = cls > 0
    n = int(valid.sum())
    frac = {name: round(float((cls == i).sum() / n), 3)
            for i, name in enumerate(["nodata", "unburned", "low", "moderate", "high"])
            if i > 0}
    d = dnbr[valid] - offset
    return {
        "n_px": n,
        "class_fraction": frac,
        "median_dnbr": round(float(np.median(d))),
        "p95_dnbr": round(float(np.percentile(d, 95))),
        "offset_applied": round(offset, 1),
    }


rep = severity_report(cls, dnbr, offset)
assert rep["class_fraction"]["unburned"] < 0.5, \
    f"most of the perimeter classed unburned — check the offset sign: {rep}"
assert rep["class_fraction"]["high"] < 0.8, \
    f"almost everything high severity — check the pre-fire date: {rep}"

The unburned-fraction assertion catches the sign error, which is the single most common mistake: subtracting post from pre gives positive dNBR for a burn, and reversing it produces a map where the fire looks like vegetation growth.

Feed the classifier something it must reject as a standing test:

rng = np.random.default_rng(20260809)
# An unburned scene differenced against itself: every pixel must land in class 1.
noise = rng.normal(0.0, 15.0, size=20000)
cls_null = classify_severity(noise, offset=0.0)
assert (cls_null == 1).mean() > 0.99, \
    "noise around zero produced burn classes — the low break is too tight"

Common pitfalls

• Skipping the offset. Every published threshold assumes an offset-corrected input. Without it the classes are shifted by whatever the phenological difference happened to be.
• Taking the offset from immediately outside the perimeter. Edge pixels are partly burned and bias the offset upward.
• Using an initial-assessment threshold on an extended assessment. Immediate post-fire imagery and next-season imagery need different breaks, because delayed mortality has not yet expressed in the first.
• Applying forest thresholds to shrubland. Shrub systems lose far less NBR in an equally lethal fire, so forest breaks under-classify them severely.
• Reporting class areas without stating the breaks. Two agencies mapping the same fire will disagree by tens of percent, and the difference is always in the thresholds.

Frequently Asked Questions

Related

• Burn Severity Mapping with Python — the wider severity workflow
• Calculating dNBR from Sentinel-2 in Python — produces the continuous index this classifies
• Classifying Recovery Classes from Multi-Year dNBR — the same threshold discipline applied to recovery
• Masking Clouds with the Sentinel-2 Scene Classification Layer — an unmasked cloud reads as high severity

Up: Burn Severity Mapping with Python · Fire Risk & Fuel Assessment