Classifying Point Clouds with PDAL and Python

Every bare-earth product a survey delivers rests on one decision made millions of times: is this point ground, or is it something standing on the ground. The classifier that makes that decision is a handful of geometric rules with four or five parameters, and the difference between a defensible digital terrain model and a plausible-looking fiction is whether those parameters matched the terrain and whether anybody checked.

This page covers classification as an automated stage in a Python pipeline: how the progressive morphological family of filters actually works, which parameter to reach for when a specific artefact appears, how to use the colour photogrammetry gives you for free, and how to prove the result rather than admire it. The classes produced here feed the surface generation in generating DSM and DTM from point clouds with PDAL and every volume computed in computing volumes and stockpiles in Python.

Audience and prerequisites. Python 3.10+, PDAL 2.5+ with its Python bindings, and a denoised dense cloud in a projected CRS. Classification on a cloud still in geographic coordinates will not work: every window and threshold parameter is in the horizontal unit of the file, and degrees are not metres.

# Prerequisites

# Conceptual architecture

The progressive morphological approach, and the Simple Morphological Filter (SMRF) that refines it, both rest on the same idea. Start from the lowest point in each cell of a coarse grid. Those points are provisionally ground. Then grow a window outward, repeatedly, and at each scale ask whether a point sits too far above the surface implied by its neighbourhood. Points that do are non-ground; the surface is re-fitted; the window grows again.

The consequence worth internalising is that the window parameter is a size of object, not a tuning knob. It is the largest thing the filter can recognise as standing on the ground. A window of 10 m over a site with a 30 m warehouse cannot see the warehouse as an object — from the filter’s perspective at that scale the roof is the terrain, and it will be classified as ground. Almost every “the building is in my DTM” report resolves to a window smaller than the building.

The threshold parameter is the complementary one: it is the height above the provisional surface at which a point stops being ground. Too small and real terrain roughness is stripped away; too large and low vegetation survives into the ground class, producing exactly the systematic bias described in validating classification against manual samples.

Figure 1 — Why a building ends up in the terrain model. The filter has no concept of “building”; it has a scale, and anything wider than that scale is ground.

# Step 1: Classify with SMRF and keep the parameters

The minimum useful pipeline is three stages: assert a projected CRS, run the filter, write the result. The parameters go into the manifest alongside the output, because a cloud classified with different parameters is a different cloud and six months later nobody will remember which was used.

import json
import subprocess
from pathlib import Path


