Python LiDAR

Machine Learning Point Classification for LiDAR

Rule-based classification — height bands, planarity thresholds, per-segment tests — gets a project a long way, but it plateaus. Each new rule fixes one confusion and introduces another, the thresholds drift between flights, and by the time you are separating low vegetation from ground clutter from car roofs from rooftop plant, the rule set has become a program nobody wants to maintain. A supervised classifier replaces those hand-set thresholds with ones learned from labelled examples. The features are the same ones the rules used — height above ground, planarity, linearity, scattering, intensity, return structure — and the model’s job is to combine them better than a person can. This topic covers the practical loop in the classification and feature extraction section: features from PDAL, a random forest in scikit-learn, evaluation that does not lie, and predictions written back as ASPRS codes.

Train once, predict everywhere Two rows. The upper training row runs from labelled tiles through PDAL feature computation into a feature table, then to model training and a saved model file. The lower inference row runs from new tiles through the same PDAL feature pipeline, then prediction with the saved model, then writing classes to LAS. An arrow shows the saved model feeding prediction, and a note stresses that both rows must use the identical feature pipeline. training labelled tiles feature pipeline feature table fit model model.joblib inference new tiles same pipeline predict write classes if the two feature pipelines differ by one option, the model silently degrades — version them together

# Prerequisites

  • PDAL 2.5+ with filters.covariancefeatures, filters.hag_nn and filters.normal, plus the Python bindings.
  • Python 3.10+ with NumPy, pandas, scikit-learn 1.3+ and joblib.
  • Labelled training data. Tiles whose classification you trust — manually edited, or a vendor delivery that passed QA — covering every class you want to predict. A few square kilometres from varied parts of the project is usually enough for a random forest.
  • Ground already classified. Train the model on above-ground classes and leave ground to a dedicated filter such as SMRF; ground filters are better at ground than any per-point model.
  • Consistent sensor and density between training and target tiles, or at least training data that spans the variation. A model trained on 8 pts/m² data will not transfer cleanly to 40 pts/m².

# Core Workflow Architecture

  1. Compute features. Run one PDAL feature pipeline — HAG, covariance features at one or more scales, normals, intensity, return number ratio — on training and target tiles alike.
  2. Assemble a training table. Stack features and labels from training tiles into a data frame, recording the tile each row came from.
  3. Balance the classes. Down-sample dominant classes (usually high vegetation) so the model does not ignore rare ones such as wires or bridges.
  4. Split spatially. Hold out whole tiles or blocks for testing, never random points, because neighbouring points share features and would leak.
  5. Train and evaluate. Fit a random forest, compute a confusion matrix and per-class precision and recall on held-out tiles.
  6. Predict and write. Apply the model to new tiles, optionally smooth labels with a neighbour vote, and write ASPRS codes back to the LAS.

# Full Implementation

