Pipeline Filtering Logic in PDAL: Attribute, Spatial, and Statistical Filters

Pipeline filtering logic determines which points survive each stage of a processing graph and which are discarded — making it the primary quality-control mechanism in any LiDAR workflow. In production environments that depend on the PDAL Pipeline Architecture & Execution model, filters run as discrete buffer-transforming stages arranged in a directed acyclic graph: each stage reads a structured point array from its predecessor, applies a predicate or transformation, and emits a reduced or annotated array downstream. Getting filter order, parameter values, and dimension contracts right separates deterministic, auditable pipelines from brittle scripts that silently drop valid data.

This guide targets LiDAR analysts, Python GIS developers, and surveying teams who need a rigorous, production-ready reference for attribute filtering, spatial subsetting, statistical cleaning, and dimension management inside PDAL pipelines.

# Prerequisites

Before implementing filtering stages, confirm your environment meets these requirements:

• PDAL 2.5+ with the filters.* stage library compiled in (verify with pdal --version).
• Python 3.10+ with python-pdal installed via conda install -c conda-forge python-pdal or pip install pdal.
• NumPy for post-execution buffer inspection; pyproj if you need to validate CRS alignment programmatically.
• Input data: LAS 1.2–1.4 or LAZ files with at minimum the X, Y, Z, Classification, and ReturnNumber dimensions populated. Statistical filters additionally require enough neighbouring points to form meaningful density statistics — sparse test tiles with fewer than 1,000 points per square metre may produce unreliable results.
• CRS metadata present in the file header. Spatial filters (filters.crop) silently return zero points when bounds are expressed in a different CRS than the input. Run pdal info --metadata input.laz to confirm spatialreference is populated before writing any spatial predicate.
• Familiarity with PDAL stage chaining — filter stages obey the same dimension propagation and buffer-passing contracts as all other stages.

# Core Workflow Architecture

Filtering execution unfolds in five deterministic phases:

1. Objective definition — Specify the exact reduction goal: attribute bounds, spatial extent, noise threshold, or dimension set.
2. Stage sequencing — Order stages so that coordinate transformations (e.g., spatial reprojection) precede spatial filters, and attribute filters precede computationally expensive statistical passes.
3. Pipeline construction — Encode the stage sequence as a Python list of dicts or a JSON string and instantiate pdal.Pipeline.
4. Lazy execution — Call pipeline.execute(); PDAL evaluates the full DAG synchronously, allocating new point buffers at each stage rather than modifying data in-place.
5. Validation and audit — Assert output point counts, inspect pipeline.metadata for dimension lists, and compare CRS round-trip values.

Stage order has hard dependencies. filters.range or filters.expression should run before filters.outlier to reduce the neighbour-search space. filters.crop must run after any reprojection stage. filters.assign for dimension injection must precede any stage that reads the injected dimension.

# Full Implementation

The following self-contained Python module applies a representative filtering sequence: classification and height bounding, spatial subsetting, statistical noise removal, and metadata capture.

import logging
import pdal
from pathlib import Path

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


def run_filter_pipeline(
    input_path: str,
    output_path: str,
    classification: int = 2,
    z_min: float = 0.5,
    z_max: float = 200.0,
    bbox: tuple[float, float, float, float] | None = None,
    outlier_mean_k: int = 12,
    outlier_multiplier: float = 2.5,
) -> dict:
    """
    Apply attribute, spatial, and statistical filters to a LiDAR point cloud.

    Parameters
    ----------
    input_path        : Path to input LAS/LAZ file.
    output_path       : Destination LAS/LAZ path.
    classification    : ASPRS class code to retain (default 2 = ground).
    z_min / z_max     : Elevation bounds in the native CRS vertical unit.
    bbox              : Optional (xmin, ymin, xmax, ymax) crop in the input CRS.
    outlier_mean_k    : Neighbour count for statistical outlier evaluation.
    outlier_multiplier: Standard-deviation multiplier; lower = more aggressive.

    Returns
    -------
    dict with 'point_count', 'dimensions', and 'metadata' keys.
    """
    stages: list[dict] = [
        {
            "type": "readers.las",
            "filename": str(input_path),
        },
        # 1. Attribute filter: keep only the target classification within Z bounds.
        # filters.range is evaluated per-point before any spatial or statistical
        # stage, minimising the point set that expensive neighbour searches must handle.
        {
            "type": "filters.range",
            "limits": f"Classification[{classification}:{classification}],Z[{z_min}:{z_max}]",
        },
    ]

    if bbox is not None:
        xmin, ymin, xmax, ymax = bbox
        # 2. Spatial crop — applied after attribute filtering so the neighbour
        #    graph in the outlier stage is not polluted by out-of-area points.
        stages.append(
            {
                "type": "filters.crop",
                "bounds": f"([{xmin},{xmax}],[{ymin},{ymax}])",
            }
        )

    stages += [
        # 3. Statistical outlier removal — computes mean distance to the
        #    outlier_mean_k nearest neighbours for each point and discards those
        #    whose distance exceeds outlier_multiplier * stddev of the local mean.
        {
            "type": "filters.outlier",
            "method": "statistical",
            "mean_k": outlier_mean_k,
            "multiplier": outlier_multiplier,
        },
        # 4. Remove points flagged as noise by the outlier stage.
        {
            "type": "filters.range",
            "limits": "Classification![7:7]",
        },
        {
            "type": "writers.las",
            "filename": str(output_path),
            "extra_dims": "all",
            "compression": "laszip",
        },
    ]

    pipeline = pdal.Pipeline(stages)
    try:
        pipeline.validate()
    except Exception as exc:
        log.error("Pipeline validation failed: %s", exc)
        raise

    count = pipeline.execute()
    log.info("Filtered %d points → %s", count, output_path)

    arrays = pipeline.arrays
    dims = list(arrays[0].dtype.names) if arrays else []

    return {
        "point_count": count,
        "dimensions": dims,
        "metadata": pipeline.metadata,
    }


