Detecting Reburn Risk from Recovery and Fuel Rasters in Python

A burned stand is not automatically safe. Within three to eight years the fine surface fuel that a fire consumed has usually rebuilt as grass, shrub and standing dead material — often more continuous and more flammable than what was there before — while the canopy that would have shaded and moderated it has not returned. That window is when a scar reburns, and reburns are disproportionately severe because there is no mature overstory left to lose slowly. This guide builds a hazard surface that finds those stands, inside Post-Fire Vegetation Recovery Monitoring, part of Fire Risk & Fuel Assessment.

When to use this approach

Reburn hazard is a susceptibility surface, not a probability. It says which parts of a scar would carry fire well; turning that into an annual likelihood needs ignition and weather, which belong to the risk model.

The three components

The hazard is a conjunction, not a sum of independent effects. A stand is reburn-prone when all three hold: recovery has stalled, surface fuel has rebuilt, and topography or wind exposure would carry fire.

$$H={\left(1-R\right)}^{w_r}\cdot F^{w_f}\cdot E^{w_e}$$

with $R$ the relative recovery, $F$ the normalised surface fuel index, $E$ the topographic exposure, and the exponents setting how strongly each term dominates. A multiplicative form is deliberate: a stand with abundant fuel but a closed recovering canopy should score low, and an additive index would let the fuel term carry it.

Minimal reproducible example

Each component is normalised to 0–1 before combining, and each normalisation is a modelling decision worth stating.

import numpy as np
import rasterio


def normalise(arr: np.ndarray, lo: float, hi: float) -> np.ndarray:
    """Linear scale to 0-1 between stated breakpoints, clipped at both ends."""
    out = (arr - lo) / (hi - lo)
    return np.clip(out, 0.0, 1.0).astype("float32")


def reburn_hazard(rri: np.ndarray, surface_fuel_t_ha: np.ndarray,
                  slope_deg: np.ndarray, exposure: np.ndarray,
                  w_r: float = 1.0, w_f: float = 1.0, w_e: float = 0.6) -> np.ndarray:
    """Multiplicative reburn susceptibility in 0-1."""
    stalled = 1.0 - normalise(rri, 0.0, 0.8)              # 0 recovery -> 1.0
    fuel = normalise(surface_fuel_t_ha, 2.0, 15.0)        # t/ha of fine surface fuel
    slope_term = normalise(slope_deg, 0.0, 35.0)
    expose = np.clip(0.5 * slope_term + 0.5 * normalise(exposure, 0.0, 1.0), 0.0, 1.0)

    h = (stalled ** w_r) * (fuel ** w_f) * (expose ** w_e)
    h[~np.isfinite(rri) | ~np.isfinite(surface_fuel_t_ha)] = np.nan
    return h.astype("float32")

The fuel breakpoints carry most of the model’s behaviour. Two tonnes per hectare of fine surface fuel is around the threshold below which fire will not carry continuously in most temperate systems; 15 t/ha is a heavily loaded shrub or slash bed. Those numbers are system-specific and should be replaced with local fuel-model values rather than inherited.

The exposure input is whatever wind-exposure proxy you have — a topographic exposure index, a distance-to-ridge measure, or a modelled wind multiplier. Where none exists, slope alone is a reasonable stand-in, which is why it is weighted more lightly than the other two terms.

Parameter reference

Test the weights before trusting the map

A weighted index looks authoritative and is usually more sensitive to its weights than to its data. Vary each weight across its plausible range and measure how much the high-hazard area moves; if it swings by a factor of two, the weights are the model.

import itertools

import numpy as np


def weight_sensitivity(rri, fuel, slope, exposure, threshold: float = 0.5) -> list[dict]:
    """How much does the high-hazard area move as the weights vary?"""
    rows = []
    for w_r, w_f, w_e in itertools.product((0.8, 1.0, 1.4), repeat=3):
        h = reburn_hazard(rri, fuel, slope, exposure, w_r=w_r, w_f=w_f, w_e=w_e)
        rows.append({
            "w_r": w_r, "w_f": w_f, "w_e": w_e,
            "high_frac": round(float(np.nanmean(h >= threshold)), 4),
        })
    fr = [r["high_frac"] for r in rows]
    print(f"high-hazard fraction ranges {min(fr):.3f}-{max(fr):.3f} "
          f"({max(fr) / max(min(fr), 1e-6):.1f}x) across weight combinations")
    return rows

