Measuring Ground Point Density Under Canopy

TL;DR: Classify ground first, then compute density over the ground returns only, on the same grid you intend to rasterize. Total point density over forest is dominated by canopy hits and tells you nothing about whether a terrain model is possible.

# Context and Motivation

This guide is part of Point Density Metrics, which covers pulse density, point density and their per-cell distribution. This page is about the one measurement that actually predicts whether a bare-earth product will succeed: how many ground returns landed in each cell.

The distinction matters because acquisition specifications are almost always written in total points per square metre, and total density is met comfortably by a canopy. A forested block delivered at a contractual eight points per square metre may carry twelve returns per square metre in the crowns and 0.4 on the forest floor. A DTM at one metre needs the second number, and no amount of filter tuning invents returns that were never recorded — a point made in more detail under tuning SMRF for forested terrain.

# Prerequisites and Assumptions

# Step-by-Step Implementation

# Step 1 — Classify before measuring

{"type": "filters.smrf", "window": 33, "slope": 0.2, "threshold": 0.6, "cell": 1.0,
 "returns": "last, only"}

# Step 2 — Keep only ground

{"type": "filters.range", "limits": "Classification[2:2]"}

# Step 3 — Rasterize a count at the target cell size

{"type": "writers.gdal", "filename": "ground_count.tif",
 "output_type": "count", "resolution": 1.0, "nodata": 0}

The count reducer writes points per cell, which at a one-metre cell is points per square metre directly.

# Step 4 — Report the distribution, not the mean

The mean is flattered by the open parts of the block. The fraction of cells with zero ground returns is the number that predicts void area in the DTM.

# Complete Working Example

"""Measure ground-return density per cell and report the distribution."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pdal


def ground_counts(src: Path, cell: float = 1.0) -> np.ndarray:
    """Points per cell, ground only, on the target grid."""
    spec = json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.range", "limits": "Classification[2:2]"},
    ]})
    p = pdal.Pipeline(spec)
    p.execute()
    arr = p.arrays[0]
    if len(arr) == 0:
        raise ValueError("no ground-classified points — classify before measuring")

    ix = np.floor((arr["X"] - arr["X"].min()) / cell).astype(np.int64)
    iy = np.floor((arr["Y"] - arr["Y"].min()) / cell).astype(np.int64)
    nx, ny = ix.max() + 1, iy.max() + 1
    flat = np.bincount(iy * nx + ix, minlength=int(nx * ny))
    return flat.reshape(int(ny), int(nx))


def report(counts: np.ndarray, cell: float) -> dict:
    total_cells = counts.size
    empty = int((counts == 0).sum())
    occupied = counts[counts > 0]
    return {
        "cell_size_m": cell,
        "cells": total_cells,
        "empty_cells": empty,
        "empty_fraction": round(empty / total_cells, 4),
        "mean_over_occupied": round(float(occupied.mean()), 2),
        "p5": float(np.percentile(counts, 5)),
        "p50": float(np.percentile(counts, 50)),
        "supportable": bool(empty / total_cells < 0.05),
    }


if __name__ == "__main__":
    counts = ground_counts(Path("forest_tile.laz"), cell=1.0)
    summary = report(counts, cell=1.0)
    print(json.dumps(summary, indent=2))
    if not summary["supportable"]:
        print("coarsen the DTM cell size or accept interpolated voids")

# Key Parameter Table

# Verification

The classification actually ran. The example raises rather than reporting zero density, because an unclassified cloud and a treeless one produce the same number otherwise.

The grid matches the product. Measuring at two metres and building at one metre understates voids fourfold.

Open ground looks right. Ground density over an open field should be close to total density. If it is far below, the classifier is rejecting real ground and the density measurement is really a classifier problem.

# Gotchas and Edge Cases

Water returns nothing and is not a void. Mask water before computing the empty-cell fraction, or a lake makes a good block look unsupportable.

Overlap inflates density in strips. Flight-line overlap doubles the count where swaths meet, which raises the mean and leaves the fifth percentile unmoved — one reason to quote the percentile.

Density is a property of a cell size. Quoting “0.4 points per square metre” without a cell size is meaningless; the same cloud yields different numbers at different grids.

# Frequently Asked Questions

Why is total point density misleading over forest?

Because it is dominated by canopy returns. A block delivered at a contractual eight points per square metre can carry twelve per square metre in the crowns and 0.4 on the forest floor. The second number is what a bare-earth product depends on, and it is not the one in the contract.

Which statistic should I quote?

The empty-cell fraction and the fifth percentile, both at the cell size the product will use. The mean is raised by the open parts of a block and hides exactly the areas where a terrain model will fail.

Can better filter tuning fix low ground density?

No. Tuning recovers ground where sparse returns exist; it cannot invent returns in cells no pulse reached. Once ground density falls below roughly one return per cell, the only real remedies are a coarser grid, interpolation you declare as such, or a reflight.

Does density depend on the cell size I choose?

Yes, and quoting a density without one is meaningless. The same cloud yields 71 percent empty cells at half a metre and 12 percent at two metres, so the metric only means something alongside the grid it was measured on.

# Related

• Point Density Metrics — the parent guide to the three things called density
• Calculating Point Density for Drone Surveys — the acquisition-side view of the same measurement
• Tuning SMRF for Forested Terrain — recovering the ground returns that do exist
• Filling NoData Voids in DTM Rasters — what to do with the cells that stay empty
• Point Cloud Data Standards and Fundamentals — the section overview