if __name__ == "__main__":
    result = run_filter_pipeline(
        input_path="survey_raw.laz",
        output_path="survey_ground_clean.laz",
        classification=2,
        z_min=0.5,
        z_max=150.0,
        bbox=(363000.0, 5621000.0, 364000.0, 5622000.0),
        outlier_mean_k=12,
        outlier_multiplier=2.5,
    )
    print(f"Output points : {result['point_count']:,}")
    print(f"Dimensions    : {result['dimensions']}")

# Code Breakdown

Stage 1 — filters.range (attribute + height bound) filters.range uses PDAL’s compact interval notation Dimension[min:max]. Chaining Classification[2:2] and Z[0.5:200.0] in a single limits string is equivalent to a logical AND: both conditions must hold. Placing this stage first reduces the point count that all downstream stages must process, which is especially important before the O(n log n) neighbour search in filters.outlier. Running range filtering last would not only be slower but would produce incorrect statistical baselines.

Stage 2 — filters.crop (spatial subsetting) filters.crop evaluates bounds in the current coordinate reference system of the buffer — after any upstream reprojection. The bounds string ([xmin,xmax],[ymin,ymax]) uses projected metre values here; mixing these with degree-based coordinates causes silent empty output. If your pipeline includes a spatial reprojection stage, place filters.crop after it. For a systematic approach to diagnosing and fixing CRS conflicts before they reach the filter chain, see fixing CRS mismatches in point clouds.

Stage 3 — filters.outlier (statistical noise removal) method=statistical computes the mean Euclidean distance to each point’s mean_k nearest neighbours, then computes the global mean and standard deviation of those distances. Any point whose mean neighbour distance exceeds mean_k_mean + multiplier * stddev is assigned Classification = 7 (noise per ASPRS classification codes). It does not delete points — it labels them, which is why the next stage is needed.

Stage 4 — filters.range (noise removal) The notation Classification![7:7] is PDAL’s negation syntax: retain all points not in the range [7, 7]. This two-stage pattern (outlier detection → range exclusion) is preferred over in-place deletion because it preserves the full audit trail in the metadata and allows inspection of flagged points before final write.

Stage 5 — writers.las extra_dims=all forwards every dimension in the buffer — including any custom attributes added by upstream attribute mapping stages — to the output file. compression=laszip halves typical file size with no point-data loss.

# Parameter Reference

For compound predicates that exceed what filters.range can express — such as (Classification == 2 OR Classification == 6) AND Intensity > 200 — replace it with filters.expression:

{
    "type": "filters.expression",
    "expression": "(Classification == 2 || Classification == 6) && Intensity > 200"
}

# Validation and Data Integrity Checks

After every filter pipeline, run the following assertions before treating output as production-ready:

import numpy as np

result = run_filter_pipeline("input.laz", "output.laz", classification=2)

# 1. Point count must be non-zero and less than input
assert result["point_count"] > 0, "Filter removed all points — check bounds and class code"

# 2. Verify the classification dimension survived
assert "Classification" in result["dimensions"], "Classification dimension missing from output"

# 3. No noise-class points should remain
arrays = pipeline.arrays  # re-use the pipeline object from run_filter_pipeline
if arrays:
    classes = arrays[0]["Classification"]
    assert 7 not in np.unique(classes), "Noise points (class 7) present in output"

# 4. CRS round-trip check via metadata
import json
meta = json.loads(result["metadata"])
srs = meta.get("metadata", {}).get("readers.las", [{}])[0].get("spatialreference", "")
assert srs, "No CRS in output metadata — verify input file has spatial reference set"

