Building a Point Density Raster with PDAL

TL;DR: Run writers.gdal with "output_type": "count", a resolution of 5–10 m and "radius" equal to half the cell diagonal, then divide the raster by the cell area to get points per square metre. Filter to ReturnNumber == 1 first for pulse density, or to class 2 for ground density. Cells far below the median are density voids worth reporting.

# Context and Motivation

This guide is part of Point Density Metrics. A single number — “the project averages 12 pts/m²” — hides everything that matters about density: the gaps between flightlines, the stripes where overlap doubles it, the lakes where it drops to nothing, the forests where ground density collapses under canopy. A density raster shows all of that at a glance, and it is the evidence that specifications ask for when they require density “uniform” or “without voids larger than” some size. PDAL produces one in a single stage, because writers.gdal can write a per-cell count instead of an interpolated elevation.

# Prerequisites and Assumptions

• PDAL 2.x with writers.gdal.
• A projected CRS in metres, so cell area is meaningful.
• Clarity about what is being counted: all points, pulses (first returns), or ground returns.

# Step-by-Step Implementation

# Step 1 — Decide what to count

All returns measure point density; first returns (ReturnNumber == 1) measure pulse density, which specifications usually mean; ground returns (class 2) measure bare-earth density, the one that limits DTM quality.

# Step 2 — Choose a cell size

Big enough that every cell has many points (5–10 m for 2–20 pts/m²), small enough to reveal gaps you care about. A specification about voids larger than a given size implies a cell no larger than that size.

# Step 3 — Count per cell

writers.gdal with output_type: "count" counts points within radius of each cell centre. Set the radius to half the cell diagonal (0.707 × resolution) so that circles cover the cell without big overlaps.

# Step 4 — Convert to density

Divide counts by the cell area. Because circular windows overlap slightly, the result is a close estimate, not an exact partition; for exact per-cell counts, bin in NumPy.

# Step 5 — Summarize and find voids

Median density, the share of cells below a threshold, and connected regions of low density give the numbers a report needs.

# Complete Working Example

{
  "pipeline": [
    "tiles/t_0431.laz",
    { "type": "filters.range", "limits": "Classification![7:7],Classification![18:18]" },
    { "type": "filters.expression", "expression": "ReturnNumber == 1" },
    { "type": "writers.gdal", "filename": "density/t_0431_pulse_count.tif",
      "resolution": 5.0, "radius": 3.54, "output_type": "count",
      "data_type": "uint32", "nodata": 0, "gdalopts": "COMPRESS=DEFLATE,TILED=YES" }
  ]
}

Then in Python, density and void statistics:

"""Pulse density statistics and low-density regions from a count raster."""
import numpy as np
import rasterio
from scipy import ndimage as ndi

with rasterio.open("density/t_0431_pulse_count.tif") as ds:
    counts = ds.read(1).astype(float)
    cell_area = abs(ds.res[0] * ds.res[1])

density = counts / cell_area
valid = density[density > 0]
median = float(np.median(valid))
low = density < 0.5 * median
labels, n = ndi.label(low)
sizes = ndi.sum(low, labels, range(1, n + 1)) * cell_area
print(f"median pulse density {median:.1f} /m²; {low.mean():.1%} of cells below half the median")
print(f"{int((sizes > 400).sum())} low-density regions larger than 400 m²")

# Key Parameter Table

# Verification

• Totals. The sum of counts is close to — slightly above, because windows overlap — the number of points passed to the writer.
• Median against the specification. Compare the median pulse density with the nominal pulse density of the collection.
• Void inspection. Overlay low-density regions on imagery: water and dark roofs are expected voids; a gap between flightlines over land is a coverage problem.

# Gotchas and Edge Cases

Counts versus density. Forgetting to divide by cell area makes a 5 m raster look 25 times denser than a 1 m raster of the same data. Always convert before comparing with specifications.

Edge cells. Cells at tile edges are only partly covered by data, so they read low. Exclude a one-cell border when summarizing, or compute density on a merged mosaic.

Overlap inflation. Pulse density doubles in sidelap. A project that meets its density target only in overlap stripes does not meet it; evaluate the distribution, not just the mean, as in checking pulse spacing against USGS quality levels.

Withheld and overlap-flagged points. Some specifications exclude withheld points and overlap-flagged points from density. Filter them out first if yours does.

# Frequently Asked Questions

How do I make a point density map with PDAL?

Use writers.gdal with output_type set to count, a resolution such as 5 or 10 metres and a radius of about 0.7 times the resolution. Divide the resulting counts by the cell area to get points per square metre.

What is the difference between point density and pulse density?

Point density counts every return; pulse density counts emitted pulses, estimated from first returns. Vegetation produces several returns per pulse, so point density can be much higher than pulse density in forests.

What cell size should a density raster use?

Large enough that typical cells hold dozens of points, commonly 5 to 10 metres for airborne data. If a specification limits void size, use a cell no larger than that size so voids are visible.

How do I measure ground point density?

Filter to classification 2 before writers.gdal. Ground density under canopy is often a small fraction of overall density and is the figure that limits DTM quality there.

# Related

• Point Density Metrics — definitions and methods
• Checking Pulse Spacing Against USGS Quality Levels — compliance testing
• Measuring Ground Point Density Under Canopy — the bare-earth case
• Calculating Point Density for Drone Surveys — high-density UAV data
• Classifying Water from Intensity and Returns — using density voids as evidence