Clip and Tile Point Clouds in PyQGIS

A LiDAR delivery arrives tiled the way the supplier's flight lines happened to fall, which is almost never the way your work is organised. You want the survey cut to a catchment boundary, or split into pieces small enough to process in parallel, or merged into a single logical layer so a script does not have to know how many files there are. All three are one algorithm call each — the interesting part is the order you do them in and what you keep.

This recipe belongs to Point Cloud & LiDAR Workflows in PyQGIS. It covers clipping to a polygon, retiling with a buffer, building a virtual point cloud manifest, and structuring the loop so a four-hour job survives being interrupted.

From the supplier's tiling to yoursThe delivery is a set of tiles laid out along flight lines. Clipping to the catchment boundary discards the tiles and parts of tiles outside it. Retiling lays the remaining points onto a regular grid of equal-sized pieces. A virtual point cloud manifest then presents those pieces to QGIS as a single layer.Clip first, retile second, manifest lastas deliveredstrips along the flight linesclipped to the boundaryeverything outside is goneretiled on a gridequal pieces, parallel-friendlysurvey.vpc — one layer over all six piecesclipping before retiling means retiling has far less to move

Prerequisites

  • QGIS 3.34 LTR or newer with the PDAL provider (virtual point clouds need 3.34).
  • Indexed clouds — see loading a point cloud layer.
  • A boundary polygon in the same CRS as the cloud. These algorithms do not reproject the overlay for you.

Clip to a boundary

import processing

clipped = processing.run("pdal:clip", {
    "INPUT": "/data/indexed/survey.vpc",
    "OVERLAY": "/data/vector/catchment.gpkg",
    "OUTPUT": "/data/clip/catchment.copc.laz",
})["OUTPUT"]

Breakdown: OVERLAY takes a polygon layer or a path to one; every polygon in it contributes, so a multi-feature layer clips to the union rather than producing one output per feature. Passing the .vpc as input is what makes this fast on a large survey: the manifest carries per-tile extents, so tiles that do not intersect the catchment are never opened. Writing to .copc.laz keeps the result indexed. If the boundary and the cloud are in different projections the result is silently empty, which is the single most common failure here — check layer.crs().authid() on both before running.

Where the area of interest is a rectangle rather than a polygon, the EXTENT parameter on pdal:filter does the same job in one pass alongside a class filter, which saves a full read.

from qgis.core import QgsPointCloudLayer

boundary = QgsVectorLayer("/data/vector/catchment.gpkg", "catchment", "ogr")
cloud = QgsPointCloudLayer("/data/indexed/survey.vpc", "survey", "pdal")

if boundary.crs() != cloud.crs():
    raise SystemExit(
        f"CRS mismatch: boundary {boundary.crs().authid()} vs cloud {cloud.crs().authid()}"
    )

Breakdown: A three-line guard at the top of a script converts a silent empty output into an immediate, readable failure. Comparing QgsCoordinateReferenceSystem objects directly is correct — the equality operator compares the definitions, not the authids, so it also catches a cloud whose CRS is a custom definition equivalent to the boundary's but labelled differently. Where they genuinely differ, reproject the boundary, which is a handful of vertices, rather than the cloud, which is a billion points.

Retile onto a regular grid

Flight-line tiles vary wildly in size — a strip over open country may hold ten times the points of one over water. Retiling gives you pieces of predictable cost, which is what a parallel run needs.

processing.run("pdal:tile", {
    "INPUT": clipped,
    "LENGTH": 500,
    "OUTPUT": "/data/tiles",
})

Breakdown: LENGTH is the tile edge in map units, so 500 gives 500 m squares. The OUTPUT here is a directory, not a file, which is unusual among Processing algorithms and easy to miss — pointing it at a filename produces a confusing error. Choosing the length is a memory decision: aim for tiles that hold a few tens of millions of points, which on a typical 4 pts/m² survey means 500 m to 1 km squares. Smaller tiles mean more files and more edge effects; larger ones mean a single worker holds more in memory.

Tile edges need an overlapA neighbourhood operation applied tile by tile sees a truncated neighbourhood at each edge, producing a visible seam. Reading each tile with a buffer of surrounding points, computing over the buffered area and then trimming the result back to the tile boundary removes the seam entirely.Anything that looks at neighbours needs points beyond the edgeno buffera seam along every shared edgeread buffered, write trimmeddashed: what is read — solid: what is kepta buffer wider than the largest neighbourhood radius you will use

Simple per-point operations — a class filter, a format conversion — need no buffer. Anything that looks at neighbouring points does: gridding with an interpolated statistic, a TIN export, a density surface. The practical approach is to grid each tile with a buffered read and trim the resulting raster back to the tile boundary before mosaicking, exactly as with any windowed raster operation.

