Estimating Canopy Bulk Density from LiDAR with the Running-Mean Method

Canopy bulk density (CBD) is the single most consequential input to a crown-fire spread model, and estimating it from LiDAR comes down to one narrow algorithm: build a vertical bulk-density profile from a normalized point cloud, smooth it with a 1 m running mean, and take the maximum. This page covers that exact sub-task — the running-mean, vertical-bin estimate that operational fuel mapping adopted from Reinhardt and Andersen. It is a focused piece of the Fuel Load Mapping from LiDAR workflow, which itself is the first stage of the broader Fire Risk & Fuel Assessment pipeline. If you have not yet normalized heights or classified ground, do that first; everything here assumes a cloud whose z is already height above ground.

When to use the running-mean method

CBD can be estimated from a vertical profile in several ways, and the running-mean maximum is the operational default for good reasons. The table compares it with the common alternatives.

Use the running-mean maximum whenever the CBD feeds a crown-fire or rate-of-spread model, which is nearly always. It represents the densest continuous fuel layer — the physically relevant quantity for sustained crown fire — and it is robust to the single-layer noise that makes the raw peak untrustworthy.

Minimal reproducible example

The snippet below takes a normalized tile, isolates the canopy, bins returns into 1 m layers, distributes an allometric canopy fuel load across those layers, converts to a bulk-density profile, and returns the running-mean CBD. It is complete and copy-pasteable; the only external input is a height-normalized LAS/LAZ file.

import numpy as np
import laspy

# --- Parameters (see the reference table below) ---
BIN_HEIGHT = 1.0        # metres per vertical layer
HEIGHT_CUTOFF = 1.5     # metres; below this is surface fuel, excluded
MAX_HEIGHT = 60.0       # metres; clip implausible returns
RUN_WINDOW_M = 1.0      # metres; running-mean window
W_TOTAL = 0.85          # kg/m^2 available canopy fuel load (allometric anchor)


def estimate_cbd(path: str) -> float:
    """Estimate canopy bulk density (kg/m^3) via the running-mean method."""
    las = laspy.read(path)
    z = np.asarray(las.z)                      # height above ground, metres
    classification = np.asarray(las.classification)

    # Canopy returns only: exclude ground (2) and noise (7, 18).
    veg = (classification != 2) & (classification != 7) & (classification != 18)
    canopy = veg & (z >= HEIGHT_CUTOFF) & (z <= MAX_HEIGHT)
    z_canopy = z[canopy]
    if z_canopy.size == 0:
        return 0.0

    # 1 m vertical bins -> return profile.
    edges = np.arange(0.0, MAX_HEIGHT + BIN_HEIGHT, BIN_HEIGHT)
    counts, _ = np.histogram(z_canopy, bins=edges)

    # Distribute the allometric load by return fraction, then per unit area.
    fraction = counts / counts.sum()
    layer_load = fraction * W_TOTAL            # kg/m^2 per layer

    # kg/m^2 over a 1 m layer -> kg/m^3 bulk-density profile.
    bulk_density = layer_load / BIN_HEIGHT

    # 1 m running mean, then its maximum.
    window = max(int(round(RUN_WINDOW_M / BIN_HEIGHT)), 1)
    kernel = np.ones(window) / window
    smoothed = np.convolve(bulk_density, kernel, mode="same")
    return float(smoothed.max())


if __name__ == "__main__":
    cbd = estimate_cbd("normalized_tile.laz")
    print(f"CBD = {cbd:.4f} kg/m^3")

The running mean divides the profile’s noise while preserving the location and height of the densest continuous layer. Because the profile is already per-unit-area (kg/m²) before the division by BIN_HEIGHT, the horizontal cell area cancels and CBD comes out in kg/m³ independent of the grid footprint — which is why the allometric anchor W_TOTAL is expressed per square metre. Where W_TOTAL varies across a tile, replace the constant with a per-cell lookup from the fitted allometric raster described in the parent fuel-mapping guide.

Parameter reference

These are the parameters that change the CBD estimate. The defaults suit temperate conifer stands at typical airborne point density; the edge cases in the rationale column are where they need tuning.

The allometric anchor itself carries coefficients you calibrate to the forest type — typically an equation of the form $W=a\cdot H^b\cdot C^c$ relating available canopy fuel load to LiDAR-derived height $H$ and cover $C$. Those a, b, c coefficients belong to the allometric model, not this estimator, and should come from a regional fit against destructively sampled plots.

Expected output and verification

A single conifer stand typically returns a CBD between roughly 0.05 and 0.20 kg/m³; values above about 0.3 kg/m³ are physically rare and usually signal a bug — an un-normalized cloud, noise returns, or a mis-scaled allometric anchor. Guard the estimate:

def check_cbd(cbd: float) -> float:
    """Range and plausibility check on a CBD estimate."""
    assert cbd >= 0.0, "CBD cannot be negative"
    assert cbd <= 0.5, f"CBD = {cbd:.3f} kg/m^3 is implausibly high — check inputs"
    if cbd > 0.0:
        # A real canopy that carries crown fire clears the CBH threshold.
        assert cbd >= 0.011, "CBD below the crown-fire threshold — likely sparse or bad data"
    return cbd

The definitive verification is a regression of LiDAR-derived CBD against field-measured CBD on calibration plots: report the residual standard error and bias, and confirm the estimator does not systematically under- or over-predict across the CBD range. A round-trip check also helps — the height at which the smoothed profile hits its maximum should sit in the upper-middle canopy, not at the ground or the very top.

Common pitfalls

• Forgetting the running mean. Taking the raw profile maximum instead of the running-mean maximum lets a single noisy 1 m layer set the CBD, which is exactly the instability the method exists to remove. Always smooth before taking the maximum.
• Mismatched window and bin. RUN_WINDOW_M is expressed in metres and converted to a bin count; if you hard-code a window of 1 while binning at 0.5 m, you are averaging over 0.5 m, not 1 m. Keep the window in metres and derive the bin count.
• Scaling by cell area. Dividing the profile by both the cell area and the layer thickness double-counts and drives CBD toward zero. The profile is already per-unit-area after the load distribution; only divide by BIN_HEIGHT.
• Trusting sparse cells. At low point density a 1 m layer may hold only a few returns, so the running-mean maximum swings wildly. Flag cells below about 4 returns/m² rather than reporting their CBD as if it were reliable.

Frequently Asked Questions

Related

• Fuel Load Mapping from LiDAR — the parent workflow that derives all four canopy fuel metrics
• Fire Risk & Fuel Assessment — the section this CBD estimate feeds into
• LiDAR Point Cloud Preprocessing — the ground classification and outlier removal this method assumes
• Canopy Height Model Creation — the height normalization that makes the vertical profile valid