Compute Point Cloud Statistics and Density in PyQGIS

A LiDAR delivery is billions of points spread over hundreds of files, and its quality is invisible until it is measured. Was the specified density of 8 points per square metre achieved everywhere, or are there thin strips between flight lines? Are all the classes present — ground, vegetation, buildings — and in plausible proportions? Do elevations fall in a sensible range, or are there birds and multipath errors hundreds of metres in the air? Answering these questions takes minutes with QGIS's point cloud statistics and the PDAL algorithms, and saves days of trouble later.

This recipe belongs to Point Cloud & LiDAR Workflows. It reads summary statistics from a point cloud layer, gets per-file metadata with pdal:info, computes class distributions, builds a density raster and a coverage boundary, and checks density against a specification.

Four checks on a LiDAR deliveryFour checks: counts and attribute ranges, which catch outliers such as points hundreds of metres above the terrain; class distribution, which shows whether ground, vegetation and buildings were classified; density, which shows whether points per square metre meet the specification everywhere, including between flight lines; and coverage, which shows gaps and the exact boundary of the data.Measure before you trustrangesZ min / maxintensityoutliersclassesground, veg,buildingsproportionsdensitypoints per m²per cellflight-line gapscoverageboundarypolygonholes, edges

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series, with the PDAL Processing provider (included since 3.32).
  • Point cloud files — LAS, LAZ or COPC. Loading them is covered in loading a point cloud layer.
  • The delivery specification: required density, classes, vertical accuracy.

Read summary statistics from the layer

QGIS computes statistics when it indexes a point cloud and stores them with the layer. They give point count and per-attribute minimum, maximum, mean and standard deviation without reading all points again.

from qgis.core import QgsPointCloudLayer, QgsProject

pc = QgsPointCloudLayer("/data/lidar/tile_571_5934.copc.laz", "tile 571/5934", "pdal")
if not pc.isValid():
    raise RuntimeError(pc.error().summary())
QgsProject.instance().addMapLayer(pc)

print("points:", pc.pointCount())
stats = pc.statistics()
for attr in ("Z", "Intensity", "ReturnNumber", "Classification"):
    s = stats.statisticsOf(attr)
    print(f"{attr:<15} min {s.minimum:10.2f}  max {s.maximum:10.2f}  mean {s.mean:10.2f}")

Breakdown: pointCount() is the total in the file or index. statistics() returns the precomputed summary; for a freshly loaded layer it may still be computing in the background, in which case wait for the layer's statistics to finish or read them from pdal:info instead. A maximum Z hundreds of metres above the terrain points to noise — birds, clouds, multipath — which should be classified as noise (class 7 or 18) before building surfaces. Unexpected intensity ranges can reveal mixed sensors or uncalibrated strips.

Get per-file metadata with pdal:info

For a folder of tiles, pdal:info reports each file's header — point count, bounds, CRS, point format — which is enough to check a delivery's completeness and consistency.

import json
import processing
from pathlib import Path

rows = []
for laz in sorted(Path("/data/lidar/delivery").glob("*.laz")):
    out = processing.run("pdal:info", {"INPUT": str(laz), "OUTPUT": "TEMPORARY_OUTPUT"})
    info = json.loads(Path(out["OUTPUT"]).read_text()) if Path(str(out.get("OUTPUT", ""))).exists() else {}
    md = info.get("metadata", info)
    rows.append((laz.name, md.get("count"), md.get("minz"), md.get("maxz"),
                 md.get("srs", {}).get("horizontal", "")[:30] if isinstance(md.get("srs"), dict) else ""))
for r in rows[:5]:
    print(r)

Breakdown: Each file's header gives the facts needed for a delivery check: are all expected tiles present, do counts look similar for similar areas, is every file in the same CRS, are elevation ranges plausible. A tile with a fraction of the typical count is usually partially covered — at the project edge, over water — or truncated in transfer. The exact structure of the JSON varies with the PDAL version, so the code reads keys defensively; print one file's full output once to see what your version provides. Run in the background for large deliveries; reading headers is fast, but hundreds of files add up.

Class distribution

The share of points in each class shows whether classification was done and whether it is plausible: a rural tile with 2 % ground and 90 % "unclassified" was not classified, an urban tile without building points was classified with the wrong routine.

A plausible class mixBars for the share of points per ASPRS class in a sample rural tile: ground about 40 percent, low, medium and high vegetation together about 45 percent, buildings 8 percent, water 2 percent, noise under 1 percent, unclassified 4 percent. A tile with mostly unclassified points or no ground signals a classification problem.Share of points per class2 ground3–5 vegetation6 building9 water1 unclassifiedillustrative shares for a rural tile

from collections import Counter

cls_stats = stats.statisticsOf("Classification")
counts = getattr(cls_stats, "classCount", {}) or {}
total = sum(counts.values()) or pc.pointCount()
names = {1: "unclassified", 2: "ground", 3: "low veg", 4: "medium veg", 5: "high veg",
         6: "building", 7: "low noise", 9: "water", 17: "bridge", 18: "high noise"}
for cls, n in sorted(counts.items(), key=lambda kv: -kv[1]):
    print(f"{int(cls):>3} {names.get(int(cls), '?'):<14} {n / total:6.1%}")

Breakdown: For categorical attributes such as Classification, the statistics include a count per class. Mapping codes to ASPRS names makes the table readable. Plausibility depends on the landscape: a forest tile with little high vegetation, or a city tile with few building points, deserves a visual check. Noise classes should be small; a high share suggests a sensor or processing problem. Reclassifying with filtering and classifying point clouds fixes many issues.

