Computing Vertical Canopy Complexity Metrics from LiDAR in Python

Two stands can share a mean canopy height of 24 m and be entirely different forests: one an even-aged plantation with a single dense layer, the other an old-growth remnant with a broken overstory, a mid-storey, and regeneration beneath. Structural complexity metrics are how that difference becomes a number a model can use — and they are among the strongest predictors available for bird and bat occupancy, for microclimate, and for crown fire potential. This guide computes them from a normalized point cloud as part of Forest Gap & Understory Analysis, inside Canopy Height Modeling & Terrain Extraction.

When to use each metric

FHD and VCI describe the arrangement of foliage; rumple describes the shape of the top; CRR describes where the mass sits. They are complementary rather than interchangeable, and a habitat model usually wants at least one from each family.

Foliage height diversity and the Vertical Complexity Index

FHD is Shannon entropy applied to the proportion of returns in each height bin. With $p_i$ the proportion of returns in bin $i$ across $n$ bins:

$$\mathrm{FHD}=-\sum _{i=1}^n{p}_i\ln p_i,\mathrm{VCI}=\frac{\mathrm{FHD}}{\ln n}$$

FHD is maximised when returns are spread evenly through every layer and near zero when they all sit in one. Because its maximum is $\ln n$, comparing two stands binned differently is meaningless — VCI exists to remove that dependence, giving a 0–1 value that is comparable as long as the bin width is fixed.

Minimal reproducible example

Every metric below reads one array of normalized heights for a plot or a raster cell neighbourhood.

import numpy as np


def height_profile(heights: np.ndarray, bin_m: float = 1.0,
                   min_h: float = 0.5, max_h: float | None = None):
    """Proportion of returns per fixed-width height bin, above the ground layer."""
    h = heights[np.isfinite(heights) & (heights >= min_h)]
    if h.size == 0:
        return np.array([]), np.array([])
    top = max_h if max_h is not None else float(np.ceil(h.max()))
    edges = np.arange(min_h, top + bin_m, bin_m)
    counts, _ = np.histogram(h, bins=edges)
    total = counts.sum()
    return (counts / total if total else counts.astype("float64")), edges


def fhd_vci(heights: np.ndarray, bin_m: float = 1.0) -> dict:
    """Foliage height diversity and the Vertical Complexity Index."""
    p, edges = height_profile(heights, bin_m=bin_m)
    if p.size == 0:
        return {"fhd": np.nan, "vci": np.nan, "n_bins": 0}
    nz = p[p > 0]
    fhd = float(-(nz * np.log(nz)).sum())
    n_bins = int(len(p))
    return {"fhd": fhd, "vci": fhd / np.log(n_bins) if n_bins > 1 else 0.0,
            "n_bins": n_bins}

Excluding returns below min_h is not a detail. Ground and near-ground returns usually dominate the count, and including them makes FHD a measure of how much bare ground the laser saw rather than how the foliage is arranged.

Canopy relief ratio and the height percentiles are one-liners over the same array:

def crr(heights: np.ndarray, min_h: float = 0.5) -> float:
    """Canopy relief ratio: (mean - min) / (max - min). Low = top-heavy."""
    h = heights[np.isfinite(heights) & (heights >= min_h)]
    if h.size < 2:
        return float("nan")
    lo, hi = float(h.min()), float(h.max())
    return float((h.mean() - lo) / (hi - lo)) if hi > lo else float("nan")


def structure_metrics(heights: np.ndarray, bin_m: float = 1.0) -> dict:
    """The standard plot-level bundle in one pass."""
    h = heights[np.isfinite(heights) & (heights >= 0.5)]
    out = {"n_returns": int(h.size)}
    if h.size < 10:
        return out | {"p95": np.nan, "crr": np.nan, "vci": np.nan}
    out |= {
        "p25": float(np.percentile(h, 25)),
        "p50": float(np.percentile(h, 50)),
        "p95": float(np.percentile(h, 95)),
        "height_sd": float(h.std(ddof=1)),
        "crr": crr(h),
        "cover_above_2m": float((h >= 2.0).mean()),
    }
    return out | fhd_vci(h, bin_m=bin_m)

A canopy relief ratio near 0.3 says most foliage sits low relative to the tallest tree — a young stand or one with emergent veterans over a short matrix. Near 0.7 says the mass is concentrated near the top, which is the classic even-aged closed plantation.

Rumple: the roughness of the canopy top

