Wildfire Risk Modeling in Python

A district fire planner asking “where should the fuel-treatment budget go?” is really asking a spatial question: which stands combine dry, abundant fuel, steep aspect-driven exposure, and a realistic chance of ignition. Answering it means fusing several raster layers — fuel, terrain, weather, and ignition history — into a single, defensible risk surface. This guide builds that surface in Python, and it is one workflow inside the Fire Risk & Fuel Assessment section, drawing terrain inputs from the Canopy Height Modeling & Terrain Extraction work upstream. Along the way we separate three words that are used loosely in the field but mean different things numerically — hazard, risk, and susceptibility — because conflating them is the fastest route to a map that misdirects money.

The output is a georeferenced GeoTIFF where each pixel carries a 0–1 index, aligned to a common grid so it can be tiled, differenced between years, and overlaid on ownership polygons. We cover both the transparent weighted-overlay approach that agencies can audit, and a data-driven susceptibility model (logistic regression or random forest) trained on historical ignitions, and we say plainly when each is appropriate.

Prerequisites checklist

Every layer below must resolve to the same CRS, resolution, and pixel grid before any arithmetic. Misalignment is the dominant source of silent error in raster risk work — a half-pixel shift between the fuel and slope layers biases every stand near a boundary.

• Python 3.10+ with rasterio >= 1.3, numpy >= 1.24, scikit-learn >= 1.3, and geopandas >= 0.14
• A fuel layer — fuel model class, fuel load (t/ha), or canopy bulk density from the fuel load mapping workflow
• A digital terrain model at the analysis resolution, from Digital Terrain Model Generation, to derive slope and aspect
• A fire-weather summary raster (or a single design-day value) — see computing the Fire Weather Index
• Historical ignition points (agency fire-occurrence database) for the data-driven path, in a known CRS
• One project-wide equal-area or UTM CRS in metres, chosen before resampling anything
• Enough disk for intermediate aligned rasters; do not hold a national mosaic in RAM

Concept background: hazard vs risk vs susceptibility

These three terms are not synonyms, and the formula you write depends on which one you actually mean.

• Susceptibility (or likelihood) is the probability that a pixel ignites and/or burns given its intrinsic conditions — fuel, terrain, and long-run climate. It is dimensionless, bounded 0–1, and is what a statistical model estimates from historical fire locations. It deliberately excludes the value of what is exposed.
• Hazard describes the potential fire behavior at a pixel — intensity, rate of spread, flame length — if a fire were present. It is a physical prediction (see fire behavior modeling) driven by fuel, weather, and slope, and it says nothing about how often fire occurs.
• Risk, in the wildfire-management sense, couples likelihood with consequence: $\text{Risk}=P(\text{fire})\times \text{Intensity}\times \text{Value-at-risk}$. A remote scree slope can be high-hazard yet low-risk because nothing of value is exposed.

For a first, auditable district map, a weighted linear overlay of normalized criteria is the workhorse. Each input layer $r_i$ is rescaled to $[0,1]$, assigned a weight $w_i$ (weights summing to one), and combined:

$$R=\sum _{i=1}^n{w}_i{r}_i,\sum _{i=1}^n{w}_i=1$$

The transparency is the point: a reviewer can see that fuel got weight 0.4 and ignition history 0.15. Its weakness is that the weights are expert judgement, not fitted to observed fire. The data-driven alternative replaces those weights with coefficients learned from where fires actually occurred. A logistic-regression susceptibility model estimates

$$P(\text{fire}\mid \mathbf{x})=\frac{1}{1+e^{-({\beta }_0+{\sum }_i{\beta }_i{x}_i)}}$$

where $\mathbf{x}$ are the same predictor layers, now fitted against presence/absence of historical ignition. A random forest relaxes the linear-and-additive assumption and captures interactions (steep and south-facing and dry) that a weighted sum cannot.

