Point Density Explained: What Resolution a Point Cloud Supports
On this page
Problem statement
Point density decides what a lidar survey can produce. It also decides what it cannot, and the honest answer is usually coarser than people expect.
Two traps make the number itself unreliable.
The CRS. A real 3DEP survey stored in EPSG:3857 reports 5.0 points per square metre. Measured on the ground it is 18.3 — the same points over an area 3.68 times smaller, because Web Mercator inflates distance by 1/cos(latitude) and area by the square of that.
The class. Overall density is not the density that matters. That same survey is 22.9% ground, so its ground density is 4.2 points per square metre, and a bare-earth model can only resolve what the ground points support.
The area. Dividing by a bounding box is a third, quieter error. Reprojected to UTM, the survey's square footprint sits 2.4° off the grid, and its bounding box is 8% larger than the ground it covers.
Quick answer
Reproject, measure the area the points actually cover, then compute density per class:
import laspy
import numpy as np
from pyproj import Transformer
from scipy.spatial import ConvexHull
las = laspy.read("tile.laz")
transformer = Transformer.from_crs(3857, 32605, always_xy=True) # Web Mercator -> UTM 5N
x, y = transformer.transform(las.x, las.y)
area = ConvexHull(np.c_[x, y]).volume # a 2-D hull's "volume" is its area
ground = np.asarray(las.classification) == 2
print(f"all points {len(x) / area:6.1f} per m²")
print(f"ground points {ground.sum() / area:6.1f} per m²")
print(f"nominal spacing {np.sqrt(area / len(x)):.2f} m")all points 18.3 per m²
ground points 4.2 per m²
nominal spacing 0.23 m
Reproject first, density per class, a usable resolution
1. Reproject before you count
Web Mercator's scale factor is 1/cos(latitude) — 1.92 at 58.6° north, and worse further from the equator. Area is inflated by its square.
extent in Web Mercator: 1,613 x 1,613 "m"
UTM bounding box: 877 x 876 m 76.8 ha
footprint on the ground: 841 x 841 m 70.9 ha
density: 5.0 reported, 16.9 per bounding box, 18.3 real
(Example 1 prints these.) The bounding box is the trap after the reprojection: the patch is square in Web Mercator, and UTM zone 5's grid is turned 2.4° against it there, so a square footprint gets a box 8% larger than itself.
Any projected CRS suited to the area works — UTM is the usual choice. Geographic coordinates are worse still, since degrees are not a unit of length at all.
2. Compute density per class, not overall
The density that constrains a product is the density of the points that product uses:
| product | points it uses | density here |
|---|---|---|
| DSM | all first returns | 10.9 per m² |
| DTM | ground only | 4.2 per m² |
| CHM | both | limited by the ground |
| building footprints | building class | none in this survey |
A survey measured at 18.3 points per square metre supports a bare-earth model at 4.2 — under a quarter of the headline.
3. Turn density into a usable cell size
The real test is what fraction of cells contain a point — counting only cells inside the survey's footprint (Example 2):
cell size ground-empty cells any-return empty
0.25 m 76.53% 33.14%
0.50 m 32.99% 1.52%
1.00 m 3.67% 0.00%
2.00 m 0.13% 0.00%
5.00 m 0.00% 0.00%
At 0.25 m, three in four cells have no ground point. At 1 m, one in 27. Below about 1 m this survey cannot produce a gap-free bare-earth grid; at 2 m it very nearly can.
4. Count only the cells the survey covered
Run the same sweep over the bounding box and the "any-return empty" column stops falling at about 7.3%, whatever the cell size. That floor looks like a genuine hole in the data — water absorbing the laser, say — and it is not: it is the 7.7% of the box outside the survey's footprint, the corners left by a square turned 2.4°. Inside the footprint, at 1 m, every cell had a return.
A floor in an empty-cell sweep means some cells can never be filled; before explaining it, check whether they were ever inside the survey.
5. Report density with its context
"18.3 points per square metre" is incomplete. The useful form states the CRS, the classes, and the fraction of cells occupied at the intended resolution.
A density report that will not mislead
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, 12,968,770 points — joined into tile.laz. The archive stores it in Web Mercator.
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")287 nodes, 12,968,770 points, 82.8 MB of nodes
Example 1 — a density report that will not mislead
import laspy
import numpy as np
from pyproj import CRS, Transformer
from scipy.spatial import ConvexHull
CLASS_NAMES = {1: "unclassified", 2: "ground", 3: "low vegetation",
4: "medium vegetation", 5: "high vegetation", 6: "building",
7: "low noise", 9: "water"}
def density_report(x, y, classification, source_crs, target_crs=None):
"""Density per class, in a CRS where a square metre is a square metre."""
src = CRS.from_user_input(source_crs)
if not src.is_projected or src.to_epsg() == 3857:
target = target_crs or CRS.from_epsg(32605)
transformer = Transformer.from_crs(src, target, always_xy=True)
x, y = transformer.transform(x, y)
print(f" reprojected {src.to_string()} -> {target.to_string()}")
elif target_crs:
transformer = Transformer.from_crs(src, target_crs, always_xy=True)
x, y = transformer.transform(x, y)
box = (x.max() - x.min()) * (y.max() - y.min())
area = ConvexHull(np.c_[x, y]).volume # the ground the points cover
print(f" bounding box {x.max() - x.min():.0f} x {y.max() - y.min():.0f} m "
f"= {box / 1e4:.1f} ha; footprint {area / 1e4:.1f} ha")
print(f" {len(x):,} points, {len(x) / area:.1f} per m² "
f"({len(x) / box:.1f} per bounding-box m²), "
f"nominal spacing {np.sqrt(area / len(x)):.2f} m")
codes, counts = np.unique(classification, return_counts=True)
for code, count in sorted(zip(codes.tolist(), counts.tolist()),
key=lambda t: -t[1]):
print(f" {code:2d} {CLASS_NAMES.get(code, '?'):18} "
f"{count / area:7.2f} per m² ({count / len(x):6.2%})")
return {"area_m2": float(area), "density": float(len(x) / area)}
las = laspy.read("tile.laz")
wx, wy = np.asarray(las.x), np.asarray(las.y)
classification = np.asarray(las.classification)
print(f" in Web Mercator: {wx.max() - wx.min():,.0f} x {wy.max() - wy.min():,.0f}, "
f"{len(wx) / ((wx.max() - wx.min()) * (wy.max() - wy.min())):.1f} per square unit")
report = density_report(wx, wy, classification, "EPSG:3857")
first = np.asarray(las.return_number) == 1
print(f" first returns {first.sum() / report['area_m2']:.1f} per m²") in Web Mercator: 1,613 x 1,613, 5.0 per square unit
reprojected EPSG:3857 -> EPSG:32605
bounding box 877 x 876 m = 76.8 ha; footprint 70.9 ha
12,968,770 points, 18.3 per m² (16.9 per bounding-box m²), nominal spacing 0.23 m
1 unclassified 14.10 per m² (77.13%)
2 ground 4.18 per m² (22.87%)
7 low noise 0.00 per m² ( 0.00%)
first returns 10.9 per m²
Reprojecting automatically when the source is Web Mercator is the defensive touch that matters. It is the CRS most likely to be attached to a cloud served for the web, and the one that breaks density silently.
Example 2 — the cell-size sweep that chooses your resolution
import numpy as np
import shapely
def cell_size_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?
footprint: a polygon of the surveyed area; cells outside it are not counted.
"""
px, py = (x[mask], y[mask]) if mask is not None else (x, y)
left, top = x.min(), y.max()
print(f" {'cell':>6} {'cells':>12} {'empty':>8} {'median pts/cell':>16}")
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)
if footprint is not None:
cx, cy = np.meshgrid(left + (np.arange(width) + 0.5) * cell,
top - (np.arange(height) + 0.5) * cell)
counts = counts[shapely.contains_xy(footprint, cx, cy).ravel()]
rows.append({"cell": cell, "cells": counts.size,
"empty": float((counts == 0).mean()),
"median": float(np.median(counts))})
print(f" {cell:6.2f} {counts.size:12,} {(counts == 0).mean():7.2%} "
f"{np.median(counts):16.0f}")
return rows
x, y = Transformer.from_crs(3857, 32605, always_xy=True).transform(wx, wy)
hull = ConvexHull(np.c_[x, y])
footprint = shapely.Polygon(np.c_[x, y][hull.vertices])
ground = classification == 2
for label, kwargs in [("ground points, bounding box", {}),
("ground points, footprint", {"footprint": footprint}),
("any return, footprint", {"footprint": footprint, "mask": None})]:
print(label)
cell_size_sweep(x, y, mask=kwargs.get("mask", ground), footprint=kwargs.get("footprint"))ground points, bounding box
cell cells empty median pts/cell
0.25 12,288,524 78.32% 0
0.50 3,073,008 38.11% 1
1.00 768,252 10.93% 4
2.00 192,282 7.51% 16
5.00 30,976 7.33% 100
ground points, footprint
cell cells empty median pts/cell
0.25 11,349,660 76.53% 0
0.50 2,837,417 32.99% 1
1.00 709,354 3.67% 4
2.00 177,338 0.13% 17
5.00 28,376 0.00% 103
any return, footprint
cell cells empty median pts/cell
0.25 11,349,660 33.14% 1
0.50 2,837,417 1.52% 4
1.00 709,354 0.00% 17
2.00 177,338 0.00% 66
5.00 28,376 0.00% 413
Choose the finest cell size whose empty fraction is close to the floor. Finer than that and you are producing a grid whose holes you will then have to interpolate — which is a model, not a measurement.
Example 3 — separating resolution holes from real ones
import numpy as np
def hole_diagnosis(x, y, classification, footprint, cell=1.0):
"""Which empty cells are sampling, and which are genuine absence?"""
left, top = x.min(), y.max()
width = int(np.ceil((x.max() - left) / cell))
height = int(np.ceil((top - y.min()) / cell))
cx, cy = np.meshgrid(left + (np.arange(width) + 0.5) * cell,
top - (np.arange(height) + 0.5) * cell)
surveyed = shapely.contains_xy(footprint, cx, cy)
def occupancy(mask):
px, py = x[mask], y[mask]
col = np.clip(((px - left) / cell).astype(int), 0, width - 1)
row = np.clip(((top - py) / cell).astype(int), 0, height - 1)
grid = np.zeros(width * height, bool)
grid[row * width + col] = True
return grid.reshape(height, width)
any_return = occupancy(np.ones(len(x), bool))
ground = occupancy(classification == 2)
outside = ~surveyed
no_data = surveyed & ~any_return
no_ground = any_return & ~ground
print(f" cells outside the footprint : {outside.mean():6.2%} (never surveyed)")
print(f" surveyed, no return at all : {no_data.mean():6.2%} (water, dark roofs, gaps)")
print(f" returns but no ground point : {no_ground.mean():6.2%} "
f"(canopy or buildings blocking the ground)")
print(f" ground observed : {ground.mean():6.2%}")
return {"no_data": float(no_data.mean()), "no_ground": float(no_ground.mean())}
holes = hole_diagnosis(x, y, classification, footprint) cells outside the footprint : 7.67% (never surveyed)
surveyed, no return at all : 0.00% (water, dark roofs, gaps)
returns but no ground point : 3.46% (canopy or buildings blocking the ground)
ground observed : 89.07%
The split is the useful output. Every 1 m cell inside the footprint had a return; the 7.67% without one were never surveyed, and 3.46% had returns but no ground point. Those can sometimes be improved with a better ground filter; cells with no returns cannot, and cells outside the survey should not be counted at all.
Why the CRS error is so large
Why the CRS error is so large and so easy to miss
Web Mercator preserves angles and inflates distances by 1/cos(latitude). At the equator that is 1; at 58.6° north it is 1.92; at 70° it is 2.92.
The inflation applies to both axes, so areas are inflated by the square: 3.68 in the measurement above. Density, being points per area, is deflated by the same factor.
Nothing about the numbers looks wrong. A survey reporting 5 points per square metre is entirely plausible, and the file opens, plots and analyses normally. Only comparing the stated extent against a known ground distance reveals it.
Why nominal spacing is more intuitive than density
Density in points per square metre is hard to reason about. Nominal spacing — sqrt(area / n) — is the average distance between neighbouring points, and it maps directly onto resolution.
At 18.3 points per square metre the spacing is 0.23 m, which suggests a 0.25 m grid should work. It does not, because the points are not evenly spread: at 0.25 m, 33% of the surveyed cells had no return at all and 77% had no ground point.
Spacing gives an optimistic bound. The occupancy sweep gives the real one.
Why ground density is the binding constraint
Vegetation intercepts most pulses, so the ground gets whatever penetrates. In this survey 22.9% of returns were classified ground — under closed canopy the local figure will be far lower still.
Every bare-earth product is limited by that fraction, not by the headline density. A survey specified as "8 points per square metre" over forest may deliver 1 ground point per square metre in the parts that matter.
Why more density stops helping
Density improves a surface until it exceeds what the surface has structure at. Beyond that, extra points measure the same facet repeatedly, and the limiting factors become the ground classifier and the vertical accuracy rather than the sampling.
For bare-earth mapping, that point arrives around a few ground points per square metre. For canopy structure it is much higher, because the quantity of interest — the vertical distribution of returns — genuinely needs many pulses per square metre.
Geographic coordinates and reporting
- Reproject before computing density. Web Mercator understated it by a factor of 3.68 here.
- Never compute density in geographic coordinates. Degrees are not lengths.
- Report density per class. Ground density is what limits a DTM.
- Bounding-box area overstates coverage for irregular or partial tiles; use a concave hull if the shape matters.
- Flight-line overlap doubles local density, producing stripes. Check
point_source_id. - An empty-cell floor that never falls is usually cells outside the survey. Count only cells inside the footprint before calling it water.
- Water returns nothing to a near-infrared laser. Expect holes over it.
- Nominal spacing is optimistic because points are not evenly spread.
Internal links
- LiDAR point clouds explained: returns, classes and intensity — the classes density is computed over
- laspy returns the wrong coordinates: scale, offset and CRS — the Web Mercator trap
- DSM, DTM and CHM from LiDAR: how each surface is derived — what the density supports
- How to create a DTM from LiDAR ground points in Python — choosing a cell size
- How to rasterise a point cloud to a grid in Python — the occupancy grid in practice
- My LiDAR DTM has holes, spikes or terraces — what happens when the cell is too small
- How to choose a cell size for an interpolated surface — the same question for point interpolation
- LAS, LAZ and COPC explained — where the CRS is recorded
- Density depends on the boundary — the denominator problem for population density
FAQ
How do I calculate lidar point density?
Points divided by ground area, after reprojecting to a projected CRS suited to the area. In Web Mercator at 58.6° north the answer is understated by a factor of 3.68.
What point density do I need?
For a bare-earth DTM, a few ground points per square metre. For canopy structure, much more, because the quantity of interest is the vertical distribution of returns.
Why is my density lower than the survey specification?
Often because you computed it in the wrong CRS, or because you are looking at overall density where the specification meant a different class or a different area.
What cell size can my point cloud support?
Run an occupancy sweep over the surveyed cells. On an 18.3 points per square metre survey, 0.25 m cells left 77% of cells without a ground point, 1 m left 3.7% and 2 m 0.1%.
Why do some cells have no points at all?
Water can absorb the laser and return nothing, and flight lines have edges — but check the simplest cause first. Here the 7.3% that never filled at any cell size were the corners of a bounding box outside the survey; inside the footprint no 1 m cell was empty.
Does higher density always give a better surface?
No. Beyond a few ground points per square metre, the limiting factors become the ground classifier and vertical accuracy rather than the sampling.
Why is ground density so much lower than total density?
Because vegetation intercepts most pulses. This survey was 22.9% ground overall, and under closed canopy the local fraction is far lower.