Calculating dNBR from Sentinel-2 in Python

Computing the delta Normalized Burn Ratio (dNBR) from a pre-fire and a post-fire Sentinel-2 scene is a narrow, repeatable task with three specific ways to get it wrong: using the wrong band pair, mishandling the reflectance scale and offset, and skipping the phenological offset correction that makes two dates comparable. This page is the exact rasterio + NumPy recipe for the Sentinel-2 case — narrow-NIR B8A and SWIR2 B12, both native 20 m — so the severity map you hand to a burned-area response team is real signal. It is one concrete step within the parent workflow Burn Severity Mapping with Python, which sits inside the wider Fire Risk & Fuel Assessment section. The hard parts are not the subtraction but the band choice, the scaling, and the offset discipline.

Why B8A and B12, and when to use this approach

The Normalized Burn Ratio contrasts near-infrared against shortwave-infrared: $\text{NBR}=({\rho }_{\text{NIR}}-{\rho }_{\text{SWIR}})/({\rho }_{\text{NIR}}+{\rho }_{\text{SWIR}})$. Sentinel-2 offers two NIR options, and the choice matters for burn severity:

For dNBR, B8A + B12 is the standard pairing: both bands are natively 20 m, so the two arrays already overlay pixel-for-pixel and no interpolation is baked into the index. Reach for this direct recipe when you control the band files and want an auditable, offset-corrected dNBR with explicit scaling. If your two scenes are not yet on a common grid — different CRS, tile, or resolution — align them first through Coordinate Reference Systems for Forestry; dNBR assumes the pre and post arrays already match.

Reflectance scaling and the ×1000 convention

Sentinel-2 L2A delivers B8A and B12 as 16-bit unsigned integers scaled by 10000. Products on processing baseline 04.00 or later (acquisitions from 25 January 2022 onward) also carry a radiometric BOA_ADD_OFFSET of -1000, so the correct conversion is reflectance = (DN - 1000) / 10000; pre-baseline products use the same scale with no offset. Direct uint16 arithmetic truncates the ratio, so cast to float32 before any math. Once NBR is a normalized difference in $[-1,1]$, dNBR is also small — typically a few tenths — so by the Key & Benson (USGS FIREMON) convention dNBR is multiplied by 1000 and stored as an integer, which is why the published severity thresholds read as values like 100, 270, 440, and 660 rather than 0.10, 0.27, 0.44, and 0.66.

Minimal reproducible example

The function below ingests the pre- and post-fire B8A and B12 GeoTIFFs, validates the grid, applies the scale and offset, masks zero denominators, differences the two NBRs, removes the phenological offset using an unburned reference mask, and returns scaled dNBR ready to threshold.

import numpy as np
import rasterio


def _reflectance(path, offset=-1000.0, scale=10000.0):
    """Read a Sentinel-2 L2A band as float32 surface reflectance."""
    with rasterio.open(path) as src:
        dn = src.read(1).astype("float32")
        profile = src.profile
    return np.clip((dn + offset) / scale, 0.0, 1.0), profile


def _nbr(nir, swir):
    """NBR = (NIR - SWIR) / (NIR + SWIR) with a guarded denominator."""
    den = nir + swir
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(den == 0, np.nan, (nir - swir) / den)


def dnbr_sentinel2(nir_pre, swir_pre, nir_post, swir_post, unburned_mask):
    """Offset-corrected, ×1000 scaled dNBR from pre/post B8A and B12 paths.

    unburned_mask: boolean array, True over stable vegetation outside the
    fire perimeter, used to estimate the phenological offset.
    """
    nbr_p, prof = _reflectance(nir_pre)
    swir_p, _ = _reflectance(swir_pre)
    nbr_q, _ = _reflectance(nir_post)
    swir_q, _ = _reflectance(swir_post)

    # Grid check: all four bands must share shape and transform.
    if not (nbr_p.shape == swir_p.shape == nbr_q.shape == swir_q.shape):
        raise ValueError("Pre/post B8A and B12 bands have mismatched shapes.")

    nbr_pre = _nbr(nbr_p, swir_p)
    nbr_post = _nbr(nbr_q, swir_q)

    dnbr = nbr_pre - nbr_post
    offset = np.nanmean(dnbr[unburned_mask])  # phenological/illumination bias
    dnbr_corrected = dnbr - offset

    return (dnbr_corrected * 1000.0).astype("float32"), offset, prof

Handling the SCL cloud mask with the standard approach makes the offset estimate honest — a cloud left in the unburned reference region would corrupt the mean:

from rasterio.enums import Resampling

SCL_DROP = {0, 1, 3, 8, 9, 10}  # no-data, sat, shadow, cloud med/high, cirrus


def scl_keep(scl_path, shape):
    """Nearest-neighbour keep-mask (True = valid) from a 20 m SCL band."""
    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

Parameter reference

Key & Benson dNBR severity thresholds (scaled ×1000)

These breaks were derived in western-US conifer forests; treat them as a starting point and recalibrate against local Composite Burn Index plots for other biomes.

Expected output and verification

A correct run returns a float32 array of scaled dNBR on the same 20 m grid as the inputs, centered near zero over unburned land and strongly positive over the burn. Assert the physical bounds and the offset residual before trusting it:

import numpy as np


def verify_dnbr(dnbr_scaled, unburned_mask):
    finite = dnbr_scaled[np.isfinite(dnbr_scaled)]
    assert finite.size > 0, "dNBR is entirely NaN — check masks and scaling"
    # Scaled dNBR is bounded by ±2000 (two normalized differences × 1000).
    assert finite.min() >= -2000 and finite.max() <= 2000, "out of range: bad scaling"
    resid = np.nanmedian(dnbr_scaled[unburned_mask])
    assert abs(resid) < 20, f"unburned dNBR not centered ({resid:.1f}): re-check offset"
    print(f"unburned median scaled dNBR: {resid:.1f}")
    return float(resid)

An unburned median far from zero means the phenological offset was estimated over a contaminated or partly-burned reference region; an out-of-range extreme usually means the reflectance offset was wrong or the band order was swapped.

Common pitfalls

• Unequal scene extents or resolution. B8A and B12 are both 20 m, so keep every band at 20 m; do not silently resample one to 10 m, which would sample the wrong pixels. Assert identical shape and transform across all four inputs.
• Skipping the offset correction. dNBR from two dates carries a cross-date phenology and illumination bias. Subtract the mean dNBR of an unburned area before scaling, or the whole scene can shift a severity class.
• Cloud masking with the wrong resample. The SCL is categorical — resample it nearest-neighbour to 20 m. Bilinear invents fractional class codes and leaks cloud pixels into both the burn and the unburned reference.
• Forgetting or double-applying BOA_ADD_OFFSET. Applying -1000 to a pre-baseline product, or omitting it on a newer one, biases reflectance and drifts every NBR. Read the baseline from MTD_MSIL2A.xml.

Frequently Asked Questions

Related

• Burn Severity Mapping with Python — the parent workflow this dNBR step belongs to
• Fire Risk & Fuel Assessment — the full fire section this recipe sits within
• Vegetation Index Calculation in Python — the shared normalized-difference band math
• Coordinate Reference Systems for Forestry — align the pre- and post-fire scenes before differencing