For pipeline structure validation before any point data is read, pipeline.validate() checks JSON syntax, stage compatibility, and dimension dependency graphs without loading the point cloud. Use it in CI/CD checks against pipeline configuration files before deploying to production. The dedicated pipeline validation guide covers the full validation lifecycle, including dry-run execution patterns and schema version gating.

# Performance Tuning

Filtering performance degrades predictably at known bottlenecks. The table below shows the primary factors and the remedies for each.

For tile-parallel workloads, keep each worker’s pipeline self-contained (no shared state) and bound OMP_NUM_THREADS so worker threads do not saturate the CPU:

import os
from concurrent.futures import ProcessPoolExecutor
from pathlib import Path


def _worker(args: tuple) -> dict:
    tile, out_dir = args
    os.environ["OMP_NUM_THREADS"] = "2"
    return run_filter_pipeline(
        input_path=str(tile),
        output_path=str(Path(out_dir) / (tile.stem + "_clean.laz")),
    )


def filter_tiles(tile_dir: str, out_dir: str, workers: int = 4) -> list[dict]:
    tiles = list(Path(tile_dir).glob("*.laz"))
    with ProcessPoolExecutor(max_workers=workers) as pool:
        return list(pool.map(_worker, [(t, out_dir) for t in tiles]))

For memory management on very large single files, insert filters.splitter with length=500 (500 m grid cells) before the filter chain and process each cell independently; see the dedicated memory management guide for a full tiling strategy.

# Common Errors and Troubleshooting

1. filters.crop returns 0 points

Root cause: bounds expressed in a different CRS than the point cloud. Check with pdal info --metadata input.laz | grep spatialreference. Fix: restate bounds in the native CRS, or insert filters.reprojection with out_srs=EPSG:32632 (or your target EPSG) before filters.crop. This is the single most common silent failure in spatial subsetting.

2. KeyError: 'Classification' in post-execution array access

Root cause: the output LAS file was written without extra_dims=all and the Classification dimension exceeded the standard LAS 1.2 byte-width constraint, causing it to be silently dropped. Fix: set extra_dims=all in the writers.las stage, or ensure the output format version is LAS 1.4 with Point Data Record Format 6 or higher.

3. filters.outlier removes valid edge points in sparse datasets

Root cause: mean_k set too low (e.g., 4) on tiles with low point density — edge points have fewer than mean_k neighbours available, so their distances are artificially high. Fix: increase mean_k to 12–16 for sparse data, or switch to method=radius with a radius matched to the average point spacing.

4. RuntimeError: schema violation: dimension 'ScanAngleRank' not found

Root cause: an intermediate filter (e.g., filters.smrf or filters.pmf) dropped ScanAngleRank from the buffer, but a downstream stage declared a dependency on it. Fix: add filters.assign after the offending stage to re-introduce the dimension with a default value, or remove the downstream dependency.

5. json.dumps raises TypeError on pipeline serialization

Root cause: a Python dict built from NumPy scalars (e.g., np.float64(100.0)) fails standard JSON serialization. Fix: cast all values to native Python types before constructing the pipeline list: float(z_max), int(classification), str(output_path).

# Frequently Asked Questions

What is the difference between filters.range and filters.expression?

filters.range uses PDAL’s compact interval syntax (Dim[min:max]) and is optimised for simple numeric bounds on one or more dimensions chained with AND. filters.expression accepts a full C-style boolean expression — including OR, parentheses, and arithmetic — and is the right choice when your predicate cannot be expressed as a flat list of intervals.

Why does filters.outlier not delete points directly?

PDAL filters operate on the buffer without shrinking it in-place; instead, filters.outlier marks noise points by setting Classification = 7. A subsequent filters.range with limits=Classification![7:7] performs the actual removal. This preserves reversibility: you can inspect flagged points before discarding them, which is essential for pipeline validation in production workflows.

How do I keep custom attributes through a filter chain?

Declare extra_dims=all in writers.las and verify the dimension list with pipeline.arrays[0].dtype.names after execution. If an intermediate filter drops a custom dimension, use filters.assign to re-inject it with a default value after the offending stage. The attribute mapping guide covers this pattern in depth.

# Related

• Applying Statistical Outlier Filters in PDAL — parameter tuning, radius vs statistical methods, and benchmark comparisons
• PDAL Stage Chaining — buffer-passing model, dimension propagation, and stage ordering rules
• Spatial Reprojection — CRS transformation before spatial filtering
• Pipeline Validation — schema checks, dry-run execution, and CI/CD integration
• PDAL Pipeline Architecture & Execution — parent overview: DAG execution model, stage categories, and production deployment