Training a Random Forest Point Classifier

TL;DR: Sample up to a fixed number of points per class from every labelled tile into one table with a tile column, evaluate with GroupKFold so each fold holds out whole tiles, tune only max_depth, min_samples_leaf and max_features, then refit on all tiles and save a joblib bundle containing the model, the ordered feature list, the class list and the cross-validated scores.

# Context and Motivation

This guide is part of Machine Learning Point Classification for LiDAR. A random forest is the pragmatic default for point classification: it trains quickly on hundreds of thousands of rows, handles features on different scales without normalization, copes with correlated features, and exposes feature importances you can reason about. What decides whether it works in production is not the algorithm but three pieces of discipline around it — how the training set is sampled, how the model is evaluated, and how it is packaged so that prediction uses exactly the features training used.

The procedure below assumes features have already been computed and cached per tile, as in computing geometric features for classification.

# Prerequisites and Assumptions

• Feature-enriched LAZ files for labelled tiles, each with trustworthy Classification.
• Python with scikit-learn 1.3+, pandas, NumPy, joblib and PDAL bindings.
• At least eight to ten labelled tiles from different parts of the project, so grouped cross-validation has something to hold out.

# Step-by-Step Implementation

# Step 1 — Sample a balanced table from every tile

Read each tile, keep above-ground classes, and sample up to N points per class per tile. Sampling per tile keeps every landscape represented; capping per class keeps rare classes visible.

# Step 2 — Set up grouped cross-validation

GroupKFold(n_splits=5) with the tile column as groups guarantees that no tile contributes to both training and testing within a fold.

# Step 3 — Search a small parameter grid

Forests are forgiving. Tuning max_depth, min_samples_leaf and max_features over a dozen combinations captures nearly all the available gain; more search mostly fits noise in the validation scores.

# Step 4 — Score with a class-balanced metric

Use macro-averaged F1 as the selection metric, so a model that ignores wires cannot win by being good at vegetation.

# Step 5 — Refit on all tiles and save a bundle

Refit the best configuration on every labelled tile and save the model with its feature order, class list, parameters and cross-validated scores.

# Complete Working Example

"""Balanced sampling, grouped CV and a versioned model bundle."""
from __future__ import annotations

import json
from datetime import date
from pathlib import Path

import joblib
import numpy as np
import pandas as pd
import pdal
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV, GroupKFold

FEATURES = ["HeightAboveGround", "ReturnRatio", "NormalZ", "Curvature",
            "Linearity", "Planarity", "Scattering", "Verticality",
            "Linearity_s", "Planarity_s", "Scattering_s", "Verticality_s"]
CLASSES = [3, 4, 5, 6, 9, 14, 17]


def sample_tile(path: Path, per_class: int, rng: np.random.Generator) -> pd.DataFrame:
    p = pdal.Pipeline(json.dumps({"pipeline": [str(path)]}))
    p.execute()
    a = p.arrays[0]
    df = pd.DataFrame({k: a[k] for k in FEATURES})
    df["y"] = a["Classification"]
    df = df[df.y.isin(CLASSES)]
    parts = [g.sample(min(len(g), per_class), random_state=int(rng.integers(1e9)))
             for _, g in df.groupby("y")]
    out = pd.concat(parts)
    out["tile"] = path.stem
    return out


def build_table(tiles: list[Path], per_class: int = 15_000) -> pd.DataFrame:
    rng = np.random.default_rng(42)
    table = pd.concat([sample_tile(t, per_class, rng) for t in tiles], ignore_index=True)
    print(table.groupby("y").size().rename("rows"))
    return table


def train(table: pd.DataFrame, out: Path) -> dict:
    X, y, groups = table[FEATURES], table["y"], table["tile"]
    grid = GridSearchCV(
        RandomForestClassifier(n_estimators=200, class_weight="balanced_subsample",
                               n_jobs=-1, random_state=0),
        param_grid={"max_depth": [14, 18, 24],
                    "min_samples_leaf": [2, 5, 10],
                    "max_features": ["sqrt", 0.5]},
        scoring="f1_macro", cv=GroupKFold(n_splits=5), n_jobs=1, refit=True, verbose=1,
    )
    grid.fit(X, y, groups=groups)
    bundle = {
        "model": grid.best_estimator_,
        "features": FEATURES,
        "classes": CLASSES,
        "params": grid.best_params_,
        "cv_f1_macro": float(grid.best_score_),
        "tiles": sorted(groups.unique()),
        "trained": date.today().isoformat(),
    }
    joblib.dump(bundle, out, compress=3)
    print(f"best {grid.best_params_}  grouped-CV macro F1 {grid.best_score_:.3f}")
    return bundle


if __name__ == "__main__":
    tiles = sorted(Path("features").glob("*.laz"))
    train(build_table(tiles), Path("models/rf_2026-09.joblib"))

refit=True means best_estimator_ is already fitted on all rows after the search, so no separate final fit is needed.

# Key Parameter Table

# Verification

• Grouped score versus random-split score. Run one random KFold for comparison. A gap of five or more points of macro F1 shows how much a naive evaluation would have overstated performance.
• Per-fold spread. Look at grid.cv_results_ for the best configuration: a fold far below the others is a tile unlike the rest, and a sign you need more labelled data from that kind of landscape.
• Bundle round trip. Load the saved bundle and predict on a few rows to make sure the feature order is honoured.

bundle = joblib.load("models/rf_2026-09.joblib")
row = table[bundle["features"]].head(5)
assert list(row.columns) == bundle["features"]
print(bundle["model"].predict(row), bundle["cv_f1_macro"], bundle["params"])

# Gotchas and Edge Cases

Class present in only one tile. Grouped CV will hold that tile out in one fold, leaving the class absent from training there and producing a very low score for it. That is honest — it says the model cannot generalize that class — but it means you need labels for the class in more tiles.

Label noise. Vendor classifications used as labels contain errors, especially low versus medium vegetation. A larger min_samples_leaf helps the forest average over them; do not chase perfect training accuracy.

Model size. Deep forests on many rows produce joblib files of hundreds of megabytes. compress=3 helps; so does capping depth. Check the size before shipping the model into a container image.

Data leakage through overlap. Two tiles that share a flightline overlap strip contain near-identical points. Group by a spatial block larger than a tile if your tiling has overlap, or remove overlap points before sampling.

# Frequently Asked Questions

Why use GroupKFold instead of ordinary cross-validation?

Neighbouring LiDAR points are nearly identical, so an ordinary split puts near-duplicates on both sides and measures memory rather than generalization. Grouping by tile ensures each test fold contains only areas the model has never seen.

How many trees should the forest have?

Two hundred is a sensible default. Beyond that, predictions barely change while prediction time and model size grow linearly. Tune depth and leaf size instead.

Should I scale or normalize features for a random forest?

No. Trees split on thresholds, so monotonic rescaling does not change them. Consistency matters more than scale: the prediction pipeline must compute features exactly as training did.

What should go into the saved model bundle?

The fitted model, the ordered feature list, the class list, the chosen parameters, the cross-validated score, the training tiles and a date. Those are what you need to reproduce a prediction or explain a result months later.

# Related

• Machine Learning Point Classification — the end-to-end loop
• Computing Geometric Features for Classification — the inputs to this model
• Evaluating Point Classification with a Confusion Matrix — reading the results per class
• Extracting Objects from Segment Labels — voting predictions per object
• Pinning PDAL Versions with conda-lock — keeping feature computation reproducible