Post-Fire Vegetation Recovery Monitoring in Python

Burn severity tells you what the fire did; recovery monitoring tells you what happened next, and it is the part that determines management. A stand that regains 80% of its pre-fire greenness within three years needs nothing; one that is still bare after five is a reforestation liability, a reburn risk, and — where it is a legal obligation — a compliance problem. The concrete problem this page solves: you have a fire perimeter, a pre-fire baseline, and a decade of Sentinel-2 or Landsat imagery, and you need a per-pixel recovery trajectory plus a defensible recovery class for every stand. It belongs to Fire Risk & Fuel Assessment and picks up where Burn Severity Mapping with Python leaves off.

Recovery analysis is fundamentally a time-series problem rather than a two-date one, and that changes what can go wrong. A single cloudy summer, a change in the satellite’s processing baseline, or a drought year all look like recovery signals. The defences are a reference-based normalisation, a dense enough series that a median absorbs contamination, and a trajectory model whose parameters have physical meaning.

Prerequisites

• A fire perimeter with a known ignition or containment date, in a projected metre-based CRS.
• A severity raster from Calculating dNBR from Sentinel-2 in Python, so recovery can be stratified by how hard each pixel burned.
• At least five post-fire growing seasons of imagery, plus three pre-fire seasons for the baseline.
• Python ≥ 3.10 with rasterio, numpy, scipy, xarray, geopandas, and pandas.
• A cloud mask applied per scene — Masking Clouds with the Sentinel-2 Scene Classification Layer is the cheap defensible option.
• A forest type or stand layer, so reference pixels are matched to what actually burned.
• One consistent seasonal window — a fixed anniversary window, typically peak growing season, for every year.

The anniversary window is the item most often skipped and the one that most often ruins a series. Comparing a June composite in one year with an August composite in the next measures phenology, not recovery.

Concept: Which Indicator, and Why Normalise

Three indicators dominate, and they answer different questions.

NBR — the normalised burn ratio — is built from near-infrared and shortwave-infrared reflectance and is the standard for both severity and recovery, because SWIR responds strongly to the loss and return of canopy moisture:

$$\mathrm{NBR}=\frac{{\rho }_{\mathrm{NIR}}-{\rho }_{SWIR2}}{{\rho }_{\mathrm{NIR}}+{\rho }_{SWIR2}}$$

NDVI tracks green biomass and recovers faster than NBR, because grasses and forbs green up long before canopy structure returns. That difference is informative rather than a nuisance: a stand where NDVI has recovered and NBR has not is a grass-dominated scar, which is a different management case from a regenerating conifer stand.

Relative recovery turns either into a proportion of what was lost:

$${\mathrm{RRI}}_t=\frac{X_t-X_{\text{post}}}{X_{\text{pre}}-X_{\text{post}}}$$

where $X_{\text{pre}}$ is the pre-fire baseline, $X_{\text{post}}$ the first post-fire measurement, and $X_t$ the value in year $t$. An RRI of 1.0 means the indicator has returned to its pre-fire level; values above 1.0 are common in productive sites where fast-growing regeneration overshoots the mature canopy’s signal.

The normalisation step is separate and equally important. Climate variability moves every pixel in a region together, so subtracting the concurrent value of unburned reference pixels of the same forest type removes the shared signal and leaves the recovery.

Step-by-Step Python Pipeline

Step 1 — Build one composite per year in a fixed window

import numpy as np
import xarray as xr


def annual_composites(stack: xr.DataArray, month_start: int = 7,
                      month_end: int = 8) -> xr.DataArray:
    """Median composite per year over a fixed anniversary window."""
    sel = stack.sel(time=stack["time"].dt.month.isin(range(month_start, month_end + 1)))
    out = sel.groupby("time.year").median(dim="time", skipna=True)
    # Carry the clear-observation count so thin years can be flagged downstream.
    n = sel.notnull().groupby("time.year").sum(dim="time")
    return out.assign_coords(n_clear=("year", n.min(dim=("y", "x")).values))

Keeping the observation count as a coordinate is what lets you distinguish a genuine recovery plateau from a year where two cloudy scenes produced a noisy median.

Step 2 — Establish the pre-fire baseline and the immediate post-fire floor

def pre_post_levels(annual: xr.DataArray, fire_year: int, n_pre: int = 3):
    """Pre-fire baseline (median of n_pre years) and the first post-fire value."""
    pre_years = [fire_year - k for k in range(1, n_pre + 1)]
    pre = annual.sel(year=[y for y in pre_years if y in annual["year"]]).median(dim="year")
    post = annual.sel(year=fire_year + 1)
    return pre, post

