⚙️ How-To LiDAR & point clouds Batch processing LiDAR & Point Clouds · step 8 of 14 Updated

How to Clip and Tile a Point Cloud in Python

Problem statement

Point clouds are too large to process whole and have no spatial index, so "give me the points in this box" is the expensive operation the format was not designed for.

Two consequences shape every workflow. Clipping a plain LAZ file means decompressing all of it. And processing a large area means tiling it, which introduces edge effects at every tile boundary.

A real 3DEP delivery arrives already tiled: 287 LAZ files, 82.8 MB, 12,968,770 points over 877 × 876 m. That tiling is what makes a bounding-box read cheap — you can reject most tiles from their headers alone.

Quick answer

Filter on tile headers before decompressing anything:

import glob
import laspy
import numpy as np


def read_bbox(pattern, bounds):
    """Read only the tiles whose header bounds intersect the box."""
    left, bottom, right, top = bounds
    parts, read, skipped = [], 0, 0

    for path in sorted(glob.glob(pattern)):
        with laspy.open(path) as reader:          # header only
            h = reader.header
            if (h.maxs[0] < left or h.mins[0] > right or
                    h.maxs[1] < bottom or h.mins[1] > top):
                skipped += 1
                continue
        read += 1
        las = laspy.read(path)
        x, y = np.asarray(las.x), np.asarray(las.y)
        inside = (x >= left) & (x <= right) & (y >= bottom) & (y <= top)
        if inside.any():
            parts.append((x[inside], y[inside], np.asarray(las.z)[inside],
                          np.asarray(las.classification)[inside]))

    print(f"  {read} tiles read, {skipped} skipped from headers")
    return [np.concatenate(a) for a in zip(*parts)]

Opening a header is microseconds; decompressing a tile is a second or more. For a 287-tile delivery that ratio is the whole optimisation.

Tile headers tested against a query box so most tiles are rejected without decompression, with only the intersecting few read.
The header is the poor man's spatial index, and for a tiled delivery it is usually enough.

Step-by-step solution

1. Clip by geometry, not just by box 

from shapely.prepared import prep
from shapely.geometry import Point

prepared = prep(polygon)
inside = np.fromiter(
    (prepared.contains(Point(px, py)) for px, py in zip(x, y)),
    dtype=bool, count=len(x))

Point-in-polygon over millions of points is slow even prepared. Always filter by the polygon's bounding box first — that is a vectorised comparison — and run the exact test only on the survivors.

2. Tile with a buffer, then trim

Every neighbourhood operation reads beyond the cell it writes. A DTM interpolation, a ground filter, a smoothing pass — each needs points outside the tile to produce correct values at the edge.

BUFFER = 20.0     # metres, larger than any neighbourhood you will use

read_bounds = (left - BUFFER, bottom - BUFFER, right + BUFFER, top + BUFFER)
write_bounds = (left, bottom, right, top)

Process the buffered points, write only the unbuffered extent. Without this, every tile boundary shows as a seam.

3. Size the tiles by memory, not by aesthetics

A tile of n points costs roughly 24 bytes for coordinates as float64 plus a few more for attributes — call it 40 bytes per point. At 16.9 points per square metre, a 500 m tile is 4.2 million points, about 170 MB.

Pick the tile size so the buffered tile fits comfortably, with room for the intermediate arrays your processing creates.

4. Write tiles that align to a shared grid 

left = np.floor(bounds[0] / tile_size) * tile_size

Snapping tile origins to a multiple of the tile size makes the outputs mosaic without resampling, and makes it trivial to work out which tile contains a coordinate.

5. Use COPC or EPT when you control the format

If you are producing the data, or converting it once, an indexed format turns every future bounding-box query into a range request rather than a full read. That is the difference between an interactive tool and a batch job.

A tile read with a buffer so neighbourhood operations are correct at the edge, then written without the buffer.
Read buffered, write unbuffered. Skipping this produces a seam at every boundary.

Code examples

Example 1 — clipping to a polygon, efficiently 

