Point Cloud Filtering for Multibeam Sonar #

Multibeam echosounder (MBES) surveys generate dense, high-dimensional acoustic returns contaminated by water-column noise, multipath artifacts, and biological backscatter. Isolating valid seabed returns requires deterministic filtering before any terrain modeling step. Within the Bathymetric Processing & Terrain Modeling pipeline, this workflow transforms raw acoustic point clouds into georeferenced, validated LAZ files optimized for gridding and interpolation. The operational objective is statistical outlier rejection, geometric consistency checks, and depth-range validation while maintaining sub-gigabyte memory footprints across multi-terabyte survey blocks.

---

Reference Configuration and Specification Table #

---

Filtering Pipeline Architecture #

The four-stage sequence below shows how raw MBES returns are progressively reduced to validated seabed points. Each stage produces a deterministic intermediate state; the pipeline is designed to restart from any checkpoint without reprocessing earlier stages.

---

Data Ingestion and CRS Enforcement #

Raw MBES data originates in vendor-specific formats (Kongsberg .all, Teledyne .kmall, R2Sonic .s7k) or standardized exchange formats (LAS/LAZ, ASCII XYZ). Production ingestion mandates conversion to LAZ with embedded EPSG codes and explicit vertical datum metadata before any filter stage executes. Coordinate-agnostic files — common in legacy XYZ exports — must be rejected at the ingestion gate; silent passthrough corrupts every downstream calculation.

Converting LAS to XYZ for bathymetry covers baseline schema mapping for format interchange. Production environments should enforce stricter validation: laspy for header inspection and pyproj for datum verification. Vertical references must resolve to MLLW, NAVD88, or ellipsoidal heights aligned with the survey’s tidal reduction model — a process documented in detail in the tidal datum transformation pipeline. Missing LASF_Projection VLRs or undefined Z-units trigger pipeline aborts, not warnings.

import logging
from pathlib import Path
import laspy
import pyproj

logger = logging.getLogger(__name__)


def validate_crs_and_vertical_datum(input_path: Path, expected_epsg: int) -> None:
    """Enforce strict CRS and vertical datum validation before pipeline execution.

    Raises ValueError if the file lacks a projection VLR, uses an unsupported
    point format, or references an invalid EPSG code.  Never silently continues.
    """
    with laspy.open(str(input_path)) as reader:
        has_projection_vlr = any(
            getattr(vlr, "user_id", "") == "LASF_Projection"
            for vlr in reader.header.vlrs
        )
        if not has_projection_vlr:
            raise ValueError(
                f"{input_path.name}: lacks embedded LASF_Projection VLR — "
                "convert with PDAL before ingestion."
            )
        # Point formats 0-5 lack GPS time and full return-number fields
        if reader.header.point_format.id not in (3, 4, 6, 7, 8, 10):
            raise ValueError(
                f"{input_path.name}: point format {reader.header.point_format.id} "
                "unsupported; requires Z, intensity, and classification dimensions."
            )
    # Validate the EPSG code itself — raises CRSError for unknown codes
    try:
        crs = pyproj.CRS.from_epsg(expected_epsg)
    except pyproj.exceptions.CRSError as exc:
        raise ValueError(f"Invalid EPSG {expected_epsg}: {exc}") from exc

    logger.info(
        "CRS validation passed: %s → EPSG:%d (%s)",
        input_path.name,
        expected_epsg,
        crs.name,
    )

---

Memory-Constrained Python Implementation #

The pipeline below constructs, validates, and executes a PDAL JSON pipeline with explicit chunking via filters.splitter. It avoids monolithic file reads: PDAL’s pipeline engine processes tiles in parallel across spatial partitions.

import json
import logging
import subprocess
import tempfile
from pathlib import Path
from typing import Any

logger = logging.getLogger(__name__)