A practical consequence of this distinction: the inputs you assemble are shared, but the combination rule is not. Weighted overlay asks you to defend a set of weights; a susceptibility model asks you to defend a training dataset and a validation scheme. Many agencies run both — the overlay as a communicable baseline, the fitted model as a check on whether the expert weights actually match where fire has occurred. When the two disagree sharply over a stand, that disagreement is itself informative and worth investigating rather than averaging away.

Step-by-step Python pipeline

The pipeline has five ordered stages: align all layers to a reference grid, derive slope and aspect from the terrain model, normalize each criterion, combine by weighted overlay, and (optionally) fit a data-driven susceptibility model. Each step is runnable.

Step 1 — Align every layer to one reference grid

Pick the layer that defines your target resolution — usually the fuel raster — and reproject everything else onto its exact transform. This guarantees that array_a[i, j] and array_b[i, j] describe the same ground.

import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling


def align_to_reference(src_path: str, ref_profile: dict, categorical: bool = False) -> np.ndarray:
    """Resample a raster onto the reference grid, returning a float32 array."""
    resampling = Resampling.nearest if categorical else Resampling.bilinear
    dst = np.empty((ref_profile["height"], ref_profile["width"]), dtype="float32")
    with rasterio.open(src_path) as src:
        reproject(
            source=rasterio.band(src, 1),
            destination=dst,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=ref_profile["transform"],
            dst_crs=ref_profile["crs"],
            resampling=resampling,
            dst_nodata=np.nan,
        )
    return dst


with rasterio.open("fuel_load.tif") as ref:
    ref_profile = ref.profile
    fuel = ref.read(1).astype("float32")

dtm = align_to_reference("dtm.tif", ref_profile)
fwi = align_to_reference("fire_weather_fwi.tif", ref_profile)

Use Resampling.nearest for categorical inputs (a fuel-model class raster) and Resampling.bilinear for continuous surfaces (elevation, FWI) — averaging class codes invents fuel models that do not exist. This is the same discipline the Raster-Vector Overlay Techniques workflow enforces before any zonal step.

Step 2 — Derive slope and aspect from the terrain model

Fire spreads faster upslope, and aspect controls solar drying, so both are first-class risk criteria. Compute them from the aligned DTM with numpy gradients; pixel spacing comes from the reference transform.

def slope_aspect(dtm: np.ndarray, pixel_size: float):
    """Return slope (degrees) and aspect (degrees from north) via a 3x3 gradient."""
    dz_dy, dz_dx = np.gradient(dtm, pixel_size)          # rise over run per axis
    slope_rad = np.arctan(np.hypot(dz_dx, dz_dy))
    slope_deg = np.degrees(slope_rad)
    aspect_rad = np.arctan2(dz_dy, -dz_dx)
    aspect_deg = (np.degrees(aspect_rad) + 360.0) % 360.0
    return slope_deg, aspect_deg


pixel_size = ref_profile["transform"].a                  # metres per pixel (x)
slope_deg, aspect_deg = slope_aspect(dtm, pixel_size)

# South/southwest aspects dry fastest in the northern hemisphere: score exposure
# as the cosine distance from 202.5 degrees (SSW), mapped to 0-1.
aspect_exposure = (1.0 + np.cos(np.radians(aspect_deg - 202.5))) / 2.0

For the physics behind slope-driven acceleration — the ${\varphi }_s$ term — see the rate-of-spread modeling guide; here slope enters only as a normalized exposure criterion.

Step 3 — Normalize each criterion to a common scale

A weighted sum only makes sense if the inputs share a scale. Rescale each to $[0,1]$ with a robust min–max on percentiles, so a single extreme pixel does not compress everything else.

def normalize(a: np.ndarray, invert: bool = False, p_low=2, p_high=98) -> np.ndarray:
    """Percentile min-max scale to 0-1, ignoring NaNs. invert=True flips the sense."""
    lo, hi = np.nanpercentile(a, [p_low, p_high])
    if hi == lo:
        return np.zeros_like(a, dtype="float32")
    scaled = np.clip((a - lo) / (hi - lo), 0.0, 1.0)
    return (1.0 - scaled) if invert else scaled


