Fixing Misclassified Ground Under Dense Canopy

A woodland edge runs across the survey. The digital terrain model is smooth and plausible everywhere, including under the trees, and a client using it to design a drainage run finds the levels are two metres out. Nothing in the file says so. The contours are as clean under the canopy as over the adjacent field.

What happened is not a classification error in the ordinary sense. The classifier did its job on the points it was given. Under a closed canopy a photogrammetric cloud contains no ground points to classify, so the lowest returns — leaves, branches, the occasional glimpse through a gap — became the ground class, and the interpolation between them produced a surface that looks exactly like terrain.

Why photogrammetry has no data under trees

The distinction that matters is between a sensor that can see through gaps and one that cannot. Airborne lidar emits a pulse that penetrates a canopy through small openings and returns several times; the last return is usually the ground, which is why lidar produces bare-earth models in forest at all. Photogrammetry reconstructs a point only where the same piece of the world is visible and matchable in at least two images. Under a closed canopy there is no such piece of ground: it is occluded in every frame, or visible in one frame from one angle and never matched.

The result is not sparse ground. It is absent ground. And a filter such as SMRF, described in classifying point clouds with PDAL and Python, has no way to express “there was nothing here”. It selects the lowest points in each window, and under a canopy those are canopy.

Why a canopy produces a terrain surface that is entirely invented A cross-section through a field running into woodland. Over the open field, reconstructed points sit densely on the true ground and the classified terrain follows it exactly. Under the closed canopy, the only reconstructed points are on the leaf surface eighteen metres up, with nothing beneath; the classified terrain therefore follows the underside of the canopy rather than the ground, sitting about two metres above the true surface at the woodland edge and far higher inside. A note records that the interpolated surface is visually indistinguishable from the measured one. true ground open field: classified = measured every reconstructed point is canopy classified "ground" — invented error grows inward Nothing in the output distinguishes the left half from the right half. That is the whole problem.

Figure 1 — The surface under the canopy is not badly measured. It was never measured.

Minimal reproducible solution

The fix is not a better filter. It is a mask: find where the ground class is really ground and where it is a guess, and carry that distinction into the deliverable. The signal is the vertical spread of points above each ground candidate together with the density of the ground class itself.

import json
import numpy as np
import pdal


def ground_confidence_grid(las_path: str, cell: float = 5.0,
                           min_ground_per_cell: int = 3,
                           canopy_cover_max: float = 0.6) -> dict:
    """Per-cell verdict: measured ground, sparse ground, or no ground at all.

    Two independent signals. Ground point density says whether the classifier
    had anything to work with. Canopy cover — the share of points in the cell
    that sit well above the local ground — says whether an absence is because
    the area is occluded or simply because it was outside the flight.
    """
    pipe = pdal.Pipeline(json.dumps({"pipeline": [
        las_path, {"type": "filters.hag_nn"}]}))
    pipe.execute()
    a = pipe.arrays[0]

    ix = np.floor(a["X"] / cell).astype(np.int64)
    iy = np.floor(a["Y"] / cell).astype(np.int64)
    keys = ix * 1_000_003 + iy

    is_ground = a["Classification"] == 2
    is_canopy = a["HeightAboveGround"] > 2.0

    verdict = {}
    for k in np.unique(keys):
        m = keys == k
        n_ground = int(np.count_nonzero(is_ground & m))
        cover = float(np.count_nonzero(is_canopy & m) / max(np.count_nonzero(m), 1))
        if n_ground >= min_ground_per_cell and cover < canopy_cover_max:
            verdict[int(k)] = "measured"
        elif n_ground > 0:
            verdict[int(k)] = "sparse"
        else:
            verdict[int(k)] = "absent"
    return verdict

Three states, not two. “Sparse” is the interesting one: a cell with one or two ground points under 80 % canopy cover has ground of a kind, but a surface interpolated from it carries an uncertainty an order of magnitude larger than the open-field value, and it should be delivered with that uncertainty attached rather than silently averaged into the same product.

