Chaining PDAL Stages for Data Cleaning

TL;DR: Chain readers.las → filters.range → filters.outlier → filters.assign → writers.las, call pdal.Pipeline(config).execute(), and PDAL passes a shared in-memory PointView buffer through each stage — no intermediate files, no redundant I/O.

# Context and Motivation

Raw LiDAR surveys contain acquisition artifacts that make direct analysis unreliable: atmospheric scatter above the sensor ceiling, sub-terrain noise from ground-penetrating returns, scan-angle-dependent intensity distortion, and unclassified or misclassified returns from water bodies and low vegetation. Cleaning these before downstream terrain modelling or feature extraction is non-negotiable in production workflows.

This guide is part of PDAL Stage Chaining, which explains the full buffer-passing execution model. The cleaning-specific techniques here build directly on how PDAL Pipeline Architecture & Execution allocates memory and propagates dimensions between stages. Understanding that broader execution context helps you reason about why stage order matters, not just what order to use.

# Stage-Flow Diagram

The diagram below shows the five-stage cleaning chain, what each stage removes or transforms, and the shared PointView buffer that eliminates intermediate disk writes.

# Prerequisites and Assumptions

• PDAL 2.5 or later with Python bindings installed (pip install pdal)
• Python 3.10+ with numpy available
• Input LAS/LAZ file with at least X, Y, Z, ScanAngleRank (LAS 1.2/1.3) or ScanAngle (LAS 1.4), and Classification dimensions
• Known input CRS — set spatialreference in the reader rather than relying on embedded VLR data, which is sometimes absent in older files
• Test dataset recommendation: any USGS 3DEP tile from the 3D Elevation Program or the PDAL sample data at https://github.com/PDAL/data

# Step-by-Step Implementation

# Step 1 — Declare the reader with a CRS anchor

{
  "type": "readers.las",
  "filename": "raw_survey.las",
  "spatialreference": "EPSG:32618"
}

Setting spatialreference here overrides any embedded projection record and prevents silent datum mismatches when a downstream spatial reprojection stage or range-based spatial filter assumes a specific unit system. Always verify the EPSG code matches the sensor’s recording datum — not the delivery datum — before adding a reprojection stage.

# Step 2 — Clamp to physically valid bounds with filters.range

{
  "type": "filters.range",
  "limits": "Z[10:450],ScanAngleRank[-20:20]"
}

The limits string is a comma-separated list of Dimension[min:max] expressions. Place filters.range before filters.outlier so the statistical algorithm only sees points that are geometrically plausible — running outlier detection on atmospheric noise wastes CPU and skews the neighbourhood statistics. The Pipeline Filtering Logic section covers filter ordering rules in depth.

# Step 3 — Remove statistical noise with filters.outlier

{
  "type": "filters.outlier",
  "method": "statistical",
  "mean_k": 12,
  "multiplier": 2.5
}

filters.outlier marks isolated points as noise (Classification 7) rather than dropping them immediately, giving downstream stages a chance to inspect or preserve them. Set mean_k to roughly 10–16 for airborne surveys and 6–10 for terrestrial scans where point density is higher but spatial extent is smaller. See the dedicated Applying Statistical Outlier Filters in PDAL guide for multiplier tuning guidance.

# Step 4 — Assign or correct classification codes with filters.assign

{
  "type": "filters.assign",
  "value": "Classification = 2 WHERE Classification == 0"
}

The value parameter is a PDAL expression string. The WHERE clause restricts the assignment to unclassified points (code 0), leaving ground (2), vegetation (3–5), and building (6) labels from the original data intact. Classification codes follow the ASPRS classification scheme — consult that reference when you need to preserve or recategorise specific return types.

# Step 5 — Write compressed output

{
  "type": "writers.las",
  "filename": "cleaned_survey.laz",
  "compression": "true",
  "extra_dims": "all"
}

extra_dims: all forwards any non-standard dimensions (HAG, NDVI-derived bands, custom sensor fields) that earlier stages may have added, preventing data loss when chaining this pipeline into a larger workflow.

# Complete Working Example

The following self-contained script wires all five stages together and extracts metadata for audit logging.

import pdal
import json
import sys
from pathlib import Path


def build_cleaning_pipeline(input_file: str, output_file: str) -> list[dict]:
    """Return a PDAL stage list for LiDAR data cleaning."""
    return [
        {
            "type": "readers.las",
            "filename": input_file,
            "spatialreference": "EPSG:32618"
        },
        {
            "type": "filters.range",
            "limits": "Z[10:450],ScanAngleRank[-20:20]"
        },
        {
            "type": "filters.outlier",
            "method": "statistical",
            "mean_k": 12,
            "multiplier": 2.5
        },
        {
            "type": "filters.assign",
            "value": "Classification = 2 WHERE Classification == 0"
        },
        {
            "type": "writers.las",
            "filename": output_file,
            "compression": "true",
            "extra_dims": "all"
        }
    ]


