Computing Point Density and Coverage Gaps from LAZ with PDAL

This page measures point density (points per square metre) and coverage uniformity from a .laz file with PDAL — using filters.hexbin to estimate density and trace the survey boundary, reading the EPSG code straight from the LAS header (EPSG:32618 in the examples), and reporting a median density plus the fraction of the extent that falls into gaps. The pipeline is a PDAL JSON stage list driven from the pdal Python bindings, with numpy only for the final reduction, so the same logic runs unchanged in a batch job over a whole delivery.

You hit this when a survey lands and you need one authoritative density figure before it enters the twin — and, more importantly, a map of where the cloud is too thin to trust. A headline mean hides dropped flight lines and occlusion shadows; PDAL’s filters.hexbin is purpose-built for both jobs because it estimates density and emits a boundary polygon whose interior holes are exactly the voids you care about. This complements the standalone laspy/numpy binning in the point cloud density standards reference and the per-asset targets in LiDAR point density best practices for infrastructure; here the tool is PDAL end to end.

PDAL hexbin density and boundary with a coverage gap Scattered LAZ points in EPSG:32618 are aggregated by filters.hexbin into hexagonal cells; occupied hexagons yield a density in points per square metre while an unfilled hexagon in the interior is reported as a hole in the boundary polygon. LAZ points · EPSG:32618 hexbin aggregate density (ppsm)+ boundarywith holes gap = hole
filters.hexbin aggregates points into hexagons, reports density per square metre, and returns a boundary polygon whose interior holes are the coverage gaps.

Prerequisites

  • PDAL 2.5+ with the Python bindings and numpy>=1.24. Install the whole stack from conda-forge so the LAZ backend and GEOS (which filters.hexbin needs for the boundary polygon) are matched: conda install -c conda-forge pdal python-pdal numpy.
  • A .laz or .las tile whose header declares a projected metric CRS. The examples assume EPSG:32618 (UTM 18N, metres); density computed in a geographic CRS such as EPSG:4326 is wrong by orders of magnitude because a degree is not a metre.
  • A density target to check against — for example a QL1 floor of 8 ppsm for urban work. The Quality Level table lives in the point cloud density standards reference.
  • The authoritative filters.hexbin and pipeline documentation at pdal.io if you need the full option list.

Step-by-Step

1. Read the CRS from the LAS header with PDAL

Run a metadata-only pipeline and pull the horizontal CRS from the reader stage. Refuse to proceed unless it is projected and metric — every density figure downstream depends on it.

python
import pdal
import json

LAZ = "survey_18N.laz"
info = pdal.Pipeline(json.dumps([{"type": "readers.las", "filename": LAZ}]))
info.execute()
meta = info.metadata["metadata"]["readers.las"]

srs_wkt = meta["srs"]["horizontal"]
print("declared CRS:", srs_wkt[:70])
assert "32618" in srs_wkt or "UTM zone 18N" in srs_wkt, \
    "expected EPSG:32618; reproject to a projected metric CRS before density"
print("header point count:", meta["count"])

2. Estimate density and boundary with filters.hexbin

filters.hexbin lays a hexagonal grid over the cloud, keeps hexagons holding at least threshold points, and reports an average density plus a boundary polygon. Set edge_size to the resolution you care about (1 m gives density directly comparable to a ppsm target).

python
import pdal
import json

pipeline = pdal.Pipeline(json.dumps([
    {"type": "readers.las", "filename": "survey_18N.laz"},
    {"type": "filters.hexbin", "edge_size": 1.0, "threshold": 1},
]))
pipeline.execute()

hb = pipeline.metadata["metadata"]["filters.hexbin"]
print("avg density (pts/m^2):", round(hb["density"], 2))
print("covered area (m^2):", round(hb["area"], 1))
print("avg point spacing (m):", round(hb.get("avg_pt_spacing", 0.0), 3))
Bounding-box area against the hexbin hull area The same curved flight swath measured two ways. Its bounding box encloses a large area of empty ground, so dividing the point count by that area understates density. The concave hull that filters.hexbin derives follows the swath, giving a much smaller area and a density figure that reflects where the points actually are. bounding box hexbin concave hull area 1.00 km² → 6.1 points/m² area 0.62 km² → 9.8 points/m² Same file, same points — a 38% difference in the number you would report against the specification
Area is the denominator, so the hull choice is not cosmetic. A bounding box over a curved or L-shaped block reports a density the survey never had.

3. Turn the boundary into a coverage gap fraction

The hexbin boundary is a WKT polygon; interior rings (holes) are voids the survey never covered. Compare the polygon’s net area against its outer-ring area with shapely to get the fraction of the footprint lost to gaps.

python
from shapely import wkt

poly = wkt.loads(hb["boundary"])              # POLYGON, holes = coverage gaps
from shapely.geometry import Polygon
outer = Polygon(poly.exterior)                 # footprint ignoring holes

gap_fraction = 1.0 - (poly.area / outer.area)
print(f"coverage gaps: {gap_fraction:.2%} of the footprint")
print("gap (hole) count:", len(list(poly.interiors)))

4. Compute a median density and empty-cell fraction over a grid

The hexbin average is a single number; acceptance also needs the median and the worst served regions. Read the point coordinates PDAL already streamed, bin them into a metric grid, and reduce over occupied cells. This uses PDAL as the reader and keeps the reduction to a compact numpy step — the standalone binning routine is detailed in the density standards reference.

python
import numpy as np