def build_pdal_pipeline(
    input_path: Path,
    output_path: Path,
    epsg: int,
    sor_mean_k: int = 10,
    sor_std_multiplier: float = 2.0,
    radius: float = 1.5,
    min_neighbors: int = 5,
    max_depth: float = -200.0,
    min_depth: float = 0.0,
    tile_length: float = 300.0,
) -> dict[str, Any]:
    """Construct a deterministic PDAL pipeline for MBES point cloud filtering.

    Stages:
      1. readers.las          — load with explicit CRS override
      2. filters.range        — depth-range clip (Z positive-up)
      3. filters.splitter     — tile-based I/O for memory bounding
      4. filters.outlier      — statistical outlier removal (SOR)
      5. filters.range        — drop Class 7 (noise) tags from SOR
      6. filters.radiusoutlier — remove isolated multipath singletons
      7. filters.normal       — compute local surface normals for slope check
      8. filters.range        — drop returns violating slope continuity
      9. filters.assign       — reclassify surviving returns to Class 2 (Ground)
      10. writers.las         — LAS 1.4, point format 6, full VLR passthrough
    """
    pipeline: list[dict[str, Any]] = [
        {
            "type": "readers.las",
            "filename": str(input_path),
            "spatialreference": f"EPSG:{epsg}",
        },
        {
            # Z is positive-up in PDAL; negate survey depth for lower bound
            "type": "filters.range",
            "limits": f"Z[{max_depth}:{min_depth}]",
        },
        {
            # Tile-based I/O: prevents OOM on surveys > 500 M points
            "type": "filters.splitter",
            "length": tile_length,
            "origin_x": 0.0,
            "origin_y": 0.0,
        },
        {
            # SOR: tags outliers as Classification 7 (Low Point / noise)
            "type": "filters.outlier",
            "method": "statistical",
            "mean_k": sor_mean_k,
            "multiplier": sor_std_multiplier,
        },
        {
            # Remove SOR-tagged noise points before density filter
            "type": "filters.range",
            "limits": "Classification![7:7]",
        },
        {
            # Radius density filter: removes isolated acoustic returns
            "type": "filters.radiusoutlier",
            "radius": radius,
            "min_points": min_neighbors,
        },
        {
            # Surface normal computation for slope validation
            # knn=8 matches SOR mean_k; provides stable normals at this density
            "type": "filters.normal",
            "knn": 8,
        },
        {
            # Slope gate: NormalZ < cos(60°) ≈ 0.5 indicates >60° gradient — reject
            # Adjust threshold for survey area expected morphology
            "type": "filters.range",
            "limits": "NormalZ[0.5:]",
        },
        {
            # Promote surviving valid returns to ASPRS Class 2 (Ground).
            # PDAL >= 2.5 takes a single assignment expression with an inline
            # WHERE clause; a separate "condition" key is silently ignored.
            "type": "filters.assign",
            "value": "Classification = 2 WHERE Classification == 1",
        },
        {
            "type": "writers.las",
            "filename": str(output_path),
            "minor_version": 4,
            "dataformat_id": 6,
            # Forward all VLRs so CRS and metadata survive the round-trip
            "forward": "all",
        },
    ]
    return {"pipeline": pipeline}


def execute_filtering_pipeline(
    pipeline_json: dict[str, Any],
    timeout_seconds: int = 3600,
) -> None:
    """Execute a PDAL pipeline via CLI subprocess.

    Raises RuntimeError with the full PDAL stderr on non-zero exit.
    timeout_seconds guards against hung I/O on network-mounted storage.
    """
    with tempfile.NamedTemporaryFile(
        mode="w", suffix=".json", delete=False
    ) as tmp:
        json.dump(pipeline_json, tmp, indent=2)
        tmp_path = tmp.name

    logger.info("Executing PDAL pipeline from %s", tmp_path)
    result = subprocess.run(
        ["pdal", "pipeline", tmp_path],
        capture_output=True,
        text=True,
        timeout=timeout_seconds,
    )
    if result.returncode != 0:
        raise RuntimeError(
            f"PDAL pipeline failed (exit {result.returncode}):\n{result.stderr}"
        )
    logger.info("PDAL pipeline completed successfully:\n%s", result.stdout)