def classify_ground(src: str, dst: str, *, window: float = 18.0,
                    slope: float = 0.15, threshold: float = 0.45,
                    scalar: float = 1.25, cell: float = 1.0) -> dict:
    """Ground-classify a denoised cloud with SMRF and record the parameters.

    Every value here is in the file's horizontal units. Running this on a
    cloud in EPSG:4326 silently produces nonsense, so the caller must have
    reprojected first.
    """
    params = {"window": window, "slope": slope, "threshold": threshold,
              "scalar": scalar, "cell": cell}
    pipeline = {"pipeline": [
        src,
        # assign everything to 'never classified' first so a re-run is
        # idempotent rather than layering on a previous classification
        {"type": "filters.assign", "value": "Classification = 0"},
        {"type": "filters.smrf", **params},
        {"type": "writers.las", "filename": dst, "compression": "laszip",
         "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

    Path(dst).with_suffix(".params.json").write_text(
        json.dumps({"filter": "smrf", **params}, indent=2))
    return params

filters.assign resetting the classification first is what makes the function safe to re-run. Without it a second pass over an already-classified file layers new decisions on top of old ones, and the result depends on how many times the pipeline has been run — a class of non-determinism that is miserable to debug because the file looks fine.

"forward": "all" on the writer carries the source header through: scale, offset, CRS and the vendor metadata. Omit it and PDAL writes defaults, which is the most common route to the precision loss described in fixing LAS scale and offset precision loss.

# Step 2: Use the colour before using the geometry

A photogrammetric point carries the RGB of the pixels that produced it, and vegetation is spectrally distinct from soil, gravel and concrete even in plain visible light. An excess-green index separates most of it in one arithmetic pass, and doing so before the morphological filter runs removes the class of points most likely to confuse it.

import json
import subprocess


EXCESS_GREEN = "(2.0 * Green - Red - Blue) / (Red + Green + Blue + 1)"


def pre_mask_vegetation(src: str, dst: str, *, egi_threshold: float = 0.08) -> None:
    """Mark strongly green points as high vegetation (class 5) before SMRF.

    The +1 in the denominator avoids a divide-by-zero on pure-black points,
    which a shadowed cloud has plenty of. The threshold is deliberately
    conservative: this stage should only catch points that are obviously
    vegetation, and leave the ambiguous ones to the geometric filter.
    """
    pipeline = {"pipeline": [
        src,
        {"type": "filters.ferry", "dimensions": "=>ExcessGreen"},
        {"type": "filters.assign",
         "value": [f"ExcessGreen = {EXCESS_GREEN}",
                   f"Classification = 5 WHERE ExcessGreen > {egi_threshold}"]},
        {"type": "writers.las", "filename": dst, "compression": "laszip",
         "forward": "all", "extra_dims": "ExcessGreen=float"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

Keeping ExcessGreen as an extra dimension rather than discarding it is worth the file size. When a classification result is questioned, being able to colour the cloud by the index that drove the decision turns an argument into an inspection. The threshold itself is site-dependent — dry grass in August is not green — and separating buildings from vegetation in point clouds covers how to calibrate it against a known patch.

# Step 3: Separate buildings from the rest of the non-ground

SMRF answers one question: ground or not. Splitting the “not” into vegetation, buildings and everything else needs a second pass, and the cheapest discriminator is local planarity. A roof is flat over several metres; a tree canopy is not flat at any scale.

import json
import subprocess


def classify_buildings(src: str, dst: str, *, knn: int = 24,
                       planarity_min: float = 0.92,
                       min_height: float = 2.0) -> None:
    """Split non-ground points into buildings (6) and vegetation (5).

    filters.covariancefeatures computes eigenvalue-derived shape descriptors
    per point; Planarity near 1 means the neighbourhood lies in a plane.
    Combined with a height above ground it is a serviceable roof detector.
    """
    pipeline = {"pipeline": [
        src,
        {"type": "filters.hag_nn"},                       # HeightAboveGround
        {"type": "filters.covariancefeatures", "knn": knn,
         "feature_set": "Dimensionality", "threads": 4},
        {"type": "filters.assign", "value": [
            f"Classification = 6 WHERE Classification != 2 "
            f"&& Planarity > {planarity_min} && HeightAboveGround > {min_height}",
            f"Classification = 5 WHERE Classification != 2 "
            f"&& Classification != 6 && HeightAboveGround > {min_height}",
            f"Classification = 3 WHERE Classification != 2 "
            f"&& Classification != 6 && HeightAboveGround <= {min_height} "
            f"&& HeightAboveGround > 0.3",
        ]},
        {"type": "writers.las", "filename": dst, "compression": "laszip",
         "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

filters.hag_nn must come first because every rule after it is expressed in height above ground rather than absolute elevation — the only frame in which “two metres up” means the same thing on a hillside as on a flat apron. The ordering constraint is easy to miss and produces a classification that is correct on flat sites and wrong on sloped ones, which is the worst kind of bug to find in production.

Figure 2 — Two cheap features do most of the work. The disagreements live in one band, and that band is where sampled validation should concentrate.

# Step 4: Run it over a survey without loading it

A site cloud is routinely 300 million points. Classification is a neighbourhood operation, so it cannot be done point-by-point in a stream, but it can be done tile-by-tile with a buffer — and the buffer is the part that is usually wrong.

import json
import subprocess
from pathlib import Path


def classify_tiled(src: str, out_dir: str, *, tile: float = 500.0,
                   buffer: float = 40.0, **smrf) -> None:
    """Tile, classify each tile with a halo, then trim the halo away.

    The buffer must exceed the SMRF window, or points near a tile edge are
    classified against a truncated neighbourhood and a grid of seams appears
    in the ground class exactly one tile apart.
    """
    if buffer <= smrf.get("window", 18.0):
        raise ValueError("buffer must exceed the morphological window")

    Path(out_dir).mkdir(parents=True, exist_ok=True)
    pipeline = {"pipeline": [
        src,
        {"type": "filters.splitter", "length": tile, "buffer": buffer},
        {"type": "filters.smrf", **smrf},
        # drop the halo so tiles can be concatenated without duplicates
        {"type": "filters.crop", "a_srs": "", "bounds": ""},
        {"type": "writers.las", "filename": f"{out_dir}/tile_#.laz",
         "compression": "laszip", "forward": "all"},
    ]}
    subprocess.run(["pdal", "pipeline", "--stdin"],
                   input=json.dumps(pipeline), text=True, check=True)

The visible symptom of an inadequate buffer is a regular grid of artefacts in the DTM at exactly the tile spacing. It is unmistakable once seen and completely invisible in any per-tile check, because each tile is internally consistent. The tiling discipline is the same one applied to imagery in structuring drone imagery for batch processing.

# Parameter deep-dive

The three that change the answer materially are window, threshold and scalar. The rest change the runtime or the edge cases. A useful discipline is to hold everything but those three fixed across a fleet of sites and record the three in the manifest, so a year of runs is comparable.

# Verification and output inspection

Classification is verified in three ways, all cheap, and a production pipeline should do all three on every run.

import json
import numpy as np
import pdal


def classification_report(path: str) -> dict:
    """Class histogram, ground coverage and height sanity, straight from the file."""
    pipe = pdal.Pipeline(json.dumps({"pipeline": [
        path, {"type": "filters.hag_nn"}]}))
    pipe.execute()
    arr = pipe.arrays[0]

    classes, counts = np.unique(arr["Classification"], return_counts=True)
    hist = {int(c): int(n) for c, n in zip(classes, counts)}
    total = int(arr.size)

    ground = arr[arr["Classification"] == 2]
    report = {
        "total_points": total,
        "class_histogram": hist,
        "ground_fraction": len(ground) / total if total else 0.0,
        # Ground points must sit at height-above-ground zero by construction;
        # any spread here means hag_nn and the classification disagree.
        "ground_hag_p99": float(np.percentile(ground["HeightAboveGround"], 99))
        if len(ground) else float("nan"),
        "max_hag": float(arr["HeightAboveGround"].max()) if total else float("nan"),
    }
    return report


def assert_classification_sane(report: dict, *, min_ground: float = 0.15,
                               max_feature_height: float = 60.0) -> None:
    """Deterministic gates a classified survey cloud must clear."""
    if report["ground_fraction"] < min_ground:
        raise ValueError(f"ground is only {report['ground_fraction']:.1%} of points")
    if report["max_hag"] > max_feature_height:
        raise ValueError(
            f"something stands {report['max_hag']:.1f} m above ground — "
            "almost certainly un-removed noise, not a structure")
    if not (0 in report["class_histogram"] or 1 in report["class_histogram"]) \
            and len(report["class_histogram"]) < 2:
        raise ValueError("everything landed in one class — the filter did not run")

The max_hag check is the one that earns its place most often. Surviving outliers hundreds of metres above the site are common in photogrammetric clouds, they defeat every subsequent auto-scaled visualisation, and they are a single comparison away from being caught.

The third verification is the sampled residual against measured ground described in the point cloud processing section overview. Run it whenever checkpoints exist — which, on any project with ground control point optimization, they do.

Figure 3 — Recognising which of the three you have is most of the fix.

# Troubleshooting

The building roof is in my terrain model. The window is smaller than the building. Set window above the widest structure on site and re-run; no threshold adjustment will fix this, because at that scale the roof genuinely is the lowest surface across the window.

Steep slopes are being classified as non-ground. slope is too low for the terrain. It is a rise-over-run tolerance, so a 30° hillside needs at least 0.58. Raising scalar alongside it lets the threshold grow with the local gradient instead of applying a flat-ground value on a bank.

The ground class follows the top of the grass. threshold is above the vegetation height, so low growth never exceeds it. Lower it toward 0.2 m and pre-mask with the excess-green filter. Confirm with sampled residuals rather than by eye — this is a bias, and biases are invisible in renderings.

A regular grid of ridges appears in the DTM. The tile buffer is at or below the morphological window, so edge points were classified against a truncated neighbourhood. Set the buffer to at least twice the window and re-run the tiled job.

Everything is class 1, unclassified. The pipeline ran but the filter did not: either filters.smrf was never reached because an earlier stage emptied the cloud, or a filters.range upstream excluded everything. Print the point count between stages; a stage that outputs zero points is silent otherwise.

Classification differs between two runs on the same file. The classification was not reset before filtering, so the second run built on the first. Add filters.assign setting Classification = 0 at the head of the pipeline.

# Related

• Filtering Noise and Outliers from Dense Clouds
• Separating Buildings from Vegetation in Point Clouds
• Fixing Misclassified Ground Under Dense Canopy
• Validating Classification Against Manual Samples
• PDAL vs laspy for Point Cloud Automation

← Point Cloud Processing & 3D Deliverables