Rumple is the ratio of canopy surface area to the ground area beneath it — the three-dimensional analogue of a coastline’s crenulation. It is computed from the canopy height model rather than the point cloud, by triangulating each cell quad and summing triangle areas.

import numpy as np
import rasterio


def rumple(chm_path: str) -> float:
    """Canopy surface area divided by ground area, from a CHM."""
    with rasterio.open(chm_path) as src:
        z = src.read(1, masked=True).astype("float64").filled(np.nan)
        c = abs(src.transform.a)

    # Split every cell quad into two triangles and sum their 3-D areas.
    z00, z10 = z[:-1, :-1], z[1:, :-1]
    z01, z11 = z[:-1, 1:], z[1:, 1:]

    def tri_area(a, b, cc):
        # Two edge vectors per triangle; the cross-product magnitude is twice the area.
        u = np.stack([np.full_like(a, c), np.zeros_like(a), b - a], axis=-1)
        v = np.stack([np.zeros_like(a), np.full_like(a, c), cc - a], axis=-1)
        return 0.5 * np.linalg.norm(np.cross(u, v), axis=-1)

    surface = np.nansum(tri_area(z00, z01, z10)) + np.nansum(tri_area(z11, z10, z01))
    ground = np.isfinite(z[:-1, :-1]).sum() * c * c
    return float(surface / ground) if ground else float("nan")

A rumple of 1.0 is a perfectly flat canopy; a closed even-aged plantation sits near 1.5–2.5; complex old growth reaches 4 and beyond. Because it is a ratio of areas it is strongly resolution-dependent — a 0.5 m CHM resolves individual branch clusters that a 2 m CHM smooths away — so record the cell size beside every rumple value and never compare across resolutions.

Parameter reference

Expected output and verification

The metrics are cheap; the discipline is in confirming they mean what you think. Three checks catch the common failures:

import numpy as np

m = structure_metrics(plot_heights, bin_m=1.0)

assert 0.0 <= m["vci"] <= 1.0, "VCI outside 0-1: bin count or normalisation is wrong"
assert m["p95"] >= m["p50"] >= m["p25"], "percentiles out of order"
assert m["n_returns"] >= 100, "too few returns for a stable entropy estimate"

# A single-layer stand must score lower than a multi-layer one on the same bins.
single = np.random.normal(24.0, 1.0, 5000)
layered = np.concatenate([np.random.normal(24.0, 1.5, 2500),
                          np.random.normal(8.0, 2.0, 1500),
                          np.random.uniform(1.0, 30.0, 1000)])
assert fhd_vci(single)["vci"] < fhd_vci(layered)["vci"], "VCI is not tracking layering"

The last assertion is the one worth keeping in a test suite. It feeds the machinery a case it must reject — a synthetic single-layer stand that is supposed to score lower — and so it catches a binning or normalisation bug that the real data would hide.

Validating against the field is a matter of correlation rather than agreement, because no crew measures FHD directly. The usual anchors are a visual structural-stage classification, understory light measured with hemispherical photography, or a species count from point transects; a Spearman correlation above about 0.6 against any of them is a good result at plot scale.

Common pitfalls

• Comparing FHD across different bin widths. Entropy scales with bin count. Report VCI, state the bin width, and freeze both for the whole study.
• Including ground returns. Ground usually contributes the largest single bin, so an unfiltered profile measures canopy openness rather than foliage arrangement.
• Computing rumple from a pit-filled CHM without saying so. Pit filling smooths the canopy top and lowers rumple by a visible margin. Whichever surface you use, use it consistently.
• Applying plot metrics to raster cells without checking sample size. A 10 m cell at 5 returns/m² holds about 500 returns, which is fine; a 5 m cell holds 125 and the entropy estimate starts to wobble.
• Ignoring pulse density differences between acquisitions. FHD is sensitive to how deeply the laser penetrated. Comparing a 4 returns/m² survey with a 20 returns/m² survey measures the sensor, not the forest.

Frequently Asked Questions

Related

• Forest Gap & Understory Analysis — the wider treatment of what happens beneath the canopy top
• Estimating Gap Fraction and LAI from LiDAR — the light-transmission view of the same return profile
• Normalizing LiDAR Point Clouds with PDAL — produces the height-above-ground field every metric here needs
• Estimating Canopy Bulk Density from LiDAR — the fuel-modelling use of the same vertical profile

Up: Forest Gap & Understory Analysis · Canopy Height Modeling & Terrain Extraction