Cropping a Point Cloud to a Polygon Boundary

TL;DR: Read the boundary with GeoPandas, reproject it to the point cloud’s CRS, pass its WKT to filters.crop as "polygon" with a matching "a_srs", and write the result. For many polygons in one pass, pass a list of WKT strings; each produces its own output view. For COPC inputs, give the polygon to readers.copc instead so only intersecting nodes are fetched.

# Context and Motivation

This guide is part of Pipeline Filtering Logic. Clipping to a boundary is one of the most frequent LiDAR operations: deliver only the client’s project area, extract a corridor around a road or pipeline, cut out a parcel for a site survey, or remove a neighbouring jurisdiction’s data before publishing. PDAL’s filters.crop does it in a streaming, per-point test, so it works on files of any size. What goes wrong is almost never the crop itself; it is the polygon — in the wrong CRS, invalid, multipart when a single part was expected, or so detailed that the test is slow.

# Prerequisites and Assumptions

• PDAL 2.x; GeoPandas and Shapely 2.x for reading and preparing polygons.
• A boundary layer (GeoPackage, Shapefile, GeoJSON) and the point cloud’s CRS, known from pdal info --metadata.
• Valid polygons. Self-intersections make point-in-polygon tests undefined; fix them with shapely.make_valid first.

# Step-by-Step Implementation

# Step 1 — Load and validate the boundary

Read the layer, make geometries valid, and dissolve to one geometry if you want a single output.

# Step 2 — Reproject the polygon, not the points

Transforming a polygon is instant; transforming millions of points to match a polygon’s CRS is not. Reproject the boundary to the point cloud’s CRS with to_crs.

# Step 3 — Simplify if the boundary is very detailed

A coastline with a hundred thousand vertices makes every point test slow. Simplify to a tolerance well below the point spacing — 0.1 m is usually invisible in results.

# Step 4 — Crop

Pass the WKT to filters.crop’s polygon option with a_srs set to the same CRS as the points, so PDAL does not guess.

# Step 5 — Crop many polygons in one pass

A list of polygons produces one output view per polygon; a writer with # in its filename writes each to its own file.

# Complete Working Example

"""Crop a LAZ tile to each parcel in a GeoPackage, one output per parcel."""
from __future__ import annotations

import json
from pathlib import Path

import geopandas as gpd
import pdal
import shapely

SRC = Path("tiles/t_0431.laz")
CLOUD_CRS = "EPSG:6347"


def parcels_wkt(gpkg: Path, layer: str, simplify_m: float = 0.1) -> list[tuple[str, str]]:
    gdf = gpd.read_file(gpkg, layer=layer).to_crs(CLOUD_CRS)
    gdf["geometry"] = shapely.make_valid(gdf.geometry.values).simplify(simplify_m)
    info = pdal.Pipeline(json.dumps({"pipeline": [str(SRC)]})).quickinfo["readers.las"]["bounds"]
    tile = shapely.box(info["minx"], info["miny"], info["maxx"], info["maxy"])
    gdf = gdf[gdf.intersects(tile)]
    return [(str(r.parcel_id), r.geometry.wkt) for r in gdf.itertuples()]


def crop(parcels: list[tuple[str, str]], out_dir: Path) -> int:
    out_dir.mkdir(parents=True, exist_ok=True)
    total = 0
    for pid, wkt in parcels:
        spec = {"pipeline": [
            str(SRC),
            {"type": "filters.crop", "polygon": wkt, "a_srs": CLOUD_CRS},
            {"type": "writers.las", "filename": str(out_dir / f"parcel_{pid}.laz"),
             "minor_version": 4, "dataformat_id": 6, "forward": "all"},
        ]}
        n = pdal.Pipeline(json.dumps(spec)).execute()
        print(f"parcel {pid}: {n:,} points")
        total += n
    return total


if __name__ == "__main__":
    crop(parcels_wkt(Path("parcels.gpkg"), "parcels"), Path("out/parcels"))

When the parcels are many and small, one pipeline with a list of polygons avoids reading the tile once per parcel:

{
  "pipeline": [
    "tiles/t_0431.laz",
    { "type": "filters.crop", "a_srs": "EPSG:6347",
      "polygon": ["POLYGON ((431120 4471200, 431180 4471200, 431180 4471260, 431120 4471260, 431120 4471200))",
                  "POLYGON ((431300 4471420, 431370 4471420, 431370 4471490, 431300 4471490, 431300 4471420))"] },
    { "type": "writers.las", "filename": "out/parcels/parcel_#.laz" }
  ]
}

# Key Parameter Table

# Verification

• Every point inside. Test output points against the polygon with Shapely on a sample: all should be contained or on the boundary.
• No points lost inside. Crop with outside: true too; the two outputs’ counts should sum to the input count.
• Bounds shrink. The header bounds of the output should fit within the polygon’s bounds.

import numpy as np
from shapely import points, contains_xy

p = pdal.Pipeline(json.dumps({"pipeline": ["out/parcels/parcel_1042.laz"]})); p.execute()
a = p.arrays[0]
poly = shapely.from_wkt(dict(parcels_wkt(Path("parcels.gpkg"), "parcels"))["1042"])
inside = contains_xy(poly.buffer(0.01), a["X"], a["Y"])
assert inside.all(), f"{(~inside).sum()} points outside the parcel"

# Gotchas and Edge Cases

CRS mismatch produces empty output. A polygon in EPSG:4326 degrees tested against UTM metres selects nothing, and PDAL may not warn. Always reproject the polygon and set a_srs.

Multipart and holes. MULTIPOLYGON and polygons with interior rings work; points inside holes are excluded. Dissolving a layer can create holes where parcels do not touch — check that is intended.

COPC and EPT sources. Cropping after a full read of a remote COPC file downloads everything. Put the polygon on readers.copc ("polygon": wkt) so only intersecting octree nodes are fetched; see querying a COPC file by bounds and resolution.

Boundary points. Points exactly on the polygon edge may fall either way. When adjacent parcels must not share or lose points, crop with one polygon layer that tiles the area exactly and check that per-parcel counts sum to the total.

# Frequently Asked Questions

How do I crop a LAS file to a shapefile polygon with PDAL?

Read the shapefile with GeoPandas, reproject it to the LAS file’s CRS, and pass the polygon’s WKT to filters.crop with the polygon option and a matching a_srs. Write the result with writers.las.

Why is my cropped output empty?

Almost always a CRS mismatch: the polygon’s coordinates are in a different system from the points. Reproject the polygon to the point cloud’s CRS and set a_srs explicitly.

Can I crop to many polygons at once?

Yes. Pass a list of WKT polygons to filters.crop. Each polygon produces its own output view, and a writer with a # in its filename writes each view to a separate file.

How do I keep points outside the polygon instead?

Set outside to true on filters.crop. The filter then keeps everything except the points inside the polygon.

# Related

• Pipeline Filtering Logic — filtering in general
• Filtering Points with filters.expression — attribute filters
• Querying a COPC File by Bounds and Resolution — spatial queries on cloud-native data
• Buffered Tiling to Avoid Edge Artefacts — cropping back after buffered processing
• Fixing CRS Mismatches in Point Clouds — the usual cause of empty crops