python
"""Train a random-forest point classifier on PDAL features and apply it to new tiles."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import joblib
import numpy as np
import pandas as pd
import pdal
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix

log = logging.getLogger("rf")

FEATURES = ["HeightAboveGround", "Intensity", "ReturnRatio",
            "Linearity", "Planarity", "Scattering", "Verticality", "NormalZ"]
CLASSES = [1, 3, 4, 5, 6, 9, 14, 17]      # above-ground classes the model predicts


def feature_stages(src: Path, knn: int = 16) -> list[dict]:
    above = "Classification != 2 && Classification != 7 && Classification != 18"
    return [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.covariancefeatures", "knn": knn, "threads": 4,
         "feature_set": "Dimensionality", "where": above},
        {"type": "filters.normal", "knn": 12, "always_up": True, "where": above},
    ]


def features(src: Path) -> tuple[np.ndarray, pd.DataFrame]:
    p = pdal.Pipeline(json.dumps({"pipeline": feature_stages(src)}))
    p.execute()
    arr = p.arrays[0]
    df = pd.DataFrame({k: arr[k] for k in FEATURES if k in arr.dtype.names})
    df["ReturnRatio"] = arr["ReturnNumber"] / np.maximum(arr["NumberOfReturns"], 1)
    df["Classification"] = arr["Classification"]
    return arr, df[FEATURES + ["Classification"]]


def training_table(tiles: list[Path], per_class: int = 60_000) -> pd.DataFrame:
    parts = []
    for tile in tiles:
        _, df = features(tile)
        df = df[df.Classification.isin(CLASSES)].copy()
        df["tile"] = tile.stem
        parts.append(df)
    table = pd.concat(parts, ignore_index=True)
    balanced = (table.groupby("Classification", group_keys=False)
                .apply(lambda g: g.sample(min(len(g), per_class), random_state=0)))
    log.info("training rows: %d (from %d)", len(balanced), len(table))
    return balanced


def train(table: pd.DataFrame, test_tiles: set[str], out: Path) -> RandomForestClassifier:
    test = table.tile.isin(test_tiles)
    X_tr, y_tr = table.loc[~test, FEATURES], table.loc[~test, "Classification"]
    X_te, y_te = table.loc[test, FEATURES], table.loc[test, "Classification"]
    model = RandomForestClassifier(n_estimators=200, max_depth=18, min_samples_leaf=5,
                                   n_jobs=-1, class_weight="balanced_subsample",
                                   random_state=0)
    model.fit(X_tr, y_tr)
    pred = model.predict(X_te)
    log.info("\n%s", classification_report(y_te, pred, digits=3, zero_division=0))
    log.info("confusion matrix (rows = truth):\n%s",
             confusion_matrix(y_te, pred, labels=CLASSES))
    joblib.dump({"model": model, "features": FEATURES, "classes": CLASSES}, out)
    return model


def predict(src: Path, dst: Path, model_path: Path) -> None:
    bundle = joblib.load(model_path)
    arr, df = features(src)
    above = ~np.isin(arr["Classification"], [2, 7, 18])
    arr["Classification"][above] = bundle["model"].predict(df.loc[above, bundle["features"]])
    pdal.Writer.las(filename=str(dst), minor_version=4, dataformat_id=6,
                    forward="all").pipeline(arr).execute()
    log.info("%s: %d points classified", dst.name, int(above.sum()))


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    tiles = sorted(Path("labelled").glob("*.laz"))
    table = training_table(tiles)
    train(table, test_tiles={tiles[0].stem, tiles[5].stem}, out=Path("rf_v1.joblib"))
    predict(Path("new/tile_6021_4402.laz"), Path("out/tile_6021_4402.laz"), Path("rf_v1.joblib"))

# Code Breakdown

One feature function for training and prediction. features is called by both paths, so there is no way for the options to drift apart. Store the model with its feature list, as joblib.dump does here, and a mismatch becomes a loud error instead of a quiet accuracy loss.

ReturnRatio as a derived feature. Return number alone is ambiguous — a second return means different things in a two-return and a five-return pulse. The ratio is close to 1 for last returns (ground, roofs) and small for early returns inside canopy.

Above-ground only. Ground, low noise and high noise are excluded from both training and prediction. The model never has to learn the easy separation from the terrain, and ground points are never overwritten by a model that is worse at them than SMRF.

Balancing by down-sampling. A suburban tile can hold a hundred times more high-vegetation points than wire points. Capping each class at 60,000 rows and using class_weight="balanced_subsample" keeps rare classes visible to the trees. Down-sampling also keeps training fast.

Holding out whole tiles. Adjacent points have nearly identical features. A random point split puts near-duplicates on both sides and reports inflated accuracy — often above 0.98 — that collapses on new data. Tile-level hold-out is the minimum honest test; the confusion matrix guide shows the gap on real numbers.

Forest settings. Two hundred trees, depth 18 and a minimum leaf of five points are sensible defaults for a few hundred thousand rows; deeper trees memorize neighbourhoods and shallower ones underfit rare classes. n_jobs=-1 uses every core for both training and prediction.

Why the split must be spatial Left: a grid of tiles where training and test points are interleaved at random within every tile, reporting an optimistic overall accuracy of 0.98. Right: the same grid where two whole tiles are held out for testing, reporting a realistic accuracy of 0.89. A note explains that neighbouring points are near-duplicates, so a random split tests memory rather than generalization. random point split whole-tile hold-out reported accuracy 0.98 reported accuracy 0.89 test points sit next to training points shaded tiles never seen in training illustrative figures; the gap between the two is the leak

# Which Features Earn Their Place

Feature engineering matters more than model choice for point classification, and the useful features fall into four families. Knowing what each contributes makes it easier to decide what to add when a particular confusion persists.

Height. HeightAboveGround is almost always the most important single feature: it separates low from medium from high vegetation outright, and it puts roofs and wires in a band where ground clutter never appears. Adding local height statistics — the range or standard deviation of HAG among the k neighbours — helps separate roof edges from the crowns next to them.

Shape. Linearity, planarity, scattering and verticality describe the neighbourhood’s geometry and are what separate a roof from a crown at the same height, or a wire from a branch. Normals add orientation, which distinguishes flat roofs from walls and sloping roofs from hedges.

Return structure. The ratio of return number to number of returns, and whether a point is a single return, capture how the pulse interacted with the target. Solid surfaces give single or last returns; vegetation gives multiple.

Radiometry. Intensity helps with water, road markings and some roof materials, but it is the least transferable family because it depends on range, incidence angle, sensor and calibration. Use it last, and normalized.

What the forest actually used A horizontal bar chart of feature importances from a random forest trained on eight features. Height above ground is highest at 0.31, followed by scattering at 0.17, planarity at 0.14, return ratio at 0.11, verticality at 0.09, NormalZ at 0.08, linearity at 0.06 and intensity lowest at 0.04. Bars are coloured by feature family. HeightAboveGround0.31 Scattering0.17 Planarity0.14 ReturnRatio0.11 Verticality0.09 NormalZ0.08 Linearity0.06 Intensity0.04 illustrative suburban model — height, shape, return structure, radiometry

Low importance for linearity is expected in a suburban model and does not mean it is useless: wires are rare, so the feature that detects them contributes little to the average but everything to that class’s recall. Judge features by per-class results, not by importance alone.

# Parameter Reference Table

Setting Type Default here Typical range Effect
knn (covariance) int 16 10–40 Neighbourhood scale; adding a second scale often helps more than tuning one
per_class int 60,000 10k–200k Training rows per class after balancing
n_estimators int 200 100–500 More trees give smoother probabilities, slower prediction
max_depth int 18 12–30 Deeper fits more detail and memorizes more noise
min_samples_leaf int 5 1–20 Larger regularizes; raise when labels are noisy
class_weight str balanced_subsample or None Counteracts class imbalance inside each tree
held-out tiles set 2 of 10 15–25 % Area never seen in training, used only for evaluation

# Validation and Integrity Checks

Per-class, not overall. Overall accuracy is dominated by the largest class. Report precision and recall for every class, and pay most attention to the classes the client cares about — buildings and wires usually matter more than the split between medium and high vegetation.

Confusion structure. Some confusions are cheap (medium versus high vegetation) and some expensive (vegetation labelled as building). Read the off-diagonal cells, not just the diagonal.

Class histogram sanity. After prediction on a new tile, compare the class proportions with those of similar training tiles. A tile that comes out 30 percent building in a rural area has a problem, usually a density or intensity mismatch.

Feature importance. model.feature_importances_ should put height above ground and one or two shape features at the top. If intensity dominates, the model may be learning sensor calibration rather than geometry and will not transfer to another flight — see normalizing intensity across flightlines.

# Performance Tuning

Feature computation dominates both training and prediction; the forest itself is quick. On a 1 km² tile at 20 pts/m², expect roughly two to four minutes of PDAL work and well under a minute of prediction with n_jobs=-1.

  • Cache features. Write the feature-enriched tile once with extra_dims and re-read it for every training experiment. Retraining then takes seconds.
  • Predict in chunks. For very large tiles, predict on slices of a few million rows to keep memory flat; random forests are embarrassingly parallel over rows.
  • Multi-scale features, selectively. A second covariance scale (for example knn: 40) often improves wire and building edges, but doubles the dominant cost. Add it only if the per-class report shows a specific confusion it would address.
  • Smaller models for deployment. Capping max_depth or reducing n_estimators shrinks the saved model and speeds prediction with little accuracy cost; measure on the held-out tiles before and after.

# Common Errors and Troubleshooting

ValueError: X has 7 features, but RandomForestClassifier is expecting 8. A feature is missing on the prediction tile — usually because a where clause excluded every point from the stage that creates it, so the dimension was never added. Guard with an explicit check of arr.dtype.names and fail with the missing name.

Near-perfect test scores. The split is leaking. Check that no tile appears in both sets, and that overlapping flightlines do not put the same ground in two differently named tiles.

Rare classes never predicted. Balancing was skipped or the class has too few training examples. Wires and bridges need several thousand labelled points each, drawn from several tiles.

Salt-and-pepper labels. Per-point prediction produces isolated misclassifications. Smooth with filters.neighborclassifier (majority vote over k neighbours) or decide per segment by majority, as in extracting objects from segment labels.

Model degrades on a new flight. Density, sensor or season changed. Add a few labelled tiles from the new flight to training rather than tuning the model; representation matters more than hyperparameters.

# Frequently Asked Questions

Why a random forest rather than a deep network?

Random forests train in minutes on a laptop, need little tuning, handle mixed features well and are easy to inspect. Point-based deep networks can be more accurate on complex scenes but need far more labelled data, GPUs and engineering. Start with a forest and move on only if its per-class report shows a ceiling you must break.

How much labelled data do I need?

For a random forest on geometric features, a few square kilometres of trustworthy classification spread across the project’s landscapes is usually enough, with at least several thousand points for each rare class. Diversity of scenes matters more than volume.

Should the model also classify ground?

Usually not. Dedicated ground filters handle terrain better than per-point models and are easier to tune. Classify ground first, then let the model separate the above-ground classes.

Can I reuse a model across projects?

Only when the new project resembles the training data in density, sensor, season and landscape. A model trained on leaf-off suburban data will misread leaf-on forest. Treat every model as project-specific until a held-out tile from the new project shows per-class results close to the original evaluation, and keep the evaluation numbers alongside the model file.

How do I stop the model learning sensor quirks?

Prefer geometric features over intensity, normalize intensity if you do use it, and train on tiles from several flights. Feature importances that put intensity first are a warning sign.

Hands-on guides

Computing Geometric Features for Classification --- title: "Computing Geometric Features for Classification" description: "Build a multi-scale per-point feature set in PDAL for... Evaluating Point Classification with a Confusion Matrix --- title: "Evaluating Point Classification with a Confusion Matrix" description: "Measure LiDAR classification quality honestly: a... Training a Random Forest Point Classifier --- title: "Training a Random Forest Point Classifier" description: "Train a scikit-learn random forest on LiDAR point features: sampling a...

More in LiDAR Classification and Feature Extraction in Python

Building Extraction from LiDAR Point Clouds --- title: "Building Extraction from LiDAR Point Clouds" description: "Classify roof points to ASPRS class 6 and turn... Individual Tree Segmentation from LiDAR --- title: "Individual Tree Segmentation from LiDAR" description: "Split a forest canopy into individual trees:... Power Line Detection in LiDAR Point Clouds --- title: "Power Line Detection in LiDAR Point Clouds" description: "Find conductor points with linearity and height,... Water and Bridge Classification in LiDAR --- title: "Water and Bridge Classification in LiDAR" description: "Classify water to ASPRS class 9 and bridge decks to... Point Cloud Segmentation with PDAL --- title: "Point Cloud Segmentation with PDAL" description: "Group points into objects with filters.cluster and...