Counting Points per Class with PDAL

TL;DR: For a quick look, pdal info --stats --enumerate Classification tile.laz lists the classes present. For counts, run a reader-only pipeline and np.bincount(arr["Classification"], minlength=256); for files that do not fit in memory, sum bincount over pipeline.iterator(chunk_size=...). Join the counts to ASPRS class names and compute shares for a QA table.

# Context and Motivation

This guide is part of ASPRS Classification Codes. A per-class point count is the most basic classification QA there is, and it answers more questions than it looks like it should. Is anything left unclassified? Did the vendor use classes the specification does not allow? Did ground classification fail on this tile — ground share of 8 percent where neighbours show 45? Did a reclassification step move the points it was meant to, and only those? Every one of those is a glance at a histogram, provided the histogram is easy to produce for one tile or a thousand.

A tile's classification at a glance Horizontal bars of point share by class for one suburban tile: ground 38 percent, high vegetation 27 percent, building 14 percent, medium vegetation 9 percent, low vegetation 7 percent, unclassified 3 percent, water 1.5 percent, and noise under 0.5 percent. Class codes and names label each bar. 2 ground5 high vegetation6 building4 medium vegetation3 low vegetation1 unclassified9 water7 + 18 noise 38 %27 %14 %9 %7 %3 %1.5 %0.4 %

# Prerequisites and Assumptions

  • PDAL 2.x with Python bindings; NumPy and pandas.
  • LAS files with classification populated.
  • A list of classes your specification allows, to flag unexpected ones.

# Step-by-Step Implementation

# Step 1 — List classes present

pdal info --stats --enumerate Classification computes statistics and reports the distinct classification values. It reads every point but needs no code.

# Step 2 — Count with bincount

A reader-only pipeline returns an array; np.bincount(a["Classification"], minlength=256) counts every class in one vectorized call.

# Step 3 — Stream for large files

pipeline.iterator(chunk_size=2_000_000) yields chunks; summing bincount across them gives exact counts in bounded memory.

# Step 4 — Add names and shares

Join counts to a dictionary of ASPRS class names and divide by the total.

# Step 5 — Aggregate a project

Run over every tile, keep per-tile rows for QA and sum for the project total.

# Complete Working Example

