Python LiDAR

Computing Canopy Cover from LiDAR

TL;DR: Canopy cover in a cell ≈ first returns above a height threshold (commonly 2 m) ÷ all first returns in the cell. Compute it at 10–30 m cells by rasterizing two counts with writers.gdal (first returns above 2 m, and all first returns) and dividing, or bin both counts in NumPy. Report the threshold and cell size with every figure: cover estimates are only comparable when both match.

# Context and Motivation

This guide is part of Canopy Height Models. Canopy cover — the fraction of ground area covered by the vertical projection of tree crowns — is one of the most used forest structure variables: in habitat models, fire behaviour, carbon estimation, urban tree-cover targets and forest definitions themselves. LiDAR estimates it well because a first return from above the threshold is, very nearly, a pulse intercepted by canopy. Counting such returns per cell gives a physically meaningful cover fraction without segmenting a single tree.

The details that change the answer are the height threshold, which returns are counted, and the cell size. None of them is standardized, so a cover map without them stated is not reproducible.

Cover as a ratio of first returns A side view of one grid cell with trees and gaps. Ten vertical pulses fall across the cell. Seven produce first returns in crowns above the 2 metre threshold line; three reach low vegetation or the ground first. Canopy cover for the cell is seven over ten, or 70 percent. 2 m threshold cover = 7 canopy first returns ÷ 10 first returns = 70 %

# Prerequisites and Assumptions

  • A ground-classified tile for height normalization; noise removed.
  • PDAL with writers.gdal; NumPy and rasterio.
  • A definition to follow: the height threshold (2 m is common; some national forest definitions use 5 m) and the reporting cell or stand boundaries.

# Step-by-Step Implementation

# Step 1 — Normalize heights

Compute HeightAboveGround and ferry it into Z for rasterization.

# Step 2 — Keep first returns

First returns sample what a vertical view sees; later returns come from inside or below the canopy and would inflate cover.

# Step 3 — Count per cell twice

Rasterize the count of all first returns, and the count of first returns above the threshold, on the same grid.

# Step 4 — Divide

Cover = above ÷ all, where all is at least a minimum count (for example 20) to avoid noisy ratios in sparse cells.

# Step 5 — Summarize per stand

Aggregate counts (not ratios) within stand polygons and divide once, so large and small cells are weighted correctly.

# Complete Working Example

python
"""Canopy cover raster from first returns above 2 m, plus a stand summary."""
from __future__ import annotations

import json
import math

import geopandas as gpd
import numpy as np
import pdal
import rasterio
from rasterio.features import geometry_mask

SRC, RES, THRESH, MIN_COUNT = "tiles/forest_0822.laz", 20.0, 2.0, 20


def grid(src: str) -> dict:
    b = pdal.Pipeline(json.dumps({"pipeline": [src]})).quickinfo["readers.las"]["bounds"]
    ox, oy = math.floor(b["minx"] / RES) * RES, math.floor(b["miny"] / RES) * RES
    return {"origin_x": ox, "origin_y": oy,
            "width": math.ceil((b["maxx"] - ox) / RES), "height": math.ceil((b["maxy"] - oy) / RES)}


g = grid(SRC)
common = {"resolution": RES, "radius": RES * 0.71, "output_type": "count",
          "data_type": "uint32", "nodata": 0, **g}
