Module 2: Optical Radiometry and Burn Severity

Phase: I, Foundations of the data, physics, and radar
Level: Intermediate
Estimated time: ~12 hours (about 5 h concepts, 5 h lab, 2 h write-up)
Prerequisites: Module 1 (cube.py and the cloud-native stack), Phase 0.
Portfolio thread: section 2 of Portfolio Project A. You produce a validated burn-severity map and an honest accuracy table against MTBS.

Why this module

dNBR is the lingua franca of post-fire assessment, and plenty of candidates can compute it. Far fewer can explain why near-infrared falls and shortwave-infrared rises when vegetation burns, why a top-of-atmosphere dNBR is a different product from a surface-reflectance one, or how to validate a burn map honestly instead of quoting a misleading pixel-count accuracy. This module gives you the radiometric physics and the accuracy discipline that make a burn-severity map credible. The honest-accuracy skill in particular (area-weighted estimates with confidence intervals) shows up in interviews and in every serious mapping deliverable.

Learning outcomes

After this module you can:

Explain the chain from radiance to top-of-atmosphere reflectance to surface reflectance, and say what a Level-2A product has and has not corrected.
Explain the spectral physics of burning: why healthy vegetation is bright in NIR, and why char and exposed soil raise SWIR and drop NIR.
Compute NBR, dNBR, and RdNBR correctly, including the Sentinel-2 baseline offset and QA masking, and classify to Key and Benson severity.
Validate a map against MTBS with a confusion matrix and area-weighted accuracy with confidence intervals, following Olofsson and colleagues.
Explain what HLS is and why it exists.
Ship a tested indices.py and a burn-severity notebook with an honest accuracy table.

Concept lessons (about 5 hours)

Lesson 1: From radiance to surface reflectance

A sensor measures radiance, the energy arriving at the detector. That is converted to top-of-atmosphere (TOA) reflectance, which normalizes for Sun geometry but still includes the atmosphere's scattering and absorption. Atmospheric correction removes the atmosphere to give surface reflectance (Sen2Cor for Sentinel-2, LaSRC for Landsat and HLS). A Level-2A product is surface reflectance with a quality and cloud layer, and it is what you compare across dates, because two TOA images on different days differ partly because of the atmosphere, not the ground.

One concrete gotcha to internalize now: since processing baseline 04.00 (January 2022), Sentinel-2 Level-2A digital numbers carry a radiometric offset. Surface reflectance is (DN + BOA_ADD_OFFSET) / 10000, and the offset has been -1000, so reflectance = (DN - 1000) / 10000. Some catalogs expose "harmonized" collections that undo this so older and newer scenes line up. If you forget the offset, every index you compute is quietly wrong. Always check the processing baseline and whether your source is harmonized.

Lesson 2: The spectral physics of burning

Healthy vegetation has a distinctive signature: high reflectance in the near-infrared (NIR), because leaf internal structure scatters it strongly, and moderate shortwave-infrared (SWIR), which is suppressed by leaf water content. The Normalized Burn Ratio exploits exactly this contrast:

NBR = (NIR - SWIR2) / (NIR + SWIR2)

Healthy vegetation gives a high NBR. After a fire, chlorophyll and leaf structure are destroyed, so NIR drops, while char, ash, and exposed dry soil raise SWIR reflectance. NBR therefore falls sharply. The difference between a pre-fire and post-fire image captures the change:

dNBR = NBR_pre - NBR_post

A larger dNBR means a bigger drop, which means more severe burning. RdNBR (relative dNBR) divides by a function of the pre-fire NBR to avoid underestimating severity where pre-fire vegetation was already sparse:

RdNBR = dNBR / sqrt(abs(NBR_pre))

The habit to keep from Module 0.2: name the physics. When you compute dNBR you are relying on the NIR-and-SWIR contrast, so anything that mimics that contrast without being fire (water, cloud shadow, bare sandbars) will produce false severity unless masked.

Lesson 3: Severity thresholds and their limits

Key and Benson defined dNBR ranges that map to severity classes (roughly: enhanced regrowth, unburned, low, moderate-low, moderate-high, high). Indicative breakpoints in reflectance units are about: unburned around -0.1 to 0.1, low 0.1 to 0.27, moderate-low 0.27 to 0.44, moderate-high 0.44 to 0.66, and high above 0.66. Two cautions that separate careful practitioners from careless ones:

These thresholds are region and vegetation dependent. A fixed global threshold is a red flag; real work calibrates them to the fire's setting, ideally with field data. MTBS publishes the exact thresholds per fire in its metadata for this reason.
Water, cloud, and cloud shadow create spurious severe values (the classic "severe burn in a lake"). Mask them before classifying, using the product QA (Fmask for HLS, the Scene Classification Layer for Sentinel-2).

