Fuel Load Mapping from LiDAR: Canopy Fuel Metrics in Python

A fuels planner staring at a 15,000-hectare mixed-conifer district needs one thing before any treatment dollar is committed: a spatially explicit map of how much canopy fuel there is and how it is arranged vertically. Field crews can measure a handful of plots, but only airborne LiDAR resolves the vertical fuel structure across the whole landscape — and it is that vertical arrangement, not just the total load, that decides whether a surface fire stays on the ground or climbs into the crowns. This guide derives the four canopy fuel metrics that crown-fire and spread models consume — canopy fuel load, canopy bulk density (CBD), canopy base height (CBH), and canopy cover — from a normalized LiDAR point cloud, and maps them onto a standard fuel model such as one of the Scott & Burgan 40. It is the first stage of the Fire Risk & Fuel Assessment workflow, and it feeds directly into the wildfire risk modeling stage that ranks the landscape by hazard.

Prerequisites checklist

Confirm each item before running the pipeline. Most fuel-metric defects trace back to a skipped check here, and a bad canopy bulk density propagates into every downstream crown-fire and rate-of-spread estimate.

• Python 3.10+ with laspy >= 2.5 (plus a lazrs/laszip backend), numpy >= 1.24, rasterio >= 1.3, and PDAL >= 2.6 available
• A LiDAR tile that has already passed LiDAR Point Cloud Preprocessing — outliers removed and ground returns classified
• Heights normalized to height-above-ground, so a return’s z is metres above the terrain, not an absolute elevation (the same normalization a Canopy Height Model Creation step validates)
• A single projected CRS in metres, so a grid cell has a real area and bulk density comes out in kg/m³
• A stand-level or plot-level canopy fuel load reference (from allometric equations on a tree list, or an existing biomass raster) to anchor the LiDAR-distributed profile
• A target grid resolution chosen for the fuel product — 20 m or 30 m cells are typical, matching operational fire-modeling grids

Concept background: from a vertical profile to canopy fuel metrics

Canopy fuel is not a single number; it is a vertical distribution. The physically meaningful quantity for crown fire is canopy bulk density, the mass of available canopy fuel per unit volume of the canopy, defined for a layer as

$$\text{CBD}(z)=\frac{m_{\text{foliage}}(z)}{A\cdot \Delta z}$$

where $m_{\text{foliage}}(z)$ is the available canopy fuel mass in the layer centred at height $z$, $A$ is the horizontal cell area, and $\Delta z$ is the layer thickness (conventionally 1 m). Evaluating this across every layer produces the vertical bulk-density profile — a curve of kg/m³ against height.

Two of the four fuel metrics are read directly off that profile. Following the Reinhardt and Andersen approach adopted by operational fuel mapping, the stand’s single canopy bulk density is the maximum of a 1 m running mean over the profile, which smooths sensor noise and represents the densest continuous canopy layer that carries a crown fire:

$$\text{CBD}=\underset{z}{\max }\left(\frac{1}{w}\sum _i\text{CBD}(z_i)\right)$$

with the running-mean window $w$ spanning 1 m. Canopy base height is the lowest height at which the bulk-density profile first exceeds a fuel threshold — conventionally $0.011{\text{kg/m}}^3$ — marking where there is enough fuel for a surface fire to transition into the crown. Canopy fuel load is the vertical integral of the profile (total kg/m² of available canopy fuel), and canopy cover is the fraction of first returns that fall in the canopy.

The mass in each layer comes from distributing a plot- or stand-level canopy fuel load across the height bins in proportion to the vegetation returns in each bin. LiDAR does not weigh foliage directly; it tells you where the foliage is, and an allometric reference tells you how much there is in total. The profile-to-metrics reduction is shown below.

Step-by-step Python pipeline

The pipeline runs per grid cell: read normalized returns, isolate the canopy, bin into vertical layers, convert the return profile into a bulk-density profile using an allometric anchor, then reduce to metrics. The steps below process a single cell’s returns; wrapping them in a windowed raster loop scales the same logic to a tile.

Step 1 — Load normalized returns and confirm the height convention

Read the LAS/LAZ tile with laspy and confirm that z is height above ground, not elevation. A cloud that was never normalized will have z values in the hundreds or thousands of metres, which quietly ruins every layer assignment.

import numpy as np
import laspy

las = laspy.read("normalized_tile.laz")

x = np.asarray(las.x)
y = np.asarray(las.y)
z = np.asarray(las.z)                      # height above ground, metres
return_number = np.asarray(las.return_number)
classification = np.asarray(las.classification)

# Sanity check: a normalized cloud sits near zero at the ground.
if np.nanpercentile(z, 1) < -2.0 or np.nanpercentile(z, 99) > 120.0:
    raise ValueError(
        "z does not look height-normalized. Run ground classification and "
        "height normalization before fuel mapping."
    )

