Evaluating Point Classification with a Confusion Matrix

TL;DR: Align predicted and reference classifications point by point, build sklearn.metrics.confusion_matrix with an explicit class list, normalize each row by its reference count, and report per-class precision, recall, F1 and IoU plus Cohen’s kappa — never overall accuracy alone, which a large vegetation class will inflate.

# Context and Motivation

This guide is part of Machine Learning Point Classification for LiDAR. Whether the classes came from a model, from rules, or from a vendor, someone eventually asks how good they are. Overall accuracy is the number people reach for, and it is nearly useless on LiDAR: in a typical suburban tile, ground and high vegetation make up 80 percent or more of points, so a classifier that never finds a single building or wire can still report 90 percent accuracy.

The confusion matrix answers the real question — which classes are confused with which, and how often — and every useful summary metric can be read from it. The same procedure works for model evaluation on held-out tiles and for acceptance testing a vendor delivery against a manually checked reference area.

# Prerequisites and Assumptions

• Two classifications of the same points: predicted and reference. Either two LAS files with identical point order, or one file with the reference stored in an extra dimension.
• Python with scikit-learn, NumPy, pandas and PDAL bindings.
• A reference you trust: a manually edited area, or a held-out tile whose labels passed QA.

# Step-by-Step Implementation

# Step 1 — Align the two classifications

If both files came from the same source tile and neither pipeline reordered points, compare arrays index by index after checking X, Y, Z and GpsTime match. Otherwise join on those coordinates.

# Step 2 — Fix the class list

Pass an explicit labels list to confusion_matrix. Without it, a class absent from one side silently disappears and the matrix changes shape between tiles.

# Step 3 — Row-normalize

Divide each row by its sum. Each diagonal cell is then that class’s recall, and each off-diagonal cell is the share of that reference class lost to another class.

# Step 4 — Compute per-class metrics

Precision (column-wise), recall (row-wise), F1 and intersection-over-union for every class, plus Cohen’s kappa for overall agreement beyond chance.

# Step 5 — Report and archive

Save the raw counts, not just the normalized matrix, so results can be pooled across tiles later.

# Complete Working Example

"""Confusion matrix and per-class metrics for a predicted vs reference classification."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pandas as pd
import pdal
from sklearn.metrics import cohen_kappa_score, confusion_matrix

NAMES = {2: "ground", 3: "low veg", 4: "med veg", 5: "high veg", 6: "building",
         9: "water", 14: "wire", 17: "bridge"}


def read(path: Path) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [str(path)]}))
    p.execute()
    return p.arrays[0]


def aligned(pred_path: Path, ref_path: Path) -> tuple[np.ndarray, np.ndarray]:
    pred, ref = read(pred_path), read(ref_path)
    if len(pred) != len(ref):
        raise ValueError(f"point counts differ: {len(pred)} vs {len(ref)}")
    for dim in ("X", "Y", "Z", "GpsTime"):
        if not np.array_equal(pred[dim], ref[dim]):
            raise ValueError(f"point order differs on {dim}; join on coordinates instead")
    return pred["Classification"], ref["Classification"]


def report(y_pred: np.ndarray, y_ref: np.ndarray, classes: list[int]) -> pd.DataFrame:
    keep = np.isin(y_ref, classes)
    cm = confusion_matrix(y_ref[keep], y_pred[keep], labels=classes)
    tp = np.diag(cm).astype(float)
    ref_n, pred_n = cm.sum(axis=1), cm.sum(axis=0)
    recall = np.divide(tp, ref_n, out=np.zeros_like(tp), where=ref_n > 0)
    precision = np.divide(tp, pred_n, out=np.zeros_like(tp), where=pred_n > 0)
    f1 = np.divide(2 * precision * recall, precision + recall,
                   out=np.zeros_like(tp), where=(precision + recall) > 0)
    iou = np.divide(tp, ref_n + pred_n - tp, out=np.zeros_like(tp), where=(ref_n + pred_n) > 0)
    table = pd.DataFrame({"class": [NAMES.get(c, c) for c in classes], "reference_pts": ref_n,
                          "precision": precision, "recall": recall, "f1": f1, "iou": iou}).round(3)
    kappa = cohen_kappa_score(y_ref[keep], y_pred[keep], labels=classes)
    overall = tp.sum() / cm.sum()
    print(f"overall accuracy {overall:.3f}  kappa {kappa:.3f}  macro F1 {f1.mean():.3f}")
    norm = pd.DataFrame(cm / np.maximum(ref_n[:, None], 1),
                        index=[NAMES.get(c, c) for c in classes],
                        columns=[NAMES.get(c, c) for c in classes]).round(2)
    print(norm.to_string())
    np.save("confusion_counts.npy", cm)
    return table


if __name__ == "__main__":
    y_pred, y_ref = aligned(Path("out/tile_6021_4402.laz"), Path("reference/tile_6021_4402.laz"))
    print(report(y_pred, y_ref, [2, 3, 5, 6, 14]).to_string(index=False))

# Key Parameter Table

# Verification

• Row sums. Each row of the normalized matrix sums to 1 (within rounding) for every class present in the reference.
• Counts add up. The sum of the raw matrix equals the number of reference points in the listed classes.
• Reference quality. Spot-check twenty disagreement points in a viewer. If the reference is wrong in a meaningful share of them, your metrics are measuring the reference, not the classifier.

cm = np.load("confusion_counts.npy")
norm = cm / cm.sum(axis=1, keepdims=True)
assert np.allclose(norm.sum(axis=1)[cm.sum(axis=1) > 0], 1.0)

# Gotchas and Edge Cases

Class 1 in the reference. Unclassified points in the reference have no true class. Exclude them from evaluation, or you will penalize the classifier for labelling points the reference never did.

Pooling across tiles. Average per-tile metrics weight a sparse rural tile as much as a dense urban one. Sum the raw count matrices across tiles and compute metrics once from the total.

Boundary points. Most disagreement concentrates at object edges — eaves, crown fringes — where even two human editors disagree. Consider reporting metrics with a one-point buffer around class boundaries excluded, alongside the full numbers.

Point order changed. A pipeline with a sort or a merge reorders points, and index-wise comparison becomes meaningless. The aligned function checks coordinates for exactly this reason.

# Frequently Asked Questions

Why is overall accuracy misleading for LiDAR classification?

Because class sizes are extremely unequal. Ground and vegetation dominate the point count, so a classifier can miss every wire and most buildings and still score above 90 percent. Per-class recall and precision show those failures immediately.

Should I normalize the confusion matrix by rows or columns?

By rows to read recall — what share of each true class was found and where the rest went. By columns to read precision — what share of each predicted class is correct. Keep the raw counts as well so both can be derived.

What is a good F1 score for building classification?

On clean airborne data with a well-tuned classifier, building F1 above 0.9 is common. Wires and bridges usually score lower because they are small and rare. Compare against a baseline on the same data rather than a universal number.

How large should the reference area be?

Large enough that every class you report has several thousand reference points from several distinct objects, and taken from a location not used for training or tuning.

# Related

• Machine Learning Point Classification — where evaluation fits in the loop
• Training a Random Forest Point Classifier — grouped cross-validation during training
• Computing Geometric Features for Classification — features whose value this measures
• Counting Points per Class with PDAL — quick class histograms before evaluation
• ASPRS Classification Codes — the class definitions being compared