And remember a TOA dNBR is a different product from a surface-reflectance dNBR; do not mix them.

Lesson 4: Honest accuracy assessment

Reporting the percentage of pixels that matched a reference map is misleading, because map class areas are themselves biased by classification error. The defensible approach, from Olofsson and colleagues' good-practice guidance, is: draw a probability sample, build a confusion matrix of your map against reference labels, then convert counts to area-weighted proportions using each mapped class's area, and report overall, user's, and producer's accuracy with confidence intervals, plus error-adjusted class areas. This is the number you cite. You will apply it by comparing your dNBR classes against the matching MTBS product, which is itself a Landsat-based burn-severity reference for US fires.

Why HLS exists

Harmonized Landsat and Sentinel-2 (HLS) resamples and cross-calibrates Landsat (L30) and Sentinel-2 (S30) to a common 30 m grid and consistent surface reflectance, so you get a denser, sensor-consistent optical time series than either alone. It is the natural input for burn-severity work when you want frequent, comparable observations.

Guided lab (about 5 hours): a validated burn-severity map

Step 1: Pre-fire and post-fire composites

Reuse cube.py from Module 1 to build cloud-masked surface-reflectance composites just before and just after your study-area fire, over the same AOI. Use HLS (L30 and S30) for a dense, harmonized series, or Sentinel-2 L2A if you prefer, applying the baseline offset and the SCL mask. Keep the NIR and SWIR2 bands.

Step 2: `indices.py` with tests

"""indices.py: spectral indices for burn severity."""
from __future__ import annotations
import numpy as np

def nbr(nir, swir2):
    """Normalized Burn Ratio."""
    return (nir - swir2) / (nir + swir2)

def dnbr(nbr_pre, nbr_post):
    """Pre minus post; higher means more severe."""
    return nbr_pre - nbr_post

def rdnbr(nbr_pre, nbr_post, eps=1e-6):
    """Relative dNBR (normalizes for pre-fire condition)."""
    return (nbr_pre - nbr_post) / np.sqrt(np.abs(nbr_pre) + eps)

# Key & Benson indicative classes on dNBR (reflectance units). Calibrate per fire.
# Code 255 is reserved for nodata; do not reuse it as a severity class.
_NODATA = 255
_BREAKS = [(-np.inf, -0.1, 0),   # 0 enhanced regrowth / unburned-low
           (-0.1, 0.1, 1),        # 1 unburned
           (0.1, 0.27, 2),        # 2 low
           (0.27, 0.44, 3),       # 3 moderate-low
           (0.44, 0.66, 4),       # 4 moderate-high
           (0.66, np.inf, 5)]     # 5 high

def classify_severity(d, nodata=_NODATA):
    """Map a dNBR array to Key & Benson severity class codes 0..5.

    NaN or masked pixels map to `nodata` (default 255), never to class 0, so a
    genuine enhanced-regrowth pixel (dNBR < -0.1, class 0) is not confused with
    missing data. This matters because you mask water and cloud before classifying,
    which introduces NaN.
    """
    d = np.asarray(d, dtype="float64")
    out = np.full(np.shape(d), nodata, dtype="uint8")
    valid = ~np.isnan(d)
    for lo, hi, code in _BREAKS:
        out = np.where(valid & (d >= lo) & (d < hi), code, out)
    return out

# tests/test_indices.py
import numpy as np
from eo_portfolio.indices import nbr, dnbr, rdnbr, classify_severity

def test_nbr_known_value():
    assert abs(nbr(np.array(0.4), np.array(0.1)) - 0.6) < 1e-9

def test_burn_lowers_nbr():
    pre = nbr(np.array(0.45), np.array(0.10))   # healthy: high NBR
    post = nbr(np.array(0.20), np.array(0.30))  # burned: low/negative NBR
    assert dnbr(pre, post) > 0

def test_classify_bands():
    d = np.array([0.05, 0.2, 0.5, 0.8])
    assert list(classify_severity(d)) == [1, 2, 4, 5]

def test_classify_nan_is_nodata():
    # NaN must become nodata (255), not class 0; -0.2 is genuine class 0 (regrowth).
    d = np.array([np.nan, -0.2, 0.05])
    assert list(classify_severity(d)) == [255, 0, 1]

Step 3: Map and classify

Compute NBR pre and post, then dNBR and RdNBR, mask water and cloud, and classify to severity. Plot the dNBR map and the classified map over the AOI.

Step 4: Validate against MTBS (area-weighted)

Pull the matching MTBS fire (its dNBR and thematic severity) from mtbs.gov or a STAC catalog that hosts it. Build a confusion matrix of your classes against MTBS, then compute area-weighted accuracy:

# accuracy.py
import numpy as np

