Ground Classification with PDAL Pipelines #

Separating ground from everything else is the single most consequential decision in point-cloud processing for archaeology. Every bare-earth model, every microtopographic hillshade that reveals a plough-levelled enclosure, every volume calculation depends on which returns you label as class 2. Get the ground surface wrong and a subtle earthwork either disappears into the vegetation you deleted or bloats into a mound that was never there. This guide tunes PDAL’s Simple Morphological Filter (filters.smrf) for heritage terrain and exports a bare-earth Digital Terrain Model with writers.gdal, extending the cleanup covered in point cloud processing and classification.

SMRF works by simulating a morphological opening at progressively larger window sizes, treating a point as ground when it stays within a slope-dependent elevation threshold of the opened surface. Four parameters govern it — slope, window, threshold, and scalar — and understanding what each does is the difference between a clean DTM and a smeared one.

Context & When to Use #

SMRF is the default choice for the open, rolling, and terraced landscapes that dominate heritage survey: barrows, hillforts, field systems, deserted settlements. It is fast, its parameters map onto physical terrain properties, and it degrades gracefully. Use it whenever the site is broadly a surface with objects sitting on it.

Its parameters trade off directly against each other. The window sets the diameter of the largest non-ground object the filter will remove — set it to comfortably exceed your largest building footprint or spoil heap, but no larger, because an oversized window starts eating genuine terrain relief. The slope parameter is the maximum ground gradient the filter tolerates before it starts rejecting steep points; a low value on a hillfort rampart will strip the rampart. The threshold is the vertical distance a point may sit above the provisional surface and still count as ground, and scalar scales that threshold with slope so that on steep terrain the tolerance widens and terrace edges are preserved. The performance cost is modest — SMRF is near-linear in point count — so the real cost of getting it wrong is archaeological, not computational: a mis-tuned filter silently rewrites the microtopography you are trying to record. When SMRF cannot cope with dense canopy, fall back to the Cloth Simulation Filter, which is more forgiving under trees but slower; the denoising that precedes either is covered in filtering vegetation and noise from site point clouds.

Implementation #

The pipeline classifies ground, writes the ASPRS classification codes into the LAZ, and rasterises the ground-only points into an inverse-distance-weighted DTM. The Python runner exposes the four SMRF parameters as arguments so a site-specific configuration can be version-controlled.

# requirements.txt
# pdal==3.4.5
# python-pdal==3.4.5
# gdal==3.8.4

import json
import logging
import pdal

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
log = logging.getLogger("ground_classify")

def ground_pipeline(src: str, laz_out: str, dtm_out: str,
                    slope: float = 0.15, window: float = 18.0,
                    threshold: float = 0.45, scalar: float = 1.2,
                    resolution: float = 0.1) -> dict:
    """Classify ground with SMRF, persist class codes, and rasterise a bare-earth DTM.

    slope     : max ground gradient (rise/run) tolerated. Lower = stricter on slopes.
    window    : diameter (m) of the largest object removed. Exceed your largest structure.
    threshold : vertical tolerance (m) for a point to count as ground.
    scalar    : widens 'threshold' proportionally to slope so terraces survive.
    resolution: DTM cell size (m).
    """
    return {
        "pipeline": [
            src,
            # Write ground into ASPRS Classification=2; everything else stays untouched.
            {"type": "filters.smrf",
             "slope": slope, "window": window,
             "threshold": threshold, "scalar": scalar},
            # Persist the classified cloud (ground = 2) for the archive.
            {"type": "writers.las", "filename": laz_out,
             "a_srs": "EPSG:27700", "compression": "laszip"},  # substitute your site's EPSG
            # Keep only ground for rasterisation.
            {"type": "filters.range", "limits": "Classification[2:2]"},
            # IDW bare-earth DTM. window_size fills small gaps between ground points.
            {"type": "writers.gdal", "filename": dtm_out,
             "resolution": resolution, "output_type": "idw",
             "window_size": 3, "gdaldriver": "GTiff",
             "override_srs": "EPSG:27700"},  # substitute your site's EPSG
        ]
    }

def run(src, laz_out, dtm_out, **kw) -> dict:
    p = pdal.Pipeline(json.dumps(ground_pipeline(src, laz_out, dtm_out, **kw)))
    p.execute()
    arr = p.arrays[0]  # arrays reflect the last stage's view (ground only)
    log.info("Ground points rasterised: %d at %.2f m resolution",
             arr.size, kw.get("resolution", 0.1))
    return {"ground_points": int(arr.size)}

if __name__ == "__main__":
    run("data/site_ground_candidate.laz",
        "archive/site_classified.laz",
        "archive/site_dtm.tif",
        slope=0.15, window=18.0, threshold=0.45, scalar=1.2, resolution=0.1)

The relationship between the provisional morphological surface and the ground decision can be written as: a point at elevation $z$ over a provisional surface value $s$ is ground when $z - s \le t + k\,\theta$, where $t$ is threshold, $k$ is scalar, and $\theta$ is the local slope. That is why raising scalar rescues steep terrace edges without loosening the filter on flat ground.

Verification #

Check that the ground fraction is plausible and that the DTM is continuous with no unfilled NoData holes across the site interior.

import pdal, json

# Ground fraction from the classified LAZ
p = pdal.Pipeline(json.dumps({"pipeline": ["archive/site_classified.laz"]}))
n = p.execute()
cls = p.arrays[0]["Classification"]
ground = int((cls == 2).sum())
frac = ground / n
print(f"total={n}  ground={ground}  fraction={frac:.3f}")
assert 0.2 <= frac <= 0.75, f"Ground fraction {frac:.3f} implausible — retune SMRF"

Then confirm the raster is sound with GDAL:

# GDAL 3.8.4
gdalinfo -stats archive/site_dtm.tif | grep -E "Size is|NoData|Minimum|Maximum"

Expected output resembles:

Size is 3120, 2680
  NoData Value=-9999
  Minimum=41.223, Maximum=58.907

A minimum near -9999 leaking into the statistics means NoData is being counted as elevation — a sign of large unfilled gaps; raise window_size in writers.gdal or lower the DTM resolution so cells capture enough ground points.

Common Errors & Fixes #

filters.smrf: Point cloud is too small to process — the input view has fewer points than the window can span, usually because an upstream filter removed almost everything or the tile is a sliver at a survey edge. Confirm the denoising stage retained a sane point count, and merge edge slivers into their neighbours before classifying.

writers.gdal: Grid width and height...exceed... maximum — the combination of site extent and a fine resolution produced a raster larger than GDAL will allocate. Either coarsen resolution (0.25 m is ample for landscape-scale DTMs) or tile the cloud first and rasterise per tile, then build a VRT mosaic with gdalbuildvrt.

Bare-earth DTM shows ridged, corduroy artefacts along flight lines — not an error message but a classic symptom of scalar set too low relative to slope, so ground points on gentle undulations are rejected in stripes. Raise scalar toward 1.5 and re-inspect the hillshade. This same striping, if left in a DTM used for registering multitemporal scans for change detection, masquerades as real surface change between epochs.

Related #

• Point Cloud Processing and Classification — the section overview placing ground classification in the full workflow.
• Filtering Vegetation and Noise from Site Point Clouds — the denoising that must precede a reliable SMRF pass.
• Registering Multitemporal Scans for Change Detection — comparing bare-earth surfaces across field seasons.
• Mesh Generation & Optimization for Ruins — using the classified ground surface as a reconstruction input.

Part of Photogrammetry & 3D Site Mapping Pipelines.