Python LiDAR

Calculating Point Density for Drone Surveys

TL;DR: Divide the count of valid LiDAR returns by the survey area in square metres — then replace the global average with a 1 m² grid map and use the 10th-percentile cell density as your compliance threshold.

# Context and Motivation

This guide is part of Point Density Metrics, which covers the broader methodology for computing and validating spatial density across LiDAR datasets.

UAV LiDAR acquisitions rarely distribute points uniformly. Multi-rotor and fixed-wing platforms generate overlapping flight swaths, creating density spikes along track centrelines and troughs in the overlap margins. A single global average — total_points / survey_area_m² — hides exactly the zones where surface reconstruction will fail. For operational drone workflows, target densities generally sit between 50 and 500 pts/m², with infrastructure inspection, corridor mapping, and urban asset modelling requiring the upper end to resolve fine features like utility poles, curb returns, and understory vegetation.

The diagram below shows the five-stage computation path from raw LAS input to a compliance decision.

Drone Survey Point Density Pipeline Five sequential stages for computing point density from a UAV LiDAR file: 1) LAS Input with CRS check, 2) Filter by ASPRS class, 3) Grid-bin XY at 1 m², 4) Compute p10/mean density, 5) Pass or fail threshold gate. LAS Input verify projected CRS (UTM) Class Filter keep codes 1,2,5 drop noise Grid Bin XY cell_size = 1.0 m np.bincount Stats: p10/mean pts/m² per cell 10th percentile Threshold Gate p10 < limit → fail export GeoTIFF chunked I/O prevents OOM on large tiles — no full in-memory load required

# Prerequisites and Assumptions

  • Python 3.10+ with a virtual environment (venv or conda)
  • laspy 2.4+ and numpy 1.24+ (install: pip install laspy[lazrs] numpy)
  • Input: LAS 1.2–1.4 or LAZ file with complete coordinate records and valid scale/offset headers; understanding of the LAS/LAZ file structure helps interpret header fields
  • Projected CRS: the file must already be in a metric coordinate system (UTM or state plane); use a spatial reprojection stage upstream if the source data is in geographic coordinates
  • ASPRS classification knowledge: knowing which ASPRS classification codes represent noise, ground, and vegetation allows you to pre-filter correctly before computing density

# Step-by-Step Implementation

# Step 1 — Validate the projected coordinate system

Geographic coordinates (latitude/longitude) yield mathematically valid but physically meaningless density values because one degree of longitude shrinks near the poles. Enforce a metric projection before any spatial calculation:

python
import laspy
import pyproj
from pathlib import Path

def assert_metric_crs(las_path: Path) -> pyproj.CRS:
    """Raise ValueError if the LAS file lacks a metric projected CRS."""
    with laspy.open(las_path) as f:
        crs = f.header.parse_crs()
        if crs is None:
            raise ValueError(f"{las_path.name}: no CRS in VLRs — cannot compute spatial density.")
        if not crs.is_projected:
            raise ValueError(
                f"{las_path.name}: geographic CRS detected ({crs.name}). "
                "Reproject to UTM or state plane before density calculation."
            )
    return crs

# Step 2 — Filter by ASPRS classification

Exclude low-noise (Class 7), high-noise (Class 18), and unclassified (Class 0) returns that would artificially inflate local counts. Keep Classes 1 (unclassified processed), 2 (ground), 5 (high vegetation), and 6 (building) for standard mapping deliverables:

python
import numpy as np

KEEP_CLASSES = [1, 2, 5, 6]  # adjust per project specification

def load_filtered_xy(las_path: Path, keep_classes: list[int] = KEEP_CLASSES) -> tuple[np.ndarray, np.ndarray]:
    """Load X, Y arrays retaining only the specified ASPRS classification codes."""
    with laspy.open(las_path) as f:
        x_chunks, y_chunks = [], []
        for chunk in f.chunk_iterator(5_000_000):
            mask = np.isin(chunk.classification, keep_classes)
            x_chunks.append(chunk.x[mask])
            y_chunks.append(chunk.y[mask])
    x = np.concatenate(x_chunks)
    y = np.concatenate(y_chunks)
    if len(x) == 0:
        raise ValueError("No points remain after classification filtering.")
    return x, y

# Step 3 — Bin XY coordinates into a density grid

