🧭 Concept LiDAR & point clouds LiDAR & Point Clouds · step 6 of 14 Updated

Point Density Explained: What Resolution a Point Cloud Supports

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. Reprojected to UTM it is 16.9 — 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 about 3.9 points per square metre, and a bare-earth model can only resolve what the ground points support.

Quick answer

Reproject, then compute density per class:

import numpy as np
from pyproj import Transformer

transformer = Transformer.from_crs(3857, 32605, always_xy=True)
x, y = transformer.transform(las.x, las.y)

area = (x.max() - x.min()) * (y.max() - y.min())
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      16.9 per m²
ground points    3.9 per m²
nominal spacing 0.24 m
The same lidar patch measured as 1,613 metres across in Web Mercator and 877 metres on the ground, giving densities of 5.0 and 16.9 points per square metre.
Web Mercator at 58.6° north inflates area by a factor of 3.68. Density inherits the error exactly.

Step-by-step solution

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
true ground extent:       877 x   876 m
area: 260.1 ha reported,   70.7 ha real
density: 5.0 reported,     16.9 real

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 per m²
DTM ground only ~3.9 per m²
CHM both limited by the ground
building footprints building class none in this survey

A survey advertised at 16.9 points per square metre supports a bare-earth model at about 3.9 — a quarter of the headline.

3. Turn density into a usable cell size

The real test is what fraction of cells contain a point. Measured on the ground points:

cell size    ground-empty cells    any-return empty
   0.25 m           78.32%              38.23%
   0.50 m           38.11%               9.00%
   1.00 m           10.93%               7.47%
   2.00 m            7.51%               7.34%
   5.00 m            7.33%               7.29%

At 0.25 m, four in five cells have no ground point. At 1 m, one in nine. Below about 1 m this survey cannot produce a gap-free bare-earth grid.

4. Notice the floor, because it is not a resolution problem

The "any-return empty" column stops falling at about 7.3%. That is not sampling — it is a genuine hole in the data: water, which absorbs the near-infrared laser and returns nothing.

Any cell size below the floor is limited by resolution; at or above it, you are looking at real absence. Distinguishing the two is what stops you chasing a cell size that cannot exist.

5. Report density with its context

"16.9 points per square metre" is incomplete. The useful form states the CRS, the classes, and the fraction of cells occupied at the intended resolution.

Empty-cell fraction falling from 78% at 0.25 m to 11% at 1 m for ground points, with an irreducible floor at 7.3% from water.
The floor at 7.3% is water, not resolution. No cell size removes it.

Code examples

Example 1 — a density report that will not mislead 

import numpy as np
from pyproj import CRS, Transformer

CLASS_NAMES = {1: "unclassified", 2: "ground", 3: "low vegetation",
               4: "medium vegetation", 5: "high vegetation", 6: "building",
               9: "water", 11: "snow"}


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)

    area = (x.max() - x.min()) * (y.max() - y.min())
    print(f"  extent {x.max() - x.min():.0f} x {y.max() - y.min():.0f} m "
          f"= {area / 1e4:.1f} ha")
    print(f"  {len(x):,} points, {len(x) / area:.1f} per 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)}

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


def cell_size_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()

    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)

        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

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, 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))

    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)

    no_data = ~any_return
    no_ground = any_return & ~ground

    print(f"  cells with no return at all : {no_data.mean():6.2%} "
          f"(water, or outside the flight line)")
    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())}
  cells with no return at all :  7.47% (water, or outside the flight line)
  returns but no ground point :  3.46% (canopy or buildings blocking the ground)
  ground observed             : 89.07%

The split is the useful output. Cells with returns but no ground can sometimes be improved with a better ground filter; cells with no returns at all cannot be improved at all.

Explanation

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 16.9 points per square metre the spacing is 0.24 m, which suggests a 0.25 m grid should work. It does not, because the points are not evenly spread: at 0.25 m, 38% of cells had no return at all and 78% 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.

Overall density of 16.9 points per square metre against 10.0 first returns and 3.9 ground points.
Every product is limited by the density of the class it uses, not by the headline figure.

Edge cases or notes

  • 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.
  • The empty-cell floor is real absence, usually water. Below it is resolution; at it is data.
  • Water returns nothing to a near-infrared laser. Expect holes over it.
  • Nominal spacing is optimistic because points are not evenly spread.

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. On a 16.9 points per square metre survey, 0.25 m cells left 78% of cells without a ground point; 1 m left 11%.

Why do some cells have no points at all?

Water absorbs the near-infrared laser and returns nothing, and flight lines have edges. That floor — 7.3% here — does not improve with a coarser cell.

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.