Forestry GIS Python · Forestry · Ecology

Classifying Recovery Classes from Multi-Year dNBR in Python

A continuous recovery surface is the right analytical product and the wrong management one. Decisions are categorical — replant, monitor, do nothing — so at some point the trajectory has to become a class, and where the thresholds fall determines how many hectares get a planting budget. This guide makes that conversion explicit and auditable, inside Post-Fire Vegetation Recovery Monitoring, part of Fire Risk & Fuel Assessment.

When to use this approach

Product Form Audience
Relative recovery raster Continuous 0–1 Analysts; the input to everything else
Recovery class raster 4–5 ordinal classes Operations; drives the field programme
Per-stand class One class per polygon Planning and budgeting
Recovery rate raster Continuous, years⁻¹ Modellers and monitoring design
Reburn hazard flag Boolean Fire management, as in Detecting Reburn Risk from Recovery and Fuel Rasters

Produce the continuous surface first and keep it. Classes are a view of it, and when the thresholds are revisited — and they always are — you want to reclassify rather than reprocess.

Relative recovery from a dNBR series

The severity index measures the drop; the recovery index measures how much of that drop has closed. With the pre-fire baseline and the first post-fire year:

The second form is the useful one operationally: if you already produce a dNBR raster each year against the same pre-fire baseline, the recovery index is one division. It also makes the failure mode obvious — where the initial dNBR is small, the denominator is small, and the ratio explodes. Lightly burned pixels must be excluded rather than divided.

Why lightly burned pixels must be excluded from the recovery ratio Two panels. On the left, a high severity pixel with a large initial dNBR: each year's dNBR shrinks steadily and the relative recovery index rises smoothly from zero towards one. On the right, a low severity pixel whose initial dNBR is barely above the noise: small year to year fluctuations in the numerator produce a recovery index that swings between minus one and plus two, because the denominator is close to zero. High severity — dNBR₀ = 0.72 Low severity — dNBR₀ = 0.06 RRI = 1 0 1 years since fire RRI = 1 −1 1 2 years since fire The numerator noise is identical in both panels. Only the denominator differs. Mask pixels whose initial dNBR is below the severity detection threshold before dividing.

Minimal reproducible example

import numpy as np
import rasterio

DNBR_MIN = 0.10          # below this the burn is not reliably detected


def relative_recovery(dnbr_initial: np.ndarray, dnbr_year: np.ndarray) -> np.ndarray:
    """RRI = 1 - dNBR_t / dNBR_initial, masked where the denominator is unreliable."""
    out = np.full(dnbr_initial.shape, np.nan, dtype="float32")
    usable = np.isfinite(dnbr_initial) & np.isfinite(dnbr_year) & (dnbr_initial >= DNBR_MIN)
    out[usable] = 1.0 - (dnbr_year[usable] / dnbr_initial[usable])
    # RRI above 1 is real (overshoot); below -0.5 is noise or a reburn.
    out[out < -0.5] = np.nan
    return np.clip(out, -0.5, 2.0)


def classify(rri: np.ndarray, rate: np.ndarray) -> np.ndarray:
    """Ordinal recovery classes. 0 = unclassified, 1 = none … 4 = rapid."""
    cls = np.zeros(rri.shape, dtype="uint8")
    valid = np.isfinite(rri)
    cls[valid] = 1
    cls[valid & (rri >= 0.25)] = 2
    cls[valid & (rri >= 0.50)] = 3
    cls[valid & (rri >= 0.75)] = 4
    # Demote anything whose fitted rate says it has stalled, whatever its level.
    stalled = valid & np.isfinite(rate) & (rate < 0.08) & (rri < 0.75)
    cls[stalled] = np.minimum(cls[stalled], 2)
    return cls

The demotion rule is what stops the class map from being a repackaged snapshot. A pixel that reached 60% in the first two years and has not moved since is a different management case from one climbing steadily through 60%, and only the rate distinguishes them.

Clean up to a minimum mapping unit

A raw class raster is speckled, and a field crew cannot act on a 400 m² island of “replant” inside a recovering stand. A majority filter followed by a small-region merge produces something operational.

from scipy import ndimage as ndi


def apply_mmu(cls: np.ndarray, cell_m: float = 20.0, mmu_ha: float = 0.5) -> np.ndarray:
    """Merge regions smaller than the minimum mapping unit into their neighbours."""
    min_px = int(np.ceil(mmu_ha * 10_000.0 / (cell_m ** 2)))
    out = ndi.median_filter(cls, size=3, mode="nearest")

    for value in np.unique(out):
        if value == 0:
            continue
        labels, n = ndi.label(out == value)
        if n == 0:
            continue
        sizes = ndi.sum(np.ones_like(labels), labels, index=range(1, n + 1))
        small = np.isin(labels, [i + 1 for i, s in enumerate(sizes) if s < min_px])
        if small.any():
            # Replace each small island with the modal class of its dilated ring.
            ring = ndi.binary_dilation(small, np.ones((3, 3), bool)) & ~small
            fill = np.bincount(out[ring], minlength=int(out.max()) + 1)
            fill[0] = 0
            out[small] = int(fill.argmax())
    return out

State the minimum mapping unit alongside the map. Half a hectare at 20 m is 13 pixels, which removes speckle without erasing genuine small openings; a 2 ha unit produces a much tidier map that has silently deleted every seed-tree gap.