numpy.bincount with linear indices is the fastest single-threaded approach for grids up to ~20,000 × 20,000 cells. The two-pass strategy (bounds scan then populate) keeps peak memory to one chunk at a time:

python
def compute_density_grid(
    x: np.ndarray,
    y: np.ndarray,
    cell_size: float = 1.0,
) -> tuple[np.ndarray, tuple[float, float, float, float]]:
    """
    Bin XY points into a density grid and return pts/m² per cell.

    Returns:
        density  : 2-D float32 array (rows × cols) in pts/m²
        bounds   : (x_min, y_min, x_max, y_max) in projected units
    """
    x_min, x_max = float(x.min()), float(x.max())
    y_min, y_max = float(y.min()), float(y.max())

    cols = max(1, int(np.ceil((x_max - x_min) / cell_size)))
    rows = max(1, int(np.ceil((y_max - y_min) / cell_size)))

    # Grid indices (Y flipped so row 0 = north, matching raster convention)
    col_idx = np.clip(((x - x_min) / cell_size).astype(np.int32), 0, cols - 1)
    row_idx = np.clip(((y_max - y) / cell_size).astype(np.int32), 0, rows - 1)

    flat_idx = row_idx * cols + col_idx
    counts = np.bincount(flat_idx, minlength=rows * cols).reshape(rows, cols)
    density = counts.astype(np.float32) / (cell_size ** 2)

    return density, (x_min, y_min, x_max, y_max)

# Step 4 — Derive percentile statistics

The 10th-percentile cell density (p10_density) is the industry-standard compliance metric because it reflects worst-case coverage, not average coverage:

python
def density_stats(density: np.ndarray) -> dict:
    """Compute summary statistics across non-empty grid cells."""
    populated = density[density > 0]
    if len(populated) == 0:
        return {}
    return {
        "mean_density":   float(np.mean(populated)),
        "median_density": float(np.median(populated)),
        "p10_density":    float(np.percentile(populated, 10)),
        "p90_density":    float(np.percentile(populated, 90)),
        "max_density":    float(np.max(populated)),
        "coverage_pct":   float(len(populated) / density.size * 100),
    }

# Step 5 — Assert threshold and surface results

python
def check_density_threshold(stats: dict, min_p10: float) -> None:
    """Raise AssertionError if p10 density falls below the project minimum."""
    actual = stats.get("p10_density", 0.0)
    assert actual >= min_p10, (
        f"Density QA failed: p10 = {actual:.1f} pts/m² < required {min_p10} pts/m²"
    )

# Complete Working Example

Paste this into a script alongside the helper functions above. It validates, filters, bins, and reports density for any LAS/LAZ file:

python
from pathlib import Path

def run_drone_density_check(
    las_path: str | Path,
    cell_size: float = 1.0,
    keep_classes: list[int] | None = None,
    min_p10_density: float = 50.0,
) -> dict:
    """
    End-to-end drone survey density check.

    Args:
        las_path         : Path to LAS 1.2–1.4 or LAZ file.
        cell_size        : Grid cell dimension in projected units (metres). Default 1.0.
        keep_classes     : ASPRS class codes to retain. Defaults to [1, 2, 5, 6].
        min_p10_density  : Minimum acceptable p10 density (pts/m²). Raises on failure.

    Returns:
        dict with density array, bounds, cell_size, shape, and stats sub-dict.
    """
    path = Path(las_path)
    if keep_classes is None:
        keep_classes = [1, 2, 5, 6]

    # 1. CRS guard
    crs = assert_metric_crs(path)
    print(f"CRS OK: {crs.name}")

    # 2. Filtered XY load
    x, y = load_filtered_xy(path, keep_classes)
    print(f"Retained {len(x):,} points after classification filter.")

    # 3. Grid density
    density, bounds = compute_density_grid(x, y, cell_size=cell_size)

    # 4. Statistics
    stats = density_stats(density)
    print(
        f"Density stats: mean={stats['mean_density']:.1f}, "
        f"p10={stats['p10_density']:.1f}, p90={stats['p90_density']:.1f} pts/m²  "
        f"(coverage {stats['coverage_pct']:.1f} %)"
    )

    # 5. Threshold assertion
    check_density_threshold(stats, min_p10_density)
    print(f"QA PASS: p10 = {stats['p10_density']:.1f} pts/m² ≥ {min_p10_density} pts/m²")

    return {
        "density_array": density,
        "bounds": {"x_min": bounds[0], "y_min": bounds[1],
                   "x_max": bounds[2], "y_max": bounds[3]},
        "cell_size": cell_size,
        "shape": density.shape,
        "stats": stats,
    }