n_fuel = normalize(fuel)                 # more fuel -> higher risk
n_slope = normalize(slope_deg)           # steeper -> higher risk
n_aspect = aspect_exposure               # already 0-1
n_weather = normalize(fwi)               # drier/windier -> higher risk

Step 4 — Combine by weighted overlay

Declare weights explicitly and assert they sum to one — this single line is what makes the map auditable. Keep the weight dictionary in version control alongside the output.

weights = {
    "fuel": 0.35,
    "slope": 0.20,
    "aspect": 0.10,
    "weather": 0.25,
    "ignition": 0.10,
}
assert abs(sum(weights.values()) - 1.0) < 1e-9, "weights must sum to 1"

n_ignition = normalize(align_to_reference("ignition_density.tif", ref_profile))

risk = (
    weights["fuel"] * n_fuel
    + weights["slope"] * n_slope
    + weights["aspect"] * n_aspect
    + weights["weather"] * n_weather
    + weights["ignition"] * n_ignition
).astype("float32")

out_profile = ref_profile.copy()
out_profile.update(dtype="float32", count=1, nodata=np.nan)
with rasterio.open("wildfire_risk_index.tif", "w", **out_profile) as dst:
    dst.write(risk, 1)

Step 5 — Fit a data-driven susceptibility model (optional)

When you have a credible historical ignition dataset, let the data set the weights. Sample predictor values at ignition points (presence) and at random background pixels (absence), then fit a classifier. This is the same presence/background logic used in species distribution modeling, applied to fire.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

# Stack aligned predictors into a (pixels, features) design matrix.
predictors = np.dstack([n_fuel, n_slope, n_aspect, n_weather]).reshape(-1, 4)

# y: 1 at ignition pixels, 0 at background pixels (rows indexed as above).
# ignition_rows / background_rows are flat pixel indices you sampled earlier.
X = np.vstack([predictors[ignition_rows], predictors[background_rows]])
y = np.concatenate([np.ones(len(ignition_rows)), np.zeros(len(background_rows))])
valid = ~np.isnan(X).any(axis=1)
X, y = X[valid], y[valid]

logit = LogisticRegression(max_iter=1000, class_weight="balanced")
print("logit AUC:", cross_val_score(logit, X, y, cv=5, scoring="roc_auc").mean())

rf = RandomForestClassifier(n_estimators=400, min_samples_leaf=5,
                            class_weight="balanced", random_state=42)
rf.fit(X, y)

# Predict susceptibility for every valid pixel, reshape back to the grid.
full = predictors.copy()
mask = ~np.isnan(full).any(axis=1)
susceptibility = np.full(full.shape[0], np.nan, dtype="float32")
susceptibility[mask] = rf.predict_proba(full[mask])[:, 1]
susceptibility = susceptibility.reshape(fuel.shape)

Configuration: the weight and predictor manifest

Keep model configuration in a small, version-controlled file rather than scattered through code. A JSON manifest makes the map reproducible and lets a reviewer see the whole specification at a glance.

{
  "target_crs": "EPSG:5070",
  "resolution_m": 30,
  "reference_layer": "fuel_load.tif",
  "method": "weighted_overlay",
  "weights": {
    "fuel": 0.35,
    "slope": 0.20,
    "aspect": 0.10,
    "weather": 0.25,
    "ignition": 0.10
  },
  "normalization": { "type": "percentile_minmax", "p_low": 2, "p_high": 98 },
  "resampling": { "continuous": "bilinear", "categorical": "nearest" }
}

Validation and verification

Never ship a risk map without checking that it is finite, aligned, and calibrated. Two cheap tests catch most defects: a structural check on the array, and — for the data-driven path — an AUC against held-out fires.

