Filtering Vegetation and Noise from Site Point Clouds #

A summer survey of an earthwork or a lynchet almost always captures more grass, bracken, and scrub than archaeology. The dense cloud that comes back from reconstruction is a mixture of the surface you want and a green haze of vegetation returns sitting centimetres to metres above it, plus a scatter of gross outliers — insects caught mid-flight, reflective wet stone, reconstruction ghosts floating in empty air. Left in place, that contamination corrupts every terrain surface and every volume you later compute. This guide covers the two-stage cleanup that turns a raw cloud into a defensible ground candidate, as part of the wider work on point cloud processing and classification.

The distinction that matters is between noise — isolated points with no real-world referent — and vegetation — dense, coherent returns from something that is genuinely there but is not the archaeological surface. They need different tools. Noise is removed statistically; vegetation is removed by reasoning about height above the ground.

Context & When to Use #

Reach for statistical outlier removal on every cloud without exception — it is cheap and the failure mode of skipping it (a single stray point kilometres away) wrecks bounding-box-driven downstream tools. The decision that needs judgement is how aggressively to strip vegetation.

There are two viable strategies. The first, used here, is a height-above-ground gate: run a ground filter to estimate the terrain, compute each point’s height above that surface with filters.hag_nn, then drop everything above a threshold. This is fast, deterministic, and works well on grassland and low scrub. The trade-off is that a low threshold on a genuinely rough masonry surface will clip real stone; a high threshold leaves tall grass behind. The second strategy, the Cloth Simulation Filter, is more robust under dense canopy but costs more compute and is covered in ground classification with PDAL pipelines. For most open-air heritage landscapes the height-gate approach in this guide is the right first pass; escalate to CSF only when the ground filter visibly fails under trees.

The compliance cost of over-filtering is real: strip too hard and you delete the low upstanding archaeology you were surveying. Always keep the unfiltered cloud as the master and treat the filtered product as a derivative.

Implementation #

The pipeline below chains four PDAL stages into one pass: statistical outlier detection, removal of the flagged noise, a SMRF ground pass used purely to seed a height model, and a height-above-ground range gate. It is driven from Python so the point counts before and after each stage can be captured for the audit log.

# requirements.txt
# pdal==3.4.5
# python-pdal==3.4.5
# numpy==1.26.4

import json
import logging
import pdal

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

def build_pipeline(src: str, dst: str, hag_max: float = 0.5) -> dict:
    """Two-stage cleanup: statistical noise removal + height-above-ground vegetation gate.

    hag_max is the maximum height (metres) a point may sit above the local ground
    before it is treated as vegetation. 0.5 m suits grazed grassland; raise to ~1.0 m
    over rough masonry so upstanding stone is not clipped.
    """
    return {
        "pipeline": [
            src,
            # 1. Flag isolated points as ASPRS class 7 (low noise). mean_k neighbours,
            #    anything beyond multiplier * std of mean distance is an outlier.
            {"type": "filters.outlier", "method": "statistical",
             "mean_k": 12, "multiplier": 2.5},
            # 2. Physically drop the flagged noise (keep everything NOT in [7:7]).
            {"type": "filters.range", "limits": "Classification![7:7]"},
            # 3. Seed a ground estimate so height-above-ground can be computed.
            {"type": "filters.smrf", "slope": 0.15, "window": 18.0,
             "threshold": 0.45, "scalar": 1.2},
            # 4. Nearest-neighbour height above the ground class, written to HeightAboveGround.
            {"type": "filters.hag_nn", "count": 6},
            # 5. Keep only points at or below hag_max -> vegetation and canopy fall away.
            {"type": "filters.expression",
             "expression": f"HeightAboveGround <= {hag_max}"},
            {"type": "writers.las", "filename": dst,
             "a_srs": "EPSG:27700", "compression": "laszip"},  # substitute your site's EPSG
        ]
    }

def run(src: str, dst: str, hag_max: float = 0.5) -> dict:
    p = pdal.Pipeline(json.dumps(build_pipeline(src, dst, hag_max)))
    n_out = p.execute()
    # count_get gives per-stage metadata; total input is on the reader
    meta = p.metadata["metadata"]
    n_in = meta["readers.las"]["count"]
    log.info("Input %d points -> output %d points (%.1f%% retained)",
             n_in, n_out, 100.0 * n_out / n_in)
    return {"points_in": n_in, "points_out": n_out}

if __name__ == "__main__":
    run("data/site_dense.laz", "data/site_ground_candidate.laz", hag_max=0.5)

The filters.expression stage (PDAL 2.7+) replaces the older filters.range string syntax for the height gate and reads far more clearly for anyone auditing the pipeline later. If you are pinned to an older PDAL, substitute {"type": "filters.range", "limits": "HeightAboveGround[:0.5]"}.

Verification #

Confirm the cleanup worked by comparing point counts and inspecting the height distribution that remains. A correct filter removes a small single-digit percentage as noise and a larger, site-dependent fraction as vegetation, and leaves a residual cloud whose HeightAboveGround is tightly clustered near zero.

import pdal, json, numpy as np

p = pdal.Pipeline(json.dumps({"pipeline": [
    "data/site_ground_candidate.laz",
    {"type": "filters.hag_nn", "count": 6},
]}))
p.execute()
hag = p.arrays[0]["HeightAboveGround"]

print(f"points: {hag.size}")
print(f"HAG max: {hag.max():.3f} m   HAG 99th pct: {np.percentile(hag, 99):.3f} m")
assert hag.max() <= 0.51, "Vegetation gate leaked points above the threshold"
assert np.percentile(hag, 99) < 0.4, "Residual canopy suspected — tighten hag_max"

Expected output resembles:

points: 4193776
HAG max: 0.499 m   HAG 99th pct: 0.271 m

If the retained count is a tiny fraction of the input, the ground seed failed and the height gate measured against a bad surface — re-check the SMRF stage before trusting the result.

Common Errors & Fixes #

filters.hag_nn: input PointView does not have any points classified as ground — the SMRF stage found no ground, so there is nothing to measure height against. Root cause is almost always that SMRF ran on a cloud still full of noise, or the window is far too small for the point spacing. Fix by confirming the outlier stage ran first and enlarging window (try 24.0) so the morphological opening spans real terrain.

PDAL: filters.expression: expression references dimension 'HeightAboveGround' which does not exist — the expression stage ran before filters.hag_nn created the dimension, or the two stages were transposed in the JSON. PDAL pipelines execute in array order, so filters.hag_nn must appear before any stage that reads HeightAboveGround.

readers.las: Global encoding WKT bit is not set... ignoring SRS with points landing far from site — the source LAZ has no coordinate system and PDAL declined to guess. Set the reader’s SRS explicitly with {"type": "readers.las", "filename": "...", "override_srs": "EPSG:32633"} so the reprojection and height computation operate in real metres. This mirrors the datum-reconciliation discipline described for registering multitemporal scans for change detection, where an unset SRS makes two epochs impossible to align.

Related #

• Point Cloud Processing and Classification — the section overview and where this cleanup sits in the full pipeline.
• Ground Classification with PDAL Pipelines — deeper SMRF and CSF tuning once the cloud is denoised.
• Registering Multitemporal Scans for Change Detection — why consistent filtering matters when comparing epochs.
• Mesh Generation & Optimization for Ruins — the surface reconstruction a clean ground candidate feeds.

Part of Photogrammetry & 3D Site Mapping Pipelines.