Burn Severity Mapping with Python

A fire has burned 40,000 hectares across a mixed-conifer landscape, and the burned-area emergency response team needs a defensible severity map within days: where did the fire consume the canopy and sterilize soils, where did it creep through the understory, and where did it barely touch. Field crews can walk a few dozen plots; they cannot walk the whole perimeter. Satellite burn severity mapping bridges that gap by differencing a spectral index before and after the fire, calibrating the result against the field plots, and classifying every pixel into an interpretable severity class. This page builds that workflow as a reproducible Python pipeline within the broader Fire Risk & Fuel Assessment section, turning a pair of scenes into a projection-correct, field-anchored severity raster suitable for rehabilitation planning and long-term monitoring.

The core signal is the Normalized Burn Ratio (NBR). Live, healthy vegetation is bright in the near-infrared and dark in the shortwave-infrared; fire inverts that relationship by removing chlorophyll and exposing char, ash, and soil, which are bright in the shortwave-infrared. Differencing NBR between a pre-fire and a post-fire scene isolates that change from the static background, and the magnitude of the change scales with severity. The discipline is in the details: pairing the right scenes, aligning them to the same grid, masking cloud and cloud shadow, and correcting for the phenological offset between two dates that were never meant to be compared pixel-for-pixel.

Prerequisites checklist

Confirm each item before differencing anything. Every one of these is a real precondition — skipping it is the usual root cause of a severity map that looks plausible but is wrong.

• Python 3.10+ with rasterio >= 1.3, numpy >= 1.24 (built against GDAL >= 3.6)
• A pre-fire and a post-fire scene from the same sensor, ideally the same relative orbit, with minimal cloud over the burn perimeter
• Both scenes converted to surface reflectance with the documented scale and offset (Sentinel-2 L2A: scale 10000, BOA_ADD_OFFSET -1000 for baseline 04.00+)
• NIR and SWIR bands resampled to a common ground resolution — for Sentinel-2 use B8A (narrow NIR) and B12 (SWIR2), both native 20 m, so no resampling mismatch
• A single project-wide UTM or equal-area target CRS chosen via Coordinate Reference Systems for Forestry so the two dates overlay exactly
• A per-scene cloud/shadow mask — the Sentinel-2 Scene Classification Layer (SCL) or a Landsat QA band — aligned pixel-for-pixel to the bands
• Field-observed severity plots (Composite Burn Index or similar) for calibrating thresholds, ideally laid out through Spatial Plot Sampling Design

Concept background: NBR, dNBR, and RdNBR

The Normalized Burn Ratio is a normalized-difference index built from the same band-math family as the vegetation indices in Vegetation Index Calculation in Python, but it swaps the red band for the shortwave-infrared to target the char-and-moisture signal rather than chlorophyll:

$$\text{NBR}=\frac{{\rho }_{\text{NIR}}-{\rho }_{\text{SWIR}}}{{\rho }_{\text{NIR}}+{\rho }_{\text{SWIR}}}$$

Healthy canopy drives NBR strongly positive (high NIR, low SWIR); recently burned ground drives it toward zero or negative (NIR collapses, SWIR rises). Because NBR is a normalized difference bounded to $[-1,1]$, a single-date value is portable across scenes — but a single date cannot separate a genuinely low-NBR burn scar from a naturally low-NBR surface such as rock or water. That separation requires the difference of two dates.

The delta Normalized Burn Ratio (dNBR) subtracts the post-fire NBR from the pre-fire NBR, so that a large positive dNBR means NBR dropped sharply — the signature of severe burning:

$$\text{dNBR}={\text{NBR}}_{\text{pre}}-{\text{NBR}}_{\text{post}}$$

By convention (Key & Benson, USGS FIREMON), dNBR is scaled by 1000 and stored as an integer, so the published severity thresholds are expressed in scaled units (for example, an unburned pixel falls in roughly $-100$ to $+99$). dNBR is an absolute measure: a stand that starts at NBR 0.7 and burns to 0.1 yields the same dNBR as a stand that starts at 0.4 and burns to $-0.2$, even though the first lost proportionally less of its canopy.

To make severity relative to pre-fire condition, Miller & Thode introduced the Relativized dNBR (RdNBR), which normalizes dNBR by the square root of the pre-fire NBR magnitude:

$$\text{RdNBR}=\frac{\text{dNBR}}{\sqrt{|{\text{NBR}}_{\text{pre}}|}}$$

RdNBR reduces the tendency of dNBR to under-report severity in sparse, pre-disturbed, or low-biomass vegetation, at the cost of instability where ${\text{NBR}}_{\text{pre}}$ approaches zero (the denominator shrinks and RdNBR explodes). The engineering rule that follows from all three formulas is constant: cast to float before dividing, guard every denominator, carry masked pixels as NaN, and correct the phenological offset before you threshold.

Step-by-step Python pipeline

The following steps build a burn-severity mapper that loads paired scenes, computes NBR for each date, differences them, corrects the offset, and classifies the result. The narrow, band-level mechanics — scaling B8A and B12, the *1000 convention, and the exact threshold table — are documented step-by-step in the child guide Calculating dNBR from Sentinel-2 in Python.

