Forestry GIS Python · Forestry · Ecology

Measuring Post-Fire NDVI Recovery Trajectories in Python

A recovery percentage is a snapshot; a trajectory is a prediction. Fitting a curve to a post-fire NDVI series converts a column of annual numbers into two parameters — how fast the stand is recovering and how far it will get — which is what a silviculturist actually needs in order to decide whether to intervene. This guide fits those curves, as the modelling step inside Post-Fire Vegetation Recovery Monitoring, part of Fire Risk & Fuel Assessment.

When to use each model form

Model Parameters Shape Use when
Saturating exponential rate, asymptote Fast early, flattening The default; matches how regeneration behaves
Linear slope Constant Fewer than five observations, or as a null model
Logistic rate, asymptote, inflection S-curved, slow start Delayed regeneration after site preparation or planting
Piecewise linear breakpoint, two slopes Two phases A known intervention — planting, salvage, a second fire
Non-parametric (Theil–Sen) slope Robust monotone trend Outlier-heavy series where a fit will not converge

Start with the saturating exponential. Move to logistic only when residuals show a systematic early under-prediction, and to piecewise only when you know the date of the event that split the trajectory — fitting a free breakpoint to seven noisy points finds a breakpoint every time, whether or not one exists.

The model and its parameters

is the immediate post-fire floor, held fixed from the data rather than fitted; is the total amount the index will eventually regain; is the rate, and is the half-life — the number of years to close half the remaining gap. Reporting the half-life rather than is worth the extra line, because “seven years to halfway” is a sentence a manager can act on and “k = 0.099” is not.

Recovery rate and what the half-life means Four curves of NDVI against years since fire, all starting from the same post-fire floor. A rapid curve with a rate of 0.5 reaches half of its eventual gain in under one and a half years. A moderate curve at 0.25 takes about three years. A slow curve at 0.1 takes about seven. A fourth, nearly flat curve represents a stalled stand whose fitted rate is near zero. Vertical marks show each curve's half-life on the time axis. 0 4 8 12 16 20 floor pre-fire years since fire NDVI k = 0.50 · half-life 1.4 yr k = 0.25 · half-life 2.8 yr k = 0.10 · half-life 6.9 yr stalled — k ≈ 0 All four reach the same asymptote in principle; only the rate says whether it happens within a rotation.

Minimal reproducible example

import numpy as np
from scipy.optimize import curve_fit


def saturating(t, A, k, floor):
    """NDVI(t) = A(1 - exp(-k t)) + floor."""
    return A * (1.0 - np.exp(-k * t)) + floor


def fit_trajectory(years: np.ndarray, ndvi: np.ndarray, floor: float,
                   min_obs: int = 5) -> dict:
    """Fit a saturating recovery curve to one pixel's or stand's NDVI series."""
    ok = np.isfinite(ndvi)
    if ok.sum() < min_obs:
        return {"status": "too_few_obs", "n": int(ok.sum())}
    t, y = years[ok].astype("float64"), ndvi[ok].astype("float64")

    span = max(float(y.max() - floor), 1e-3)
    try:
        popt, pcov = curve_fit(
            lambda tt, AA, kk: saturating(tt, AA, kk, floor),
            t, y,
            p0=[span, 0.25],
            bounds=([0.0, 0.02], [1.5, 2.0]),      # NDVI gain <= 1.5; rate 0.02-2.0
            maxfev=4000,
        )
    except (RuntimeError, ValueError) as exc:
        return {"status": f"no_fit: {type(exc).__name__}", "n": int(ok.sum())}

    A, k = float(popt[0]), float(popt[1])
    pred = saturating(t, A, k, floor)
    ss_res = float(((y - pred) ** 2).sum())
    ss_tot = float(((y - y.mean()) ** 2).sum())
    se = np.sqrt(np.diag(pcov))
    return {
        "status": "ok",
        "n": int(ok.sum()),
        "asymptote": A,
        "rate": k,
        "half_life_yr": float(np.log(2.0) / k),
        "years_to_80pct": float(-np.log(0.2) / k),
        "rate_se": float(se[1]),
        "r2": 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan"),
        "rmse": float(np.sqrt(ss_res / len(y))),
    }

The bounds are the difference between a usable result and confident nonsense. Without a lower bound on k, a noisy series is fitted by an almost-straight line with a huge asymptote — mathematically a fine local optimum, and a prediction of full recovery in four centuries. Without an upper bound on A, the optimiser will happily propose an NDVI gain of 12.

Returning rate_se matters as much as the rate. A rate of 0.3 with a standard error of 0.04 is a measurement; the same 0.3 with a standard error of 0.9 is noise wearing a number.

Parameter reference

Parameter Type Default Range Rationale
floor float measured — The first full post-fire year; fixing it from data avoids a third free parameter
min_obs int 5 4–8 Two free parameters need at least four points; five leaves one degree of freedom
p0 rate float 0.25 0.1–0.5 A ~3-year half-life is a reasonable start for most temperate systems
bounds on k tuple (0.02, 2.0) — Lower bound rejects the flat-line fit; upper bound rejects instant recovery
bounds on A tuple (0.0, 1.5) — NDVI is bounded, so its recoverable gain is too
maxfev int 4000 1000–10000 Raise it only if convergence failures are common and the bounds are correct