There is one more reason to buffer that has nothing to do with seams. Several PDAL stages estimate a local property — a plane fit, a nearest-neighbour distance, a local density — and their estimate degrades near a boundary even when the output raster is later trimmed, because the estimate itself was computed from a truncated neighbourhood. Trimming hides the symptom in the middle of the tile and leaves a band of subtly wrong values just inside the trim line. The fix is the same but the buffer has to be generous: at least twice the radius the stage uses, not just enough to cover the output cells.

Build a manifest over the result

import glob

pieces = sorted(glob.glob("/data/tiles/*.copc.laz"))

processing.run("pdal:virtualpointcloud", {
    "LAYERS": pieces,
    "OUTPUT": "/data/tiles/catchment.vpc",
})

Breakdown: LAYERS takes the list of member paths, and the algorithm reads each header to record extents and counts. The manifest is a few kilobytes regardless of how much data it describes, so rebuilding it is cheap and there is no reason to keep a stale one. Open the result with the pdal provider — a .vpc is not itself an octree, so copc will refuse it. From that point on, every recipe on this site that takes a cloud takes the manifest instead, and reads only the pieces it needs.

Keeping a long job restartable

A survey-wide clip and retile is measured in hours, and hours of work should not be lost to a full disk at the ninetieth tile.

import os

log_path = "/data/tiles/done.txt"
done = set()
if os.path.exists(log_path):
    done = set(open(log_path).read().split())

for tile in sorted(glob.glob("/data/indexed/*.copc.laz")):
    name = os.path.basename(tile)
    if name in done:
        continue
    try:
        processing.run("pdal:clip", {
            "INPUT": tile,
            "OVERLAY": "/data/vector/catchment.gpkg",
            "OUTPUT": f"/data/clip/{name}",
        })
    except Exception as error:
        print("failed", name, error)
        continue
    with open(log_path, "a") as handle:
        handle.write(name + "\n")

Breakdown: Appending to the log after the algorithm returns is the whole design — a tile is only marked done once its output exists, so an interrupted run leaves the in-flight tile unmarked and the next run redoes it. Opening the file per line rather than holding a handle is deliberate: an append-and-close is durable across a kill, while a buffered handle loses whatever had not been flushed. A tile that produces no output because it falls entirely outside the catchment is a legitimate result rather than a failure, so check the point count rather than the file's existence if you care about the distinction. The same restart pattern appears in running an algorithm over a folder of files.

QGIS version compatibility

pdal:clip, pdal:tile and pdal:merge arrived in 3.32. Virtual point cloud support — both the pdal:virtualpointcloud builder and reading a .vpc — needs 3.34, and some 3.32 builds expose the builder under a different identifier, so enumerate the registry rather than assuming. Behaviour is unchanged through 3.44.

Troubleshooting

  • The clip output is empty. CRS mismatch between the overlay and the cloud, or the boundary genuinely does not overlap.
  • pdal:tile errors on the output. It wants a directory, not a filename.
  • The manifest will not open. Opened with the copc provider instead of pdal, or a member path moved after the manifest was built — paths inside a .vpc are resolved as written.
  • Seams in a gridded mosaic. Tiles were processed without a buffer for a neighbourhood operation.
  • Disk fills during retiling. Retiling writes a full copy of everything it keeps. Clip first so there is less of it.
  • Two tiles both contain the same points. pdal:tile assigns each point to exactly one tile, but a clip run per delivery tile can duplicate points in the overlap zone if the delivery tiles themselves overlap. Merge and re-tile rather than concatenating clips.

Conclusion

Clip before you retile so the expensive step moves less data, keep tiles in the tens of millions of points, buffer any neighbourhood operation and trim afterwards, build a manifest so downstream code sees one layer, and log completions after the fact so the job can be resumed. That is the whole of multi-tile LiDAR handling, and it is the same shape as every other batch problem in QGIS.

Frequently Asked Questions

Can I clip to several polygons and get one output each? Not from pdal:clip, which unions the overlay. Loop over the features, writing each to its own extent-filtered output, or use the algorithm's iteration support in the batch interface.

Does merging tiles back together lose anything?pdal:merge concatenates points and takes the union of the extents. It does not deduplicate, so merging overlapping inputs gives you duplicate points, which quietly doubles the weight of the overlap zone in any statistic.

How large should a tile be? Large enough that per-tile overheads are negligible and small enough to fit comfortably in memory — for most work, a few tens of millions of points, which is 500 m to 1 km on a typical survey.

Can a virtual point cloud reference files on a remote server? Yes, where the members are HTTP-accessible COPC files, since the manifest just holds URIs. Expect throughput to be dominated by how well the server honours range requests, as with reading a cloud optimized GeoTIFF.