def run_mbes_filter(
    input_path: Path,
    output_path: Path,
    epsg: int,
    **kwargs: Any,
) -> None:
    """Orchestrate CRS validation, pipeline build, and execution."""
    validate_crs_and_vertical_datum(input_path, epsg)
    pipeline = build_pdal_pipeline(input_path, output_path, epsg, **kwargs)
    execute_filtering_pipeline(pipeline)
    logger.info("Filtered output written to %s", output_path)

---

Validation Gates and Quality Control #

Post-filtering validation must be automated and threshold-gated — manual spot-checks are insufficient at multi-terabyte survey scale. Three mandatory checkpoints apply before the output LAZ is promoted to the gridding stage.

Gate 1 — Point retention ratio. Extract pre- and post-filter counts with pdal info --summary and assert that fewer than 15% of points were removed. Aggressive rejection (>25%) usually indicates a CRS mismatch or an incorrect depth-range limit rather than genuine noise.

# Capture summary JSON and extract point count
pdal info --summary cleaned_output.laz | python3 -c "
import json, sys
data = json.load(sys.stdin)
n = data['summary']['num_points']
print(f'Retained: {n:,} points')
assert n > 0, 'Output is empty — check depth range and CRS'
"

Gate 2 — IHO S-44 vertical accuracy. RMS error against check-line crossings must meet the survey order threshold. For Special Order surveys:

$$\text{RMSE}_Z \leq 0.25\,\text{m} + 0.0075 \times d$$

where $d$ is depth in metres. Order 1a relaxes this to $0.5\,\text{m} + 0.013 \times d$. Compute crossline RMS from the filtered cloud before gridding — residuals are easier to trace at point level than in a rasterized DEM.

Gate 3 — Classification distribution. Verify that surviving points carry Classification 2 (Ground) and that no Classification 7 (noise) tags remain:

import logging

import laspy
import numpy as np

logger = logging.getLogger(__name__)


def assert_classification_distribution(output_path: str) -> None:
    """Raise AssertionError if noise points remain or no ground points exist."""
    with laspy.open(output_path) as reader:
        las = reader.read()
    codes, counts = np.unique(np.asarray(las.classification), return_counts=True)
    class_map: dict[int, int] = {int(c): int(n) for c, n in zip(codes, counts)}
    logger.info("Classification distribution: %s", class_map)
    assert class_map.get(7, 0) == 0, (
        f"Classification 7 (noise) still present: {class_map.get(7, 0)} points — "
        "re-run filters.range to drop them."
    )
    assert class_map.get(2, 0) > 0, (
        "No Class 2 (Ground) points in output — check filters.assign condition."
    )

For deeper debugging of outlier removal edge cases, see Using PDAL for Bathymetric Point Cloud Cleaning, which covers SOR parameter tuning and radius selection for variable-density survey lines.

---

Common Failure Modes and Diagnosis #

1. Silent Z-unit inversion (depth reported as elevation) #

Symptom: pdal info --summary shows Z range of +0.0 to +200.0 where negative values are expected, and the depth-range clip removes all points.

Root cause: Source data was exported with Z positive-down (depth convention) rather than PDAL’s positive-up (elevation convention). Teledyne .s7k and some .all exports apply the depth sign before writing XYZ.

Fix: Negate Z during readers.las via a PDAL expression, or apply filters.ferry to remap the Z dimension before range clipping:

{
  "type": "filters.assign",
  "value": "Z = Z * (-1.0)"
}

Insert this stage immediately after readers.las and before filters.range.

---

2. Missing LASF_Projection VLR after format conversion #

Symptom: validate_crs_and_vertical_datum raises ValueError: lacks embedded LASF_Projection VLR.

Root cause: The conversion tool (e.g., las2las, older MB-System export) wrote a LAS 1.0–1.2 file without WKT or GeoTIFF VLR records. The CRS alignment workflow covers how to diagnose and repair CRS records in GIS files generally; for LAZ specifically, PDAL’s writers.las with an explicit spatialreference is the authoritative fix.

