Programmable Python Filters in PDAL

PDAL ships with well over a hundred stages, and sooner or later your workflow needs the one it does not have. A local quality flag computed from three existing dimensions; a vendor’s intensity normalisation formula; a rejection rule that depends on the acquisition’s flight-line geometry. filters.python is the escape hatch: it hands your function a dictionary of NumPy arrays — one entry per dimension, one element per point in the current buffer — and takes back whatever you put in it. This topic belongs to PDAL Pipeline Architecture and Execution, and it is where a pipeline stops being configuration and starts being code.

The trade is real and worth stating plainly. A programmable filter is the most flexible stage in PDAL and the slowest one available, because every point crosses the C++/Python boundary as part of an array and your function runs in the interpreter. Written with vectorised NumPy it costs perhaps twenty percent over a native stage. Written with a for loop over points it costs two orders of magnitude, and no amount of tuning elsewhere in the pipeline will hide that.

# Prerequisites

The interpreter mismatch is the single most common setup failure. A conda-forge PDAL uses the conda Python; a system PDAL uses the system one. Installing NumPy into the wrong one produces an import error inside the stage that mentions neither.

# Core Workflow Architecture

1. Declare the stage. filters.python takes either a script path or an inline source, plus a function name and a module label used in error messages.
2. Declare added dimensions. Any dimension your function creates must be listed in add_dimension, with an explicit type. PDAL allocates it before the function runs; a key written into outs that was never declared is discarded silently.
3. Receive the buffer. PDAL calls your function once per buffer — once for the whole cloud in standard mode, once per chunk if the surrounding pipeline streams.
4. Mutate or add. Assign into outs to change a dimension or fill a new one. The arrays are the right length by construction; producing a different length is an error.
5. Return True. Returning False tells PDAL the stage failed. There is no partial success.
6. Hand on. The next stage sees the modified buffer with its new dimension, indistinguishable from one a native stage produced.

# Full Implementation

The example computes a per-point quality flag from intensity and return geometry, adds it as a new dimension, and demotes suspect points to a review class rather than deleting them.

"""filters.python stage: derive a QA flag from intensity and return geometry."""
import numpy as np


def flag_quality(ins, outs):
    intensity = ins["Intensity"].astype(np.float64)
    ret = ins["ReturnNumber"].astype(np.int16)
    n_ret = ins["NumberOfReturns"].astype(np.int16)
    z = ins["Z"]

    # Robust intensity bounds from this buffer, not a hard-coded threshold.
    lo, hi = np.percentile(intensity, [1.0, 99.0])
    weak = intensity < lo
    saturated = intensity > hi

    # A last return that is also the only return over a rough surface is the
    # most trustworthy ground candidate; an intermediate return is the least.
    intermediate = (ret > 1) & (ret < n_ret)

    # Elevation blunders: more than five robust deviations from the median.
    med = np.median(z)
    mad = np.median(np.abs(z - med)) * 1.4826
    blunder = np.abs(z - med) > (5.0 * mad) if mad > 0 else np.zeros_like(z, dtype=bool)

    flag = np.zeros(len(z), dtype=np.uint8)
    flag[weak] |= 1
    flag[saturated] |= 2
    flag[intermediate] |= 4
    flag[blunder] |= 8

    outs["QaFlag"] = flag

    # Anything with a blunder bit goes to class 12 (overlap/reserved for review)
    # rather than being deleted, so the decision stays auditable.
    cls = ins["Classification"].copy()
    cls[blunder] = 12
    outs["Classification"] = cls

    return True

Wired into a pipeline:

{
  "pipeline": [
    {"type": "readers.las", "filename": "tile_0431.laz"},
    {
      "type": "filters.python",
      "script": "qa_flag.py",
      "function": "flag_quality",
      "module": "qa",
      "add_dimension": "QaFlag=uint8"
    },
    {
      "type": "writers.las",
      "filename": "tile_0431_qa.laz",
      "compression": "laszip",
      "minor_version": 4,
      "dataformat_id": 6,
      "extra_dims": "QaFlag=uint8",
      "forward": "all"
    }
  ]
}

# Code Breakdown

Every operation is vectorised. There is no for loop over points anywhere. np.percentile, boolean masks and in-place bitwise assignment all run at C speed over the whole array; the Python interpreter executes about twenty statements regardless of whether the buffer holds ten points or ten million.

Thresholds are derived from the buffer, not hard-coded. Intensity scaling varies between sensors and even between flight lines, so a fixed cutoff transfers badly. Percentiles adapt. The cost is that in streaming mode each chunk computes its own percentiles, which is a real difference in behaviour — see the gotcha below.

A median absolute deviation is used instead of a standard deviation. A cloud with a handful of returns at 4,000 m has a standard deviation dominated by exactly the points you are trying to find. The MAD does not move.

Suspect points are reclassified, not dropped. Deleting evidence inside a QA stage makes the result unauditable. A downstream filters.range can remove class 12 when the pipeline actually wants it gone.

add_dimension and extra_dims are both required. The first makes the dimension exist inside the pipeline; the second makes the writer put it in the file. Omitting the second is the commonest reason a computed dimension “disappears”, as attribute mapping explains at length.

# Parameter Reference Table

# Validation and Integrity Checks

The dimension exists and is populated. After the run, read it back and confirm it is not uniformly zero:

p = pdal.Pipeline(json.dumps({"pipeline": [{"type": "readers.las", "filename": "tile_0431_qa.laz"}]}))
p.execute()
arr = p.arrays[0]
assert "QaFlag" in arr.dtype.names, "dimension was declared but never written to the file"
assert arr["QaFlag"].any(), "every point has flag zero — the rule matched nothing"