# Structural: no all-NaN output, values in range, grid matches the reference.
assert np.isfinite(risk[~np.isnan(risk)]).all()
assert np.nanmin(risk) >= 0.0 and np.nanmax(risk) <= 1.0
with rasterio.open("wildfire_risk_index.tif") as chk:
    assert chk.transform == ref_profile["transform"]
    assert chk.crs == ref_profile["crs"]

# Skill: an ROC-AUC below ~0.7 on spatially independent folds means the model
# is not separating burned from unburned better than chance-plus-a-little.
auc = cross_val_score(logit, X, y, cv=5, scoring="roc_auc").mean()
assert auc > 0.65, f"weak susceptibility skill (AUC={auc:.2f})"

Cross-validation for spatial models must use spatially separated folds, not random ones — nearby pixels are correlated, and random folds leak information across the train/test split and inflate AUC. Block the study area into tiles and hold out whole tiles.

Beyond a single AUC number, inspect where the model succeeds and fails. Map the residuals — predicted susceptibility minus observed burn — and look for spatial structure. A ring of over-prediction around towns usually signals sampling bias in the ignition record (fires are recorded where people are), while systematic under-prediction on north aspects may mean an important predictor, such as fuel moisture or canopy shading, is missing. A calibration plot (binned predicted probability against observed burn frequency) is the other cheap diagnostic: a well-calibrated susceptibility surface should track the 1:1 line, so that a pixel scored 0.3 burns about 30% as often as a pixel scored 1.0. Weighted-overlay maps cannot be calibrated this way, which is one more reason to fit the data-driven model even when the overlay is the published deliverable.

Failure modes and gotchas

• Unaligned grids. A one-pixel offset between fuel and slope silently biases every boundary stand. Always reproject onto one reference transform (Step 1) and assert the transform on the output.
• Averaging categorical fuel models. Bilinear-resampling a fuel-class raster fabricates non-existent classes. Use Resampling.nearest for any categorical layer.
• Weights that do not sum to one. An un-normalized weight vector rescales the whole index and breaks comparability between map versions. Assert the sum.
• Random cross-validation on spatial data. Spatial autocorrelation leaks between random folds and inflates AUC by 0.05–0.15. Use spatial block folds.
• Presence-only ignition data with sampling bias. Fire databases over-record near roads and towns; unmodelled, the classifier learns detectability, not flammability. Add road distance as a predictor or use target-group background.
• Confusing hazard with risk. A high-hazard, zero-value pixel is not high-risk. If the deliverable drives budget, multiply the likelihood surface by a value-at-risk layer before calling it risk.

Performance and scale notes

• Tile and stream. Process the raster in windows (src.block_windows()) so a state-wide 30 m mosaic never has to fit in RAM; apply the weighted sum window by window.
• Persist aligned intermediates. Warping is the expensive step — reproject each input once to the reference grid, write a Cloud-Optimized GeoTIFF, and reuse it across weight experiments.
• Vectorize, do not loop pixels. All the normalization and the weighted sum are numpy array operations; a Python per-pixel loop is orders of magnitude slower.
• Predict in chunks. For a random forest over tens of millions of pixels, call predict_proba on row blocks to bound memory, and set n_jobs=-1 to use all cores.
• Cache the fitted model. Serialize the trained classifier (joblib.dump) so the annual re-run only re-scores pixels rather than re-fitting.

Frequently Asked Questions

Related

• Fire Risk & Fuel Assessment — the parent section this risk workflow belongs to
• Fuel Load Mapping from LiDAR — produces the fuel criterion that drives the index
• Computing the Fire Weather Index in Python — builds the fire-weather input layer
• Fire Behavior and Rate-of-Spread Modeling — the hazard side: how intense a fire would be
• Digital Terrain Model Generation — the DTM that slope and aspect are derived from
• Raster-Vector Overlay Techniques — the grid-alignment discipline every layer depends on