How to Rasterise a Point Cloud to a Grid in Python
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 becomes unusably slow.
A Python loop over 12,968,770 points takes minutes. The vectorised equivalent takes under a second, and the difference is one NumPy idiom.
The other half of the problem is choosing a reduction and a cell size that the data supports. Measured on that survey:
0.25 m cells: 38.23% of cells contain no point at all
1.00 m cells: 7.47%
5.00 m cells: 7.29%
The floor at about 7.3% is water. Below it you are looking at resolution; at it, at real absence.
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)
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)
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)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.
Step-by-step solution
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
out[flat] += 1 # WRONG: repeated indices are applied once
np.add.at(out, flat, 1) # correctThis 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.
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.
Code examples
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)
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, transformThe 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, 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)
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 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
The empty fraction stops falling between 1 m and 2 m. That plateau is the honest resolution limit for this cloud.
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)
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)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.
Explanation
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.
Edge cases or notes
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.
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 one where the empty-cell fraction stops falling. On this survey that was between 1 m and 2 m, with a floor of about 7.3% from water.
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.