Applying Statistical Outlier Filters in PDAL

TL;DR: Insert filters.outlier with "method": "statistical" into your PDAL pipeline immediately after the reader to remove noise points; set mean_k between 10 and 30 and multiplier between 2.0 and 3.0 for typical airborne LiDAR.

This guide is part of Pipeline Filtering Logic, the parent reference for sequencing and combining filter stages in a PDAL execution model. Statistical outlier removal (SOR) is the most common first cleaning pass applied to raw LiDAR returns—before ground classification, rasterization, or any measurement that depends on surface continuity.

# Context and Motivation

Raw LiDAR acquisitions contain returns that no classification code can meaningfully represent: atmospheric scatter at apogee, bird strikes crossing the scan plane, multipath reflections off glass facades, and solar noise in single-photon systems. These isolated points have no spatial neighbors that resemble them; they appear as lone spikes well above or below the true surface.

Fixed-threshold Z-range filters catch gross spikes but miss noise at mid-elevation because they have no knowledge of local point density. Statistical outlier removal solves this by measuring each point’s relationship to its k nearest neighbors rather than checking an absolute bound. A point is an outlier only if it is unusually far from its local neighborhood—not because it exceeds a fixed threshold. This density-awareness makes SOR essential for mixed-resolution datasets where a single Z cutoff would destroy valid sparse returns in rural areas while leaving dense urban noise untouched.

Within a broader PDAL Pipeline Architecture & Execution workflow, SOR sits at the boundary between raw acquisition data and analytically meaningful point clouds. Everything downstream—ground models, tree height extraction, building footprints—benefits from a clean input. The pipeline validation stage you run before deployment will flag a missing or misnamed filter immediately, catching configuration errors before they propagate to expensive downstream steps.

# Prerequisites and Assumptions

• PDAL 2.4 or later — earlier releases have different parameter defaults for filters.outlier.
• Python 3.10+ with pdal bindings installed via conda install -c conda-forge python-pdal or pip install pdal.
• Input file format: LAS 1.2–1.4 or LAZ. The stage reads the X, Y, Z dimensions; it does not require Classification or any custom dimension.
• numpy for in-memory array access after execution.
• A test tile of 5–20 million points is recommended for initial parameter calibration before processing full flight lines.

Verify the setup:

import pdal
print(pdal.__version__)   # should be 2.4.x or higher

# Step-by-Step Implementation

# Step 1 — Build the minimal pipeline

The bare-minimum pipeline to apply statistical outlier removal reads a file, filters it, and writes a cleaned output:

{
  "pipeline": [
    {
      "type": "readers.las",
      "filename": "input_raw.laz"
    },
    {
      "type": "filters.outlier",
      "method": "statistical",
      "mean_k": 20,
      "multiplier": 3.0
    },
    {
      "type": "writers.las",
      "filename": "output_clean.laz",
      "forward": "all"
    }
  ]
}

"forward": "all" on the writer preserves every original point attribute (intensity, return number, scan angle, any custom dimensions) that the filter stage did not modify. Omitting it drops non-standard dimensions silently.

# Step 2 — Choose parameters for your acquisition type

The two parameters that control how aggressively points are removed:

For a first pass on an unknown airborne dataset, mean_k=20 and multiplier=3.0 are safe starting values. Validate on a representative 5 km² tile before committing to full-project parameters.

# Step 3 — Position the stage correctly in a real pipeline

In any production pipeline that continues to PDAL stage chaining for classification or surface modelling, place filters.outlier immediately after the reader and any CRS-setting step, but before any classification or rasterization stage:

pipeline_def = [
    {"type": "readers.las", "filename": "survey.laz"},
    # Assign CRS if not already embedded in the header
    {"type": "filters.reprojection", "in_srs": "EPSG:32618", "out_srs": "EPSG:32618"},
    # Statistical outlier removal first — before SMRF or any feature extractor
    {"type": "filters.outlier", "method": "statistical", "mean_k": 20, "multiplier": 3.0},
    # Ground classification after cleaning
    {"type": "filters.smrf", "slope": 0.15, "window": 18.0, "threshold": 0.5, "scalar": 1.2},
    {"type": "writers.las", "filename": "classified_clean.laz", "forward": "all"},
]

Running SOR after filters.smrf risks removing valid ground points that SMRF has already classified as class 2, which corrupts DTM generation downstream.

# Step 4 — Execute and capture results

import pdal

pipeline = pdal.Pipeline(pipeline_def)
count = pipeline.execute()
print(f"Points retained: {count}")

pipeline.execute() returns the total number of points that passed through to the final stage. For typical airborne acquisitions you should see 0.1–2 % of points removed by SOR; more than 5 % warrants investigating whether multiplier is too aggressive for your data density.

# Complete Working Example

The following self-contained module can be copy-pasted and run against any LAS/LAZ file. It handles both file-output and in-memory array modes, logs progress, and validates results before returning:

import json
import sys
from typing import Any, Dict, Optional

import numpy as np
import pdal