Use the year after the fire as the floor, not the fire year itself: a mid-season fire leaves the fire-year composite as a mixture of pre- and post-fire observations, which biases the floor upward and every subsequent recovery ratio downward.

Step 3 — Normalise against unburned reference pixels

import numpy as np


def reference_normalise(annual: xr.DataArray, burned_mask: np.ndarray,
                        forest_type: np.ndarray, buffer_mask: np.ndarray) -> xr.DataArray:
    """Subtract the concurrent median of unburned pixels of the same forest type."""
    out = annual.copy()
    for t in np.unique(forest_type[burned_mask]):
        ref_sel = (~burned_mask) & buffer_mask & (forest_type == t)
        if ref_sel.sum() < 100:
            continue                       # too few reference pixels to trust
        ref_series = annual.where(xr.DataArray(ref_sel, dims=("y", "x"))).median(
            dim=("y", "x"), skipna=True)
        target = xr.DataArray(burned_mask & (forest_type == t), dims=("y", "x"))
        out = xr.where(target, annual - ref_series, out)
    return out

The buffer_mask restricts reference pixels to a ring around the fire — typically 1–5 km — so they share the same weather, soils and elevation band as the burned area. Reference pixels drawn from a hundred kilometres away introduce a climate gradient that looks exactly like differential recovery.

Step 4 — Fit a recovery trajectory per pixel

A saturating exponential is the standard form: it has two parameters, both of which mean something, and it cannot predict unbounded growth.

$$X_t=A\left(1-e^{-kt}\right)+X_{\text{post}}$$

$A$ is the asymptote — how much of the lost signal will eventually return — and $k$ is the rate, whose reciprocal is the recovery time constant in years.

import numpy as np
from scipy.optimize import curve_fit


def _recovery_model(t, A, k, floor):
    return A * (1.0 - np.exp(-k * t)) + floor


def fit_recovery(years_since: np.ndarray, values: np.ndarray, floor: float) -> dict:
    """Fit the saturating recovery curve for one pixel's time series."""
    ok = np.isfinite(values)
    if ok.sum() < 4:
        return {"A": np.nan, "k": np.nan, "r2": np.nan, "n": int(ok.sum())}
    t, y = years_since[ok], values[ok]
    try:
        (A, k, _), _ = curve_fit(
            lambda tt, AA, kk: _recovery_model(tt, AA, kk, floor),
            t, y, p0=[max(y.max() - floor, 1e-3), 0.3],
            bounds=([0.0, 0.01], [np.inf, 3.0]), maxfev=2000)
    except (RuntimeError, ValueError):
        return {"A": np.nan, "k": np.nan, "r2": np.nan, "n": int(ok.sum())}
    pred = _recovery_model(t, A, k, floor)
    ss_res = float(((y - pred) ** 2).sum())
    ss_tot = float(((y - y.mean()) ** 2).sum())
    return {"A": float(A), "k": float(k),
            "half_life_yr": float(np.log(2) / k),
            "r2": 1.0 - ss_res / ss_tot if ss_tot else np.nan,
            "n": int(ok.sum())}

The bounds are doing real work. Without a lower bound on $k$ the optimiser will happily return a rate near zero with an enormous asymptote, which fits a straight line through five noisy points and predicts recovery in three centuries.

Step 5 — Classify and summarise per stand

import numpy as np


def classify_recovery(rri_year5: np.ndarray, k: np.ndarray) -> np.ndarray:
    """Recovery classes from year-5 relative recovery and the fitted rate."""
    cls = np.zeros(rri_year5.shape, dtype="uint8")     # 0 = no data
    cls[np.isfinite(rri_year5)] = 1                    # 1 = not recovering
    cls[(rri_year5 >= 0.25) & (k >= 0.10)] = 2         # 2 = slow
    cls[(rri_year5 >= 0.50) & (k >= 0.20)] = 3         # 3 = moderate
    cls[(rri_year5 >= 0.75) & (k >= 0.35)] = 4         # 4 = rapid
    return cls

Using both the level and the rate is deliberate. A pixel at 60% of pre-fire with a high rate is on a healthy trajectory; the same 60% with a rate near zero has stalled, and only the second needs intervention.

Validation & Verification

Recovery maps are hard to validate directly because nobody measures “percent recovered” in the field. Three checks make the product defensible anyway:

import numpy as np


