PDAL vs Open3D for Point Cloud Filtering: Choosing the Right Tool

Both PDAL and Open3D can strip outliers, thin density, and classify ground on a LiDAR scan, but they are built on opposite philosophies: PDAL is a declarative, CRS-aware pipeline engine that reads .laz headers and streams tiles through a JSON stage graph, while Open3D is an in-memory research library that hands you a numpy-backed PointCloud and expects you to manage coordinates and metadata yourself. This page decides between them for a filtering step — comparing outlier removal, ground classification, and voxel downsampling in each, their behaviour on large .laz, and how each treats a coordinate reference system such as EPSG:32618 (UTM zone 18N) — and ends with a verdict on when to reach for one, the other, or both.

You hit this decision the moment your filtering step outgrows a one-off script. If you already have both installed after reading the broader point cloud filtering techniques guide, the question is no longer how to run statistical outlier removal but which engine owns the stage — because the choice dictates whether your CRS is preserved automatically, whether a billion-point survey fits in RAM, and whether the run is reproducible from a committed JSON file or a Python notebook.

Prerequisites

• Python 3.10+ with both stacks: conda install -c conda-forge pdal python-pdal (PDAL’s native GDAL/PROJ chain installs far more reliably through conda) and pip install "open3d>=0.18" "laspy[lazrs]>=2.5" "numpy>=1.24" "pyproj>=3.6".
• A test scan in LAS/LAZ (ASPRS point format 6 or 7) carrying a valid CRS in its header. The worked examples use EPSG:32618 (UTM zone 18N, metres) throughout; every distance parameter below is in metres because of it.
• Confirm the PDAL stage registry: pdal --drivers | grep filters.outlier should list the outlier filter, and pdal --version should report a PROJ build so CRS reprojection stages resolve.

The core difference: declarative pipeline vs in-memory API

PDAL models filtering as a directed graph of stages serialised to JSON. You name a reader, a chain of filters.* stages, and a writer; PDAL streams points through them, and — critically — the GDAL/PROJ machinery underneath reads the source EPSG straight from the LAS header and can reproject inside the same pipeline. Nothing about the CRS is your job unless you want to change it. The whole run is a text document you commit to version control, which is why PDAL is the natural fit for a CI-driven ingestion stage.

Open3D takes the opposite stance. o3d.io.read_point_cloud and the LAS-via-laspy path both discard georeferencing; you get a PointCloud wrapping a numpy array of raw XYZ and nothing else. That is a feature when your next step is normal estimation, registration, or reconstruction — the data is already a numpy array you can hand to any scientific library — but it means the EPSG code lives only in your head (or a sidecar you write yourself), and a whole survey must fit in memory because there is no streaming layer.

PDAL: one declarative pipeline, CRS carried end to end

The following pipeline reads a .laz in EPSG:32618, runs statistical then radius outlier removal, classifies ground with SMRF, voxel-thins, and writes a LAZ with the CRS re-stamped — all without a line of coordinate handling on your side.

import pdal
import json

pipeline = pdal.Pipeline(json.dumps({
    "pipeline": [
        {"type": "readers.las", "filename": "block_utm18n.laz"},
        {"type": "filters.outlier", "method": "statistical",
         "mean_k": 20, "multiplier": 2.0},
        {"type": "filters.outlier", "method": "radius",
         "radius": 0.25, "min_k": 8},
        {"type": "filters.range", "limits": "Classification![7:7]"},  # drop flagged noise
        {"type": "filters.smrf", "scalar": 1.2, "slope": 0.2,
         "threshold": 0.45, "window": 16.0},
        {"type": "filters.voxelcentroid", "cell": 0.10},
        {"type": "writers.las", "filename": "block_filtered.laz",
         "a_srs": "EPSG:32618", "compression": "laszip"}
    ]
}))
n_points = pipeline.execute()
meta = pipeline.metadata           # stage-by-stage counts, bounds, SRS
print(f"{n_points:,} points written in EPSG:32618")

filters.outlier flags outliers by setting Classification to 7 rather than deleting them, so the subsequent filters.range is what actually removes them — an audit-friendly separation the declarative model makes explicit. The a_srs on the writer guarantees the output header carries EPSG:32618 forward.

Open3D: numpy in memory, CRS is yours to track

Open3D cannot read .laz, so laspy loads it and you build the cloud from a numpy array, recording the EPSG yourself. The same filters are one method call each, returning the surviving cloud and the kept indices.

import laspy
import numpy as np
import open3d as o3d

SRC_EPSG = "EPSG:32618"                    # you record this — Open3D will not
las = laspy.read("block_utm18n.laz")
xyz = np.vstack((las.x, las.y, las.z)).T.astype(np.float64)

# Shift to a local origin so KD-tree math stays precise on 6-figure UTM eastings.
origin = xyz.min(axis=0)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz - origin)

pcd, keep_sor = pcd.remove_statistical_outlier(nb_neighbors=20, std_ratio=2.0)
pcd, keep_rad = pcd.remove_radius_outlier(nb_points=8, radius=0.25)   # metres, EPSG:32618
pcd = pcd.voxel_down_sample(voxel_size=0.10)

out_xyz = np.asarray(pcd.points) + origin  # restore absolute EPSG:32618 coordinates
print(f"{len(out_xyz):,} points; EPSG {SRC_EPSG} tracked out of band")

