Point Cloud Processing and Classification for Archaeological Site Recording #

A dense photogrammetric or lidar point cloud straight out of reconstruction is not an archaeological record — it is raw evidence contaminated with reconstruction noise, vegetation returns, spurious floating points, and coordinate ambiguity. Before a cloud can support a terrain model, a stratigraphic surface, or season-over-season change analysis, it has to be filtered, classified, and locked to a known datum. This topic area sits inside the wider Photogrammetry & 3D Site Mapping Pipelines framework and covers the deterministic, scriptable middle of that pipeline: taking the dense cloud produced by automated drone image processing workflows and turning it into a clean, classified, archive-ready product that feeds surface reconstruction in mesh generation and optimization for ruins.

The three guides in this section address the operations that most often go wrong in practice: filtering vegetation and noise from site point clouds, ground classification with PDAL pipelines, and registering multitemporal scans for change detection. Everything here is built on PDAL pipelines expressed as JSON and driven from Python, because reproducibility — not a one-off click in a desktop viewer — is what makes a heritage dataset defensible years later.

Prerequisites #

The whole section standardises on PDAL and its Python bindings so that pipelines can be version-controlled and run headless on a compute node. Pin these versions in an isolated environment; PDAL filter names and stage options change between minor releases, and a heritage archive needs to reproduce a result on demand.

Environment assumptions: a POSIX shell, roughly 8 GB RAM per 100M-point tile, and write access to a projected working directory. Install with conda install -c conda-forge pdal python-pdal gdal to keep the C++ engine and the Python binding on matching ABI versions — mixing a pip pdal against a system libpdal is the single most common cause of RuntimeError: Unable to convert stage output.

Implementation #

The field-to-archive order below is deliberate: never classify before you have removed gross noise, and never derive a terrain surface before you have reprojected to your archive datum.

1. Inspect and validate the incoming cloud #

Before touching a single point, read the header. pdal info reports the spatial reference, point count, and bounding box, which is where you catch a cloud that Metashape exported in a local engineering frame rather than a real EPSG.

# PDAL 2.7.2
pdal info data/site_dense.laz --summary
pdal info data/site_dense.laz --metadata | grep -i srs

2. Strip statistical outliers #

Photogrammetric clouds carry sparse floating points — birds, dust, reconstruction ghosts above the site. filters.outlier in statistical mode flags points whose mean neighbour distance exceeds a multiple of the global standard deviation.

{
  "pipeline": [
    "data/site_dense.laz",
    { "type": "filters.outlier", "method": "statistical", "mean_k": 12, "multiplier": 2.5 },
    { "type": "filters.range", "limits": "Classification![7:7]" },
    "data/site_denoised.laz"
  ]
}

The filters.range stage keeps everything the outlier filter did not mark as class 7 (low noise). This two-step idiom recurs throughout filtering vegetation and noise from site point clouds, where it is combined with height-above-ground gating to remove canopy.

3. Reproject to the archive datum #

Reproject once, early, so every subsequent metric is expressed in metres on the British National Grid.

{
  "pipeline": [
    "data/site_denoised.laz",
    { "type": "filters.reprojection", "in_srs": "EPSG:32633", "out_srs": "EPSG:27700" },
    "data/site_bng.laz"
  ]
}

EPSG:27700 is British National Grid; swap in your site EPSG. If the source cloud has no SRS in its header, set in_srs explicitly — PDAL will not guess.

4. Classify ground versus non-ground #

Ground extraction is the pivot of the whole section. PDAL ships three algorithms: Simple Morphological Filter (filters.smrf), Progressive Morphological Filter (filters.pmf), and Cloth Simulation Filter (filters.csf). SMRF is the workhorse for open, gently undulating heritage landscapes; CSF handles steep, terraced, or heavily built sites more gracefully. The full parameterisation lives in ground classification with PDAL pipelines.

{
  "pipeline": [
    "data/site_bng.laz",
    { "type": "filters.smrf", "slope": 0.15, "window": 18.0, "threshold": 0.45, "scalar": 1.2 },
    "data/site_classified.laz"
  ]
}

5. Decimate and tile for delivery #

Archive and web-delivery products rarely need full density. Thin with filters.decimation for lightweight derivatives, and tile large sites so downstream tools stream rather than load whole.

{
  "pipeline": [
    "data/site_classified.laz",
    { "type": "filters.decimation", "step": 4 },
    { "type": "writers.copc", "filename": "archive/site_thinned.copc.laz" }
  ]
}

Cloud-Optimized Point Cloud (COPC) is a single LAZ with an internal octree, so a web viewer or a spatial query can fetch only the octree nodes it needs over HTTP range requests — the point-cloud equivalent of a Cloud-Optimized GeoTIFF.

6. Derive the terrain surface #

The last step rasterises the classified ground points into a bare-earth DTM with writers.gdal, and (optionally) all first returns into a DSM. The difference between the two surfaces is a canopy/structure height model.