Expected output and verification

Fit quality alone is not enough — a saturating curve fits almost anything over five points. The check that matters is whether the fitted parameters behave the way the ecology requires across many pixels:

import numpy as np


def audit_fits(results: list[dict]) -> dict:
    """Population-level checks over many fitted trajectories."""
    ok = [r for r in results if r.get("status") == "ok"]
    if not ok:
        return {"fitted": 0}
    rates = np.array([r["rate"] for r in ok])
    r2 = np.array([r["r2"] for r in ok])
    se_ratio = np.array([r["rate_se"] / max(r["rate"], 1e-6) for r in ok])
    return {
        "fitted": len(ok),
        "failed": len(results) - len(ok),
        "median_r2": round(float(np.median(r2)), 3),
        "median_half_life": round(float(np.median(np.log(2) / rates)), 2),
        "frac_at_lower_bound": round(float((rates <= 0.021).mean()), 3),
        "frac_unreliable": round(float((se_ratio > 0.5).mean()), 3),
    }


report = audit_fits(all_results)
assert report["frac_at_lower_bound"] < 0.15, \
    f"many fits pinned at the rate floor — the series is too short or too noisy: {report}"
assert report["median_r2"] > 0.6, f"poor fits overall: {report}"

frac_at_lower_bound is the diagnostic worth watching. A fit that lands exactly on the bound has not converged to anything meaningful — the optimiser simply ran out of room — and reporting those pixels as “recovering very slowly” rather than “not fitted” turns a modelling failure into a false ecological claim.

Feed the fitter something it must refuse, as a standing test:

rng = np.random.default_rng(20260809)
years = np.arange(1, 9, dtype="float64")

flat = 0.15 + rng.normal(0, 0.01, size=years.size)          # no recovery at all
res_flat = fit_trajectory(years, flat, floor=0.15)
assert res_flat["status"] != "ok" or res_flat["rate"] <= 0.05, \
    "a flat series was fitted as recovery — the rate bound is not doing its job"

short = np.array([0.2, 0.3, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan])
assert fit_trajectory(years, short, floor=0.15)["status"] == "too_few_obs"
Reading residuals to choose the model form Three residual plots against years since fire. The first, from a saturating fit, shows residuals scattered randomly above and below the zero line with no pattern. The second, from a linear fit, shows a clear arc: residuals are positive early, negative in the middle and positive again late, which is the signature of fitting a straight line to a curve. The third, from a free-breakpoint piecewise fit, has residuals almost exactly on zero, which indicates the model has absorbed the noise rather than the signal. Saturating — well specified Linear — under specified Free breakpoint — over fitted no pattern — accept arc — wrong model form too good — suspect A residual arc means the shape is wrong, not that the data is noisy. Residuals that are all but zero over eight points mean the model has more freedom than the data. Choose the form from the residual pattern, never from R² alone. Three fits to the same eight observations; only the first is reportable. Fitting across a reburn against splitting at it A post-fire index series rises for six years, drops sharply at a second fire, and rises again. A single saturating curve fitted across the whole record is nearly flat and reports a very low recovery rate. Two curves fitted separately, one before the reburn and one after, each recover a plausible rate and asymptote, and the reburn date is marked with a vertical line. reburn one curve — rate 0.03, meaningless rate 0.31 rate 0.26 years since the first fire A reburn is a new event, not a slow recovery — and the single-curve fit hides it behind a plausible number.

Common pitfalls

  • Fitting the floor as a free parameter. With five points and three parameters the fit is barely identified, and the floor drifts to absorb noise. Measure it from the first post-fire year and hold it.
  • Reporting the asymptote from a short series. Before the curve flattens, the asymptote is pure extrapolation. Report the observed level and the rate; add the asymptote once the series is long enough to constrain it.
  • Treating a pinned rate as a slow recovery. A fit sitting on the bound is a failure, and it belongs in the “not fitted” class, not in the slowest recovery class.
  • Comparing rates across indices. NDVI rates are systematically higher than NBR rates because herbaceous cover returns first. A rate is only comparable within one indicator.
  • Fitting through a reburn. A second fire is a discontinuity. Split the series at the event rather than letting one curve average across it.

Frequently Asked Questions

What if the fit will not converge at all?

Fall back to a robust monotone trend — a Theil–Sen slope over the observations — and report it as a slope rather than a rate. It cannot give a half-life, but it is stable under outliers and it will not silently invent a parameter. A large fraction of non-convergent pixels usually means the series is too short rather than the model being wrong.

Should I fit per pixel or per stand?

Per stand for management, per pixel for maps. Aggregating the series to a stand median before fitting gives far more stable parameters because the noise averages out, and it is two orders of magnitude cheaper. Per-pixel fits are worth it when within-stand variation is the question — patchy regeneration around surviving seed trees, for instance.

How do I compare recovery between two fires?

Only through normalised quantities. Rates from different fires are comparable when both series were normalised against local unburned reference pixels, used the same indicator, and were fitted with the same floor convention. Raw NDVI levels are not comparable at all, because pre-fire forest type differs.

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