Fix:

pdal translate input_no_crs.las output_with_crs.laz \
  --writers.las.spatialreference="EPSG:32618" \
  --writers.las.minor_version=4 \
  --writers.las.dataformat_id=6

---

3. OOM crash on surveys exceeding 500 million points #

Symptom: Python process killed (MemoryError) or PDAL exits with std::bad_alloc when loading a full tile.

Root cause: filters.splitter was omitted or tile_length was set too large (>500 m). At 1 m beam spacing, a 500 m tile holds 250,000 points per layer — manageable — but 1 km tiles at 0.5 m spacing can exceed 4 million points in memory simultaneously.

Fix: Set tile_length to 200–250 m for high-density surveys. For surveys stored on network-mounted NFS or S3-mounted FUSE filesystems, also set PDAL_DRIVER_OPTIONS to increase read buffer size and reduce seek overhead.

---

4. Slope gate removes valid steep-slope morphology #

Symptom: Post-filter point density drops sharply in submarine canyon walls or coral pinnacle flanks; cross-sections show missing data at slope angles >45°.

Root cause: The NormalZ >= 0.5 threshold (cosine of 60°) is too aggressive for high-relief bathymetry. It was tuned for soft-sediment continental shelf surveys.

Fix: Relax to NormalZ[0.25:] (≈ 75° maximum slope) for reef or canyon surveys, or remove the slope gate entirely and rely solely on SOR and radius outlier removal. Log the rejection count per tile to detect anomalously high slope-gate losses during QC.

The slope gate keys on NormalZ — the vertical component of each return’s surface normal, equal to the cosine of the seabed slope angle. Flat seabed yields NormalZ ≈ 1.0; a vertical wall yields NormalZ ≈ 0.0. The threshold sets the steepest slope the filter accepts, so choosing it requires knowing the survey area’s expected relief.

---

Pipeline Integration and Downstream Handoff #

Once filtering completes, the output LAZ is staged with a sidecar metadata manifest and fed to the DEM interpolation pipeline for rasterization. The manifest records filter parameters, PDAL version, point counts, and the IHO gate results so that gridding operators can reconstruct the provenance of any DEM cell.

The manifest is written as JSON alongside the LAZ output:

import json
import datetime
from pathlib import Path


def write_filter_manifest(
    output_path: Path,
    pipeline_json: dict,
    points_retained: int,
    points_input: int,
    iho_order: str,
) -> Path:
    """Write a sidecar JSON manifest recording filter provenance."""
    manifest = {
        "generated_utc": datetime.datetime.utcnow().isoformat(),
        "output_file": output_path.name,
        "points_input": points_input,
        "points_retained": points_retained,
        "rejection_rate_pct": round(
            100.0 * (points_input - points_retained) / max(points_input, 1), 2
        ),
        "iho_order": iho_order,
        "pdal_pipeline": pipeline_json,
    }
    manifest_path = output_path.with_suffix(".filter_manifest.json")
    manifest_path.write_text(json.dumps(manifest, indent=2))
    return manifest_path

The filtered and manifested LAZ then enters surface smoothing to mitigate residual gridding artifacts while preserving genuine morphological relief — channelized features, coral crests, and scour pockets that over-smoothing would flatten.

Removing persistent systematic artifacts (banding, nadir lines, motion-induced spikes) that survive the statistical filters is addressed by the bathymetric artifact removal pipeline.

---

Related #

• Using PDAL for Bathymetric Point Cloud Cleaning — SOR parameter tuning and PDAL pipeline JSON deep-dive
• How to Convert LAS to XYZ for Bathymetry — format interchange and schema mapping
• DEM Interpolation Techniques for Seafloor Mapping — downstream gridding from validated point clouds
• Removing Bathymetric Artifacts and Noise — systematic artifact removal after statistical filtering
• Surface Smoothing Algorithms in Python — post-grid smoothing to preserve morphological features

Up: Bathymetric Processing & Terrain Modeling