Classified noise (classification 7/18) and ground (classification 2) returns are excluded from the canopy so they cannot inflate the understory layers. The upstream classification is exactly what LiDAR Point Cloud Preprocessing delivers.

Step 2 — Isolate the canopy and bin returns into 1 m layers

Drop ground and near-ground returns below a live-crown cutoff, then assign every remaining return to a 1 m height layer. The cutoff (commonly 1–2 m) removes low shrubs and surface fuel that belong to a separate surface-fuel model, not the canopy.

HEIGHT_CUTOFF = 1.5        # metres; below this is surface fuel, not canopy
BIN_HEIGHT = 1.0           # metres per vertical layer
MAX_HEIGHT = 60.0          # clip absurd returns

veg = (classification != 2) & (classification != 7) & (classification != 18)
canopy = veg & (z >= HEIGHT_CUTOFF) & (z <= MAX_HEIGHT)

z_canopy = z[canopy]
edges = np.arange(0.0, MAX_HEIGHT + BIN_HEIGHT, BIN_HEIGHT)
layer_counts, _ = np.histogram(z_canopy, bins=edges)
layer_centres = edges[:-1] + BIN_HEIGHT / 2.0

layer_counts is now the vertical return profile: how many canopy returns fall in each 1 m layer. The shape of this profile — not its raw magnitude — is what carries the vertical fuel arrangement.

Step 3 — Convert the return profile into a foliage-biomass profile

LiDAR return counts are proportional to intercepted foliage, but they are not mass. Anchor them to a real canopy fuel load: distribute a stand-level available canopy fuel load W (kg/m², from allometric equations on the plot’s tree list) across the layers in proportion to their return fraction.

CELL_AREA = 20.0 * 20.0    # m^2, the fuel-grid cell footprint
W_TOTAL = 0.85             # kg/m^2 available canopy fuel load (allometric anchor)

total_returns = layer_counts.sum()
if total_returns == 0:
    layer_fraction = np.zeros_like(layer_counts, dtype=float)
else:
    layer_fraction = layer_counts / total_returns

# Mass per layer (kg) over the whole cell, then per unit area (kg/m^2).
layer_load = layer_fraction * W_TOTAL          # kg/m^2 in each layer
layer_mass = layer_load * CELL_AREA            # kg in each layer over the cell

The allometric anchor W_TOTAL is where field data enters the pipeline. In production it is a per-cell value looked up from a biomass raster or predicted from LiDAR-derived height and cover with a fitted allometric model, not a constant — the constant here keeps the example readable.

Step 4 — Build the vertical bulk-density profile

Divide each layer’s load by the layer thickness to convert kg/m² into kg/m³, the units CBD is defined in.

bulk_density = layer_load / BIN_HEIGHT         # kg/m^3 per 1 m layer

bulk_density is the vertical bulk-density profile plotted in the diagram above: near-zero in the understory, peaking in the mid-to-upper canopy, tapering at the treetops.

Step 5 — Derive CBD (running mean) and CBH (threshold crossing)

Smooth the profile with a 1 m running mean, take its maximum for CBD, and find the lowest layer whose bulk density exceeds the fuel threshold for CBH. This is the reduction detailed further in Estimating Canopy Bulk Density from LiDAR.

CBD_THRESHOLD = 0.011      # kg/m^3, canopy base height fuel threshold
RUN_WINDOW = 1             # layers in the running mean (1 m at 1 m bins)

kernel = np.ones(RUN_WINDOW) / RUN_WINDOW
smoothed = np.convolve(bulk_density, kernel, mode="same")

cbd = float(smoothed.max())                    # kg/m^3

above = np.where(smoothed >= CBD_THRESHOLD)[0]
cbh = float(layer_centres[above[0]]) if above.size else np.nan   # metres

canopy_fuel_load = float(layer_load.sum())     # kg/m^2, profile integral

cbd, cbh, and canopy_fuel_load are three of the four metrics. The fourth, canopy cover, comes from the proportion of first returns in the canopy rather than the biomass profile:

first = return_number == 1
cover = float((z[first] >= HEIGHT_CUTOFF).sum() / max(first.sum(), 1))

Step 6 — Map metrics to a standard fuel model

The continuous metrics are assigned to a discrete fuel model — one of the Scott & Burgan 40 — through a classification lookup on vegetation type, cover, and load. The mapping is a rule table, applied per cell after the metrics are computed:

def assign_fuel_model(cover: float, cbd_val: float, veg_type: str) -> str:
    """Assign a Scott & Burgan 40 fuel model code from canopy metrics."""
    if veg_type == "timber":
        if cover >= 0.60 and cbd_val >= 0.10:
            return "TU5"    # very high load, dry-climate timber-understory
        if cover >= 0.40:
            return "TU1"    # low-load timber-grass-shrub
        return "TL1"        # low-load compact conifer litter
    if veg_type == "shrub":
        return "SH5" if cover >= 0.50 else "SH1"
    return "GR2"            # low-load dry-climate grass fallback