Parameter reference

Parameter Type Default Range Rationale
DNBR_MIN float 0.10 0.08–0.15 The severity detection floor; below it the recovery ratio is undefined
RRI class breaks tuple 0.25 / 0.50 / 0.75 — Even quartiles; adjust to local regeneration standards and say so
stalled rate float 0.08 0.05–0.12 Demotion threshold; a half-life beyond ~9 years is not a trajectory
mmu_ha float 0.5 0.25–2.0 Smallest reportable patch; larger values erase real openings
assessment year int 5 3–10 Year 3 is too early in slow systems; year 10 is too late to intervene cheaply
median filter size int 3 3–5 Removes single-pixel speckle without moving class boundaries

Expected output and verification

A class map is validated the way any classification is: against independent observations, with a confusion matrix, and reporting the errors by class rather than one overall figure.

import numpy as np


def confusion(observed: np.ndarray, predicted: np.ndarray, n_classes: int = 5) -> dict:
    """Confusion matrix plus per-class producer and user accuracy."""
    cm = np.zeros((n_classes, n_classes), dtype="int64")
    for o, p in zip(observed, predicted):
        cm[int(o), int(p)] += 1
    total = cm.sum()
    correct = np.trace(cm)
    producer = np.divide(np.diag(cm), cm.sum(axis=1),
                         out=np.zeros(n_classes), where=cm.sum(axis=1) > 0)
    user = np.divide(np.diag(cm), cm.sum(axis=0),
                     out=np.zeros(n_classes), where=cm.sum(axis=0) > 0)
    # Cohen's kappa: agreement beyond what class prevalence alone would give.
    pe = float((cm.sum(axis=0) * cm.sum(axis=1)).sum()) / (total ** 2)
    po = correct / total
    return {
        "overall": round(float(po), 3),
        "kappa": round(float((po - pe) / (1 - pe)), 3) if pe < 1 else float("nan"),
        "producer": [round(float(x), 3) for x in producer],
        "user": [round(float(x), 3) for x in user],
        "n": int(total),
    }


res = confusion(plot_class_observed, plot_class_predicted)
assert res["overall"] > 0.6, f"class map is not usable: {res}"
assert res["kappa"] > 0.4, f"agreement no better than prevalence would give: {res}"

Kappa is the assertion that matters. In a fire where 70% of the area recovered rapidly, always predicting “rapid” scores 70% overall accuracy and is worthless; kappa near zero exposes that immediately, and overall accuracy never will.

Report the per-class figures separately, because the operationally important class is usually the rarest. A map with 85% overall accuracy that finds only half the genuinely non-recovering stands has failed at the one job it was built for.

Confusion matrix for a four-class recovery map A four by four confusion matrix comparing mapped recovery class against field survey class. Counts on the diagonal are much larger than off it, and almost all disagreement falls in cells adjacent to the diagonal, meaning the map confuses neighbouring classes rather than mistaking rapid recovery for none. The rarest class, not recovering, has the lowest producer accuracy and is annotated as the class the map is weakest at. none slow moderate rapid mapped class none slow moderate rapid field class 14710 53161 18449 021160 64% producer 72% 71% 82% Overall 75%, kappa 0.65 — but the rarest and most consequential class is the weakest. What each minimum mapping unit removes Three panels of the same recovery class map. The raw classification is heavily speckled with single-pixel islands. At a half hectare minimum mapping unit the speckle is gone and the genuine small openings around surviving seed trees remain. At a two hectare unit the map is very tidy and those openings have disappeared entirely, taking with them the pattern a regeneration survey would target. raw 0.5 ha unit 2 ha unit speckled openings kept uniform — and wrong The tidiest map is not the most informative one; state the unit so a reader knows what was removed.

Common pitfalls

  • Dividing by a small initial dNBR. Lightly burned pixels produce recovery values in the hundreds. Mask below the detection floor rather than clipping afterwards.
  • Thresholds chosen after seeing the hectare totals. Set the breaks from a regeneration standard or a prior study, record them, and then compute the areas.
  • A minimum mapping unit large enough to erase the signal. In a fire whose regeneration is patchy at the 0.2 ha scale, a 2 ha unit maps a uniformity that is not there.
  • Overall accuracy as the headline. With an imbalanced class distribution it is uninformative. Lead with kappa and per-class producer accuracy.
  • Reclassifying from the class raster. Always reclassify from the continuous surface; going class-to-class compounds every earlier threshold decision.

Frequently Asked Questions

How many recovery classes should I use?

Four is the usual compromise: none, slow, moderate, rapid. Fewer loses the distinction between “watch” and “act”; more produces classes that field crews cannot separate reliably, which shows up immediately as off-diagonal mass in the confusion matrix.

Should classes be defined on the recovery level or on the rate?

On both, with the level leading and the rate demoting. The level is what has been measured; the rate is what predicts the future. A rule that uses only one of them will either flag every slow-but-healthy stand or miss every stalled one.

What if I have no field plots to validate against?

Use whatever independent evidence exists — aerial photography interpreted at points, a regeneration survey from a neighbouring fire in the same forest type, or operational records of which stands were replanted. State clearly that the accuracy figure comes from a proxy. A class map with no validation at all should be published as provisional, not as a product.

Up: Post-Fire Recovery Monitoring · Fire Risk & Fuel Assessment