def run_cleaning_pipeline(input_file: str, output_file: str) -> dict:
    """Execute the cleaning pipeline and return audit metadata."""
    if not Path(input_file).exists():
        return {"status": "error", "message": f"Input file not found: {input_file}"}

    stages = build_cleaning_pipeline(input_file, output_file)

    try:
        pipeline = pdal.Pipeline(json.dumps(stages))
        point_count = pipeline.execute()
        meta = pipeline.metadata
        arrays = pipeline.arrays

        return {
            "status": "success",
            "points_out": point_count,
            "array_shape": arrays[0].shape if arrays else (0,),
            "stages_run": len(stages),
            "bounds": meta.get("metadata", {}).get("filters.range", {}).get("bbox", {}),
        }

    except RuntimeError as exc:
        print(f"[pdal] RuntimeError: {exc}", file=sys.stderr)
        return {"status": "error", "message": str(exc)}
    except Exception as exc:
        print(f"[pdal] Unexpected error: {exc}", file=sys.stderr)
        return {"status": "error", "message": str(exc)}


if __name__ == "__main__":
    INPUT_LAS = "raw_survey.las"
    OUTPUT_LAZ = "cleaned_survey.laz"

    result = run_cleaning_pipeline(INPUT_LAS, OUTPUT_LAZ)

    if result["status"] == "success":
        print(f"Points written: {result['points_out']}")
        print(f"Array shape:    {result['array_shape']}")
    else:
        print(f"Pipeline failed: {result['message']}", file=sys.stderr)
        sys.exit(1)

# Key Parameter Table

# Verification

After execution, confirm the pipeline produced a valid, complete output:

import pdal
import json

# Re-read the output and check point count and dimension names
verify_pipeline = [
    {"type": "readers.las", "filename": "cleaned_survey.laz"}
]
pipe = pdal.Pipeline(json.dumps(verify_pipeline))
count = pipe.execute()
schema = pipe.schema

dim_names = [d["name"] for d in schema["schema"]["dimensions"]]
assert "Classification" in dim_names, "Classification dimension missing"
assert count > 0, "Output file is empty"

print(f"Verified: {count} points, dimensions: {dim_names}")

Also inspect bounding boxes in the metadata to catch silent data loss from over-aggressive range clamps:

meta = pipe.metadata
bbox = meta["metadata"]["readers.las"]["bbox"]
print(f"Output bounds: Z {bbox['minz']:.2f} – {bbox['maxz']:.2f} m")

If minz equals your lower filters.range bound, you may have clipped valid ground returns — raise the lower Z limit and re-run. The memory management page explains how to verify buffer allocation and detect truncation from out-of-memory conditions during execution.

# Gotchas and Edge Cases

1. filters.outlier flags valid thin-vegetation returns. In sparse forests or low-shrub areas, isolated single returns from sub-metre vegetation look statistically identical to noise. Lower multiplier to 1.8 and inspect a sample tile with pdal info --stats before applying fleet-wide. Alternatively, run filters.outlier only on ground-classified subsets by combining it with a filters.range on Classification.

2. filters.range on ScanAngleRank vs ScanAngle. PDAL 2.5+ uses ScanAngle (float, degrees) in LAS 1.4 files and ScanAngleRank (int8, 1/10-degree units scaled −128 to +127) in LAS 1.2/1.3. Applying a ScanAngleRank[-20:20] limit to a LAS 1.4 file that stores ScanAngle will silently pass all points because the dimension name does not match. Run pdal info --schema raw_survey.las | grep -i scan to confirm which dimension is present before writing the limits expression.

3. filters.assign without a WHERE clause overwrites existing labels. Omitting the WHERE Classification == 0 guard sets every point to code 2, destroying building, vegetation, and water classifications that the sensor vendor may have pre-computed. Always scope assignment expressions with a WHERE predicate. Cross-reference the ASPRS classification codes table to confirm you are targeting the correct numeric codes before running a bulk assignment.

4. compression: "true" requires a .laz file extension. If you set compression: "true" but write to a .las filename, writers.las raises a RuntimeError at execution time. The filename extension and the compression flag must agree — use .laz for compressed output and .las for uncompressed.

# Frequently Asked Questions

# What order should PDAL filters run in a cleaning pipeline?

Apply coordinate normalisation first (if a reprojection is needed), then geometric range clamps with filters.range, then statistical outlier removal with filters.outlier, then classification updates with filters.assign. Running outlier detection before range clamps wastes CPU cycles on points that will be discarded anyway and skews the neighbourhood density statistics that filters.outlier relies on.

# Does PDAL write intermediate files when stages are chained?

No. PDAL passes a shared PointView buffer between stages entirely in memory. Intermediate disk writes only occur when you explicitly insert a writers.* stage mid-pipeline. This is what makes the chain fast for large files — I/O happens once on read and once on write.

# How do I handle multi-gigabyte LAS files in a cleaning chain?

Insert filters.splitter early in the chain. Set the length parameter to tile the point cloud into chunks that fit in available RAM, forcing sequential processing and releasing memory between tiles. See the memory management page for recommended length values relative to available RAM and point density.

# Related

• PDAL Stage Chaining — parent guide covering the full buffer-passing execution model
• Applying Statistical Outlier Filters in PDAL — deep dive on mean_k and multiplier calibration
• Pipeline Filtering Logic — filter ordering rules and dimension propagation
• Spatial Reprojection — inserting a CRS transformation before geometric filters
• PDAL Pipeline Architecture & Execution — top-level guide to the execution model