Removing Data Pits and Spikes From a Canopy Height Model

A rasterized canopy height model straight out of a highest-hit interpolation is rarely clean: it is pockmarked with data pits — single cells that drop far below their neighbours because a laser pulse slipped past the crown apex and registered a first return deep inside the canopy — and occasionally scarred by spikes, where a spurious high return lifts one cell above the true surface. Both defects corrupt every downstream metric, from tree-top detection to biomass allometry. This guide is a focused task within Canopy Height Model Creation, part of the broader Canopy Height Modeling & Terrain Extraction workflow: it explains the two point-cloud-level algorithms that prevent pits at rasterization time — the pit-free (Khosravipour) triangulation and the spike-free TIN — and gives a runnable scipy.ndimage pit-fill you can apply to a CHM that is already on disk.

A pit is not missing data. The cell has a value; it is simply the wrong value, because the interpolator honoured a below-canopy first return as if it were the crown top. That distinction drives the whole choice of remedy: you either stop the bad return from winning the cell during rasterization, or you detect and repair the resulting depression afterwards.

When to Use Each Approach

The right method depends on whether you still hold the point cloud and how much control you have over the rasterization step. All three below target the same visual defect but act at different stages of the pipeline.

The two point-cloud algorithms are strictly better when available, because they never fabricate height — they choose among real returns. Reach for the post-hoc fallback only when the point cloud is gone or the CHM came from a third party. If your pits stem from sparse or noisy returns rather than interpolation, the fix belongs further upstream in LiDAR Point Cloud Preprocessing, not in raster patching.

Minimal Reproducible Example

The point-cloud algorithms live in dedicated tools. The pit-free surface is a one-liner in lidR:

# R / lidR — reference pit-free CHM (Khosravipour et al. 2014)
library(lidR)
las <- readLAS("normalized.laz", filter = "-keep_first")
chm <- rasterize_canopy(las, res = 0.5,
                        algorithm = pitfree(thresholds = c(0, 2, 5, 10, 15, 20),
                                            max_edge = c(0, 1.5)))
writeRaster(chm, "chm_pitfree.tif")

The spike-free surface is produced by LAStools, driven from Python via subprocess; the freeze distance is the only critical knob:

import subprocess

# Spike-free CHM (Khosravipour et al. 2016) via LAStools las2dem.
subprocess.run([
    "las2dem", "-i", "normalized.laz", "-keep_first",
    "-spike_free", "1.0",        # freeze distance in metres
    "-step", "0.5", "-o", "chm_spikefree.tif",
], check=True)

When you only hold the finished raster, repair it in pure Python. Grey-scale morphological closing (dilation then erosion) raises narrow depressions to the level of their surroundings; restricting the write to cells the closing lifted by more than a set depth means genuine canopy structure is left untouched:

import numpy as np
import rasterio
from scipy import ndimage


def fill_chm_pits(chm_path: str, out_path: str,
                  size: int = 3, pit_depth: float = 1.0) -> str:
    """Fill isolated data pits in a rasterized CHM with grey-scale closing.

    A cell is treated as a pit only where closing raises it by more than
    `pit_depth` metres, so real crowns and edges are preserved.
    """
    with rasterio.open(chm_path) as src:
        chm = src.read(1).astype(np.float32)
        profile = src.profile
        nodata = src.nodata

    valid = np.isfinite(chm) if nodata is None else (chm != nodata) & np.isfinite(chm)
    # Grey-scale closing ignores NaN badly, so fill voids with a low sentinel first.
    filled = np.where(valid, chm, -9999.0)
    closed = ndimage.grey_closing(filled, size=size)

    is_pit = valid & ((closed - chm) > pit_depth)
    out = chm.copy()
    out[is_pit] = closed[is_pit]

    profile.update(dtype="float32")
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(out, 1)
    print(f"Filled {int(is_pit.sum())} pit cells")
    return out_path


