How to Clip and Tile a Point Cloud in Python
On this page
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. The nine used here — a corner of the USGS survey of New York City after Hurricane Sandy — are 1,500 m squares, 393 MB and 105,039,668 points between them. 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)]
import time
start = time.perf_counter()
x, y, z, classification = read_bbox("tiles/*.laz", reservoir.total_bounds)
print(f" {len(x):,} points in the reservoir's box, {time.perf_counter() - start:.1f} s")
start = time.perf_counter()
for path in sorted(glob.glob("tiles/*.laz")):
with laspy.open(path) as reader:
reader.header.maxs
print(f" all nine headers: {(time.perf_counter() - start) * 1000:.1f} ms") 1 tiles read, 8 skipped from headers
1,631,125 points in the reservoir's box, 0.4 s
all nine headers: 1.6 ms
Opening a header is a fraction of a millisecond; decompressing a tile takes seconds. For a delivery of hundreds of tiles that ratio is the whole optimisation.
Clip by geometry, tile with a buffer, size by memory
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 in a Python loop even prepared — 7.2 s for 1.6 million points here. Filter by the polygon's bounding box first, a vectorised comparison, and run the exact test on the survivors with shapely.contains_xy, which is vectorised too (Example 1).
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. Measure the density from the headers rather than assuming it:
counts, area = 0, 0.0
for path in sorted(glob.glob("tiles/*.laz")):
with laspy.open(path) as reader:
h = reader.header
counts += h.point_count
area += (h.maxs[0] - h.mins[0]) * (h.maxs[1] - h.mins[1])
density = counts / area
for size in (250, 500, 1000):
n = density * size ** 2
print(f"{size:5d} m tile: {n / 1e6:5.2f} million points, about {n * 40 / 1e6:4.0f} MB")
print(f"{counts:,} points over {area / 1e6:.2f} km², {density:.2f} points per m²") 250 m tile: 0.32 million points, about 13 MB
500 m tile: 1.30 million points, about 52 MB
1000 m tile: 5.19 million points, about 207 MB
105,039,668 points over 20.25 km², 5.19 points per m²
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_sizeSnapping 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.
Clipping efficiently and a tiling scheme
The data used on this page
Nine tiles of the USGS 3DEP delivery NY_New_York_CMGP_SANDY_LiDAR_15 — the post-Hurricane-Sandy survey of New York City, flown in 2013–14 — covering 4.5 × 4.5 km of Manhattan around Central Park: USNG tiles 18TWL 850–880 E, 135–165 N, each 1,500 m square, LAS 1.2 compressed as LAZ, in NAD83(2011) / UTM zone 18N with NAVD88 heights. They come from the USGS rockyweb server, which is slow — expect well over an hour for the 393 MB. The clipping polygon is the Jacqueline Kennedy Onassis Reservoir from OpenStreetMap (way 6678417).
import os
import subprocess
import geopandas as gpd
BASE = ("https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/Projects/"
"NY_New_York_CMGP_SANDY_LiDAR_15/NY_CMPG_2013/LAZ/USGS_LPC_NY_New_York_CMGP_SANDY_LiDAR_15_18TWL")
TILES = [f"{e}{n}" for n in (135, 150, 165) for e in (850, 865, 880)]
os.makedirs("tiles", exist_ok=True)
for tile in TILES:
path = f"tiles/{tile}.laz"
for _ in range(20): # rockyweb drops long transfers: resume
if os.path.exists(path):
break
if subprocess.run(["curl", "-sfL", "-C", "-", "--retry", "3", "-o", path + ".part",
f"{BASE}{tile}.laz"]).returncode == 0:
os.replace(path + ".part", path)
if not os.path.exists("reservoir.geojson"):
subprocess.run(["curl", "-sf", "-A", "spatialworkflow.io guide", "-o", "reservoir.geojson",
"https://nominatim.openstreetmap.org/lookup?osm_ids=W6678417&format=geojson&polygon_geojson=1"],
check=True)
reservoir = gpd.read_file("reservoir.geojson").to_crs("EPSG:6347")
print(len(os.listdir("tiles")), "tiles,", f"{sum(os.path.getsize(f'tiles/{t}.laz') for t in TILES) / 1e6:.0f} MB;",
f"reservoir {reservoir.area.sum() / 1e6:.3f} km²")9 tiles, 393 MB; reservoir 0.387 km²
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]
import shapely
polygon = reservoir.geometry.iloc[0]
start = time.perf_counter()
clipped = clip_to_polygon(x, y, z, classification, polygon)
print(f" prepared contains, point by point: {time.perf_counter() - start:.1f} s")
start = time.perf_counter()
inside = shapely.contains_xy(polygon, x, y)
print(f" shapely.contains_xy on all {len(x):,}: {time.perf_counter() - start:.2f} s, "
f"same result: {inside.sum() == len(clipped[0])}")
print(f" water (class 9) inside: {(clipped[3] == 9).mean():.1%}") bounding box keeps 1,631,125 of 1,631,125 (100.0%)
polygon keeps 519,699 (31.9%)
prepared contains, point by point: 7.3 s
shapely.contains_xy on all 1,631,125: 0.08 s, same result: True
water (class 9) inside: 46.1%
The two-stage filter is not an optimisation detail when the exact test is a Python-level call per point: here the prepared-geometry loop took 7.2 s for 1.6 million points. Shapely 2 does the same test vectorised — shapely.contains_xy gave the identical answer in 0.07 s, a hundred times faster — so with a current Shapely the box filter matters less than never looping in Python.
The reservoir also shows what the clip keeps. Only 31.9% of the points in its bounding box were inside it, and 46.1% of those were classed as water: a lidar pulse that reaches open water mostly does not come back, so the reservoir is sparse compared with the park around it.
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
delivery = (585_000, 4_513_500, 589_500, 4_518_000)
for size in (1500, 500, 100):
scheme = tile_scheme(delivery, size) 9 tiles of 1500 m with a 20.0 m buffer
buffered area is 1.05x the written area
81 tiles of 500 m with a 20.0 m buffer
buffered area is 1.17x the written area
2025 tiles of 100 m with a 20.0 m buffer
buffered area is 1.96x the written area
The overhead ratio is worth printing. A 20 m buffer on the delivery's own 1,500 m tiles costs 5% extra reading, on 500 m tiles 17%, and on 100 m tiles 96% — the signal to use larger tiles.
Example 3 — processing tiles in parallel
A process pool has to send each worker the function to run, and since Python 3.14 Linux starts workers with forkserver, which can only send a function its processes can import. A worker defined in a notebook cell cannot be imported, and the fork context is no escape here: laspy's LAZ decoder runs its own thread pool, and a process forked from one with live threads hung with its workers at 0% CPU. So the worker goes in a module:
%%writefile tile_worker.py
import glob
import os
import laspy
import numpy as np
import rasterio
from rasterio.fill import fillnodata
from rasterio.transform import from_origin
def read_bbox(pattern, bounds):
"""Read only the tiles whose header bounds intersect the box."""
left, bottom, right, top = bounds
parts = []
for path in sorted(glob.glob(pattern)):
with laspy.open(path) as reader:
h = reader.header
if h.maxs[0] < left or h.mins[0] > right or h.maxs[1] < bottom or h.mins[1] > top:
continue
las = laspy.read(path)
x, y = np.asarray(las.x), np.asarray(las.y)
inside = (x >= left) & (x <= right) & (y >= bottom) & (y <= top)
parts.append((x[inside], y[inside], np.asarray(las.z)[inside],
np.asarray(las.classification)[inside]))
return [np.concatenate(a) for a in zip(*parts)]
def dtm_from_points(x, y, z, classification, cell, bounds):
"""Mean height of ground (class 2) points per cell, gaps filled from neighbours."""
left, bottom, right, top = bounds
width, height = int(round((right - left) / cell)), int(round((top - bottom) / cell))
ground = classification == 2
col = ((x[ground] - left) / cell).astype(int)
row = ((top - y[ground]) / cell).astype(int)
ok = (col >= 0) & (col < width) & (row >= 0) & (row < height)
index = row[ok] * width + col[ok]
counts = np.bincount(index, minlength=width * height).reshape(height, width)
sums = np.bincount(index, weights=z[ground][ok], minlength=width * height).reshape(height, width)
dtm = np.where(counts > 0, sums / np.maximum(counts, 1), 0).astype("float32")
dtm = fillnodata(dtm, mask=counts > 0, max_search_distance=15) # the neighbourhood operation
dtm[dtm == 0] = np.nan
return dtm, counts, from_origin(left, top, cell, cell)
def write_raster(path, array, transform, offset=(0, 0)):
profile = {"driver": "GTiff", "width": array.shape[1], "height": array.shape[0], "count": 1,
"dtype": "float32", "crs": "EPSG:6347", "nodata": np.nan,
"transform": transform * transform.translation(*offset)}
with rasterio.open(path, "w", **profile) as dst:
dst.write(array, 1)
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, spec["read"])
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, offset=(c0, r0))
return spec["id"], int(np.isfinite(trimmed).sum())from concurrent.futures import ProcessPoolExecutor
import rasterio
from tile_worker import dtm_from_points, process_tile
def run_tiles(tiles, pattern, workers=4, **kwargs):
with ProcessPoolExecutor(workers) as pool:
futures = [pool.submit(process_tile, t, pattern, **kwargs) for t in tiles]
for future in futures:
tile_id, cells = future.result()
print(f" {tile_id}: {cells:,} cells written")
around_reservoir = tile_scheme((587_000, 4_515_000, 588_000, 4_516_000), 500)
run_tiles(around_reservoir, "tiles/*.laz", workers=2) 4 tiles of 500 m with a 20.0 m buffer
buffered area is 1.17x the written area
587000_4515000: 195,665 cells written
587500_4515000: 143,772 cells written
587000_4515500: 182,816 cells written
587500_4515500: 122,014 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.
Does the buffer matter? Run the same four tiles with no buffer, and compare both mosaics with one DTM made in a single pass over the whole square:
def mosaic(out_dir, block=(587_000, 4_515_000), size=500, cell=1.0):
n = int(size / cell)
out = np.full((2 * n, 2 * n), np.nan, dtype="float32")
for i, y0 in enumerate((block[1] + size, block[1])): # top row first
for j, x0 in enumerate((block[0], block[0] + size)):
with rasterio.open(f"{out_dir}/{int(x0)}_{int(y0)}.tif") as src:
out[i * n:(i + 1) * n, j * n:(j + 1) * n] = src.read(1)
return out
unbuffered = tile_scheme((587_000, 4_515_000, 588_000, 4_516_000), 500, buffer=0.0)
run_tiles(unbuffered, "tiles/*.laz", workers=2, out_dir="dtm_nobuffer")
whole = (586_980, 4_514_980, 588_020, 4_516_020)
reference = dtm_from_points(*read_bbox("tiles/*.laz", whole), 1.0, whole)[0][20:-20, 20:-20]
near_seam = np.zeros(reference.shape, bool)
near_seam[485:515, :] = near_seam[:, 485:515] = True
for label, out_dir in [("20 m buffer", "dtm"), ("no buffer", "dtm_nobuffer")]:
diff = np.abs(mosaic(out_dir) - reference)
for where, sel in [("within 15 m of a seam", near_seam), ("elsewhere", ~near_seam)]:
d = diff[sel]
print(f" {label:12} {where:22} max {np.nanmax(d):5.2f} m "
f"cells differing by > 1 cm: {(d > 0.01).sum():6,} no value: {np.isnan(d).sum():6,}") 4 tiles of 500 m with a 0.0 m buffer
buffered area is 1.00x the written area
587000_4515000: 195,665 cells written
587500_4515000: 143,526 cells written
587000_4515500: 182,780 cells written
587500_4515500: 122,009 cells written
4 tiles read, 5 skipped from headers
20 m buffer within 15 m of a seam max 0.00 m cells differing by > 1 cm: 0 no value: 36,270
20 m buffer elsewhere max 0.00 m cells differing by > 1 cm: 0 no value: 319,463
no buffer within 15 m of a seam max 1.64 m cells differing by > 1 cm: 765 no value: 36,486
no buffer elsewhere max 2.42 m cells differing by > 1 cm: 4,026 no value: 319,534
With a 20 m buffer the four tiles reproduced the single-pass DTM exactly — not one cell off by more than a centimetre. Without it, 765 cells within 15 m of the internal seams differed, by up to 1.64 m, where hole-filling near a tile edge could only look one way. The 4,026 cells "elsewhere" are mostly along the square's outer edge, which is a seam too. About a third of the square has no value in either version: the reservoir and the buildings have no ground returns within the 15-cell fill distance.
Why LAZ makes clipping expensive
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 nine tiles used here, the reservoir query rejected eight of them from their headers — all nine headers took 1.6 ms to read — and decompressed one, in 0.4 s. A delivery of hundreds of tiles scales the same way: a small query costs its header reads plus the few tiles it touches.
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.
Header filters, bounding boxes and process pools
- Filter on tile headers first. 1.6 ms for nine headers against 0.4 s to decompress one tile here.
- Bounding box before exact polygon containment, and the exact test vectorised with
shapely.contains_xy. - A process pool needs an importable worker on Python 3.14's
forkserverdefault;forkdeadlocked here against laspy's decoder threads. - 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.
Internal links
- LAS, LAZ and COPC explained — why clipping is expensive
- How to read a LAS or LAZ file in Python with laspy — the header-only read
- How to rasterise a point cloud to a grid in Python — what each tile produces
- How to create a DTM from LiDAR ground points in Python — the operation needing a buffer
- How to split a large layer into tiles for batch processing — the vector equivalent
- How to speed up batch GIS jobs with parallel processing in Python — running the tiles
- How to process GIS files in nested subfolders and mirror the output tree — organising the output
- Point density explained — sizing tiles from density
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 the exact test on the survivors with shapely.contains_xy — a hundred times faster than a prepared geometry in a Python loop.
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.