def area_adjusted_accuracy(cm, mapped_area):
    """Olofsson-style area-weighted accuracy from a confusion matrix.

    cm[i][j]: sample count, rows = your map class i, cols = reference class j.
    mapped_area[i]: total mapped area (or pixel count) of class i.
    Returns overall, user's, and producer's accuracy on an area basis.
    """
    cm = np.asarray(cm, float); W = np.asarray(mapped_area, float)
    W = W / W.sum()
    ni = cm.sum(axis=1)
    p = W[:, None] * (cm / ni[:, None])         # area-weighted proportions
    overall = float(np.trace(p))
    users = np.diag(p) / p.sum(axis=1)
    producers = np.diag(p) / p.sum(axis=0)
    return {"overall": overall, "users": users, "producers": producers}

# tests/test_accuracy.py
import numpy as np
from eo_portfolio.accuracy import area_adjusted_accuracy

def test_perfect_map_scores_one():
    cm = [[10, 0], [0, 5]]
    r = area_adjusted_accuracy(cm, [100, 50])
    assert abs(r["overall"] - 1.0) < 1e-9

def test_known_error_case():
    r = area_adjusted_accuracy([[9, 1], [0, 5]], [100, 50])
    assert abs(r["overall"] - 0.9333) < 1e-3

Extend it with the error-adjusted area and 95 percent confidence intervals from the Olofsson paper; those are the numbers you actually quote in the write-up.

Step 5: Write-up and commit

A one-page results note comparing your dNBR classes to MTBS, with the area-weighted accuracy table and an honest discussion of where they disagree and why (threshold calibration, date offsets, water and shadow). Commit indices.py, accuracy.py, tests, and the notebook.

git add src/eo_portfolio/indices.py src/eo_portfolio/accuracy.py tests/ notebooks/04_burn_severity.ipynb
git commit -m "Module 2: burn-severity indices + MTBS-validated accuracy"
git push

Checkpoint

Self-check (answers below).

Physically, why does NBR fall after a fire?
You computed dNBR from raw Sentinel-2 L2A scenes from 2023 and the values look wrong. What did you likely forget?
Why is a pixel-count accuracy against MTBS misleading, and what do you report instead?
Your dNBR shows severe burn over a reservoir. What happened and how do you fix it?

Interview-style questions (practice out loud).

A dNBR map shows severe burn in a lake. Walk me through the failure.
Your random-pixel accuracy is 94 percent but your area-weighted accuracy is 81 percent. Which do you report and why?
Why is a top-of-atmosphere dNBR a different product from a surface-reflectance dNBR?
How would you calibrate severity thresholds for a fire in a vegetation type you have not worked in?

Answers. (1) Fire destroys chlorophyll and leaf structure so NIR reflectance drops, while char and exposed soil raise SWIR reflectance; since NBR is (NIR - SWIR2)/(NIR + SWIR2), both effects push it down.
(2) The processing baseline 04.00 reflectance offset: surface reflectance is (DN - 1000)/10000, or use a harmonized collection.
(3) Map class areas are biased by classification error, so a raw pixel-count match overstates accuracy; report area-weighted overall, user's, and producer's accuracy with confidence intervals and error-adjusted areas (Olofsson).
(4) Water mimics the NIR-low, SWIR-low pattern and was not masked; apply the QA and a water mask before classifying.

Deliverable

A tested indices.py (NBR, dNBR, RdNBR, severity classification) and an area-weighted accuracy helper, plus notebooks/04_burn_severity.ipynb producing a classified burn-severity map over your fire AOI and a one-page note comparing it to MTBS with an honest, area-weighted accuracy table. This completes section 2 of Portfolio Project A.

What a hiring manager sees

Anyone can divide two bands. Showing that you understand the radiometry (surface versus TOA, the baseline offset), that you mask water and cloud before classifying, and above all that you validate with area-weighted accuracy and confidence intervals rather than a flattering pixel count, is what makes a reviewer trust your maps. The Olofsson-style accuracy table is a small thing that signals real rigor.

Currency note

Verified August 2026. HLS version 2.0 (L30 and S30) is distributed as COGs from NASA's LP DAAC in the Earthdata Cloud and is STAC-searchable with an Earthdata login; confirm the collection names and access method at study time. MTBS remains available at mtbs.gov (and via STAC and Earth Engine mirrors), with per-fire dNBR and published thresholds in the metadata. The Sentinel-2 Level-2A baseline 04.00 offset (BOA_ADD_OFFSET, so far -1000, giving surface reflectance = (DN - 1000)/10000) applies to scenes from January 2022 onward unless you use a harmonized collection; verify the baseline of any scene you load. Key and Benson thresholds are indicative and must be calibrated per fire. Pin all versions in your Module 0.1 environment file.