pdal.Pipeline(json.dumps({"pipeline": [
    SRC,
    {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
    {"type": "filters.hag_nn", "count": 2},
    {"type": "filters.range", "limits": "ReturnNumber[1:1]", "tag": "first"},
    {"type": "writers.gdal", "inputs": ["first"], "filename": "out/first_all.tif", **common},
    {"type": "filters.range", "inputs": ["first"], "limits": f"HeightAboveGround[{THRESH}:]", "tag": "canopy"},
    {"type": "writers.gdal", "inputs": ["canopy"], "filename": "out/first_canopy.tif", **common},
]})).execute()

with rasterio.open("out/first_all.tif") as a, rasterio.open("out/first_canopy.tif") as c:
    all_n, can_n = a.read(1).astype(float), c.read(1).astype(float)
    profile, transform = a.profile, a.transform
cover = np.where(all_n >= MIN_COUNT, can_n / np.maximum(all_n, 1), np.nan)
profile.update(dtype="float32", nodata=-1)
with rasterio.open("out/canopy_cover_20m.tif", "w", **profile) as out:
    out.write(np.nan_to_num(cover, nan=-1).astype("float32"), 1)

stands = gpd.read_file("stands.gpkg").to_crs(profile["crs"])
for s in stands.itertuples():
    m = geometry_mask([s.geometry], all_n.shape, transform, invert=True)
    pct = can_n[m].sum() / max(all_n[m].sum(), 1)
    print(f"stand {s.stand_id}: canopy cover {pct:.0%}")

With a 20 m cell, radius of 0.71 × resolution counts returns in circles that overlap slightly; counts are then approximate but the ratio is unaffected, because numerator and denominator share the same windows.

The threshold changes the answer Bars of canopy cover for one stand computed with three height thresholds. At 1.3 metres, 81 percent, including tall shrubs. At 2 metres, 74 percent. At 5 metres, 62 percent, counting only tree canopy. The same data yields a 19 point range depending on the definition, so the threshold must always be reported. 81 %74 %62 % 1.3 m2 m5 m height threshold (illustrative stand)

# Key Parameter Table

Setting Common value Effect
height threshold 2 m (1.3 m, 5 m) What counts as canopy; state it always
returns counted first Vertical-view interception
cell size 10–30 m Enough first returns per cell for a stable ratio
MIN_COUNT 20 Cells with fewer returns left undefined
stand summary sum counts, divide once Correct weighting across cells

# Verification

  • Bounds. Cover values lie between 0 and 1; open fields near 0, closed forest near 0.9 or higher.
  • Compare with imagery. Crown cover digitized from orthophotos on a few plots should agree within about 10 percentage points; LiDAR usually reads slightly higher because it sees small gaps as covered at coarse cells.
  • Threshold sensitivity. Report cover at two thresholds for a few stands so users see how much the definition matters.

# Gotchas and Edge Cases

All returns instead of first. Counting all returns above the threshold divided by all returns measures something closer to vegetation density, not cover, and depends on sensor and penetration. Keep first returns for cover.

Scan angle. At large off-nadir angles, pulses travel through more canopy and are intercepted more often, inflating cover at swath edges. Restrict to near-nadir first returns where overlap allows.

Leaf-off data. Deciduous canopy measured leaf-off gives far lower cover than leaf-on. Say which season the data represents.

Why scan angle inflates cover Two pulses through a crown with gaps. The vertical pulse at nadir passes through a gap and reaches the ground. The slanted off-nadir pulse crosses the crown diagonally, traverses more foliage and is intercepted, even though it aims at the same gap. At swath edges more pulses are intercepted, raising apparent cover. nadir: through the gap off-nadir: intercepted

Tiny cells. At 1 m cells a ratio is built from a handful of returns and is nearly binary. Use a CHM thresholded at the canopy height for fine-grained cover maps instead.

# Frequently Asked Questions

How is canopy cover calculated from LiDAR?

Typically as the number of first returns above a height threshold divided by the total number of first returns within a cell or area. It approximates the fraction of the ground covered by the vertical projection of the canopy.

What height threshold should I use?

Two metres is common in forestry research; 1.3 metres is sometimes used to align with breast height, and some national forest definitions use 5 metres. Choose one that matches your definition and report it.

Why use first returns only?

First returns represent what is seen from above. Later returns come from inside or below the canopy, so including them makes the ratio depend on penetration rather than cover.

Can I compute cover from a CHM instead?

Yes: the fraction of CHM cells above the threshold. It works well at fine resolution but depends on CHM gap filling and pits; the return-ratio method is more robust at coarser cells.

Related in Canopy Height Models with filters.hag_nn

↑ Canopy Height Models with filters.hag_nn overview --- title: "Canopy Height Models with filters.hagnn" description: "Producing a canopy height model in the point domain: classify ground,... Building a Pit-Free Canopy Height Model --- title: "Building a Pit-Free Canopy Height Model" description: "Build a pit-free CHM with PDAL and NumPy using the layered approach of... Computing Height Above Ground with filters.hag_nn --- title: "Computing Height Above Ground with filters.hagnn" description: "What filters.hagnn actually computes, why count and maxdistance... Extracting Individual Tree Heights from a CHM --- title: "Extracting Individual Tree Heights from a CHM" description: "Local maxima detection with a crown-scaled window, why every... Rasterizing a Canopy Height Model from HAG --- title: "Rasterizing a Canopy Height Model from HAG" description: "The three writers.gdal choices that decide whether you get a canopy...