Report that range with the map. A hazard surface whose high class covers 8% of the scar under one weighting and 24% under another is not wrong, but it is a scenario rather than a measurement, and presenting it as the latter is how a planning process ends up defending an arbitrary exponent.

Expected output and verification

The honest validation is retrospective: build the hazard surface as it would have looked before a reburn that has since occurred, and check whether the area that actually reburned scored higher than the area that did not.

import numpy as np


def validate_against_reburn(hazard: np.ndarray, reburned: np.ndarray) -> dict:
    """Compare hazard scores inside and outside an observed reburn footprint."""
    ok = np.isfinite(hazard)
    inside = hazard[ok & reburned]
    outside = hazard[ok & ~reburned]
    if inside.size < 50 or outside.size < 50:
        return {"status": "insufficient_sample"}

    # AUC via the Mann-Whitney relationship: P(hazard_inside > hazard_outside).
    order = np.argsort(np.concatenate([inside, outside]))
    ranks = np.empty(order.size, dtype="float64")
    ranks[order] = np.arange(1, order.size + 1)
    r_in = ranks[: inside.size].sum()
    auc = (r_in - inside.size * (inside.size + 1) / 2) / (inside.size * outside.size)
    return {
        "median_inside": round(float(np.median(inside)), 3),
        "median_outside": round(float(np.median(outside)), 3),
        "auc": round(float(auc), 3),
        "n_inside": int(inside.size),
    }


res = validate_against_reburn(hazard_2019, reburn_2024_mask)
assert res["auc"] > 0.6, f"hazard surface does not discriminate reburned ground: {res}"

An AUC near 0.5 means the surface is no better than chance and the weights need calibrating against data rather than intuition. An AUC above about 0.75 is a strong result for a susceptibility index built from three coarse components.

Where no reburn has occurred, the surface can still be sanity-checked against its own logic — the machinery must reject inputs it should refuse:

n = 1000
rng = np.random.default_rng(20260809)

# Fully recovered canopy: hazard must be near zero regardless of fuel.
h = reburn_hazard(np.full(n, 0.95), np.full(n, 20.0),
                  np.full(n, 30.0), np.full(n, 0.9))
assert np.nanmax(h) < 0.1, "a fully recovered stand scored as high hazard"

# No fuel: hazard must be zero regardless of recovery and exposure.
h = reburn_hazard(np.zeros(n), np.full(n, 0.5),
                  np.full(n, 35.0), np.ones(n))
assert np.nanmax(h) < 0.05, "a fuel-free stand scored as hazardous"

Common pitfalls

• Adding the components instead of multiplying them. An additive index lets a single high term carry a pixel, which produces high hazard on bare mineral soil with no fuel at all.
• Using a fuel layer from before the fire. The whole point is the fuel that has rebuilt since. A pre-fire fuel raster describes vegetation that no longer exists.
• Applying the model outside the reburn window. Before roughly year three there is not enough fuel; after canopy closure the surface fuel is shaded and moist. Outside that band the index means little.
• Treating susceptibility as probability. Without ignition and weather there is no annual likelihood, and presenting a 0–1 index as a probability invites exactly that misreading.
• Never testing the weights. An unexamined weighted index is an opinion with a raster attached.

Frequently Asked Questions

Related

• Post-Fire Vegetation Recovery Monitoring — supplies the recovery term
• Fuel Load Mapping from LiDAR — supplies the fuel term
• Wildfire Risk Modeling in Python — turns susceptibility into risk with ignition and weather
• Computing Slope, Aspect and TRI from a Forest DTM — supplies the topographic exposure term

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