How to Rasterise a Point Cloud to a Grid in Python
On this page
Problem statement
Turning millions of scattered points into a regular grid is the operation underneath every lidar product — DTM, DSM, CHM, density, intensity. It is also where naive code gets slow, and where the fix has a slow path of its own.
A Python loop over a lidar survey's points is slow enough to matter, and the obvious vectorised replacement hides a slow path of its own (step 4 times both on the 12,968,770 points used here). The difference is one NumPy idiom, used with the right argument.
The other half of the problem is choosing a reduction and a cell size that the data supports. Measured on that survey, counting only the cells the survey covered, 33.1% of 0.25 m cells contained no point at all, 1.5% of 0.5 m cells, and none at 1 m (Example 2).
Quick answer
import numpy as np
def rasterise(x, y, values, cell, bounds, how="mean"):
"""Reduce scattered points onto a regular grid, vectorised."""
left, bottom, right, top = bounds
width = int(np.ceil((right - left) / cell))
height = int(np.ceil((top - bottom) / cell))
col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
flat = row * width + col
if how == "count":
out = np.zeros(width * height)
np.add.at(out, flat, 1.0)
elif how == "min":
out = np.full(width * height, np.inf)
np.minimum.at(out, flat, values)
out[np.isinf(out)] = np.nan
elif how == "max":
out = np.full(width * height, -np.inf)
np.maximum.at(out, flat, values)
out[np.isinf(out)] = np.nan
elif how == "mean":
total = np.zeros(width * height)
count = np.zeros(width * height)
np.add.at(total, flat, values)
np.add.at(count, flat, 1.0)
out = np.where(count > 0, total / np.maximum(count, 1), np.nan)
else:
raise ValueError(f"unknown reduction {how!r}")
return out.reshape(height, width)
bounds = (x.min(), y.min(), x.max(), y.max())
counts = rasterise(x, y, None, 1.0, bounds, how="count")
dsm = rasterise(x, y, z, 1.0, bounds, how="max")
print(f"{counts.shape[1]} x {counts.shape[0]} cells, {int(counts.sum()):,} points, "
f"highest surface {np.nanmax(dsm):.1f} m")877 x 876 cells, 12,968,770 points, highest surface 203.1 m
np.add.at, np.minimum.at and np.maximum.at are the unbuffered in-place operations. They handle repeated indices correctly, which plain fancy indexing does not.
An explicit grid, integer cell indices, flat accumulation
1. Define the grid explicitly
left, bottom, right, top = bounds
width = int(np.ceil((right - left) / cell))
height = int(np.ceil((top - bottom) / cell))
transform = from_origin(left, top, cell, cell)Deriving the grid from the data's own extent is convenient and makes two rasters from two tiles non-aligned. For anything that will be mosaicked or stacked, snap the origin to a multiple of the cell size.
2. Compute cell indices by integer division
col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
row = np.clip(((top - y) / cell).astype(int), 0, height - 1)The top - y for rows is the usual source of vertically flipped rasters: raster row 0 is the top, while y increases upward.
np.clip handles the point exactly on the right or bottom edge, which would otherwise index one past the end.
3. Flatten to one index and reduce
flat = row * width + colA single flat index turns a 2D scatter into a 1D reduction, which is what the .at operations need. Reshape at the end.
4. Do not use a Python loop, and do not use plain indexing
import time
left, top = x.min(), y.max()
width, height = int(np.ceil(x.max() - left)), int(np.ceil(top - y.min())) # 1 m cells
flat = (np.clip((top - y).astype(int), 0, height - 1) * width
+ np.clip((x - left).astype(int), 0, width - 1))
wrong = np.zeros(width * height)
wrong[flat] += 1 # WRONG: repeated indices are applied once
right = np.zeros(width * height)
np.add.at(right, flat, 1.0) # correct
print(f"plain indexing counted {int(wrong.sum()):,} points; np.add.at {int(right.sum()):,}")
def timed(fn):
start = time.perf_counter()
fn()
return time.perf_counter() - start
def python_loop():
out = np.zeros(width * height)
for i in flat:
out[i] += 1
for label, fn in [("Python loop", python_loop),
("np.add.at(out, flat, 1)", lambda: np.add.at(np.zeros(width * height), flat, 1)),
("np.add.at(out, flat, 1.0)", lambda: np.add.at(np.zeros(width * height), flat, 1.0)),
("np.bincount(flat)", lambda: np.bincount(flat, minlength=width * height))]:
print(f"{label:28} {timed(fn):6.3f} s for {len(flat):,} points")plain indexing counted 710,856 points; np.add.at 12,968,770
Python loop 1.574 s for 12,968,770 points
np.add.at(out, flat, 1) 0.424 s for 12,968,770 points
np.add.at(out, flat, 1.0) 0.015 s for 12,968,770 points
np.bincount(flat) 0.010 s for 12,968,770 points
This is the trap. out[flat] += 1 reads, adds and writes as three separate vectorised steps, so where an index appears twenty times it is incremented once. The result looks like a valid count raster and undercounts every populated cell — here it counted one point per occupied cell.
The timings hold a second trap. np.add.at(out, flat, 1) with an integer 1 on a float array took a slow casting path and was only four times faster than a Python loop; the same call with 1.0 was about thirty times faster, and np.bincount — which only counts or sums, but does both — was faster still. Match the value's type to the array's, and reach for np.bincount(flat, weights=values) for sums and counts.
5. Return the counts alongside the values
A cell built from one point and a cell built from a hundred are indistinguishable in the output. The count raster costs one extra reduction and makes the product interpretable.
A rasteriser with percentiles and counts
The data used on this page
The USGS 3DEP lidar survey of Brooks Camp, Katmai National Park, Alaska (2012), from the USGS Entwine Point Tile archive on AWS, AK_BrooksCamp_2012: every node under octree node 4-3-13-7 — 287 LAZ files in nodes/, 12,968,770 points — joined into tile.laz, reprojected from Web Mercator to UTM zone 5N and loaded as arrays. footprint is the convex hull of the points, and surveyed() marks the grid cells inside it: the patch is turned 2.4° against the UTM grid, so its bounding box has corners that were never surveyed.
import json
import os
import urllib.request
import laspy
EPT = "https://s3-us-west-2.amazonaws.com/usgs-lidar-public/AK_BrooksCamp_2012"
ROOT = (4, 3, 13, 7) # depth, x, y, z of the octree node
def under(key):
d, x, y, z = map(int, key.split("-"))
step = 2 ** (d - ROOT[0])
return d >= ROOT[0] and (x // step, y // step, z // step) == ROOT[1:]
hierarchy = json.load(urllib.request.urlopen(f"{EPT}/ept-hierarchy/0-0-0-0.json"))
nodes = sorted(k for k in hierarchy if under(k))
os.makedirs("nodes", exist_ok=True)
for key in nodes:
path = f"nodes/{key}.laz"
if not os.path.exists(path):
urllib.request.urlretrieve(f"{EPT}/ept-data/{key}.laz", path + ".part")
os.replace(path + ".part", path)
if not os.path.exists("tile.laz"): # the whole subtree as one file
with laspy.open(f"nodes/{nodes[0]}.laz") as first:
header = first.header
with laspy.open("tile.laz", mode="w", header=header) as writer:
for key in nodes:
writer.write_points(laspy.read(f"nodes/{key}.laz").points)
with laspy.open("tile.laz") as reader:
print(len(nodes), "nodes,", f"{reader.header.point_count:,} points,",
f"{sum(os.path.getsize(f'nodes/{k}.laz') for k in nodes) / 1e6:.1f} MB of nodes")
import numpy as np
from pyproj import Transformer
las = laspy.read("tile.laz")
x, y = Transformer.from_crs(3857, 32605, always_xy=True).transform(las.x, las.y) # UTM 5N metres
z = np.asarray(las.z)
classification = np.asarray(las.classification)
return_number = np.asarray(las.return_number)
number_of_returns = np.asarray(las.number_of_returns)
import shapely
from scipy.spatial import ConvexHull
hull = ConvexHull(np.c_[x, y])
footprint = shapely.Polygon(np.c_[x, y][hull.vertices]) # the ground the survey covers
def surveyed(transform, shape):
"""True for grid cells whose centre lies inside the survey footprint."""
rows, cols = np.indices(shape)
cx, cy = transform * (cols + 0.5, rows + 0.5)
return shapely.contains_xy(footprint, cx, cy)287 nodes, 12,968,770 points, 82.8 MB of nodes
Example 1 — a rasteriser with percentiles and counts
import numpy as np
import rasterio
from rasterio.transform import from_origin
def rasterise_cloud(x, y, values, cell=1.0, bounds=None, crs=None,
reductions=("min", "max", "mean", "count"),
percentiles=(), out_path=None):
"""Several reductions in one pass over the cell index."""
if bounds is None:
bounds = (x.min(), y.min(), x.max(), y.max())
left, bottom, right, top = bounds
width = int(np.ceil((right - left) / cell))
height = int(np.ceil((top - bottom) / cell))
transform = from_origin(left, top, cell, cell)
col = np.clip(((x - left) / cell).astype(int), 0, width - 1)
row = np.clip(((top - y) / cell).astype(int), 0, height - 1)
flat = row * width + col
n = width * height
bands, names = [], []
count = np.zeros(n)
np.add.at(count, flat, 1.0)
for how in reductions:
if how == "count":
bands.append(count.copy())
elif how == "min":
g = np.full(n, np.inf); np.minimum.at(g, flat, values)
g[np.isinf(g)] = np.nan; bands.append(g)
elif how == "max":
g = np.full(n, -np.inf); np.maximum.at(g, flat, values)
g[np.isinf(g)] = np.nan; bands.append(g)
elif how == "mean":
total = np.zeros(n); np.add.at(total, flat, values)
bands.append(np.where(count > 0, total / np.maximum(count, 1), np.nan))
names.append(how)
if percentiles:
order = np.argsort(flat, kind="stable")
fs, vs = flat[order], values[order]
edges = np.searchsorted(fs, np.arange(n + 1))
for q in percentiles:
g = np.full(n, np.nan)
for i in np.nonzero(count > 0)[0]:
g[i] = np.percentile(vs[edges[i]:edges[i + 1]], q)
bands.append(g); names.append(f"p{q}")
stack = np.stack([b.reshape(height, width) for b in bands]).astype("float32")
print(f" {width} x {height} at {cell}, "
f"{(count == 0).mean():.2%} empty, bands {names}")
if out_path:
profile = dict(driver="GTiff", height=height, width=width,
count=len(bands), dtype="float32", crs=crs,
transform=transform, nodata=np.nan,
compress="deflate", tiled=True)
with rasterio.open(out_path, "w", **profile) as dst:
dst.write(stack)
for i, name in enumerate(names, start=1):
dst.set_band_description(i, name)
return stack, names, transform
stack, names, transform = rasterise_cloud(x, y, z, cell=2.0, percentiles=(95,),
crs="EPSG:32605", out_path="surfaces.tif") 439 x 438 at 2.0, 7.34% empty, bands ['min', 'max', 'mean', 'count', 'p95']
The sort-and-slice approach for percentiles is the only part that loops, and it loops over cells rather than points — thousands rather than millions.
Example 2 — choosing the cell size from occupancy
import numpy as np
def occupancy_sweep(x, y, mask=None, footprint=None, cells=(0.25, 0.5, 1.0, 2.0, 5.0)):
"""What fraction of cells contain a point, at each resolution?"""
px, py = (x[mask], y[mask]) if mask is not None else (x, y)
left, top = x.min(), y.max()
rows = []
for cell in cells:
width = int(np.ceil((x.max() - left) / cell))
height = int(np.ceil((top - y.min()) / cell))
col = np.clip(((px - left) / cell).astype(int), 0, width - 1)
row = np.clip(((top - py) / cell).astype(int), 0, height - 1)
counts = np.zeros(width * height)
np.add.at(counts, row * width + col, 1.0)
if footprint is not None: # count only surveyed cells
r, c = np.indices((height, width))
counts = counts[shapely.contains_xy(footprint, left + (c.ravel() + 0.5) * cell,
top - (r.ravel() + 0.5) * cell)]
rows.append({"cell": cell, "cells": counts.size,
"empty": float((counts == 0).mean()),
"median": float(np.median(counts[counts > 0]))})
print(f" {cell:5.2f} m {counts.size:10,} cells "
f"{(counts == 0).mean():7.2%} empty "
f"median {np.median(counts[counts > 0]):5.0f} points/cell")
return rows
print("bounding box")
box_rows = occupancy_sweep(x, y)
print("inside the footprint")
surveyed_rows = occupancy_sweep(x, y, footprint=footprint)bounding box
0.25 m 12,288,524 cells 38.23% empty median 1 points/cell
0.50 m 3,073,008 cells 9.00% empty median 4 points/cell
1.00 m 768,252 cells 7.47% empty median 17 points/cell
2.00 m 192,282 cells 7.34% empty median 66 points/cell
5.00 m 30,976 cells 7.29% empty median 410 points/cell
inside the footprint
0.25 m 11,349,660 cells 33.14% empty median 1 points/cell
0.50 m 2,837,417 cells 1.52% empty median 4 points/cell
1.00 m 709,354 cells 0.00% empty median 17 points/cell
2.00 m 177,338 cells 0.00% empty median 66 points/cell
5.00 m 28,376 cells 0.00% empty median 413 points/cell
Over the bounding box the empty fraction stops falling at about 7.3% — a plateau that looks like a resolution limit and is the 7.7% of the box outside the survey. Counted over the surveyed cells, every cell is occupied from 1 m up; the resolution limit for this cloud sits between 0.5 m and 1 m.
Example 3 — rasterising a cloud too large for memory
import glob
import numpy as np
import laspy
def rasterise_tiles(pattern, cell, bounds, crs, ground_class=2):
"""Accumulate a DTM across many tiles without loading them together."""
left, bottom, right, top = bounds
width = int(np.ceil((right - left) / cell))
height = int(np.ceil((top - bottom) / cell))
grid = np.full(width * height, np.inf)
counts = np.zeros(width * height)
tiles = 0
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
tiles += 1
for points in laspy.open(path).chunk_iterator(1_000_000):
x = np.asarray(points.x); y = np.asarray(points.y)
z = np.asarray(points.z)
keep = (np.asarray(points.classification) == ground_class) & \
(x >= left) & (x < right) & (y >= bottom) & (y < top)
if not keep.any():
continue
col = ((x[keep] - left) / cell).astype(int)
row = ((top - y[keep]) / cell).astype(int)
flat = row * width + col
np.minimum.at(grid, flat, z[keep])
np.add.at(counts, flat, 1.0)
grid[np.isinf(grid)] = np.nan
print(f" {tiles} tiles, {int(counts.sum()):,} ground points, "
f"{np.isnan(grid).mean():.2%} empty")
return grid.reshape(height, width), counts.reshape(height, width)
# the 287 EPT node files, in their own Web Mercator coordinates
with laspy.open("tile.laz") as reader:
mins, maxs = reader.header.mins, reader.header.maxs
box = (mins[0], mins[1], maxs[0] + 1e-6, maxs[1] + 1e-6)
by_tile, tile_counts = rasterise_tiles("nodes/*.laz", 4.0, box, "EPSG:3857")
ground = classification == 2
at_once = rasterise(np.asarray(las.x)[ground], np.asarray(las.y)[ground], z[ground], 4.0, box, how="min")
print(f" same grid as a single pass: {np.array_equal(by_tile, at_once, equal_nan=True)}") 287 tiles, 2,965,430 ground points, 0.14% empty
same grid as a single pass: True
Because minimum and count are both associative, the grid can be accumulated across tiles and chunks without holding any of them. Mean needs the same treatment via separate sum and count arrays; median and percentiles cannot be accumulated this way at all.
np.add.at, reductions and single-point cells
Why np.add.at and not +=
out[flat] += 1 is three operations: gather the values at flat, add one, scatter them back. With a repeated index, all copies gather the same original value and the last scatter wins — so twenty points in a cell increment it once.
np.add.at is unbuffered: it applies the operation for each index in turn. It is slower per element than a buffered operation and correct, which is the trade you want.
The failure is silent and produces a count raster that looks fine and undercounts every populated cell. Every count-based product built on it — density, occupancy, weighting — is wrong.
Why the reduction choice is not cosmetic
Each reduction answers a different question, and the wrong one is wrong in a specific direction:
- min for ground: robust to misclassified low vegetation, vulnerable to low noise.
- max for surfaces: the top of the canopy, vulnerable to high noise.
- mean for continuous fields like intensity: smooth, and biased by any outlier.
- count for density: exact, and the only one that needs no values at all.
- percentiles for canopy structure: expensive, and the only ones that describe a distribution.
Why cells with one point are the real risk
At 1 m on this survey the median cell holds 17 points, which sounds comfortable. The distribution is what matters: cells at the edge of the flight line, under canopy, or at the tile boundary may hold one.
A one-point cell has no redundancy. Whatever reduction you use returns that point, so a single bad measurement becomes a cell value with no way to detect it from the raster alone.
Returning the count raster is what lets a downstream user mask thinly supported cells.
Why grid alignment matters more than it seems
Two tiles rasterised from their own extents produce grids offset by a fraction of a cell. Mosaicking them requires resampling, which blurs; stacking them for a difference produces a systematic shift.
Snapping the origin to a multiple of the cell size — left = floor(left / cell) * cell — makes every product in a project share pixel corners for free. It costs one line and removes a whole class of alignment bug.
Row order, clipping and grid origins
np.add.at, not+=. Plain indexing applies a repeated index once.- Rows count from the top. Use
top - y, or the raster is flipped. np.clipthe indices so a point exactly on the far edge does not overflow.- Snap the grid origin to a multiple of the cell size.
- Return the count raster alongside the values.
- min, max, sum and count accumulate across tiles. Median and percentiles do not.
- Choose the cell size from the occupancy plateau, not from nominal spacing.
- Percentiles need a sort per cell; loop over cells, never over points.
Internal links
- Point density explained — the occupancy sweep in detail
- How to create a DTM from LiDAR ground points in Python — the most common application
- DSM, DTM and CHM from LiDAR: how each surface is derived — which reduction each surface needs
- How to rasterize a vector layer in Python — the vector equivalent
- My LiDAR DTM has holes, spikes or terraces — what thin cells produce
- How to read a LAS or LAZ file in Python with laspy — getting the arrays
- How to clip and tile a point cloud in Python — processing large clouds in pieces
- The raster data model explained — the transform this produces
FAQ
How do I rasterise a point cloud in Python?
Compute integer row and column indices by dividing the offset coordinates by the cell size, flatten to one index, and reduce with np.add.at, np.minimum.at or np.maximum.at.
Why is my count raster too low?
You used out[flat] += 1. Plain fancy indexing applies a repeated index once. Use np.add.at(out, flat, 1.0) or np.bincount(flat).
Why is my raster upside down?
Raster row 0 is the top and y increases upward, so the row index must be (top - y) / cell, not (y - bottom) / cell.
What cell size should I use?
The finest where the surveyed cells are occupied. On this survey that was 1 m: a third of 0.25 m cells and 1.5% of 0.5 m cells were empty, and none at 1 m. Count only cells inside the survey's footprint, or the bounding box's empty corners look like a floor.
Which reduction should I use?
Minimum for ground, maximum for surfaces, mean for continuous attributes, count for density, percentiles for distributions. Each is wrong in a specific direction.
How do I rasterise a cloud too large for memory?
Accumulate across tiles and chunks. Minimum, maximum, sum and count are associative and can be built incrementally; median and percentiles cannot.
Should I output the point count?
Yes. A cell built from one point and one built from a hundred look identical otherwise, and the difference is the only clue to reliability.