The Swiss Cheese Problem: Building Pit-Free Canopy Height Models

In forestry remote sensing, generating a Canopy Height Model (CHM) sounds delightfully simple: take your LiDAR surface model (DSM), subtract the bare-earth ground model (DTM), and boom: you have tree heights:

Then you zoom in on a mature Douglas-fir stand, and the tree crowns look like Swiss cheese.

Why Conifers Break Standard Rasters

Conifers are tricky beasts. Unlike broadleaf deciduous trees with dense, umbrella-like canopies, conifer branches are tiered and needles are sparse. When an airborne or drone LiDAR fires 40 pulses per square meter into a stand, a good chunk of those laser beams punch straight between the branches and hit needles 10 or 15 meters below the crown apex.

If you interpolate those points with standard grid binning or naive Delaunay triangulation, those deep interior hits drag the surface down into narrow, artificial craters known as canopy pits.

Run an automated tree-segmentation algorithm on that raw CHM, and it will think one Douglas-fir is actually four tiny dwarf trees clinging together for warmth.

The Tempting (and Wrong) Fix

The first instinct is usually: “Let’s just slap a Gaussian blur or morphological closing filter over it to smooth out the holes!”

Please don’t. While that fills the pits, it also rounds off the needle-sharp terminal leaders at the tops of the trees. Suddenly your 38-meter dominant trees are smoothed down to 34 meters. In a commercial timber inventory or carbon modeling project, shaving four meters off the top of every tree adds up to massive volume errors.

The Trick: Multi-Scale TIN Envelopes

The elegant fix is surprisingly simple: build multiple Triangulated Irregular Networks (TINs) at ascending height strata, and drape them over each other.

By throwing away low points before building each upper mesh, deep interior pulses can never form a triangular face with the apex returns:

import numpy as np

def build_pit_free_envelope(canopy_points, height_cutoffs=[2.0, 10.0, 20.0, 35.0], resolution=0.5):
    """
    Build successive TIN surfaces above ascending height thresholds
    and take the element-wise maximum envelope.
    """
    master_chm = initialize_empty_grid(canopy_points, resolution)

    for z_min in height_cutoffs:
        stratum = canopy_points[canopy_points[:, 2] >= z_min]
        if len(stratum) < 10:
            continue
        
        # Triangulate only points at or above this stratum
        layer_surface = rasterize_tin(stratum, resolution)
        master_chm = np.maximum(master_chm, layer_surface)

    return master_chm

By taking the maximum envelope across the strata, the pits vanish, the sharp apical tips stay crisp, and your watershed segmentation algorithms stop hallucinating micro-trees.

In our lidar-canopy-extractor pipeline, we chunk this across 15-meter buffered spatial tiles so it can chew through 100M+ points in Docker without consuming your machine’s entire memory pool.