def apply_statistical_filter(
    input_path: str,
    output_path: Optional[str] = None,
    mean_k: int = 20,
    multiplier: float = 3.0,
) -> Dict[str, Any]:
    """
    Apply PDAL's statistical outlier filter to a LAS/LAZ file.

    Args:
        input_path:  Path to source LAS/LAZ file.
        output_path: Optional path for the cleaned output. If None, results
                     are only available as an in-memory NumPy array.
        mean_k:      Number of nearest neighbours for the local mean distance.
        multiplier:  Standard-deviation threshold; points beyond this multiple
                     of the global mean distance are removed.

    Returns:
        Dict with 'points_kept' (int), 'points_removed' (int),
        and 'metadata' (dict).
    """
    # Count raw points before filtering for comparison
    count_pipeline = pdal.Pipeline(
        [{"type": "readers.las", "filename": input_path}]
    )
    raw_count = count_pipeline.execute()

    pipeline_stages: list = [
        {"type": "readers.las", "filename": input_path},
        {
            "type": "filters.outlier",
            "method": "statistical",
            "mean_k": mean_k,
            "multiplier": multiplier,
        },
    ]

    if output_path:
        pipeline_stages.append(
            {
                "type": "writers.las",
                "filename": output_path,
                "forward": "all",
            }
        )

    pipeline = pdal.Pipeline(pipeline_stages)

    try:
        kept = pipeline.execute()
    except RuntimeError as exc:
        print(f"Pipeline failed: {exc}", file=sys.stderr)
        sys.exit(1)

    removed = raw_count - kept
    pct_removed = (removed / raw_count * 100) if raw_count else 0
    print(f"Input points : {raw_count:,}")
    print(f"Points kept  : {kept:,}")
    print(f"Points removed: {removed:,} ({pct_removed:.2f}%)")

    if pct_removed > 5.0:
        print(
            "WARNING: >5% of points removed. Consider raising multiplier "
            "or increasing mean_k to reduce over-filtering.",
            file=sys.stderr,
        )

    return {
        "points_kept": kept,
        "points_removed": removed,
        "metadata": pipeline.metadata,
    }


def get_filtered_array(
    input_path: str,
    mean_k: int = 20,
    multiplier: float = 3.0,
) -> np.ndarray:
    """Return the filtered cloud as a NumPy structured array (no disk write)."""
    pipeline = pdal.Pipeline(
        [
            {"type": "readers.las", "filename": input_path},
            {
                "type": "filters.outlier",
                "method": "statistical",
                "mean_k": mean_k,
                "multiplier": multiplier,
            },
        ]
    )
    pipeline.execute()
    return pipeline.arrays[0]


if __name__ == "__main__":
    result = apply_statistical_filter(
        "input_raw.laz",
        "output_clean.laz",
        mean_k=20,
        multiplier=3.0,
    )
    # Confirm CRS survived the filter pass
    meta = json.loads(result["metadata"])
    reader_meta = meta["metadata"]["readers.las"][0]
    srs_wkt = reader_meta.get("srs", {}).get("wkt", "No CRS embedded")
    print(f"Output CRS: {srs_wkt[:80]}...")

# Key Parameter Table

# Verification

After execution, confirm three things:

1. Point count decreased as expected:

result = apply_statistical_filter("raw.laz", "clean.laz")
assert result["points_kept"] > 0, "Pipeline produced empty output"
assert result["points_removed"] > 0, "No points were removed — check parameters"
# For typical airborne data, SOR removes 0.1–2% of points
assert result["points_kept"] / (result["points_kept"] + result["points_removed"]) > 0.95, \
    "More than 5% removed — likely over-filtering"

2. Output dimensions are intact:

arr = get_filtered_array("raw.laz")
expected_dims = {"X", "Y", "Z", "Intensity", "ReturnNumber", "Classification"}
actual_dims = set(arr.dtype.names)
assert expected_dims.issubset(actual_dims), f"Missing dims: {expected_dims - actual_dims}"

3. CRS is preserved:

import json
meta = json.loads(result["metadata"])
reader_meta = meta["metadata"]["readers.las"][0]
print(reader_meta.get("srs", {}).get("wkt", "No CRS embedded"))

You can also run pdal info --stats output_clean.laz from the command line and compare the Z range and point count against the raw file to spot gross over-filtering quickly.

# Gotchas and Edge Cases

Stage name confusion. The correct PDAL stage is filters.outlier with "method": "statistical". There is no filters.statistical stage; using that name causes a RuntimeError: Couldn't create stage 'filters.statistical'. The pipeline validation workflow catches this before any data is read if you run pdal --validate on your JSON definition during CI.

Applying SOR after classification removes valid sparse classes. Transmission towers, isolated trees, and bridge decks have far fewer neighbors than the surrounding ground; SOR sees them as outliers. If your project requires preserving classified returns, run SOR on the raw, unclassified cloud first, or apply it selectively via filters.expression to subset only unclassified points before the cleaning pass.

Tile-edge artifacts. The k-d tree is built per pipeline invocation. Points at tile boundaries have artificially truncated neighborhoods because their true neighbors sit in the adjacent tile. Apply a 5–10 m overlap buffer when tiling with filters.splitter, then strip the buffer after SOR and before merging outputs. Alternatively, consult the memory management strategies for processing larger in-memory chunks that span tile boundaries without splitting.

k-d tree memory scaling. Memory consumption for the neighborhood index grows with point count at roughly O(n log n). For datasets above 50 million points per tile, monitor peak RAM with tracemalloc during development and budget at least 16 GB for 100 M-point tiles. If processing fails with an std::bad_alloc error, reduce tile size before adjusting filter parameters.

Pre-filtered sensor data. Some survey-grade scanners apply their own noise rejection before writing LAS. Running SOR on already-clean data with a low multiplier can strip valid but sparse features. Check pdal info --metadata input.laz for manufacturer-provided noise flags or pre-existing classification values before adding a redundant filter stage.

# Related

• Pipeline Filtering Logic — parent reference covering filter sequencing, predicate evaluation, and spatial constraints
• PDAL Stage Chaining — how buffers propagate between stages and how to chain SOR with classification stages
• Chaining PDAL Stages for Data Cleaning — worked example combining outlier removal, range filtering, and ground classification
• Memory Management — tiling strategies and RAM budgeting for large airborne datasets
• Pipeline Validation — catch stage name errors and schema violations before execution