Python LiDAR

Detecting Planar Roofs with Covariance Features

TL;DR: Run filters.covariancefeatures with feature_set: "Dimensionality" and a knn that spans about a metre at your density (roughly 16 at 15 pts/m²), then pick planarity and scattering thresholds from histograms of a few hand-labelled roofs and crowns — typically planarity above 0.7 and scattering below 0.1, applied per segment rather than per point.

# Context and Motivation

This guide is part of Building Extraction from LiDAR. Height above ground puts roofs and tree crowns in the same band; what separates them is shape. A roof is a surface — its points spread in two directions and hardly at all in the third. A crown is a volume — points spread in all three directions. PDAL’s filters.covariancefeatures measures exactly that, from the eigenvalues of the 3×3 covariance matrix of each point’s k nearest neighbours.

With eigenvalues sorted so that λ1 ≥ λ2 ≥ λ3, the dimensionality features are ratios:

  • Linearity = (λ1 − λ2) / λ1 — high for wires and edges.
  • Planarity = (λ2 − λ3) / λ1 — high for roofs, walls and roads.
  • Scattering = λ3 / λ1 — high for vegetation and noise.

They sum to one, which is useful: a point cannot be both strongly planar and strongly scattered. The difficulty is not the maths but the parameters, because the same roof looks planar at one neighbourhood size and noisy at another.

Three features that always sum to one Three stacked horizontal bars, each divided into linearity, planarity and scattering shares. A roof point is mostly planarity at 0.86, with linearity 0.10 and scattering 0.04. A crown point is mostly scattering at 0.52, with planarity 0.30 and linearity 0.18. A wire point is mostly linearity at 0.91. roof crown wire Linearity Planarity Scattering illustrative values at knn = 16, about 15 pts/m²

# Prerequisites and Assumptions

  • PDAL 2.4 or newer; confirm the options with pdal --options filters.covariancefeatures.
  • Ground classified, and HeightAboveGround available or computable with filters.hag_nn.
  • A handful of hand-labelled examples: ten roofs and ten crowns digitized as polygons are enough to set thresholds.
  • Python with NumPy and GeoPandas for the histogram step.

# Step-by-Step Implementation

# Step 1 — Estimate the right neighbourhood size

The neighbourhood should cover about one metre of roof: large enough that three or four returns are not accidentally coplanar, small enough that it rarely straddles a ridge or a roof edge. At density ρ points per m², a one-metre-radius disc holds about πρ points, so start with knn ≈ 3ρ and round.

Density (pts/m²) Suggested knn Neighbourhood radius
4 12 ~1.0 m
8 16 ~0.8 m
15 20 ~0.65 m
30 32 ~0.6 m

# Step 2 — Compute features on elevated points only

Restrict the stage with a where clause so ground and low objects do not consume time.

json
{
  "type": "filters.covariancefeatures",
  "knn": 20,
  "threads": 4,
  "feature_set": "Dimensionality",
  "where": "HeightAboveGround > 2.5 && Classification != 2"
}

# Step 3 — Sample features inside labelled polygons

Burn the labelled roof and crown polygons into a Label dimension with filters.overlay, then read Planarity and Scattering for each group.

# Step 4 — Choose thresholds from the histograms

Plot the two distributions for each feature and set the threshold where they cross. Planarity usually separates cleanly at 0.6 to 0.75; scattering at 0.08 to 0.12.

# Step 5 — Apply thresholds per segment, not per point

Segment elevated points with filters.cluster and test each segment’s median planarity and scattering. Medians absorb the noisy edge points that would otherwise speckle the classification.

# Complete Working Example