arr = pipeline.arrays[0]                        # structured array from the executed pipeline
x, y = np.asarray(arr["X"]), np.asarray(arr["Y"])   # metres, EPSG:32618

cell = 1.0
minx, miny = x.min(), y.min()
cols = int(np.ceil((x.max() - minx) / cell)) + 1
rows = int(np.ceil((y.max() - miny) / cell)) + 1
ci = np.clip(((x - minx) / cell).astype(np.int64), 0, cols - 1)
ri = np.clip(((y - miny) / cell).astype(np.int64), 0, rows - 1)

counts = np.zeros((rows, cols), dtype=np.int64)
np.add.at(counts, (ri, ci), 1)
occupied = counts[counts > 0] / (cell * cell)   # ppsm

print("median density:", round(float(np.median(occupied)), 2), "ppsm")
print("p05 density:", round(float(np.percentile(occupied, 5)), 2), "ppsm")
print("empty-cell fraction:", round(float((counts == 0).mean()), 4))

The gap fraction and the median tell you different things, and a pipeline needs both. The median answers “is the survey dense enough where it has data” — it is robust to a handful of extreme cells and it maps directly onto the acquisition specification. The empty-cell fraction answers “does it have data everywhere”, which no density average can express, because a cell with zero points contributes nothing to a mean over occupied cells and vanishes entirely from a mean over the hull.

What turns those two numbers into a decision is the arrangement of the failures. Scattered single empty cells are normal: water bodies absorb the pulse, wet asphalt scatters it, and a specular roof returns nothing at certain incidence angles. Those are properties of the surface and reflying will not change them. Contiguous runs of empty cells are different — they are occlusions behind a structure, a gap between flight lines, or a trajectory dropout, and every one of them is fixable by another pass. So the check worth automating is not simply empty_fraction < 0.02; it is the size of the largest connected empty region, which scipy.ndimage.label will give you in two lines over the same count grid.

Per-cell counts over a one-metre grid A grid of per-cell point counts across a tile. Most cells sit comfortably above the requirement, a scattering fall between eight and twelve, and four cells hold no points at all. The empty cells are adjacent, which is the signature of an occlusion rather than of random sparsity. Points per cell — the two statistics that matter are the median and the empty fraction 12 18 17 11 14 18 16 19 18 10 18 9 16 13 17 12 12 20 16 17 17 16 15 19 11 12 0 0 19 11 17 15 20 9 19 21 10 11 0 21 18 9 13 21 9 13 16 18 20 15 20 21 0 15 15 20 21 18 16 11 14 10 9 11 16 12 median 15 points/m² · 4 of 66 cells empty (6.1%) · no cell between 1 and 8 The empty cells touch each other — an occlusion under a structure, not thin coverage
Counting is the easy part. Reading the arrangement is what separates a survey that needs reflying from one that needs a note in the metadata.

5. Gate the result against the acquisition target

Fold the figures into a pass/fail a CI job can act on: the median must clear the contracted density and the gaps must stay under an acceptance threshold.

python
TARGET_PPSM = 8.0          # QL1 urban floor
MAX_GAP = 0.02             # 2% of footprint

median = float(np.median(occupied))
status = "PASS" if (median >= TARGET_PPSM and gap_fraction <= MAX_GAP) else "FAIL"
print(f"median {median:.1f} ppsm, gaps {gap_fraction:.2%} -> {status}")

Expected Output & Verification

A healthy QL1 urban tile in EPSG:32618 should print an average and median near or above target with a small gap fraction:

text
declared CRS: PROJCS["WGS 84 / UTM zone 18N", ...
header point count: 41880233
avg density (pts/m^2): 11.4
covered area (m^2): 998840.0
avg point spacing (m): 0.29
coverage gaps: 0.8% of the footprint
median density: 11.1 ppsm
p05 density: 6.9 ppsm
empty-cell fraction: 0.0121
median 11.1 ppsm, gaps 0.80% -> PASS

Cross-check the two independent estimates: the hexbin avg_pt_spacing (~0.29 m) and the grid median (~11 ppsm) should agree through spacing ≈ 1/√density — here 1/√11.1 ≈ 0.30 m, a match that confirms the figure is real rather than an artifact of one method’s cell geometry. If the hexbin density and the grid median disagree by more than a few percent, the cloud is clustered along flight lines; if either comes back near zero, the CRS is still geographic — recheck step 1.

It is also worth persisting the count grid itself, not just the summary statistics. The grid is small — a few hundred kilobytes for a survey tile — and it is what lets a later question be answered without re-reading the LAZ: which cells were sparse, whether the sparse region moved between acquisitions, and whether a clearance failure sits inside one. Summary statistics answer the acceptance question and nothing else.

Common Errors

PDAL: filters.hexbin: Unable to compute boundary or a missing boundary key. PDAL was built without GEOS, so hexbin can estimate density but cannot trace the polygon. Install the conda-forge pdal build (which links GEOS), or drop step 3 and rely on the grid empty-cell fraction from step 4 for gap detection.

Density reported as a tiny fraction such as 0.000004 pts/m². The header CRS is geographic (EPSG:4326), so PDAL measured area in square degrees. The assertion in step 1 is meant to catch this; reproject with filters.reprojection (in_srs/out_srs set to explicit EPSG codes, e.g. to EPSG:32618) before the hexbin stage.

KeyError: 'filters.hexbin' when reading metadata. The pipeline was constructed but never executed, or the stage name was misspelled, so no metadata block exists. Call pipeline.execute() before touching pipeline.metadata, and read the block under the exact stage key filters.hexbin.

Back to Point Cloud Density Standards.