Separating Buildings from Vegetation in Point Clouds

A ground filter answers one question and leaves a second one open. After classifying point clouds with PDAL and Python has run, the cloud holds a ground class and an undifferentiated mass of everything else — roofs, trees, hedges, parked vehicles, scaffolding and the occasional crane. The moment a deliverable needs building footprints, canopy height, or a vegetation mask for a volume boundary, that mass has to be split.

The symptom that brings people here is usually one of two: a building-footprint extraction that returns tree crowns as buildings, or a canopy height model with rectangular holes where the roofs were counted as vegetation. Both come from using a single feature where two are needed.

# Why one feature is never enough

Height above ground separates tall things from short things and nothing more. A twelve-metre tree and a twelve-metre warehouse are identical to it. Planarity — how well a point’s neighbourhood fits a plane — separates flat things from rough things, and on its own it promotes every stretch of tarmac and every flat-topped hedge into the building class.

The features are complementary because their failure modes do not overlap. A pitched roof has lower planarity than a flat one but is still far more planar than any canopy at the same height. A dense conifer can present a locally smooth surface but never at the scale of a roof plane. And the colour a photogrammetric point carries from its source imagery adds a third axis that is nearly free: a roof is rarely green, and a summer canopy almost always is.

The mistake worth naming is computing planarity at the wrong scale. Planarity is derived from the eigenvalues of the local covariance matrix over knn neighbours. With knn at 8 on a 2 cm cloud, that neighbourhood spans about 10 cm — a scale at which roof gravel is rough and leaves are smooth. The feature inverts. The neighbourhood must span at least a metre for the descriptor to mean what its name suggests.

Figure 1 — The single most common cause of an inverted building mask, and its fix.

# Minimal reproducible solution

Three features, one pass, explicit thresholds.

import json
import subprocess