import numpy as np
from shapely.geometry import Point
from shapely.prepared import prep


def clip_to_polygon(x, y, z, classification, polygon):
    """Bounding box first, exact containment only on the survivors."""
    left, bottom, right, top = polygon.bounds
    in_box = (x >= left) & (x <= right) & (y >= bottom) & (y <= top)
    print(f"  bounding box keeps {in_box.sum():,} of {len(x):,} "
          f"({in_box.mean():.1%})")

    if not in_box.any():
        return tuple(np.array([]) for _ in range(4))

    prepared = prep(polygon)
    bx, by = x[in_box], y[in_box]
    exact = np.fromiter((prepared.contains(Point(a, b)) for a, b in zip(bx, by)),
                        dtype=bool, count=len(bx))

    keep = np.zeros(len(x), bool)
    keep[np.nonzero(in_box)[0][exact]] = True
    print(f"  polygon keeps {keep.sum():,} ({keep.mean():.1%})")
    return x[keep], y[keep], z[keep], classification[keep]

The two-stage filter is not an optimisation detail. The bounding-box test is a vectorised comparison over millions of points in milliseconds; the exact test is a Python-level call per point. Reducing the number of exact tests by an order of magnitude is the difference between seconds and minutes.

Example 2 — a tiling scheme with buffers 

import numpy as np


def tile_scheme(bounds, tile_size, buffer=20.0, snap=True):
    """Tile definitions with read and write extents, on a snapped grid."""
    left, bottom, right, top = bounds
    if snap:
        left = np.floor(left / tile_size) * tile_size
        bottom = np.floor(bottom / tile_size) * tile_size

    tiles = []
    y = bottom
    while y < top:
        x = left
        while x < right:
            tiles.append({
                "id": f"{int(x)}_{int(y)}",
                "write": (x, y, x + tile_size, y + tile_size),
                "read": (x - buffer, y - buffer,
                         x + tile_size + buffer, y + tile_size + buffer),
            })
            x += tile_size
        y += tile_size

    print(f"  {len(tiles)} tiles of {tile_size} m with a {buffer} m buffer")
    print(f"  buffered area is {((tile_size + 2 * buffer) / tile_size) ** 2:.2f}x "
          "the written area")
    return tiles
  4 tiles of 500 m with a 20 m buffer
  buffered area is 1.17x the written area

The overhead ratio is worth printing. A 20 m buffer on 500 m tiles costs 17% extra reading; the same buffer on 100 m tiles costs 96%, which is usually the signal to use larger tiles.

Example 3 — processing tiles in parallel 

import glob
import os
from concurrent.futures import ProcessPoolExecutor

import numpy as np


def process_tile(spec, pattern, cell=1.0, out_dir="dtm"):
    """One tile: read buffered, process, write unbuffered."""
    x, y, z, classification = read_bbox(pattern, spec["read"])
    if len(x) == 0:
        return spec["id"], 0

    dtm, counts, transform = dtm_from_points(x, y, z, classification, cell)

    wl, wb, wr, wt = spec["write"]
    left, top = spec["read"][0], spec["read"][3]
    c0 = int((wl - left) / cell)
    c1 = int((wr - left) / cell)
    r0 = int((top - wt) / cell)
    r1 = int((top - wb) / cell)
    trimmed = dtm[r0:r1, c0:c1]

    os.makedirs(out_dir, exist_ok=True)
    write_raster(f"{out_dir}/{spec['id']}.tif", trimmed, transform, cell,
                 offset=(c0, r0))
    return spec["id"], int(np.isfinite(trimmed).sum())


def run_tiles(tiles, pattern, workers=4):
    with ProcessPoolExecutor(workers) as pool:
        futures = [pool.submit(process_tile, t, pattern) for t in tiles]
        for future in futures:
            tile_id, cells = future.result()
            print(f"  {tile_id}: {cells:,} cells written")

Tiles are independent once buffered, so process-level parallelism is straightforward and scales with cores. The buffer is what makes that true: without it, tiles would need to exchange edge data and the parallelism would be gone.

