LiDAR Classification and Feature Extraction in Python

Ground classification answers one question about every point — is it the bare earth or not? — and most LiDAR deliverables need a great deal more than that. An asset owner wants building footprints with roof heights, a utility wants every conductor span and every tree within striking distance of it, a forester wants stems per hectare and crown diameters, and a hydrologist wants lakes flattened and bridges removed from the terrain before water is routed across it. This section is for the LiDAR analysts, Python GIS developers and surveying teams who have to deliver those products from the same classified tiles. It covers the per-point geometric features PDAL can compute, the segmentation stages that group points into candidate objects, the rules and models that decide what each object is, and the checks that tell you whether the classification is good enough to ship.

# Why Feature Extraction Is Its Own Discipline

A LiDAR tile that has been through SMRF ground classification is already worth something: it yields a terrain model and, with a little more work, a surface model. But everything above the ground is still lumped together as “not ground”, and the questions people actually pay for — how many buildings, how tall, how close is that oak to the 132 kV line — live in that lump. Extracting them is a different kind of problem from ground filtering, for three reasons.

First, the unit of reasoning changes from the point to the object. SMRF can decide about each point by comparing it with an interpolated surface; deciding that a set of points is a roof requires knowing that they are coplanar, contiguous, elevated and large enough, which are properties of the set, not of any member. That is why this section leans so heavily on segmentation: until points are grouped, most object rules cannot even be expressed.

Second, the evidence is geometric and local. Intensity and return number help at the margins, but the dependable signal for most above-ground classes is the shape of each point’s neighbourhood — flat for a roof, linear for a wire, volumetric and scattered for a crown. PDAL computes those shape descriptors from the eigenvalues of each neighbourhood’s covariance matrix, and a surprisingly large share of practical classification is thresholding them sensibly.

Third, the output is usually vector or per-object, not a raster. A footprint layer, a tree table, a list of conductor spans with clearance violations: these are what people consume, and the point cloud is the evidence behind them. The pages here therefore end in GeoPackages and data frames as often as in LAS files.

The ASPRS classes you will be writing are covered in ASPRS classification codes: 6 for buildings, 3 to 5 for vegetation by height, 9 for water, 13 to 16 for the wire and tower family, 17 for bridge decks. Getting the codes right matters because every downstream consumer — a DTM builder, a GIS viewer, a client QA script — keys on them.

# Conceptual Architecture: Normalize, Describe, Group, Decide

Every workflow in this section follows the same four moves, in the same order, and most failures come from skipping or reordering one.

1. Normalize heights. Compute HeightAboveGround with filters.hag_nn or filters.hag_delaunay so that “tall” means tall relative to the terrain, not relative to sea level. A two-storey house on a hillside and one on a plain then look alike.
2. Describe neighbourhoods. Run filters.covariancefeatures, filters.normal or filters.eigenvalues to attach per-point shape descriptors — linearity, planarity, scattering, verticality, normal vectors.
3. Group into candidates. Use filters.cluster (Euclidean) or filters.dbscan to give every point in a contiguous object the same ClusterID, usually after restricting the input to a height band.
4. Decide per object. Aggregate the descriptors per ClusterID in pandas — median planarity, height range, footprint area, point count — and assign an ASPRS class or reject the group. Rules first; a trained model when rules plateau.

The order is not arbitrary. Features computed before normalization mix terrain slope into verticality. Segmentation before height filtering merges a roof with the lawn around it. And deciding per point instead of per object produces the speckle — isolated “building” points in a tree crown — that makes a classification look amateur however good its aggregate accuracy.

# Core Components

The stages below do almost all of the work in this section. None of them classifies anything on its own; they produce the evidence that your rules or your model will read.

Two practical notes. knn is the parameter that matters most in every feature stage: too small and descriptors are noisy, too large and neighbourhoods straddle object boundaries. Scale it to point density — a value that works at 8 pts/m² is far too small at 60 pts/m². And every one of these stages needs a spatial index over the whole input, which means none of them streams; the implications are covered in which filters break streaming mode.

# Annotated Reference Pipeline

The pipeline below takes a ground-classified tile to the point where every candidate building has a ClusterID and every point carries the descriptors needed to decide. JSON has no comments, so the reasoning for each stage follows the block.

{
  "pipeline": [
    { "type": "readers.las", "filename": "tile_4710_5285.laz" },
    { "type": "filters.range", "limits": "Classification![7:7],Classification![18:18]" },
    { "type": "filters.hag_nn", "count": 2 },
    { "type": "filters.range", "limits": "HeightAboveGround[2.5:120]" },
    {
      "type": "filters.covariancefeatures",
      "knn": 16,
      "threads": 4,
      "feature_set": "Dimensionality"
    },
    { "type": "filters.normal", "knn": 12, "always_up": true },
    {
      "type": "filters.cluster",
      "tolerance": 1.0,
      "min_points": 60,
      "is3d": true,
      "where": "Planarity > 0.6"
    },
    {
      "type": "writers.las",
      "filename": "tile_4710_5285_features.laz",
      "minor_version": 4,
      "dataformat_id": 6,
      "extra_dims": "all"
    }
  ]
}