python
"""Per-class counts for tiles and a project, with ASPRS names and anomaly flags."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pandas as pd
import pdal

ASPRS = {0: "never classified", 1: "unclassified", 2: "ground", 3: "low vegetation",
         4: "medium vegetation", 5: "high vegetation", 6: "building", 7: "low noise",
         8: "model key point (legacy)", 9: "water", 10: "rail", 11: "road surface",
         12: "overlap (legacy)", 13: "wire guard", 14: "wire conductor",
         15: "transmission tower", 16: "wire connector", 17: "bridge deck",
         18: "high noise", 19: "overhead structure", 20: "ignored ground",
         21: "snow", 22: "temporal exclusion"}
ALLOWED = {1, 2, 3, 4, 5, 6, 7, 9, 17, 18}


def class_counts(path: Path, chunk: int = 2_000_000) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [str(path)]}))
    counts = np.zeros(256, dtype=np.int64)
    for arr in p.iterator(chunk_size=chunk):
        counts += np.bincount(arr["Classification"], minlength=256)
    return counts


def project_table(tiles: list[Path]) -> pd.DataFrame:
    rows = []
    for t in tiles:
        c = class_counts(t)
        total = c.sum()
        for cls in np.nonzero(c)[0]:
            rows.append({"tile": t.stem, "class": int(cls), "name": ASPRS.get(int(cls), "user-defined"),
                         "points": int(c[cls]), "share": c[cls] / total,
                         "allowed": int(cls) in ALLOWED})
    return pd.DataFrame(rows)


if __name__ == "__main__":
    df = project_table(sorted(Path("tiles").glob("*.laz")))
    project = (df.groupby(["class", "name"]).points.sum().reset_index()
                 .assign(share=lambda d: (d.points / d.points.sum()).round(4)))
    print(project.to_string(index=False))

    bad = df[~df.allowed].groupby("tile").points.sum()
    if len(bad):
        print("tiles with disallowed classes:\n", bad)
    ground = df[df["class"] == 2].set_index("tile").share
    median = ground.median()
    suspicious = ground[ground < 0.5 * median]
    print(f"median ground share {median:.1%}; suspicious tiles: {list(suspicious.index)}")

The streaming iterator requires a streamable pipeline — a bare reader always is — so the same function handles 10-million and 500-million-point files.

Spotting a failed tile from its histogram A bar chart of ground share for twelve tiles. Most are between 35 and 48 percent. One tile has only 9 percent ground, less than half the median, and is highlighted as suspicious: its ground classification likely failed, or it covers mostly water. 9 % ground share per tile · dashed line: half the median

# Turning Counts into a QA Report

A project-wide table of counts is useful; a short report built from it is what reviewers read. Four derived views cover most needs.

Project totals by class, with names and shares, as the example prints. This is the headline table for a delivery report and the first thing a client’s QA team compares with their expectations.

Tiles with disallowed classes. Any class outside the specification — class 0, legacy class 12, an undocumented vendor code — listed per tile with counts. These are usually quick to fix with conditional filters.assign or a vendor-code remapping, but only if someone sees them.

Outlier tiles per class. For ground, vegetation and buildings, tiles whose share is far from the median of their neighbours. A tile with 9 percent ground among neighbours at 40 percent almost always has a failed ground run, a large water body or a data gap, and each of those needs a different response.

Change between versions. When a project is reprocessed, a table of per-class differences between the old and new delivery shows at a glance what the reprocessing actually changed — and whether it changed anything it was not supposed to.

Keep the raw per-tile counts as a CSV beside the report so any figure can be recomputed.

# Key Parameter Table

Method Reads Memory Output
pdal info --stats --enumerate Classification all points low distinct values in JSON
np.bincount on pipeline.arrays all points whole tile exact counts
bincount over iterator(chunk_size) all points one chunk exact counts
filters.stats with count all points low per-value counts in metadata
laspy chunk_iterator all points one chunk exact counts without PDAL

# Verification

  • Totals match the header. The sum of all class counts equals the header point count.
  • Before/after diffs. After any reclassification step, only the classes you meant to change should differ.
  • Cross-tool agreement. On one tile, compare with laspy: np.bincount(laspy.read(p).classification) must be identical.

# Gotchas and Edge Cases

Legacy formats and 5-bit classes. In PDRF 0–5, PDAL exposes the 5-bit class; the withheld, synthetic and key-point bits are separate dimensions. Counts above class 31 are impossible there.

Class 0 versus class 1. Class 0 (never classified) means no classification was ever attempted; class 1 (unclassified) means a classifier looked and assigned nothing. Many specifications forbid class 0 in deliveries.

Shares mislead across landscapes. A tile over a lake has little ground; a tile over a city has little high vegetation. Compare shares against neighbouring tiles, not against a fixed number.

Never classified is not unclassified Two boxes. Class 0, never classified, means no classification process ran on the point, and most specifications reject it in deliveries. Class 1, unclassified, means a classifier examined the point and did not assign any specific class, which is acceptable in small amounts. class 0 — never classified no classifier ran on the point usually rejected in deliveries class 1 — unclassified classifier ran, assigned nothing acceptable in small shares

Withheld points. Decide whether withheld points belong in the counts. For QA of the delivery, count everything; for product statistics, exclude them.

# Frequently Asked Questions

How do I count points by classification in a LAS file?

Read the file with a PDAL pipeline and apply NumPy’s bincount to the Classification array, with minlength 256 so every code has a slot. For large files, sum bincount over chunks from the pipeline iterator.

How do I list the classes present with the PDAL command line?

Run pdal info with the stats flag and enumerate Classification. The output includes the distinct classification values found in the file.

What share of points should be ground?

It depends on land cover. Open agricultural land can exceed 60 percent ground; dense forest or urban centres can be below 20 percent. Compare each tile with its neighbours rather than with a fixed value.

Why does my file contain class 0?

Class 0 means the points were never classified, typically because a processing step was skipped or a tool wrote default values. Most delivery specifications do not allow it.