Step 1 — Load paired bands as surface reflectance

Read the NIR and SWIR2 bands for both dates, cast to float32, and convert digital numbers to reflectance with the documented scale and offset. Both scenes must already share a CRS and grid; if they do not, warp them first through the Coordinate Reference Systems for Forestry workflow, because dNBR arithmetic assumes the pre and post arrays overlay pixel-for-pixel.

import numpy as np
import rasterio


def read_reflectance(path, scale=10000.0, offset=-1000.0):
    """Read a single-band raster as float32 surface reflectance."""
    with rasterio.open(path) as src:
        dn = src.read(1).astype("float32")
        profile = src.profile
    refl = (dn + offset) / scale
    # Surface reflectance is physically bounded to [0, 1]; clip stray values.
    return np.clip(refl, 0.0, 1.0), profile


# Pre- and post-fire Sentinel-2 B8A (NIR) and B12 (SWIR2), already co-registered.
nir_pre, prof = read_reflectance("pre_B8A.tif")
swir_pre, _ = read_reflectance("pre_B12.tif")
nir_post, _ = read_reflectance("post_B8A.tif")
swir_post, _ = read_reflectance("post_B12.tif")

Step 2 — Compute NBR for each date with a guarded denominator

NBR is a normalized difference, so it inherits the same zero-denominator hazard as NDVI: over deep water, shadow, or no-data fill, NIR + SWIR approaches zero and an unguarded division emits inf/nan that poisons every downstream statistic. Compute both dates through one guarded helper.

def nbr(nir, swir):
    """Normalized Burn Ratio with a guarded zero denominator."""
    num = nir - swir
    den = nir + swir
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(den == 0, np.nan, num / den)


nbr_pre = nbr(nir_pre, swir_pre)
nbr_post = nbr(nir_post, swir_post)

Step 3 — Mask cloud and shadow before differencing

A cloud that drifted over the burn between the two dates will masquerade as a severity signal — cloud raises SWIR less predictably than fire, and cloud shadow lowers both bands, so the dNBR over contaminated pixels is meaningless. Build a keep-mask for each date from the Scene Classification Layer (Sentinel-2) or the QA band (Landsat), resample it nearest-neighbour so class codes are never interpolated, and set contaminated pixels to NaN in both NBR rasters.

# SCL codes to drop: no-data(0), saturated(1), shadow(3), cloud med/high(8,9),
# thin cirrus(10). Apply the union of pre- and post-date contamination.
SCL_DROP = {0, 1, 3, 8, 9, 10}


def scl_keep_mask(scl_path, shape):
    from rasterio.enums import Resampling
    with rasterio.open(scl_path) as src:
        scl = src.read(1, out_shape=shape, resampling=Resampling.nearest)
    keep = np.ones(scl.shape, dtype=bool)
    for code in SCL_DROP:
        keep &= scl != code
    return keep


shape = nbr_pre.shape
keep = scl_keep_mask("pre_SCL.tif", shape) & scl_keep_mask("post_SCL.tif", shape)
nbr_pre = np.where(keep, nbr_pre, np.nan)
nbr_post = np.where(keep, nbr_post, np.nan)

Step 4 — Difference, offset-correct, and scale

Subtract post from pre to get dNBR. Then correct the phenological offset: two scenes captured weeks or months apart differ in sun angle, soil moisture, and green-up even where nothing burned, biasing dNBR by a small constant. Estimate that bias as the mean (or median) dNBR over an unburned reference area outside the fire perimeter, subtract it, and only then scale by 1000 into integer severity units.

def offset_corrected_dnbr(nbr_pre, nbr_post, unburned_mask):
    """dNBR corrected by the mean dNBR of an unburned reference region."""
    dnbr = nbr_pre - nbr_post
    offset = np.nanmean(dnbr[unburned_mask])  # phenological/illumination bias
    return dnbr - offset, offset


# unburned_mask: boolean array, True over stable vegetation outside the perimeter.
dnbr, offset = offset_corrected_dnbr(nbr_pre, nbr_post, unburned_mask)
dnbr_scaled = dnbr * 1000.0  # Key & Benson integer convention
print(f"phenological offset removed: {offset * 1000:.1f} scaled dNBR units")

Step 5 — Classify into Key & Benson severity classes

Threshold the scaled, offset-corrected dNBR into the standard severity classes. numpy.digitize maps each pixel to a class index in one vectorized pass; a small integer coding keeps the output compact and unambiguous for downstream mapping.

def classify_dnbr(dnbr_scaled):
    """Assign USGS Key & Benson severity classes from scaled dNBR.

    Class codes: 0 unburned/regrowth, 1 low, 2 moderate-low,
    3 moderate-high, 4 high. NaN inputs map to 255 (no-data).
    """
    # Lower edges of low, moderate-low, moderate-high, high.
    edges = [100.0, 270.0, 440.0, 660.0]
    classes = np.digitize(dnbr_scaled, edges).astype("uint8")  # 0..4
    classes[np.isnan(dnbr_scaled)] = 255
    return classes