Notice what Open3D forces you to do that PDAL did for free: subtract a local origin to protect float precision, and thread SRC_EPSG through by hand. Notice also what Open3D gives you that PDAL does not — keep_sor and keep_rad are numpy index arrays you can immediately apply to intensity or RGB, and pcd is one call away from estimate_normals for reconstruction.

Ground classification: PDAL wins outright

Ground/bare-earth separation is the sharpest capability gap. PDAL ships filters.smrf (Simple Morphological Filter), filters.pmf, and filters.csf as first-class stages that write ASPRS class 2 into the point record. Open3D has no ground filter at all — you would reimplement a morphological or cloth-simulation algorithm on the numpy array yourself, or shell out to PDAL. If your filtering step needs a DTM or a vegetation strip, that alone decides it.

# PDAL: ground stays a one-liner in the graph.
ground = pdal.Pipeline(json.dumps({
    "pipeline": [
        "block_filtered.laz",
        {"type": "filters.smrf", "slope": 0.2, "window": 18.0, "threshold": 0.45},
        {"type": "filters.range", "limits": "Classification[2:2]"},
        {"type": "writers.las", "filename": "ground_dtm.laz", "a_srs": "EPSG:32618"}
    ]
}))
ground.execute()

Performance on large .laz

PDAL streams. Stages such as filters.range, filters.voxelcentroid, and the outlier filters run in a streaming mode that processes points in fixed-size buffers, and filters.splitter or filters.chipper tile a multi-billion-point survey so nothing exceeds memory. A national-grid .laz that would never fit in RAM flows through a PDAL pipeline on a modest machine.

Open3D holds everything resident. A PointCloud over 50 M points already wants several gigabytes, and remove_statistical_outlier builds a KD-tree over the whole array, so you must tile the cloud yourself — read with laspy.open(...).chunk_iterator(n), filter each chunk, buffer the tile edges, and merge — reproducing by hand the chunking PDAL does declaratively. Where PDAL is memory-bounded by design, Open3D is memory-bounded by your discipline.

Decision table

Verdict

Choose PDAL when the filtering step is part of an ingestion pipeline: it reads and preserves the CRS automatically, ground-classifies natively, streams .laz far larger than RAM, and reduces to a JSON document you gate in CI. Choose Open3D when filtering is the front of a geometry-processing chain that stays in memory — outlier removal followed by normal estimation and surface reconstruction — where having the cloud as a live numpy array and the kept-index arrays for attribute alignment is worth managing the EPSG yourself.

In practice most production twins use both, and the seam is clean: let PDAL own ingestion — read the header, reproject to EPSG:32618, classify ground, clip, and tile — then hand each tile to Open3D for the in-memory geometry work. The terrestrial LiDAR noise-removal walkthrough is a pure-Open3D example of that second half; a PDAL splitter feeding it is the pairing that scales to a city. The one rule that spans both tools is to keep every distance parameter in a metric CRS — a radius of 0.25 is 25 cm in EPSG:32618 and a meaningless 0.25 degrees in EPSG:4326.

Expected Output & Verification

Run both paths on the same scan and the surviving point counts should agree to within a couple of percent; the small gap comes from the two libraries computing the statistical cutoff slightly differently, not from a bug. Report the retained ratio from each so a parameter drift fails loudly:

def retained(before, after, tag):
    print(f"{tag:14s} {before:>12,} -> {after:>12,}  ({after/before:.1%})")

retained(len(las.points), n_points, "PDAL")
retained(len(las.points), len(out_xyz), "Open3D")

assert abs(n_points - len(out_xyz)) / len(las.points) < 0.05, \
    "PDAL and Open3D diverged >5% — check radius units and CRS"

Sample console output for a 12 M-point urban tile:

PDAL            12,140,551 ->   11,402,338  (93.9%)
Open3D          12,140,551 ->   11,318,090  (93.2%)

Confirm the CRS actually survived the PDAL side — the most common silent defect is a writer without a_srs:

pdal info block_filtered.laz --metadata | grep -i srs

Common Errors

RuntimeError: Unable to convert ... no SRS or an empty output SRS from PDAL. The readers.las found no CRS in the source header, or the writers.las omitted a_srs. PDAL cannot reproject or stamp a frame it was never given. Set "a_srs": "EPSG:32618" on the writer, and if the source header is blank, add "spatialreference": "EPSG:32618" to the reader so every downstream stage — and every distance filter — runs in metres.

RuntimeError: [Open3D ERROR] ... Unknown file extension for file (block_utm18n.laz). Open3D has no LAS/LAZ reader. Load with laspy.read and build the PointCloud from a numpy array as shown; never call o3d.io.read_point_cloud on a .laz directly.

Open3D radius filter removes almost everything (or nothing). Either the cloud is still in geographic EPSG:4326 so a radius=0.25 “metre” is really 0.25 degrees, or the local-origin shift was skipped and large UTM eastings lost float precision in the KD-tree. Reproject to EPSG:32618 and subtract origin before filtering — the step PDAL performs for you and Open3D does not.

Related Guides

• Point Cloud Filtering Techniques — the full SOR, radius, voxel, and ground toolkit both libraries implement
• Removing Noise from Terrestrial LiDAR Scans — a pure-Open3D noise-removal walkthrough for the in-memory half
• Surface Reconstruction for Geospatial Twins — where an Open3D-filtered cloud goes next
• Coordinate Reference Systems for 3D Assets — choosing the metric EPSG both tools depend on

Back to Point Cloud Filtering Techniques.