{
  "pipeline": [
    "data/site_classified.laz",
    { "type": "filters.range", "limits": "Classification[2:2]" },
    { "type": "writers.gdal", "filename": "archive/site_dtm.tif", "resolution": 0.1, "output_type": "idw", "gdaldriver": "GTiff" }
  ]
}

CRS & Geometry Validation #

Point clouds fail quietly. A cloud in the wrong datum still opens, still renders, still classifies — it is simply in the wrong place by tens or hundreds of metres. Guard against that with three checks: header SRS presence, bounding-box plausibility, and an inverse-transform round trip.

A projected cloud on EPSG:27700 should have easting/northing values roughly in the ranges $0 \le E \le 700{,}000$ and $0 \le N \le 1{,}300{,}000$ metres. If the mean height is negative thousands or the eastings are single digits, the header SRS is wrong. Cross-check by inverse-transforming the centroid to EPSG:4326 and confirming it lands inside your survey region:

# requirements.txt
# pdal==3.4.5
# pyproj==3.6.1
# numpy==1.26.4

import pdal, json
from pyproj import Transformer

pipe = pdal.Pipeline(json.dumps({"pipeline": ["data/site_bng.laz"]}))
pipe.execute()
xyz = pipe.arrays[0]
cx, cy = float(xyz["X"].mean()), float(xyz["Y"].mean())

# EPSG:27700 British National Grid -> WGS84  # substitute your site's EPSG
lon, lat = Transformer.from_crs("EPSG:27700", "EPSG:4326", always_xy=True).transform(cx, cy)
assert -8.6 <= lon <= 1.8 and 49.9 <= lat <= 60.9, f"Centroid {lon:.3f},{lat:.3f} outside GB"

Common mismatches: a Metashape export tagged as LOCAL_CS (no EPSG), a cloud reprojected twice because in_srs was set on an already-projected file, and a vertical datum confusion where ellipsoidal heights are treated as orthometric — the last drifts the DTM by the geoid separation (around 45–50 m over Great Britain).

Automated QA/QC #

Gate every pipeline on point counts and classification proportions, and route failures to a rejects folder rather than the archive. A sound gate asserts that (a) the denoised cloud retains at least 95% of input points — a bigger loss means the outlier multiplier is too aggressive — and (b) ground points make up a plausible 20–70% of a classified open-air site. Emit a structured audit record for every run so the archive can prove provenance.

import pdal, json, datetime

def qa_record(path):
    p = pdal.Pipeline(json.dumps({"pipeline": [path]}))
    n = p.execute()
    arr = p.arrays[0]
    ground = int((arr["Classification"] == 2).sum())
    return {
        "file": path,
        "timestamp": datetime.datetime.utcnow().isoformat() + "Z",
        "total_points": int(n),
        "ground_points": ground,
        "ground_fraction": round(ground / n, 4) if n else 0.0,
        "srs": p.metadata["metadata"]["readers.las"]["srs"]["horizontal"][:60],
    }

rec = qa_record("data/site_classified.laz")
assert 0.2 <= rec["ground_fraction"] <= 0.7, f"Implausible ground fraction {rec['ground_fraction']}"

Promotion criterion: a cloud reaches the archive only when its audit record shows a resolved horizontal SRS, a ground fraction inside the plausible band, and a bounding box inside the site polygon. Everything else is quarantined with its failure reason logged.

Integration with adjacent workflows #

The classified cloud is a hub, not an endpoint. The bare-earth ground class feeds surface reconstruction in mesh generation and optimization for ruins, where a clean DTM makes screened-Poisson meshing dramatically more stable. The georeferencing established here is the same National Grid frame that automated drone image processing workflows target when aligning clouds to ground control. And the tiled COPC and derived DTM/DSM rasters are the direct inputs to orthomosaic and terrain-tile delivery further along the pipeline. Because every product carries an explicit EPSG and an audit record, handoff between these stages never requires a manual reprojection guess.

Compliance notes #

Archived point-cloud products carry metadata under ISO 19115 and, for UK deposition, MIDAS Heritage. Record at minimum: the horizontal and vertical CRS (EPSG codes, not free text), the acquisition date and sensor, the processing lineage (PDAL version and the exact pipeline JSON), the point density, and the classification scheme used (ASPRS class codes — 2 for ground, 3–5 for vegetation, 7 for low noise). For FAIR compliance the pipeline JSON should be deposited alongside the cloud so the classification is reproducible, and the COPC product should be given a persistent identifier so it can be cited. Lineage that names the algorithm and its parameters is what separates an archival dataset from an undocumented blob.

Related #

• Filtering Vegetation and Noise from Site Point Clouds — statistical outlier removal plus height gating to strip canopy returns.
• Ground Classification with PDAL Pipelines — tuning SMRF and exporting a bare-earth DTM with writers.gdal.
• Registering Multitemporal Scans for Change Detection — ICP alignment and M3C2-style distance for season-over-season erosion.
• Automated Drone Image Processing Workflows — the upstream reconstruction that produces the dense cloud.
• Mesh Generation & Optimization for Ruins — turning the classified ground surface into an analysis-ready mesh.

Part of Photogrammetry & 3D Site Mapping Pipelines.