Writing the verdict out as a raster alongside the terrain model is the part that makes it useful downstream:

import numpy as np
import rasterio
from rasterio.transform import from_origin


def write_confidence_raster(verdict: dict, bounds: tuple, cell: float,
                            crs: str, out_path: str) -> None:
    """Emit the ground-confidence verdict as a companion raster to the DTM."""
    minx, miny, maxx, maxy = bounds
    w = int(np.ceil((maxx - minx) / cell))
    h = int(np.ceil((maxy - miny) / cell))
    codes = {"measured": 1, "sparse": 2, "absent": 3}
    grid = np.zeros((h, w), dtype="uint8")

    for key, state in verdict.items():
        ix, iy = divmod(key, 1_000_003)
        col = int(ix - np.floor(minx / cell))
        row = int(h - 1 - (iy - np.floor(miny / cell)))
        if 0 <= row < h and 0 <= col < w:
            grid[row, col] = codes[state]

    with rasterio.open(out_path, "w", driver="GTiff", height=h, width=w,
                       count=1, dtype="uint8", crs=crs, nodata=0,
                       transform=from_origin(minx, maxy, cell, cell)) as dst:
        dst.write(grid, 1)
        dst.update_tags(1, CLASSES="1=measured, 2=sparse, 3=no ground data")

Edge-case matrix

Input variant Naive DTM behaviour Correct handling
Closed broadleaf canopy, summer Smooth invented surface Mark absent, do not interpolate
Same woodland, leaf-off winter flight Partly real ground Re-fly in winter if terrain under trees is required
Open woodland, scattered trees Mostly correct Mark sparse; interpolation is defensible
Hedge line, 2 m wide Interpolated across, usually fine Mark sparse; the span is short
Building overhang, canopy-like Ground absent under it Same treatment as canopy
Tall crop (maize, at height) Ground absent, error ≈ crop height Mark absent; this is frequently missed
Deep shadow under canopy edge Points present but noisy Filter noise first, then re-assess
Water body No ground, no canopy Separate mask; different failure

The tall-crop row catches teams out most often, because the visual signature is not obviously “forest”. A maize field in August has closed cover at two and a half metres and behaves exactly like a canopy, which turns an apparently open agricultural site into one where half the terrain model is invented.

Verification snippet

Verification means testing the claim rather than the appearance. If checkpoints exist under and outside the canopy, the difference between their residual distributions is the whole finding.

import numpy as np


def canopy_bias_report(residuals_open: np.ndarray,
                       residuals_canopy: np.ndarray) -> dict:
    """Compare DTM residuals inside and outside canopy; fail if they differ."""
    def summary(r):
        r = r[np.isfinite(r)]
        return {"n": int(r.size), "bias": float(np.mean(r)),
                "rmse": float(np.sqrt(np.mean(r ** 2)))}

    op, cp = summary(residuals_open), summary(residuals_canopy)
    ratio = cp["rmse"] / op["rmse"] if op["rmse"] else float("inf")

    report = {"open": op, "canopy": cp, "rmse_ratio": ratio}
    if ratio > 3.0:
        report["verdict"] = (
            "canopy terrain is not a measurement — deliver it masked, "
            f"not blended ({cp['rmse']:.2f} m vs {op['rmse']:.2f} m)")
    elif cp["bias"] > 0.3:
        report["verdict"] = (
            f"systematic +{cp['bias']:.2f} m under canopy: the ground class "
            "is following vegetation")
    else:
        report["verdict"] = "canopy and open terrain agree within tolerance"
    return report
Three ground-confidence states across a survey, and what each licenses A three-column table of the confidence states. Measured ground has three or more ground points per cell and canopy cover below sixty percent, licensing full delivery with the normal accuracy statement. Sparse ground has one or two ground points per cell under heavier cover, licensing delivery with an inflated uncertainty and an explicit note. Absent ground has no ground points at all, licensing only a hole in the deliverable and a statement that the terrain there was not surveyed. A note observes that collapsing the three into one product is what produces the two-metre design error. measured ≥ 3 ground points / cell cover < 60 % deliver normally site accuracy statement applies as written sparse 1–2 ground points / cell heavier cover deliver with a wider band inflate the stated uncertainty and say where it applies absent no ground points closed canopy or tall crop deliver a hole state that terrain here was not surveyed Collapsing three states into one product is what produced the two-metre design error. A hole in a deliverable is an honest answer. A smooth invented surface is not.