Explanation

Why LAZ makes clipping expensive

LAZ compresses in independent chunks, but the points inside them are in acquisition order — roughly, the order the aircraft flew. The points in any bounding box are therefore scattered across every chunk.

So reading a small area from a large file means decompressing all of it. The chunking does allow a reader to skip a chunk whose declared bounds miss the query, but writers do not always record those bounds, and acquisition order means most chunks overlap most queries anyway.

That is precisely the gap COPC fills: it reorders the points into an octree so that spatial locality exists in the file layout, and records the byte range of each node.

Why the header filter works so well on a tiled delivery

A tiled delivery has already done the spatial sort, at the granularity of a tile. Each header records the tile's bounds, and reading a header is a few hundred bytes.

For the 287-tile survey used here, a query covering a tenth of the area rejects roughly nine tenths of the tiles for the cost of 287 header reads — microseconds each. The remaining tiles are decompressed in full and filtered in memory.

This is a spatial index with one level and no code. It is why tiled deliveries remain the norm despite indexed formats existing.

Why buffers are not optional

Any operation with a neighbourhood — interpolation, ground filtering, smoothing, slope — needs data outside the cell it writes.

At a tile edge without a buffer, that neighbourhood is truncated. A ground filter sees half its usual context and makes different decisions; an interpolation extrapolates instead of interpolating. The result is a visible discontinuity along every boundary, and it does not average out when the tiles are mosaicked.

The buffer must exceed the largest neighbourhood in the pipeline. Twenty metres covers most raster operations at metre resolution; a ground filter with a 50 m window needs 50.

Why tile size is a memory decision

Small tiles mean more boundaries, more buffer overhead and more per-tile fixed cost. Large tiles mean more memory per worker and less parallelism.

The buffer overhead makes this concrete. A 20 m buffer costs 17% extra area on 500 m tiles and 96% on 100 m tiles — so the smallest sensible tile is set by the buffer, not by preference.

Work backwards from memory: choose the largest tile whose buffered point count fits comfortably in a worker's share of RAM, and divide the area accordingly.

Buffer overhead rising from 1.08 times at 1000 metre tiles to 1.96 times at 100 metre tiles.
A 20 m buffer nearly doubles the reading on 100 m tiles. That is what sets the minimum tile size.

Edge cases or notes

  • Filter on tile headers first. Microseconds against seconds per tile.
  • Bounding box before exact polygon containment, always.
  • Buffer must exceed the largest neighbourhood in the pipeline.
  • Write the unbuffered extent or tiles overlap in the mosaic.
  • Snap tile origins to a multiple of the tile size.
  • Buffer overhead grows as tiles shrink — 17% at 500 m, 96% at 100 m for a 20 m buffer.
  • Points exactly on a boundary belong to one tile only; use a half-open interval.
  • COPC or EPT if you control the format — a bounding-box query becomes a range request.

FAQ

How do I clip a point cloud to a polygon in Python?

Filter by the polygon's bounding box first with vectorised comparisons, then run an exact containment test only on the survivors using a prepared geometry.

Why is reading a small area from a LAZ file so slow?

LAZ has no spatial index and stores points in acquisition order, so the points in any box are scattered through the whole file. Use COPC or EPT if you need cheap spatial queries.

How do I avoid seams between tiles?

Read each tile with a buffer larger than any neighbourhood your processing uses, and write only the unbuffered extent.

How big should the buffer be?

Larger than the biggest window in the pipeline. Twenty metres covers most metre-resolution raster work; a ground filter with a 50 m window needs 50.

What tile size should I use?

The largest whose buffered point count fits in a worker's memory. Small tiles pay a large buffer overhead — 96% extra area for a 20 m buffer on 100 m tiles.

Can I process tiles in parallel?

Yes, once they are buffered. Buffering is what makes tiles independent, so a process pool scales with cores.

Should I convert my tiles to COPC?

If you will query them repeatedly by area, yes. A one-off conversion turns every later bounding-box read into a range request instead of a full decompression.