severity = classify_dnbr(dnbr_scaled)

Optionally compute RdNBR alongside dNBR when the burn spans sparse or pre-disturbed vegetation, where dNBR under-reports severity:

def rdnbr(nbr_pre, dnbr):
    """Relativized dNBR (Miller & Thode). Unstable where NBR_pre -> 0."""
    denom = np.sqrt(np.abs(nbr_pre))
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(denom < 1e-3, np.nan, dnbr / denom)

Library configuration: a creation profile for the severity raster

Write the classified severity raster as a tiled, compressed uint8 GeoTIFF, and the continuous dNBR as float32. Categorical class codes must never be interpolated, so downstream consumers must read them nearest-neighbour.

# Creation profile for the classified severity raster (uint8)
driver: GTiff
dtype: uint8
nodata: 255          # reserve 255 for masked / no-data pixels
compress: deflate    # lossless
tiled: true
blockxsize: 256
blockysize: 256
BIGTIFF: IF_SAFER    # auto-promote for a large-fire mosaic

The nodata: 255 value keeps masked pixels out of every area tabulation; tiling matched to a 256-pixel block accelerates the zonal reads used when summarizing severity by ownership or watershed.

Validation and verification

A severity map earns trust only when it is checked numerically and against the ground. First, assert the continuous product is physically plausible — NBR must fall in $[-1,1]$ and offset-corrected dNBR should center near zero over unburned land.

import numpy as np


def verify_severity(nbr_pre, nbr_post, dnbr, unburned_mask):
    for name, arr in (("NBR pre", nbr_pre), ("NBR post", nbr_post)):
        finite = arr[np.isfinite(arr)]
        assert finite.min() >= -1.0 - 1e-6, f"{name} below -1: bad band order"
        assert finite.max() <= 1.0 + 1e-6, f"{name} above 1: integer division leaked"
    resid = np.nanmedian(dnbr[unburned_mask])
    assert abs(resid) < 0.02, f"unburned dNBR not centered ({resid:.3f}): re-check offset"
    return {"unburned_median_dnbr": float(resid)}

Then anchor the thresholds to field observations. Extract the modeled dNBR or RdNBR at each Composite Burn Index plot with the Raster-Vector Overlay Techniques workflow, regress the field CBI against the spectral value, and adjust the class breaks to your ecosystem — the published Key & Benson thresholds are a starting point, not a universal constant. Partition plots spatially rather than randomly when validating, or accuracy will look optimistic. Log the offset value, the class breaks, and the plot fit so the map is auditable for rehabilitation funding.

Failure modes and gotchas

• Skipping the phenological offset. Two dates always differ where nothing burned; without subtracting the unburned-area bias, an entire scene can shift a full severity class. Always offset-correct before thresholding.
• Mismatched pre/post grids. dNBR silently samples the wrong pixels if the two dates differ in CRS, extent, or resolution by even one pixel. Co-register first; assert identical shape and transform.
• Cloud read as severity. A cloud or its shadow between the dates produces a strong false dNBR. Mask the union of pre- and post-date contamination before differencing, not after.
• Integer truncation. Computing NBR on uint16 arrays floors the ratio to 0 or ±1. Cast to float32 and apply the reflectance scale/offset first — the single most common silent corruption.
• Blindly trusting default thresholds. Key & Benson breaks were derived in specific western-US forests; grassland, shrubland, and boreal systems need recalibration against local CBI plots.
• RdNBR blow-ups. Where pre-fire NBR is near zero (bare rock, water margins), the RdNBR denominator collapses and values explode. Guard the denominator and prefer dNBR in low-NBR terrain.

Performance and scale notes

• Window large mosaics. For a multi-tile fire, iterate src.block_windows(1) so each date is read one block at a time and peak memory tracks a tile, not the whole perimeter.
• Stay in float32. dNBR precision is limited by sensor radiometry, not mantissa bits; float64 doubles memory for no gain.
• Estimate the offset on a sample. Computing nanmean over a modest unburned sample is enough; you do not need every unburned pixel to fix a scene-wide constant.
• Reuse the mask. Build the combined cloud/shadow keep-mask once and apply it to both NBR rasters rather than re-reading the SCL per band.
• Landsat parity. The same pipeline runs on Landsat 8/9 with NIR = Band 5 and SWIR2 = Band 7; only the band numbers and reflectance scaling change, feeding the same wildfire risk modeling layers downstream.

Frequently Asked Questions

Related

• Fire Risk & Fuel Assessment — the parent section this severity workflow belongs to
• Calculating dNBR from Sentinel-2 in Python — the exact B8A/B12 band math, scaling, and threshold table
• Wildfire Risk Modeling in Python — where severity history feeds forward-looking risk layers
• Vegetation Index Calculation in Python — the shared NIR/SWIR normalized-difference band math
• Coordinate Reference Systems for Forestry — align the pre- and post-fire scenes before differencing