Figure 2 — What each confidence state licenses. The mask is cheap; the alternative is a client discovering the problem on site.

Ground point recovery against canopy cover for two sensing methods A curve of recovered ground point density against canopy cover from zero to one hundred percent. Photogrammetry falls steeply, reaching effectively zero ground points by about seventy percent cover. Lidar declines gradually and still returns usable ground density at ninety-five percent cover. A shaded band beyond seventy percent cover marks the region where a photogrammetric terrain model is entirely interpolated. A note states that no processing parameter moves the photogrammetry curve, because the limitation is the absence of a matchable observation. terrain entirely interpolated 0 % 25 % 50 % 75 % 100 % canopy cover ground points photogrammetry lidar No parameter moves the lower curve: the observation does not exist to be processed.

Figure 3 — The limit is physical. Recognising it early converts an argument about filter settings into a scoping decision about the survey.

What to put in the deliverable

Once the confidence grid exists, the question becomes contractual rather than technical: what does the client receive for the cells marked absent. There are three defensible answers and one indefensible one.

A hole. The terrain raster carries NoData under closed canopy and the accompanying report states the extent. This is the honest default, and it is the only option that cannot mislead. Clients occasionally object that a model with holes “looks unfinished”, which is a conversation worth having before the flight rather than after.

An interpolated surface, flagged. The holes are filled by interpolation from the surrounding measured ground, the interpolated cells are marked in the companion raster, and the report states the expected error. This is reasonable for narrow gaps — a hedge line, a few isolated trees — where the interpolation spans metres rather than tens of metres. It stops being reasonable when the gap is large enough that the interpolation is unconstrained in the middle.

A surface from another source. Where national lidar or an earlier ground survey covers the wooded area, splicing it in produces a complete model, provided the vertical datums are reconciled first and the provenance of each cell is recorded. The reconciliation is the hard part and is covered in converting ellipsoidal to orthometric heights in bulk.

The indefensible answer is the default behaviour of most software: interpolate silently, deliver a seamless surface, and say nothing. It is indefensible not because the surface is wrong — an interpolation is a reasonable estimate — but because nothing in the product distinguishes an estimate from a measurement, and the client has no way to know which they are standing on.

import numpy as np


def apply_confidence_mask(dtm: np.ndarray, confidence: np.ndarray,
                          policy: str = "hole") -> np.ndarray:
    """Apply the delivery policy for cells where ground was never observed.

    confidence: 1 measured, 2 sparse, 3 absent (0 outside the survey).
    """
    out = dtm.astype("float32").copy()
    absent = confidence == 3
    if policy == "hole":
        out[absent] = np.nan
    elif policy == "flagged":
        pass                     # values retained; the companion raster carries the flag
    else:
        raise ValueError(f"unknown delivery policy {policy!r}")
    return out

Whichever policy is chosen, write it into the run manifest next to the file. Six months later, when somebody asks whether the terrain under the trees was measured, the answer should be in the metadata rather than in somebody’s memory.

When to escalate

  • The terrain under canopy is the deliverable. Photogrammetry cannot produce it. Either fly leaf-off, commission lidar for that area, or survey it on the ground — and price accordingly before the flight rather than after the processing.
  • The absent mask covers more than a small fraction of the site. That is a scoping conversation, not a processing one. A site that is thirty percent woodland needs its terrain deliverable defined before anyone processes anything.
  • Ground appears under canopy after a parameter change. Be suspicious. Loosening the classifier until points appear does not create measurements; verify against the confidence grid and sampled checkpoints, per validating classification against manual samples.

← Classifying Point Clouds with PDAL and Python