python
"""Pick roof/crown thresholds for covariance features from labelled polygons."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pdal


def labelled_features(tile: Path, labels_gpkg: Path, knn: int) -> np.ndarray:
    stages = [
        {"type": "readers.las", "filename": str(tile)},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.range", "limits": "HeightAboveGround[2.5:80]"},
        {"type": "filters.ferry", "dimensions": "=>Label"},
        {"type": "filters.overlay", "dimension": "Label", "datasource": str(labels_gpkg),
         "layer": "samples", "column": "label"},          # 1 = roof, 2 = crown
        {"type": "filters.range", "limits": "Label[1:2]"},
        {"type": "filters.covariancefeatures", "knn": knn, "threads": 4,
         "feature_set": "Dimensionality"},
    ]
    p = pdal.Pipeline(json.dumps({"pipeline": stages}))
    p.execute()
    return p.arrays[0]


def best_threshold(roof: np.ndarray, crown: np.ndarray, roof_high: bool) -> tuple[float, float]:
    """Threshold maximizing balanced accuracy between the two samples."""
    grid = np.linspace(0.0, 1.0, 201)
    best = (0.0, 0.0)
    for t in grid:
        tpr = (roof > t).mean() if roof_high else (roof < t).mean()
        tnr = (crown <= t).mean() if roof_high else (crown >= t).mean()
        score = 0.5 * (tpr + tnr)
        if score > best[1]:
            best = (float(t), float(score))
    return best


if __name__ == "__main__":
    for knn in (12, 16, 20, 28):
        a = labelled_features(Path("tile_5840_2710.laz"), Path("samples.gpkg"), knn)
        roof, crown = a[a["Label"] == 1], a[a["Label"] == 2]
        tp, sp = best_threshold(roof["Planarity"], crown["Planarity"], roof_high=True)
        ts, ss = best_threshold(roof["Scattering"], crown["Scattering"], roof_high=False)
        print(f"knn={knn:>2}  planarity>{tp:.2f} (bal.acc {sp:.3f})  "
              f"scattering<{ts:.2f} (bal.acc {ss:.3f})")

Typical output shows the separation improving up to a point and then degrading as neighbourhoods start to cross roof edges:

text
knn=12  planarity>0.58 (bal.acc 0.884)  scattering<0.14 (bal.acc 0.902)
knn=16  planarity>0.66 (bal.acc 0.917)  scattering<0.11 (bal.acc 0.931)
knn=20  planarity>0.70 (bal.acc 0.926)  scattering<0.10 (bal.acc 0.938)
knn=28  planarity>0.71 (bal.acc 0.912)  scattering<0.09 (bal.acc 0.927)
Reading the threshold off the histograms Two overlapping histograms of planarity. Crown points peak around 0.3 and tail off by 0.7. Roof points peak around 0.85 and tail down to about 0.5. A vertical dashed line at 0.70 marks the threshold where the two distributions cross, which maximizes balanced accuracy. threshold 0.70 crowns roofs 0 0.5 1.0 Planarity

# Key Parameter Table

Option Type Default Guidance
knn int 10 About 3 × density in pts/m²; confirm with the sweep above
threads int 1 4 roughly halves run time; little gain beyond 8
feature_set string Dimensionality all adds omnivariance, eigenentropy and more at extra cost
mode string SQRT Whether eigenvalues are square-rooted before ratios; keep one setting for training and use
min_k int 3 Minimum neighbours for a valid result
radius float unset Radius neighbourhood instead of kNN; steadier across density changes

# Verification

  • Distributions. After running on a new tile, the planarity histogram of elevated points should be bimodal, with a clear roof peak above 0.7. A unimodal hump means knn is badly mismatched to density.
  • Spot check. Colour the output by Planarity in a viewer. Roof interiors should be uniformly high, ridges and edges lower, crowns mottled.
  • Consistency across tiles. Compute the median planarity of segments classified as roofs on each tile; a tile that differs by more than 0.1 from the others was flown at different density or processed with different options.

# Gotchas and Edge Cases

Ridges and hips look non-planar. A neighbourhood that straddles two roof planes has two dominant directions and a larger third eigenvalue. That is why thresholds are applied to segment medians: ridges are a small fraction of each roof.

Flat roofs with gravel or plant. Rooftop HVAC units, solar panels on racks and parapets reduce planarity locally. They rarely change a segment median enough to matter, but they will fail a per-point test.

Green roofs and very steep roofs. Vegetated roofs scatter like crowns, and near-vertical mansard faces return few points. Both are legitimate exceptions; label a few in the training polygons if the project has many.

Where the neighbourhood falls matters A gable roof profile with three circled neighbourhoods. One circle sits in the middle of a roof plane and reads as planar. One circle straddles the ridge and mixes two planes, reading as partly scattered. A third circle inside a nearby tree crown reads as scattered. Labels give example planarity values of 0.88, 0.52 and 0.27. plane: 0.88 ridge: 0.52 crown: 0.27

Density changes within a tile. Flightline overlap doubles density, so a fixed knn covers a smaller area there. If your data has strong overlap stripes, use the radius option instead, which keeps the neighbourhood physically constant. Normalizing density first with decimation is another option.

# Frequently Asked Questions

What is the difference between planarity and filters.approximatecoplanar?

Planarity is a continuous value from 0 to 1 that you threshold yourself. filters.approximatecoplanar applies fixed eigenvalue-ratio tests and writes a boolean Coplanar flag. The continuous value is more useful because you can aggregate it per segment and tune the threshold to your data.

Why do my roof points have low planarity?

Most often knn is too small for the density, so each neighbourhood holds only a few points that happen to be noisy, or so large that it crosses roof edges. Run the knn sweep on labelled samples and pick the value with the best separation.

Should I use knn or radius neighbourhoods?

kNN adapts to density, which is good when density is uniform. A radius keeps the physical scale constant, which is better when overlap or scan pattern makes density vary strongly across the tile.

Are the thresholds transferable between projects?

Approximately, when density and knn are similar, but not exactly. Re-derive them from a few labelled samples on each project; it takes minutes.

Related in Building Extraction from LiDAR Point Clouds

↑ Building Extraction from LiDAR Point Clouds overview --- title: "Building Extraction from LiDAR Point Clouds" description: "Classify roof points to ASPRS class 6 and turn them into footprint... Estimating Roof Pitch and Aspect from Normals --- title: "Estimating Roof Pitch and Aspect from Normals" description: "Compute per-plane roof pitch and aspect from PDAL filters.normal:... Extracting Building Footprints from LiDAR --- title: "Extracting Building Footprints from LiDAR" description: "Turn class-6 roof points into clean, rectilinear footprint polygons:...