Coastal & Marine Spatial Analysis Pipelines

Using PDAL for Bathymetric Point Cloud Cleaning

Raw multibeam echosounder (MBES) datasets routinely contain water-column backscatter, surface multipath returns, acoustic sidelobe artifacts, and vessel-wake noise. Manual GUI-based culling fails at survey scale and breaks reproducibility — the same dataset cleaned twice by hand yields two different surfaces. This page shows how to use PDAL (the Point Data Abstraction Library) to express the entire cleaning sequence as a versioned JSON pipeline that runs identically on a laptop, a CI runner, or an edge processing node. It is a hands-on companion to the broader multibeam point-cloud filtering workflow that it sits under; read that parent page first if you need the full filter-selection rationale and downstream handoff contract.

Root Cause: Why Raw MBES Returns Need Programmatic Cleaning

A multibeam swath records every acoustic return the sonar resolves, not just the seafloor. Within a single ping you get genuine bottom detections interleaved with refraction-bent outer-beam returns, fish and bubble scatter in the water column, and sidelobe energy that maps to a false depth. Two structural problems make this hard to fix interactively:

  1. Scale. A single survey line can hold tens of millions of soundings; a campaign holds billions. There is no manual cursor fast enough.
  2. Non-determinism. GUI editing leaves no machine-readable record of which points were removed or why, so the cleaning cannot be re-run or audited.

You can reproduce the failure that motivates this work in one command — inspect a raw tile and watch the Z range blow far past any plausible depth corridor:

# Raw survey tile: note the Z minimum is well below the true seabed,
# and the maximum sits above the waterline — both are noise.
pdal info --stats raw_survey.laz | python3 -c "import sys,json; s=json.load(sys.stdin)['stats']['statistic']; z=[d for d in s if d['name']=='Z'][0]; print('Z min/max:', z['minimum'], z['maximum'])"
# Z min/max: -812.4 14.7   <-- seabed corridor is ~ -120..0, the rest is noise

PDAL solves both problems by executing an ordered list of filters as a directed acyclic graph: each stage reads the previous stage’s points, transforms or tags them, and passes them on. Because the pipeline is JSON, it is reviewable, diffable, and committable alongside the survey metadata.

PDAL MBES Cleaning Pipeline as a Directed Acyclic Graph Left-to-right data-flow diagram of seven PDAL stages. readers.las (input_survey.laz) feeds filters.range with limit Z[-120:0], which drops out-of-corridor blunders. filters.outlier (statistical, mean_k=24) tags noise as Classification 7. A second filters.range with Classification![7:7] drops those tagged points. filters.assign sets surviving class-1 points to Ground (class 2). filters.splitter tiles into 400 m chunks, and writers.las writes laszip-compressed LAZ with forward all. Two downward dashed arrows mark where rejected points leave the graph. readers.las input_survey.laz filters.range Z[-120:0] clip filters.outlier tag noise → class 7 filters.range drop Class![7:7] filters.assign class 1 → Ground 2 filters.splitter 400 m tiles writers.las laszip, forward:all blunders dropped class-7 noise dropped rejected points leave the graph — never reach the writer JSON pipeline executes as a directed acyclic graph — reviewable, diffable, committable

Step-by-Step Fix: A Production Cleaning Pipeline

Step 1 — Define the filter sequence as JSON

The order matters. Clip to a valid depth corridor first so the outlier statistics are not skewed by gross blunders, tag and drop statistical outliers, reclassify survivors to ASPRS Ground (class 2), then tile for parallel downstream work and write compressed LAZ.

{
  "pipeline": [
    {
      "type": "readers.las",
      "filename": "input_survey.laz"
    },
    {
      "type": "filters.range",
      "limits": "Z[-120:0], Classification[1:9]"
    },
    {
      "type": "filters.outlier",
      "method": "statistical",
      "mean_k": 24,
      "multiplier": 2.0
    },
    {
      "type": "filters.range",
      "limits": "Classification![7:7]"
    },
    {
      "type": "filters.assign",
      "assignment": "Classification[:]=2",
      "condition": "Classification==1"
    },
    {
      "type": "filters.splitter",
      "length": 400,
      "origin_x": "auto",
      "origin_y": "auto"
    },
    {
      "type": "writers.las",
      "filename": "cleaned_#.laz",
      "compression": "laszip",
      "forward": "all"
    }
  ]
}

What each stage does:

  • filters.range with Z[-120:0] clips to a valid depth corridor. Adjust the bounds to local bathymetry and survey datum; PDAL treats positive Z as up, so seabed returns are negative.
  • filters.outlier (method: statistical) runs a k-nearest-neighbour density estimator and tags outliers as Classification 7. For dense MBES (>50 pts/m²), raise mean_k to 32–48 so steep-slope returns are not over-filtered.
  • The second filters.range with Classification![7:7] drops the points the outlier stage tagged.
  • filters.assign reclassifies remaining unclassified points (class 1) to Ground (class 2); the condition expression guards which points are touched.
  • filters.splitter tiles into 400 m chunks so downstream gridding can run in parallel; # in the output filename is the tile-index placeholder.