if __name__ == "__main__":
    fill_chm_pits("chm_raw.tif", "chm_filled.tif", size=3, pit_depth=1.0)

A median filter is the alternative repair. It is more aggressive — ndimage.median_filter(chm, size=3) smooths every cell, not just pits — so gate it the same way, writing the median only where a pit was detected. Closing is generally preferable because it is asymmetric: it fills pits without also clipping legitimate spikes and sharp crown tops.

Parameter Reference

The freeze distance and pit_depth are the parameters that move results most. Set the freeze distance too small and the spike-free TIN reintroduces pits; too large and it flattens genuine gaps between neighbouring trees. Record whichever value you used in the output metadata so a reviewer can reproduce the surface.

Expected Output and Verification

A pit-corrected CHM should be free of isolated deep cells while preserving crown tops and true canopy gaps. The cheapest quantitative check is the count and depth of local minima before and after: a good fill removes narrow pits without changing the 95th–99th height percentiles, which represent real dominant-tree tops.

import numpy as np
import rasterio
from scipy import ndimage


def pit_report(path: str, pit_depth: float = 1.0) -> dict:
    with rasterio.open(path) as src:
        chm = src.read(1).astype(np.float32)
        chm = np.where(chm == src.nodata, np.nan, chm)

    # A pit is a cell more than pit_depth below its 3x3 neighbourhood max.
    neigh_max = ndimage.maximum_filter(np.nan_to_num(chm, nan=-9999.0), size=3)
    pits = np.isfinite(chm) & ((neigh_max - chm) > pit_depth)
    return {
        "n_pits": int(pits.sum()),
        "p95_height_m": float(np.nanpercentile(chm, 95)),
        "p99_height_m": float(np.nanpercentile(chm, 99)),
    }


before = pit_report("chm_raw.tif")
after = pit_report("chm_filled.tif")

# The fill must remove most pits without lowering the tall-tree percentiles.
assert after["n_pits"] < before["n_pits"] * 0.2, "pit fill left too many pits"
assert abs(after["p99_height_m"] - before["p99_height_m"]) < 0.5, "crown tops altered"

If n_pits barely drops, the window size is too small to span the pits; increase it by one and re-run. If the 99th percentile falls by more than a metre, the fill is clipping real crowns — raise pit_depth so only deeper depressions qualify. As a visual confirmation, a hillshade of the corrected CHM should show smooth crown surfaces rather than a stippled, noisy texture.

Common Pitfalls

• Filling before normalization. Pit removal assumes the CHM is height-above-ground. If you patch a raw surface model that still carries terrain, closing spans real topographic dips and injects bias. Normalize first, using the ground surface from Digital Terrain Model Generation.
• NaN leaking into the morphology. scipy.ndimage treats NaN and nodata sentinels as ordinary numbers, so a void next to a pit smears a garbage value across the neighbourhood. Fill voids with a low sentinel (as in fill_chm_pits) and restore them afterward, or mask them out entirely.
• Over-aggressive median smoothing. A blanket median_filter over the whole raster erases the fine crown structure that tree-top detection depends on. Always restrict the repair to detected pit cells rather than rewriting every pixel.
• Confusing spikes with tall trees. A genuine emergent tree looks like a spike to a naive high-pass filter. Remove spikes at the point-cloud level with the spike-free freeze distance, which distinguishes an isolated erroneous return from a coherent crown, rather than clipping the raster by a fixed height cap.

Frequently Asked Questions

Related

• Canopy Height Model Creation — the parent workflow that builds the CHM these methods clean
• Calculating Canopy Cover From CHM in Python — a cover statistic that pits would otherwise bias
• LiDAR Point Cloud Preprocessing — where sparse-return and noise problems are best solved upstream
• Digital Terrain Model Generation — the ground surface that must normalize the CHM before any pit fill

Up: Canopy Height Model Creation · Canopy Height Modeling & Terrain Extraction