def recovery_diagnostics(rri, severity, k, r2) -> dict:
    """Relationships that must hold if the trajectory fits mean anything."""
    ok = np.isfinite(rri) & np.isfinite(severity) & np.isfinite(k)
    # High-severity pixels must recover more slowly than low-severity ones.
    hi = ok & (severity >= 0.66)
    lo = ok & (severity < 0.27)
    return {
        "median_rri_high_sev": float(np.median(rri[hi])),
        "median_rri_low_sev": float(np.median(rri[lo])),
        "median_k_high_sev": float(np.median(k[hi])),
        "median_k_low_sev": float(np.median(k[lo])),
        "median_r2": float(np.nanmedian(r2)),
        "unfitted_frac": float(np.mean(~np.isfinite(k))),
    }


d = recovery_diagnostics(rri5, dnbr, k_map, r2_map)
assert d["median_rri_high_sev"] < d["median_rri_low_sev"], \
    "high-severity pixels recovering faster than low — check the severity raster sign"
assert d["median_r2"] > 0.5, f"trajectory fits are poor: {d}"
assert d["unfitted_frac"] < 0.2, f"too many pixels failed to fit: {d}"

The severity-versus-recovery relationship is the strongest available internal check. It is a robust ecological expectation — hard-burned ground recovers more slowly — so a violation almost always means a sign error in dNBR or a mismatched fire year rather than a surprising forest.

Where field data exists, validate against regeneration survey plots: a Spearman correlation between fitted recovery rate and measured seedling density above about 0.5 is a good result, and the residuals are usually informative about site preparation and planting history rather than about the remote sensing.

Failure Modes & Gotchas

• A fire-year composite used as the post-fire floor. Mid-season fires make that composite a blend of burned and unburned observations, biasing every ratio.
• Reference pixels drawn from a different forest type. A conifer scar normalised against broadleaf reference recovers “too fast” because the reference greens up on a different schedule.
• Phenology mistaken for recovery. A shifting seasonal window produces a trend that is entirely calendar. Fix the anniversary window and check the composite dates.
• Sensor changes mid-series. Landsat 7 to 8, or a Sentinel-2 processing-baseline change, introduces a step. Harmonise with published coefficients or split the series at the boundary.
• Fitting on fewer than four clear years. The two-parameter model is unidentifiable, and the optimiser returns confident nonsense rather than failing.
• Ignoring reburn. A second fire inside the window is a discontinuity, not a recovery reversal. Screen for it explicitly rather than letting the fit absorb it.

Performance & Scale Notes

Per-pixel curve fitting is the expensive step: an 8,000 ha fire at 20 m is roughly 200,000 pixels, and curve_fit costs milliseconds each. Two changes make it tractable. First, fit only where it matters — mask to the burn perimeter and skip pixels with fewer than four clear years. Second, vectorise the easy cases: relative recovery at a fixed year needs no optimiser, so compute the RRI map for everything and fit trajectories only where the classification is ambiguous.

import concurrent.futures

import numpy as np


def fit_recovery_grid(cube: np.ndarray, years_since: np.ndarray,
                      floor: np.ndarray, mask: np.ndarray, workers: int = 6):
    """Fit trajectories over a masked grid. cube is (time, y, x)."""
    rows, cols = np.where(mask)
    def one(i):
        r, c = rows[i], cols[i]
        return fit_recovery(years_since, cube[:, r, c], float(floor[r, c]))
    with concurrent.futures.ThreadPoolExecutor(max_workers=workers) as pool:
        results = list(pool.map(one, range(len(rows))))
    A = np.full(mask.shape, np.nan, dtype="float32")
    k = np.full(mask.shape, np.nan, dtype="float32")
    for i, res in enumerate(results):
        A[rows[i], cols[i]] = res["A"]
        k[rows[i], cols[i]] = res["k"]
    return A, k

For a regional multi-fire analysis, aggregate first: fitting one trajectory per stand, or per severity class within each stand, is usually more useful than a per-pixel map and is two orders of magnitude cheaper.

Frequently Asked Questions

Related

• Burn Severity Mapping with Python — produces the severity raster recovery is stratified by
• Measuring Post-Fire NDVI Recovery Trajectories in Python — the trajectory fit in isolation
• Classifying Recovery Classes from Multi-Year dNBR — turning trajectories into a management map
• Detecting Reburn Risk from Recovery and Fuel Rasters — where a stalled recovery becomes a hazard
• Wildfire Risk Modeling in Python — the risk surface recovery status feeds back into

Up: Fire Risk & Fuel Assessment