Python LiDAR

Checking Pulse Spacing Against USGS Quality Levels

TL;DR: From single and first returns only, compute aggregate nominal pulse density (ANPD, pulses per m² over the land area) and aggregate nominal pulse spacing (ANPS ≈ 1/√ANPD), then grid first returns at a cell size of 2 × the required ANPS and check that at least 90 % of land cells contain a return. For QL2 that means ANPD ≥ 2 pulses/m², ANPS ≤ 0.71 m, and 1.42 m cells; QL1 and QL0 require ANPD ≥ 8 and ANPS ≤ 0.35 m.

# Context and Motivation

This guide is part of Point Density Metrics. Quality levels in the USGS Lidar Base Specification bundle accuracy and density requirements, and the density part is the easiest to get subtly wrong. It is defined on pulses, not points, so counting every return overstates it in vegetation. It is aggregated over land, so water must be excluded. And an average is not enough: a separate spatial-distribution test checks that the pulses are spread evenly, so a collection cannot pass by concentrating pulses in overlap stripes while leaving gaps between lines.

Doing the test in Python before submission avoids surprises from a reviewer and documents exactly how the numbers were produced.

Quality levels on one axis A horizontal log-scale axis of pulse density from 0.25 to 32 pulses per square metre. Threshold marks show QL3 at 0.5, QL2 at 2, and QL1 and QL0 at 8. A measured project value of 3.4 is plotted between QL2 and QL1, so the project meets QL2 density but not QL1. QL3 ≥ 0.5QL2 ≥ 2QL1/QL0 ≥ 8 0.250.512481632 measured ANPD 3.4 pulses per m², log scale; confirm thresholds against the specification edition in your contract meets QL2 density, not QL1

# Prerequisites and Assumptions

  • LAS data with ReturnNumber and NumberOfReturns populated, and water classified or supplied as polygons.
  • Tiles covering the whole project, or a representative set, in a metric projected CRS.
  • The specification edition your contract cites; values below follow the published USGS quality-level tables.
  • PDAL and Python with NumPy and rasterio.

# Step-by-Step Implementation

# Step 1 — Keep one return per pulse

Filter to first returns (ReturnNumber == 1), which includes single returns. Each pulse contributes exactly one first return.

# Step 2 — Exclude water and withheld points

Drop class 9 and any withheld points; the specification assesses density over land.

# Step 3 — Compute ANPD and ANPS

ANPD = first returns ÷ land area. ANPS ≈ 1/√ANPD for a roughly uniform pattern.

# Step 4 — Run the spatial-distribution test

Grid first returns at 2 × the required ANPS for the target quality level (1.42 m for QL2) and compute the share of land cells containing at least one return; the requirement is commonly 90 %.

# Step 5 — Report per tile and overall

Tile-level results show where a collection fails; the aggregate decides compliance.

# Complete Working Example

python
"""ANPD/ANPS and the spatial-distribution cell test for a LiDAR tile."""
from __future__ import annotations

import json

import numpy as np
import pdal

QL = {"QL0": (8.0, 0.35), "QL1": (8.0, 0.35), "QL2": (2.0, 0.71), "QL3": (0.5, 1.41)}


def first_returns_on_land(tile: str) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        tile,
        {"type": "filters.expression",
         "expression": "ReturnNumber == 1 && Classification != 9 && Classification != 7 "
                       "&& Classification != 18 && Withheld == 0"},
    ]}))
    p.execute()
    return p.arrays[0]


