Clipping Point Clouds to CAD Boundaries in Python

Q: How do I keep the classification and intensity?

Clip the point records rather than a coordinate array. laspy lets you index the chunk's point record with a boolean mask, which keeps every dimension the format carries. Extracting XYZ into a numpy array and writing that back produces a file with geometry and nothing else.

To clip a point cloud to a CAD site boundary, read the closed polyline from the drawing, bring it into the cloud’s coordinate system, then stream the cloud in chunks applying a bounding-box prefilter followed by a vectorised point-in-polygon test, writing survivors into a new file that inherits the source header. Nearly every failure here is a coordinate-system mismatch rather than a geometry problem. This page is part of Point Cloud and Reality Capture Integration.

How the Two Coordinate Worlds Meet

A site boundary in a drawing is a closed LWPOLYLINE in drawing units on a site grid. A point cloud is metres on a projected coordinate reference system. Between them sit the two transformations this site covers at length: a unit scale from the drawing header, and a reprojection from the site grid onto the projection.

Both have to be applied to the boundary, not to the cloud. Transforming a hundred million points to meet a polygon is orders of magnitude more expensive than transforming a polygon of thirty vertices to meet the points, and it loses precision in the process.

The boundary also has to be a valid ring before it can test anything. A closed polyline stores closure as a flag rather than as a repeated coordinate, and it may carry bulge arcs that must be flattened. Both are covered in Extracting LWPOLYLINE Vertices with ezdxf; the clip inherits them.

Production-Ready Script

# laspy[lazrs]>=2.5, ezdxf>=1.1.0, shapely>=2.0, pyproj>=3.5, numpy>=1.24
from __future__ import annotations

import numpy as np
import laspy
import ezdxf
from shapely.geometry import Polygon
from shapely import contains_xy
from pyproj import Transformer

INSUNITS_TO_M = {1: 0.0254, 2: 0.3048, 4: 0.001, 5: 0.01, 6: 1.0, 7: 1000.0}


def boundary_polygon(dxf_path: str, layer: str, dst_epsg: int, src_epsg: int,
                     *, sag: float = 0.02) -> Polygon:
    """The largest closed polyline on that layer, in the cloud CRS, in metres."""
    doc = ezdxf.readfile(dxf_path)
    code = doc.header.get("$INSUNITS", 0)
    if code not in INSUNITS_TO_M:
        raise ValueError(f"$INSUNITS={code} is undefined — resolve the unit explicitly")
    scale = INSUNITS_TO_M[code]

    best: Polygon | None = None
    for pl in doc.modelspace().query(f'LWPOLYLINE[layer=="{layer}"]'):
        if not pl.closed:
            continue
        pts = [(v.x * scale, v.y * scale) for v in pl.flattening(distance=sag / scale)]
        if len(pts) < 3:
            continue
        ring = Polygon(pts)
        if best is None or ring.area > best.area:
            best = ring
    if best is None:
        raise ValueError(f"no closed polyline found on layer {layer!r}")

    t = Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)
    x, y = t.transform(*np.array(best.exterior.coords).T)
    return Polygon(np.column_stack((x, y)))


def clip(las_path: str, out_path: str, poly: Polygon, *, chunk: int = 2_000_000) -> dict:
    minx, miny, maxx, maxy = poly.bounds
    kept = seen = 0
    with laspy.open(las_path) as reader:
        header = reader.header
        with laspy.open(out_path, mode="w", header=header) as writer:
            for points in reader.chunk_iterator(chunk):
                seen += len(points)
                x, y = np.asarray(points.x), np.asarray(points.y)
                # Cheap exact prefilter: the envelope never excludes an inside point.
                envelope = (x >= minx) & (x <= maxx) & (y >= miny) & (y <= maxy)
                if not envelope.any():
                    continue
                mask = np.zeros(len(points), dtype=bool)
                mask[envelope] = contains_xy(poly, x[envelope], y[envelope])
                if mask.any():
                    writer.write_points(points[mask])   # keeps every dimension
                    kept += int(mask.sum())
    return {"read": seen, "written": kept, "retained": kept / seen if seen else 0.0}


if __name__ == "__main__":
    poly = boundary_polygon("site.dxf", layer="SITE-BOUNDARY",
                            dst_epsg=27700, src_epsg=27700)
    print(clip("survey.laz", "survey_clipped.laz", poly))

Key implementation notes:

points[mask] indexes the point record, so classification, intensity, return number and colour all survive. Rebuilding from an XYZ array silently discards them.
The writer is opened with the source header, so scale, offset and the coordinate reference system carry into the output.
contains_xy is Shapely 2’s vectorised predicate — one call for an entire array rather than a Python loop over points.
The bounding-box prefilter is applied first and the polygon test only to survivors. It is exact, so this is pure saving.
The retention ratio is returned. A clip that retains 0.02% is usually a coordinate mismatch rather than a small site.

Compatibility Matrix

Component	Supported range	Notes
`laspy`	`>=2.5`	chunked read and write, point-record indexing
`ezdxf`	`>=1.1.0`	`flattening` for bulge arcs
`shapely`	`>=2.0`	`contains_xy` vectorised predicate
`pyproj`	`>=3.5`	boundary reprojection
Output format	LAS or LAZ	LAZ output needs the compression backend

Fallback Strategies

1. Zero points retained. Print poly.bounds against the cloud header bounds. A site-grid polygon near the origin against a projected cloud is the usual answer, and it is visible at a glance.

2. $INSUNITS is 0. The code raises rather than guessing. Resolve the drawing unit explicitly — the policy is set out in Autoscaling DXF Geometry from $INSUNITS in Python.

3. The boundary is not closed. A drafted boundary is frequently several open polylines that visually meet. Merge them into a ring before polygonising, and reject the result if it does not close within tolerance rather than silently taking the largest fragment.

4. Holes in the boundary. Exclusion zones drawn as separate rings inside the site boundary are not automatically holes. Build the polygon with explicit interior rings, or the clip retains points inside them.

5. Output is much larger than expected. LAZ output written without a compression backend falls back to LAS. Assert the output extension against the driver actually used.

FAQ

Why does my clip return zero points?

Almost always a coordinate mismatch. The drawing is in millimetres on a site grid and the cloud is in metres on a projected CRS, so the boundary polygon sits a few metres from the origin while the cloud is half a million metres away. Print both bounding boxes before clipping; the discrepancy is immediately obvious and never subtle.

Is a bounding-box prefilter worth it?

Substantially. Point-in-polygon is far more expensive per point than four array comparisons, and on a typical site the envelope discards most of an airborne tile before the polygon test runs at all. The prefilter is exact — it never discards a point inside the polygon — so it costs nothing in correctness.

How do I keep the classification and intensity?

Clip the point records rather than a coordinate array. laspy lets you index the chunk’s point record with a boolean mask, which keeps every dimension the format carries. Extracting XYZ into a numpy array and writing that back produces a file with geometry and nothing else.

Point Cloud and Reality Capture Integration — parent reference on coordinate metadata and density
Reading LAS and LAZ Files with laspy — the chunked read this clip is built on
Extracting LWPOLYLINE Vertices with ezdxf — reading the boundary polyline out of the drawing