Build a density raster

Average density hides the problem that matters: thin strips between flight lines, gaps over water, sparse coverage under dense canopy. A density raster shows points per cell everywhere.

dens = processing.run("pdal:density", {
    "INPUT": pc,
    "RESOLUTION": 5,                 # metres
    "TILE_SIZE": 1000,
    "FILTER_EXPRESSION": "Classification != 7 AND Classification != 18",
    "OUTPUT": "/data/lidar/qa/density_5m.tif",
})["OUTPUT"]

from qgis.core import QgsRasterLayer
import numpy as np
from osgeo import gdal

arr = gdal.Open(dens).ReadAsArray().astype("float64")
per_m2 = arr / (5 * 5)
valid = per_m2[per_m2 > 0]
print(f"density: median {np.median(valid):.1f} pts/m², 5th percentile {np.percentile(valid, 5):.1f}")

Breakdown: pdal:density counts points in each output cell; dividing by the cell area gives points per square metre. Excluding noise classes with the filter expression keeps outliers from inflating the count. A 5 m cell balances detail and noise: small enough to show flight-line gaps, large enough that random variation does not dominate. The 5th percentile is a better quality indicator than the mean — it shows how thin the thinnest parts are. Parameter names follow the PDAL provider; check processing.algorithmHelp("pdal:density") on your version.

Check against the specification

Specifications state a minimum density, usually as a requirement that a large share of cells meet it. A simple check computes that share and maps the failures.

Density against specificationThe density raster is compared with the specified minimum, for example 8 points per square metre. Cells below it are marked as failing. The share of passing cells is compared with the contractual threshold, such as 95 percent, and failing cells are mapped, usually forming strips between flight lines or patches over water.Share of cells meeting the minimumdensity rasterpts per m²≥ 8 pts/m²?per cellpass sharevs 95 % required

SPEC_MIN = 8.0
REQUIRED_SHARE = 0.95
passing = (valid >= SPEC_MIN).mean()
status = "PASS" if passing >= REQUIRED_SHARE else "FAIL"
print(f"{passing:.1%} of covered cells meet {SPEC_MIN} pts/m² — {status} "
      f"(required {REQUIRED_SHARE:.0%})")

fail_mask = processing.run("gdal:rastercalculator", {
    "INPUT_A": dens, "BAND_A": 1,
    "FORMULA": f"(A > 0) * (A / 25 < {SPEC_MIN})", "RTYPE": 0, "NO_DATA": 0,
    "OUTPUT": "/data/lidar/qa/density_fail.tif"})["OUTPUT"]

Breakdown: Only covered cells (density above zero) count, so areas outside the survey do not fail it. The pass share against the contractual threshold is the verdict; the failure mask shows where — typically systematic strips that point to a flight-planning or processing issue rather than random variation. Water bodies often fail legitimately because water absorbs the laser; exclude them with a water mask before judging. Combined with the class and range checks, this forms the core of a delivery report like the one in building a data quality report.

Map the coverage boundary

A polygon of where data exists is useful for indexing tiles, clipping products and showing coverage on overview maps. pdal:boundary traces it.

boundary = processing.run("pdal:boundary", {
    "INPUT": pc, "RESOLUTION": 10, "THRESHOLD": 4,
    "OUTPUT": "/data/lidar/qa/coverage.gpkg"})["OUTPUT"]
from qgis.core import QgsVectorLayer
cov = QgsVectorLayer(boundary, "coverage", "ogr")
area_km2 = sum(f.geometry().area() for f in cov.getFeatures()) / 1e6
print(f"covered area: {area_km2:.2f} km²")

Breakdown: The boundary is built from cells containing at least THRESHOLD points at the given resolution, which ignores scattered stray points beyond the real edge. Holes in the polygon reveal gaps inside the coverage. Comparing the covered area with the contracted area is a quick completeness check, and the polygon itself is a good clip mask for derived DEMs.

QGIS version compatibility

Point cloud statistics are available on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. The PDAL Processing provider with pdal:info, pdal:density and pdal:boundary ships with QGIS 3.32 and later. Parameter names and output formats can change with PDAL versions; check processing.algorithmHelp(...) and inspect one output before automating.

Troubleshooting

  • Statistics are empty. The layer is still indexing; wait or use pdal:info.
  • Density looks too high everywhere. Noise and duplicate points were counted; filter classes and check for overlapping tiles.
  • Strips of low density. Flight-line gaps — a real defect to report.
  • pdal:info output cannot be parsed. The output path or format differs by version; print the raw output first.

Conclusion

Read point counts and attribute ranges from layer statistics, collect headers for every file with pdal:info, check class shares for plausibility, build a density raster with noise excluded and judge it by low percentiles, compare the share of passing cells with the specification and map failures, and trace a coverage boundary for completeness and clipping.

Frequently Asked Questions

What resolution should the density raster use? Two to five times the expected point spacing — 5 m for typical 8 pts/m² data.

Should density count all returns or first returns only? Specifications usually state which; filter by ReturnNumber = 1 for first-return density.

Can I compute statistics for a whole delivery at once? Build a virtual point cloud or iterate files and aggregate the per-file results.

How do I find noise points? Statistics show the Z range; style by elevation and look for isolated high points, then classify them as noise.