def density_report(tile: str, level: str = "QL2", water_area_m2: float = 0.0) -> dict:
    a = first_returns_on_land(tile)
    need_anpd, need_anps = QL[level]
    cell = 2 * need_anps
    x0, y0 = a["X"].min(), a["Y"].min()
    cols = int(np.ceil((a["X"].max() - x0) / cell)) + 1
    rows = int(np.ceil((a["Y"].max() - y0) / cell)) + 1
    occupied = np.zeros((rows, cols), dtype=bool)
    occupied[((a["Y"] - y0) // cell).astype(int), ((a["X"] - x0) // cell).astype(int)] = True

    # Land area: cells inside the tile minus water; approximated here from the grid extent.
    tile_area = (a["X"].max() - x0) * (a["Y"].max() - y0)
    land_area = tile_area - water_area_m2
    anpd = len(a) / land_area
    anps = 1 / np.sqrt(anpd)
    share = occupied.mean()                      # water cells should be masked out, see below
    return {"level": level, "first_returns": len(a),
            "anpd": round(anpd, 2), "anps_m": round(anps, 2),
            "cell_m": cell, "occupied_share": round(float(share), 3),
            "pass_density": anpd >= need_anpd and anps <= need_anps,
            "pass_distribution": share >= 0.90}


if __name__ == "__main__":
    print(density_report("tiles/t_0431.laz", "QL2", water_area_m2=38_500))

The occupied share above uses every cell in the tile’s bounding box; cells that are water have no returns by nature and would count as failures. In production, rasterize the water polygons to the same grid and exclude those cells from both numerator and denominator.

The cell test at twice the pulse spacing A grid of cells sized at twice the required pulse spacing. Most cells contain at least one first return and are marked occupied. A band of empty cells between two flightlines is marked as failing. A block of cells over a lake is excluded from the test. The share of occupied land cells must reach 90 percent. occupied: ≥ 1 first return empty land cell: counts against water: excluded from the test cell size = 2 × required ANPS (1.42 m for QL2) pass when ≥ 90 % of land cells are occupied

# Key Parameter Table

Quality level ANPS (m) ANPD (pulses/m²) Test cell (2 × ANPS)
QL0 ≤ 0.35 ≥ 8 0.70 m
QL1 ≤ 0.35 ≥ 8 0.70 m
QL2 ≤ 0.71 ≥ 2 1.42 m
QL3 ≤ 1.41 ≥ 0.5 2.82 m

# Verification

  • First returns only. Confirm the filter by counting: the number of first returns should equal the number of pulses, and be well below the total point count in vegetated areas.
  • Water excluded. Compare the occupied share with and without the water mask on a tile with a lake; the unmasked value should be visibly lower.
  • Spot the failing tiles. Map the per-tile occupied share. Failures cluster at flightline gaps, steep terrain facing away from the scanner, and project edges.

# Gotchas and Edge Cases

All returns inflate density. Counting every return in forest can double or triple apparent density. The specification is about pulses; use first returns.

Overlap-flagged points. Some deliveries flag sidelap points as overlap (or older ones class 12). Whether to include them depends on the specification’s wording; be explicit in the report.

Tile edges. Cells cut by the tile boundary look sparsely occupied. Run the test on a merged mosaic, or ignore a border of one cell.

Points are not pulses Paired bars for open ground and forest. In open ground, point density and pulse density are both about 4 per square metre because each pulse returns once. In forest, point density is about 11 per square metre while pulse density is still about 4, because each pulse returns from several canopy layers. open ground forest 11 points pulses

ANPS from a formula. 1/√ANPD assumes a uniform pattern. Linear scanners produce anisotropic spacing — denser along scan lines than between them — which the cell test catches even when the average passes.

# Frequently Asked Questions

What is aggregate nominal pulse spacing?

A typical distance between pulses across the collection, derived from first returns over land. For a roughly uniform pattern it is approximately one over the square root of the aggregate nominal pulse density.

Why use first returns rather than all points?

Density requirements are defined per emitted pulse. Each pulse has exactly one first return, while vegetation can produce several returns per pulse, so counting all points overstates density.

What is the spatial distribution test?

A grid with cells twice the required pulse spacing is laid over the land area, and the share of cells containing at least one first return is computed. Specifications commonly require at least 90 percent, which catches gaps that an average would hide.

Does water count against density?

No. Water often returns few or no pulses, and density requirements are assessed over land. Mask water cells out of the calculation.

Related in Point Density Metrics in Python LiDAR Workflows

↑ Point Density Metrics in Python LiDAR Workflows overview --- title: "Point Density Metrics in Python LiDAR Workflows" description: "Compute, validate, and normalize point density metrics from... Building a Point Density Raster with PDAL --- title: "Building a Point Density Raster with PDAL" description: "Map LiDAR point and pulse density as a GeoTIFF with writers.gdal... Measuring Ground Point Density Under Canopy --- title: "Measuring Ground Point Density Under Canopy" description: "Why total point density says nothing about a forested block, how to... Calculating Point Density for Drone Surveys --- title: "Calculating Point Density for Drone Surveys" description: "Step-by-step Python guide to calculating point density for UAV LiDAR...