Step 2 — Enforce CRS and vertical datum before filtering

Bathymetric pipelines fail silently when the horizontal and vertical reference systems diverge, because filters.range clips on raw Z and PDAL does not auto-resolve datum shifts. Insert an explicit reprojection ahead of the range clip. This is the same discipline the parent cluster applies through its dedicated CRS alignment stage and tidal datum transformation handling — apply it before any depth-dependent filter runs.

{
  "type": "filters.reprojection",
  "in_srs": "EPSG:4326+EPSG:5703",
  "out_srs": "EPSG:32618+EPSG:5703"
}

Validate vertical datums against local tidal benchmarks (MLLW, NAVD88, or LAT). If the provider delivers raw ellipsoidal heights, apply a geoid model via filters.reprojection with the appropriate grid. Normalising Z before filters.range is non-negotiable; skip it and you guarantee datum-induced clipping or retained artifacts.

Step 3 — Run in streaming mode to bound memory

MBES surveys routinely exceed 40–80 GB and PDAL loads many filter chains fully in memory, triggering OOM kills on standard CI runners. Force streaming and cap concurrency:

# --stream avoids in-memory point buffering; --workers caps PDAL worker threads.
pdal pipeline pipeline.json --stream --workers 4

If resident memory exceeds ~75% of available RAM, shrink the filters.splitter length to 200 and confirm --stream is set. Note that some filters (filters.outlier among them) are not streamable and will force a fallback to standard mode for that stage — check the run log for a streaming-mode warning when memory spikes unexpectedly.

Step 4 — Drive the pipeline from Python with structured logging

Embed PDAL in a Python spatial stack through its bindings, wrapping execution in deterministic error handling. This is the same pattern used when converting LAS to XYZ for legacy gridding tools downstream.

import json
import logging

import pdal

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


def run_bathy_cleaner(pipeline_dict: dict) -> int:
    """Execute a PDAL cleaning pipeline and return the surviving point count.

    Raises:
        RuntimeError: if PDAL fails to execute the pipeline.
    """
    pipeline = pdal.Pipeline(json.dumps(pipeline_dict))
    try:
        count: int = pipeline.execute()
    except RuntimeError as exc:
        logging.error("PDAL execution failed: %s", exc)
        raise

    logging.info("PDAL executed. Points retained: %d", count)
    # Tile-level acceptance (ground ratio, density) is validated in Step 6 by
    # parsing `pdal info --stats` on each written tile, not from runtime metadata.
    return count

Verification: Confirm the Cleaning Actually Worked

A cleaning run is only trustworthy if a machine check confirms it. Pull summary statistics from a written tile and assert the Z range now sits inside the survey corridor:

pdal info --stats cleaned_0.laz
import json
import subprocess

result = subprocess.run(
    ["pdal", "info", "--stats", "cleaned_0.laz"],
    capture_output=True, text=True, check=True,
)
stats = {s["name"]: s for s in json.loads(result.stdout)["stats"]["statistic"]}
z_min, z_max = stats["Z"]["minimum"], stats["Z"]["maximum"]

assert -122.0 <= z_min and z_max <= 2.0, f"Z out of corridor after cleaning: {z_min}..{z_max}"
print(f"PASS — Z range {z_min:.2f}..{z_max:.2f} within corridor")

For CI/CD gating, enforce these thresholds per tile and reject failures with their bounding coordinates logged for manual review:

Gate Shallow water Deep water
Minimum point density ≥ 15 pts/m² ≥ 5 pts/m²
Ground (class 2) share ≥ 85% of retained points ≥ 85%
Z-range compliance within survey bounds ± 2 m within survey bounds ± 2 m

Rejecting tiles at this gate stops garbage propagating into DEM generation — the next stage handled by the DEM interpolation workflow.

Edge Cases and Gotchas

  • Steep slopes read as outliers. On reefs, canyon walls, and wreck flanks, legitimate returns are sparse relative to the flat-bottom assumption baked into filters.outlier. Raise mean_k and prefer multiplier ≥ 2.5, or swap to method: radius with an explicit radius so density is judged geometrically rather than by fixed neighbour count.
  • forward: all is required to preserve provenance. Without it writers.las resets scale, offset, and VLRs, silently dropping the vertical CRS and breaking any downstream datum check. Keep forward: all and verify the output SRS with pdal info --metadata.
  • LAS version and PDAL version drift. Class codes 8–18 and per-point extra dimensions only round-trip cleanly in LAS 1.4; older 1.2 writers truncate them. Pin LAS 1.4 and PDAL ≥ 2.6.0 (the first reliably streamable line) in your environment.yml or Dockerfile so the pipeline behaves the same in CI as on the survey workstation.
  • filters.splitter boundary slivers. Tiling on a fixed origin can split a single feature across tiles; if you grid each tile independently, add a small overlap buffer or merge before gridding to avoid seam artifacts — the same boundary-continuity concern that the automated spike removal page handles with chunk overlap.

Production Deployment Checklist

Reference the PDAL pipeline documentation for schema validation and the NOAA VDatum technical documentation to confirm geoid offsets match your survey epoch.

Up: Point Cloud Filtering for Multibeam Sonar