if __name__ == "__main__":
    result = run_drone_density_check(
        las_path="survey_tile_utm32n.laz",
        cell_size=1.0,
        keep_classes=[1, 2, 5, 6],
        min_p10_density=100.0,   # engineering-grade threshold
    )

# Key Parameter Table

Parameter Type Default Guidance
cell_size float 1.0 Use 0.25–0.5 m for high-density urban scans; 2–5 m for coarse topographic tiles to reduce memory.
keep_classes list[int] [1, 2, 5, 6] Match to project spec. Remove Class 1 if only ground and structure returns matter for bare-earth deliverables.
min_p10_density float 50.0 Infrastructure inspection typically requires ≥100 pts/m²; topographic mapping 8–50 pts/m² depending on specification.
chunk_size (in load_filtered_xy) int 5_000_000 Lower to 1 M on machines with < 8 GB RAM; increase to 10 M on servers for faster throughput.

# Verification

After running the script, confirm the output is sensible before handing off to downstream processes:

python
result = run_drone_density_check("survey.laz", min_p10_density=100.0)

density = result["density_array"]
stats   = result["stats"]

# Sanity checks
assert density.shape[0] > 0 and density.shape[1] > 0, "Empty density grid."
assert stats["coverage_pct"] > 80.0, f"Low coverage: only {stats['coverage_pct']:.1f} % of grid populated."
assert stats["max_density"] < 10_000, "Implausibly high max density — check for duplicate points or wrong cell_size."

# Optional: export to GeoTIFF for spatial QA
try:
    import rasterio
    from rasterio.transform import from_bounds

    b = result["bounds"]
    transform = from_bounds(b["x_min"], b["y_min"], b["x_max"], b["y_max"],
                            density.shape[1], density.shape[0])
    with rasterio.open(
        "density_map.tif", "w",
        driver="GTiff", height=density.shape[0], width=density.shape[1],
        count=1, dtype="float32", crs="EPSG:32632", transform=transform,
    ) as dst:
        dst.write(density, 1)
    print("Exported density_map.tif")
except ImportError:
    print("rasterio not installed — skipping GeoTIFF export.")

Overlay the exported density_map.tif with flight-line tracks in QGIS to correlate density troughs with swath gaps or navigation errors.

# Gotchas and Edge Cases

1. Geographic CRS produces nonsensical density values. If assert_metric_crs is skipped and the file is in EPSG:4326, the binning will treat decimal degrees as metres. A 0.0001° cell at the equator covers ~11 m × 11 m, giving counts ~121× too low. Always run the CRS guard first, then apply a coordinate reference systems correction upstream if needed.

2. Classification filter removes all points. This happens when the LAS file stores classifications in the LAS 1.0/1.1 legacy bit-packed format rather than the dedicated classification byte. With laspy, chunk.classification reads the full byte, which on legacy files may equal zero for everything. Use chunk.classification & 0x1F for files produced before LAS 1.4 if all returns appear as Class 0.

3. Edge-cell clipping inflates border density. The np.clip call assigns boundary points to the last valid cell index. For surveys with clean convex hulls this is negligible, but for irregular corridor scans the first and last cells along each axis may appear 20–40 % denser than neighbours. Crop the density grid by one cell on each edge when computing statistics for narrow strips.

4. p10_density passes but the survey has dense stripes and empty gaps. High-overlap flight lines can push p10 above the threshold while leaving systematic gaps between swaths. Supplement percentile checks with spatial autocorrelation: flag any run of more than three contiguous empty cells in a row or column, which typically signals a missed strip rather than natural point scatter.

5. Large LAZ files exceed available RAM if loaded with laspy.read(). The complete working example uses chunk_iterator(5_000_000) in load_filtered_xy to avoid this. Never replace that with laspy.read(las_path).x on tiles above ~2 GB on machines with less than 32 GB RAM.

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...