Writing the per-cell metrics and the fuel-model code back to co-registered rasters with rasterio produces the fuel product the Fire Risk & Fuel Assessment pipeline consumes downstream.

PDAL configuration for the return-conditioning stage

When the tile still needs conditioning, run it through a PDAL pipeline before the laspy metric step. The annotated JSON below classifies ground with the Simple Morphological Filter and normalizes height, so the output cloud satisfies the Step 1 precondition:

{
  "pipeline": [
    "raw_tile.laz",
    {
      "type": "filters.outlier",
      "method": "statistical",
      "mean_k": 8,
      "multiplier": 2.5
    },
    {
      "type": "filters.smrf",
      "scalar": 1.2,
      "slope": 0.2,
      "threshold": 0.45,
      "window": 18.0
    },
    {
      "type": "filters.hag_nn"
    },
    {
      "type": "filters.ferry",
      "dimensions": "HeightAboveGround=>Z"
    },
    "normalized_tile.laz"
  ]
}

The key parameters: filters.outlier removes the high-altitude noise returns that would otherwise dominate an upper layer; filters.smrf classifies ground; filters.hag_nn computes height-above-ground from the nearest ground neighbours; and filters.ferry moves that height into Z so the downstream laspy code reads normalized heights directly.

Validation and verification

Fuel metrics derived from LiDAR must be checked against field plots and against physical plausibility before they drive a treatment decision.

def validate_fuel_metrics(cbd, cbh, load, cover):
    """Range and consistency checks on derived canopy fuel metrics."""
    assert 0.0 <= cover <= 1.0, f"cover out of range: {cover}"
    assert 0.0 <= cbd <= 0.5, f"CBD implausible: {cbd} kg/m^3"
    assert load >= 0.0, f"negative fuel load: {load}"
    if not np.isnan(cbh):
        assert 0.0 <= cbh <= 40.0, f"CBH implausible: {cbh} m"
        # A stand with real crown fuel should have CBD above the threshold.
        assert cbd >= 0.011, "CBH defined but CBD below threshold — inconsistent"
    return True

Beyond range checks, regress LiDAR-derived CBD and CBH against destructively sampled or plot-cruised values on a set of calibration plots; a defensible product reports the residual standard error and the bias, not just a map. Where field CBD is unavailable, cross-check the CBD-to-load ratio against published values for the forest type.

Failure modes and gotchas

• Un-normalized heights. If z is elevation rather than height above ground, every layer assignment is wrong and CBD is meaningless. The Step 1 percentile guard catches the common case; the real fix is upstream height normalization.
• Noise returns inflating the top layers. A single un-filtered bird or cloud return at 200 m creates a phantom upper canopy layer. Run filters.outlier first and clip with MAX_HEIGHT.
• Surface fuel leaking into the canopy. Too low a HEIGHT_CUTOFF folds shrubs and slash into the canopy profile, raising CBD and lowering CBH. Match the cutoff to the live-crown base of the forest type.
• Constant allometric anchor. Distributing one W_TOTAL across a heterogeneous tile flattens real between-stand variation. Predict per-cell load from LiDAR height and cover with a fitted model instead.
• Low point density. Below roughly 4 returns/m² the vertical profile is too sparse for a stable running mean, and CBD swings with a single return. Report point density with the product and flag under-sampled cells.
• Wrong CRS units. A geographic CRS makes CELL_AREA and the kg/m³ conversion nonsense. Enforce a projected CRS in metres, as the coordinate-system workflow requires.

Performance and scale notes

• Window the raster, not the whole survey. Process the tile in rasterio block windows so a district-wide survey never has to fit in RAM; compute metrics per cell within each window and stream the outputs.
• Bin with np.histogram, not Python loops. Layer assignment for millions of returns is a single vectorized histogram call; a per-return loop is orders of magnitude slower.
• Parallelize by tile. Fuel mapping is embarrassingly parallel across tiles — dispatch tiles to a process pool, since each is independent once normalized.
• Cache the allometric raster. If per-cell W_TOTAL comes from a fitted model, evaluate it once per tile and reuse it across the metric steps rather than recomputing per cell.

Frequently Asked Questions

Related

• Fire Risk & Fuel Assessment — the parent workflow this fuel-mapping stage feeds
• Estimating Canopy Bulk Density from LiDAR — the running-mean CBD sub-task in depth
• Wildfire Risk Modeling in Python — the hazard stage that consumes these fuel metrics
• LiDAR Point Cloud Preprocessing — the outlier removal and ground classification this guide assumes
• Canopy Height Model Creation — the height normalization that makes vertical binning valid