• Noise first. Classes 7 and 18 are removed before anything computes neighbourhoods; a single high-noise return above a roof can swing its planarity and pull the normal off vertical.
• count: 2 for HAG. Averaging two ground neighbours smooths the normalization on rough terrain without the cost of a triangulation.
• Height band 2.5 to 120 m. Below 2.5 m sit cars, hedges and garden sheds; above 120 m sit only birds and towers. Restricting the band before segmentation is what stops roofs merging with the ground around them.
• knn: 16 for covariance. At roughly 10 to 20 pts/m², sixteen neighbours span about a metre — small enough to stay on one roof plane, large enough for a stable eigen-decomposition.
• where: "Planarity > 0.6" on the cluster stage. Only planar points are grouped, so tree crowns never join a building segment even when they overhang it. Non-planar points pass through with ClusterID 0.
• extra_dims: "all". Without it the features you just paid to compute are dropped at write time, which is the single most common surprise in this workflow.

# Python Integration

The same pipeline, driven from Python, becomes a function that returns a per-object table — the form in which classification decisions are easiest to write and to test.

"""Compute per-object building candidates from a ground-classified tile."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import numpy as np
import pandas as pd
import pdal

log = logging.getLogger("features")


def feature_pipeline(src: Path, dst: Path, knn: int = 16) -> pdal.Pipeline:
    stages = [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.range", "limits": "HeightAboveGround[2.5:120]"},
        {"type": "filters.covariancefeatures", "knn": knn, "threads": 4,
         "feature_set": "Dimensionality"},
        {"type": "filters.normal", "knn": 12, "always_up": True},
        {"type": "filters.cluster", "tolerance": 1.0, "min_points": 60, "is3d": True,
         "where": "Planarity > 0.6"},
        {"type": "writers.las", "filename": str(dst), "minor_version": 4,
         "dataformat_id": 6, "extra_dims": "all"},
    ]
    return pdal.Pipeline(json.dumps({"pipeline": stages}))


def object_table(points: np.ndarray) -> pd.DataFrame:
    df = pd.DataFrame({name: points[name] for name in (
        "X", "Y", "HeightAboveGround", "Planarity", "Scattering", "NormalZ", "ClusterID")})
    df = df[df["ClusterID"] > 0]
    grouped = df.groupby("ClusterID").agg(
        n=("X", "size"),
        hag_p90=("HeightAboveGround", lambda s: float(np.percentile(s, 90))),
        planarity=("Planarity", "median"),
        scattering=("Scattering", "median"),
        flat_share=("NormalZ", lambda s: float((s > 0.95).mean())),
        x_span=("X", lambda s: float(s.max() - s.min())),
        y_span=("Y", lambda s: float(s.max() - s.min())),
    )
    grouped["bbox_area"] = grouped["x_span"] * grouped["y_span"]
    return grouped


def run(src: Path, dst: Path) -> pd.DataFrame:
    pipeline = feature_pipeline(src, dst)
    n = pipeline.execute()
    log.info("%s: %d points above 2.5 m", src.name, n)
    table = object_table(pipeline.arrays[0])
    table["is_building"] = (
        (table["planarity"] > 0.75)
        & (table["scattering"] < 0.08)
        & (table["bbox_area"] > 30.0)
        & (table["hag_p90"] > 3.0)
    )
    log.info("%d candidate objects, %d kept as buildings",
             len(table), int(table["is_building"].sum()))
    return table


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    print(run(Path("tile_4710_5285.laz"), Path("tile_4710_5285_features.laz")).head())

The design choice worth copying is the split between feature_pipeline, which is pure PDAL, and object_table, which is pure pandas. The first is expensive and changes rarely; the second is cheap and is where you will iterate on thresholds for days. Caching the feature output and re-running only the aggregation turns a ten-minute tuning loop into a two-second one.

# Schema and Data-Flow Considerations

Feature extraction adds more dimensions to a cloud than any other workflow on this site, and each one has to survive every stage and writer between where it is computed and where it is read.

New dimensions are doubles. Linearity, Planarity, HeightAboveGround and the normals are created as 8-byte floats. Four covariance features, three normal components, curvature, HAG and a ClusterID add about 76 bytes per point — more than doubling the in-memory footprint of a PDRF 6 cloud. Plan memory accordingly, or drop what you no longer need before the next blocking stage.

LAS stores them as extra bytes. writers.las with extra_dims: "all" writes each as an extra-bytes VLR entry, and readers.las restores them by name. If you want smaller files, declare types explicitly — "extra_dims": "Planarity=float,HeightAboveGround=float,ClusterID=uint32" — because a 4-byte float is ample precision for a shape descriptor.

ClusterID is per-run, not global. Segment numbers restart in every tile and every execution. Anything that must be joined across tiles — a building that straddles a tile edge — needs a global key built from the tile name plus the local ID, and a buffered tiling step so that the straddling object is seen whole in at least one tile.

Classification is written last. Keep candidate labels in a scratch dimension (CandidateClass via filters.ferry) until the per-object decision is made, then copy into Classification with filters.assign and a WHERE clause. Writing into Classification early destroys the ground/non-ground evidence that later stages read.

# Performance and Scaling

Neighbourhood features dominate the run time of every pipeline in this section, and they scale with knn roughly linearly and with point count slightly worse than linearly. The table gives illustrative timings for a 1 km² tile at about 18 pts/m² (roughly 18 million points after the height band removes ground) on an 8-core worker; use it for proportions rather than absolutes and measure your own hardware.

Three levers matter. Cut the input first — the height band and a class filter routinely remove 60 to 80 percent of points before any neighbourhood is computed. Use the threads option on filters.covariancefeatures, which is one of the few PDAL stages that parallelizes internally. And parallelize across tiles rather than within them, exactly as in parallel tile processing, giving each process a thread count so that processes times threads does not exceed the core count.

# Production Deployment Patterns

Classification pipelines change more often than any other kind on this site, because the thresholds are tuned against each new project’s vegetation, roof styles and sensor. That makes versioning and regression testing unusually important.

• Version thresholds with the pipeline. Keep rule thresholds in the same repository as the pipeline JSON, ideally as parameters in a templated pipeline, so that a delivered tile can always be traced to the exact numbers that classified it.
• Keep a labelled reference tile per project. A few hundred metres of hand-edited classification is enough to compute per-class precision and recall on every change; the method is in evaluating point classification with a confusion matrix.
• Separate the expensive pass. Store the feature-enriched LAZ once, then run the cheap per-object decisions as a second job. The split mirrors splitting a blocking pipeline into two passes and makes reclassification after a threshold change nearly free.
• Run in the same container everywhere. Feature values differ slightly between PDAL releases as neighbourhood code evolves, so the image that tuned the thresholds must be the image that applies them — see pinning PDAL versions with conda-lock.

# Failure Modes and Debugging

Everything is planar. On sparse data (under 4 pts/m²) with a small knn, three or four neighbours always look planar and Planarity saturates near 1 for crowns as well as roofs. Raise knn, or switch to a radius-based neighbourhood, and check the distribution with a histogram before trusting any threshold.

Roofs merge with trees. Euclidean clustering links any two points closer than tolerance, so an overhanging branch bridges a roof and a crown into one segment. Restrict clustering with a where clause on planarity, or use DBSCAN with a min_points high enough that a thin branch is not dense enough to bridge.

ClusterID is all zero. The where clause excluded everything, usually because the feature dimension name was misspelled or the feature stage was placed after the cluster stage. Set "loglevel": 4 in the pipeline options or read pipeline.log after execution; PDAL reports how many points each stage saw.

Extra dimensions vanish. A writer without extra_dims silently drops them. Check the output with pdal info --schema out.laz and assert the names you expect before shipping the tile.

Results change between machines. Covariance features on ties and near-duplicate points are order-sensitive, and different thread counts can visit neighbourhoods in different orders. The differences are tiny per point but can flip a borderline object. Pin threads along with the PDAL version when reproducibility matters.

# Frequently Asked Questions

Do I need machine learning to classify buildings and vegetation?

Not usually to start. Height above ground plus planarity and scattering, aggregated per segment, separates buildings from vegetation well on most airborne data. A trained model earns its place when rules plateau — mixed urban canopy, unusual roof materials, or many classes at once — and it uses the same features as input.

Which neighbourhood size should I use for covariance features?

Aim for a neighbourhood about one metre across on the objects you care about. At 8 pts/m² that is roughly 10 neighbours; at 30 pts/m² it is closer to 25. Plot the feature histogram for known roofs and known crowns and pick the knn at which the two distributions separate most cleanly.

Why compute height above ground before segmentation?

Because segmentation without a height band merges objects with the terrain beneath them. Removing everything under 2 or 3 metres above ground leaves elevated objects isolated, so connected-component segmentation finds each roof or crown as its own group.

Can these stages run in streaming mode on huge tiles?

No. Every neighbourhood stage needs a spatial index over the full input, so they all block streaming. Keep tiles to a size that fits in memory, cut points early with range filters, and use buffered tiles so objects at edges are seen whole.

How do I keep the computed features in the output file?

Set extra_dims to all, or to an explicit list with types, on writers.las. LAS 1.4 stores them as extra bytes with names, and PDAL and laspy both read them back by name.

# Related

• Building Extraction from LiDAR — planar roof segments to footprint polygons with heights
• Individual Tree Segmentation — watershed and point-based methods for stems and crowns
• Power Line Detection — conductors from linearity, catenary fits and clearance checks
• Point Cloud Segmentation — Euclidean and DBSCAN grouping, the step every object workflow depends on
• Machine Learning Point Classification — features, training and honest evaluation
• Water and Bridge Classification — the two classes that make or break a hydrologic DTM