The point count did not change. A programmable filter that assigns dimensions must not alter the buffer length. If the output count differs from the input, the function returned arrays of the wrong shape.

The type survived. A uint8 dimension that comes back as float64 means extra_dims declared a different type from add_dimension, and every consumer downstream now reads eight bytes where one was intended.

# Performance Tuning

The single decision that matters is vectorisation, and the table above quantifies it. Beyond that, three smaller levers apply.

• Compute once, not per chunk. Anything that does not depend on the points — loading a lookup table, opening a raster — belongs at module scope, not inside the function, because the function runs once per buffer.
• Prefer views to copies. ins["Z"] is a view; ins["Z"].copy() is a copy of a potentially enormous array. Copy only when you need to keep the original.
• Consider whether a native stage already does it. filters.assign with an expression covers a surprising amount of what people write Python for, at native speed and with no interpreter in the pipeline at all.

# Shipping a Programmable Filter to Production

A filters.python stage turns a pipeline into a program with a dependency, and it needs the same handling any other dependency gets.

Version the script with the pipeline JSON. They are one artefact. A pipeline that references qa_flag.py by a relative path and a script that lives in someone’s home directory is a job that works on one machine. Put both in the repository, reference the script relative to the pipeline file, and record a hash of each in the output’s provenance record.

Pin the interpreter, not just the packages. The behaviour of np.percentile on ties has changed across NumPy releases, and a QA flag that shifts by a percent between runs is worse than one that is slightly wrong consistently. The container is the right unit here — the PDAL Docker containers guide covers pinning PDAL, GDAL, PROJ and the Python stack together so none of them drifts independently.

Parameterise through pdalargs rather than editing the script. Thresholds, sensor identifiers and file paths belong in the pipeline JSON, which changes per run, not in the function, which should not. Inside the function they arrive as a dictionary:

def flag_quality(ins, outs):
    args = globals().get("pdalargs", {})
    lo_pct = float(args.get("intensity_low_percentile", 1.0))
    ...

Test the function without PDAL. The signature is just two dictionaries of arrays, so a unit test can build them with NumPy and call the function directly. That test runs in milliseconds, needs no LiDAR data, and catches the majority of logic errors before a pipeline is involved at all:

def test_blunder_is_flagged():
    n = 1000
    ins = {
        "Z": np.concatenate([np.full(n - 1, 100.0), [4000.0]]),
        "Intensity": np.full(n, 3000, dtype=np.uint16),
        "ReturnNumber": np.ones(n, dtype=np.uint8),
        "NumberOfReturns": np.ones(n, dtype=np.uint8),
        "Classification": np.ones(n, dtype=np.uint8),
    }
    outs = {}
    assert flag_quality(ins, outs) is True
    assert outs["QaFlag"][-1] & 8, "the 4000 m return should carry the blunder bit"
    assert outs["Classification"][-1] == 12

Fail loudly inside the function. A stage that catches every exception and returns True produces a file full of zeros and a run that reports success. Let the exception propagate; PDAL will surface it with your module label attached, which is what the label is for.

# Common Errors and Troubleshooting

Unable to import module. PDAL is using a different interpreter from the one you installed NumPy into. pdal --debug prints the Python it links against.

The new dimension is missing from the output file. add_dimension was set but extra_dims on the writer was not, or the LAS version cannot carry it. Write LAS 1.4 with point format 6 or higher.

Invalid dimension at pipeline start. A key written into outs that PDAL does not know about. Every created dimension needs a matching add_dimension entry.

Results differ between runs on the same tile. Almost always a threshold derived from buffer statistics combined with streaming, so each chunk computes different percentiles. Either compute the statistics in a prior pass and pass them in through pdalargs, or accept that the stage is chunk-dependent and document it.

The pipeline is no longer streamable. filters.python itself can stream, but a function that needs whole-cloud statistics cannot honestly do so. That tension is the subject of streaming mode execution.

# Frequently Asked Questions

Why is my computed dimension missing from the output file?

Because add_dimension and extra_dims do different jobs. The first makes the dimension exist inside the pipeline; the second makes the writer store it. You need both, and the LAS version has to be able to carry extra bytes — point format 6 or above in a 1.4 file.

How much slower is filters.python than a native stage?

About twenty percent if the function is vectorised NumPy, and roughly a hundred times slower if it loops over points in Python. The interpreter overhead is per call, not per point, so the whole question is whether your function makes one call or ten million.

Can filters.python stream?

The stage itself can, and PDAL will call your function once per chunk. Whether that is correct depends on the function: one that computes a threshold from buffer statistics gets different thresholds per chunk. If the rule needs whole-cloud context, compute the statistics in a prior pass and pass them in through pdalargs.

Why does PDAL say it cannot import numpy?

PDAL links against a specific Python interpreter, which is often not the one first on your PATH. Run pdal --debug to see which one, then install numpy into that interpreter. This is the most common setup failure with programmable filters and the error message does not name the interpreter.

# Related

• PDAL Pipeline Architecture and Execution — the section this stage type belongs to
• Writing a filters.python Stage with NumPy — the hands-on walkthrough from empty file to working stage
• Adding a New Dimension from a Python Filter — add_dimension, extra_dims and getting the type right
• Debugging and Profiling filters.python — finding out why the stage is slow or silent
• Attribute Mapping — how dimensions move through a pipeline and into a file
• Streaming Mode Execution in PDAL — what changes when your function is called once per chunk