def split_nonground(src: str, dst: str, *, knn: int = 32,
                    planarity_min: float = 0.92, roof_min_height: float = 2.5,
                    egi_max_for_roof: float = 0.05) -> None:
    """Split non-ground points into buildings (6), high (5) and low (3) vegetation.

    The rules are evaluated in order, and each one excludes the classes already
    assigned, so a point can only be claimed once. Colour is used as a veto on
    the building class rather than as a positive vegetation test: a grey tree
    is rare, a green roof is not.
    """
    egi = "(2.0 * Green - Red - Blue) / (Red + Green + Blue + 1)"
    pipeline = {"pipeline": [
        src,
        {"type": "filters.hag_nn"},
        {"type": "filters.covariancefeatures", "knn": knn,
         "feature_set": "Dimensionality", "threads": 4},
        {"type": "filters.ferry", "dimensions": "=>ExcessGreen"},
        {"type": "filters.assign", "value": [
            f"ExcessGreen = {egi}",
            # buildings: planar, tall, and not green
            f"Classification = 6 WHERE Classification != 2"
            f" && Planarity > {planarity_min}"
            f" && HeightAboveGround > {roof_min_height}"
            f" && ExcessGreen < {egi_max_for_roof}",
            # high vegetation: tall, not already a building
            f"Classification = 5 WHERE Classification != 2 && Classification != 6"
            f" && HeightAboveGround > {roof_min_height}",
            # low vegetation: short but off the ground
            f"Classification = 3 WHERE Classification != 2 && Classification != 6"
            f" && Classification != 5 && HeightAboveGround > 0.3",
        ]},
        {"type": "writers.las", "filename": dst, "compression": "laszip",
         "forward": "all", "extra_dims": "Planarity=float,ExcessGreen=float"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

Using colour as a veto rather than a test is the design decision that matters. A positive rule — “green means vegetation” — misfires on green roofs, algae-covered flat roofs and any surface photographed under a green canopy’s reflected light. A veto — “a building must not be strongly green” — fails safe: the worst case is a green roof landing in vegetation, which is visible and rare, rather than a whole tree line landing in buildings.

Carrying Planarity and ExcessGreen out as extra dimensions costs about eight bytes per point and buys the ability to answer “why was this classified that way” without re-running anything.

# Edge-case matrix

# Verification snippet

The useful test is not “did it run” but “does the split agree with something independent”. Building footprints are usually available as a vector layer, and a spatial agreement rate against them is a one-off cost that pays for itself.

import json
import numpy as np
import pdal
from shapely.geometry import Point, shape
from shapely.strtree import STRtree


def footprint_agreement(las_path: str, footprints_geojson: str,
                        sample: int = 50_000) -> dict:
    """Agreement between classified buildings and a known footprint layer."""
    geoms = [shape(f["geometry"])
             for f in json.load(open(footprints_geojson))["features"]]
    tree = STRtree(geoms)

    pipe = pdal.Pipeline(json.dumps({"pipeline": [las_path]}))
    pipe.execute()
    arr = pipe.arrays[0]
    idx = np.random.default_rng(0).choice(arr.size, size=min(sample, arr.size),
                                          replace=False)

    inside_and_building = inside = building = 0
    for i in idx:
        p = Point(float(arr["X"][i]), float(arr["Y"][i]))
        in_fp = any(geoms[j].contains(p) for j in tree.query(p))
        is_b = int(arr["Classification"][i]) == 6
        inside += in_fp
        building += is_b
        inside_and_building += in_fp and is_b

    recall = inside_and_building / inside if inside else float("nan")
    precision = inside_and_building / building if building else float("nan")
    return {"precision": precision, "recall": recall,
            "sampled": int(idx.size), "in_footprint": int(inside)}

Interpret the two numbers separately. Low recall with high precision means the thresholds are conservative — roofs are being missed, usually pitched ones. High recall with low precision means something else is being called a building, and the usual culprit is a flat-topped hedge or a stretch of raised hardstanding that cleared the height test.

Figure 2 — Two numbers, four diagnoses, four different parameters. A single “accuracy” figure collapses this into something nobody can act on.

Figure 3 — The cascade, with a typical mixed-use site’s shares. The unclassified remainder is a feature: forcing every point into a class hides the ambiguity rather than resolving it.

# Calibrating the thresholds on a site you know

Default thresholds are a starting point, not an answer. The values that work on a suburban business park fail on a livestock farm with polytunnels, and both fail on a construction site where the “buildings” are shipping containers and partially erected steel frames. Calibration takes about twenty minutes per site type and only has to happen once.

The procedure is to pick one small area where the truth is obvious — a block of four known buildings and a stand of known trees — and sweep the two thresholds over it, reading the precision and recall at each combination rather than choosing by eye.

import itertools
import numpy as np


def sweep_thresholds(features: np.ndarray, truth: np.ndarray,
                     planarity_grid=(0.86, 0.89, 0.92, 0.95),
                     height_grid=(1.5, 2.0, 2.5, 3.5)) -> list[dict]:
    """Grid-search the two building thresholds against a hand-labelled patch.

    `features` carries Planarity and HeightAboveGround per point; `truth` is a
    boolean array that is True for points a reviewer called building.
    """
    rows = []
    for pmin, hmin in itertools.product(planarity_grid, height_grid):
        pred = (features["Planarity"] > pmin) & (features["HeightAboveGround"] > hmin)
        tp = int(np.count_nonzero(pred & truth))
        fp = int(np.count_nonzero(pred & ~truth))
        fn = int(np.count_nonzero(~pred & truth))
        rows.append({
            "planarity_min": pmin, "min_height": hmin,
            "precision": tp / (tp + fp) if tp + fp else 0.0,
            "recall": tp / (tp + fn) if tp + fn else 0.0,
        })
    return sorted(rows, key=lambda r: -(2 * r["precision"] * r["recall"]
                                        / max(r["precision"] + r["recall"], 1e-9)))

Two things usually emerge from the sweep. The first is that the surface is flat near the optimum — several combinations score within a percentage point of each other, which means the exact value does not matter and arguing about it is wasted effort. The second is that the ridge runs diagonally: a lower planarity cutoff can be compensated by a higher height floor, because both exclude the same population of flat low clutter. Knowing the shape of that trade-off is more useful than knowing the peak, because it tells you which knob to turn when a new site misbehaves.

Record the chosen pair per site type in the job definition rather than in code. A pipeline that reads thresholds from the job is one where a new site type is a configuration change; a pipeline with the numbers inline is one where it is a release.

# When to escalate

• Precision and recall are both poor after fixing the scale. The reconstruction is probably too sparse on roofs — a common outcome when flight overlap was tuned for terrain rather than structures. See calculating optimal flight overlap for Python processing.
• The site has buildings under canopy. Photogrammetry has no data there at all, and no classifier can produce it. Flag the region and state it in the deliverable.
• Footprints are the actual deliverable. Point classification is the wrong tool for a vector product; extract planes and fit polygons instead, using this classification only as the input mask.

# Related

• Classifying Point Clouds with PDAL and Python
• Fixing Misclassified Ground Under Dense Canopy
• Removing Vegetation from DTM Ground Surfaces

← Classifying Point Clouds with PDAL and Python