⚙️ How-To 3D & buildings LiDAR & point clouds 3D & Buildings · step 9 of 20 Updated

How to classify roof shapes from LiDAR in Python

Problem statement

An LoD2 model tells you the roof shape. If you do not have one, you have a point cloud and a footprint, and getting from there to "this is a gable, that is a hip, that one is flat" is a segmentation problem with a classification on top.

The useful news is that the first and most valuable distinction is nearly free. On a real LoD2 dataset, 63.5% of buildings were flat-roofed, and for those every downstream calculation — volume, solar, shadow — is the simple case. Separating flat from pitched needs one number per building; identifying gable versus hip versus mansard needs plane segmentation.

This guide does both, with the checks that stop a chimney or a tree from deciding the answer.

Quick answer

Height range within the footprint separates flat from pitched:

import numpy as np, geopandas as gpd

def roof_relief(points_in_footprint, low=10, high=90):
    z = points_in_footprint
    return float(np.percentile(z, high) - np.percentile(z, low))

b["relief_m"] = [roof_relief(z) for z in per_building_z]
b["roof_class"] = np.where(b.relief_m < 0.8, "flat", "pitched")

Percentiles rather than min and max, because a chimney, an aerial, a parapet or an overhanging branch sets the extremes. On the reference LoD2 data, buildings with ridge height equal to eaves height — genuinely flat — were 63.5% of the stock, which is the figure this threshold should roughly reproduce.

Four scenes showing flat, gable, hip and mansard roof forms with their plane counts.
Plane count and orientation are what separate the forms; height alone only finds flat.

Step-by-step solution

1. Clip the point cloud to the footprint, inset

Buffer the footprint inwards by half a metre to a metre before selecting points. Otherwise the wall returns, the guttering and the ground immediately outside contaminate every statistic.

2. Keep the right returns

First returns only, and filter out points classified as vegetation. An overhanging tree adds a second surface above the roof and turns a gable into a nonsense.

3. Require enough points

A footprint with twenty returns cannot be segmented into planes. Set a minimum point count and a minimum density, and report how many buildings fall below it rather than classifying them anyway.

4. Separate flat from pitched by relief

The 10th-to-90th percentile height range within the inset footprint. Below about 0.8 m is flat — the residual is parapets, plant and noise. This one number handles the majority of the stock.

5. Segment planes with RANSAC for the rest

Fit a plane, remove its inliers, repeat. Two roughly symmetric planes with opposite azimuths is a gable; four planes meeting at a point or a short ridge is a hip; two planes per side with different tilts is a mansard; one tilted plane is a shed.

6. Classify from the plane statistics

Plane count, the tilt of each, the azimuth spread, and whether the azimuths are opposite. That is enough for the common forms; anything more exotic should be labelled "complex" rather than forced into a category.

7. Validate against something

If any LoD2 coverage exists nearby, classify the same buildings and compare. Roof type codes in such models follow a convention — in the reference tile, code 1000 (flat) covered 1,238 of 1,990 buildings, 1030 411, 1130 138 and 1010 110 — and mapping your classes onto them makes the comparison concrete.

Flow from point cloud through footprint clipping, return filtering, relief test and plane segmentation to a roof class.
Most buildings exit at the relief test; only the pitched ones need segmenting.

Code examples

Example 1 — clip, filter and compute relief 

import numpy as np, geopandas as gpd, laspy

def roof_points(las_path, buildings, inset=0.75, min_points=30):
    las = laspy.read(las_path)
    keep = (las.return_number == 1) & (~np.isin(las.classification, [3, 4, 5]))  # no vegetation
    pts = gpd.GeoDataFrame({"z": las.z[keep]},
                           geometry=gpd.points_from_xy(las.x[keep], las.y[keep]),
                           crs=buildings.crs)

    inner = buildings.copy()
    inner["geometry"] = inner.buffer(-inset)
    inner = inner[~inner.is_empty]

    j = gpd.sjoin(pts, inner[["id", "geometry"]], predicate="within")
    stats = j.groupby("id")["z"].agg(
        n="size",
        p10=lambda s: np.percentile(s, 10),
        p50="median",
        p90=lambda s: np.percentile(s, 90))
    stats["relief_m"] = stats.p90 - stats.p10
    return stats[stats.n >= min_points], j

stats, points = roof_points("tile.laz", b)
print(f"{len(stats):,} of {len(b):,} buildings had {30}+ roof returns")
print(f"flat (relief < 0.8 m): {(stats.relief_m < 0.8).mean():.1%}")

Example 2 — plane segmentation with RANSAC 

import numpy as np
from sklearn.linear_model import RANSACRegressor, LinearRegression

def segment_planes(xyz, max_planes=4, residual_threshold=0.15, min_inliers=25):
    """Fit planes z = ax + by + c, removing inliers each round."""
    remaining = xyz.copy()
    planes = []
    for _ in range(max_planes):
        if len(remaining) < min_inliers:
            break
        X, z = remaining[:, :2], remaining[:, 2]
        model = RANSACRegressor(LinearRegression(), residual_threshold=residual_threshold,
                                max_trials=200, random_state=0).fit(X, z)
        inliers = model.inlier_mask_
        if inliers.sum() < min_inliers:
            break
        a, bcoef = model.estimator_.coef_
        normal = np.array([-a, -bcoef, 1.0])
        normal /= np.linalg.norm(normal)
        planes.append({
            "n_points": int(inliers.sum()),
            "tilt_deg": float(np.degrees(np.arccos(abs(normal[2])))),
            "azimuth_deg": float((np.degrees(np.arctan2(normal[0], normal[1])) + 360) % 360),
        })
        remaining = remaining[~inliers]
    return planes

residual_threshold is the sensitive parameter: too small and a single roof plane splits into strips, too large and a gable fits as one plane through the ridge. 0.10–0.20 m suits airborne LiDAR at 10 points/m².

Example 3 — classify from the plane statistics 

import numpy as np

def classify_roof(relief_m, planes, flat_threshold=0.8):
    if relief_m < flat_threshold:
        return "flat"
    big = [p for p in planes if p["n_points"] >= 25 and p["tilt_deg"] > 5]
    if not big:
        return "flat"
    if len(big) == 1:
        return "shed"
    azimuths = np.array([p["azimuth_deg"] for p in big])
    tilts = np.array([p["tilt_deg"] for p in big])

    def opposite(a, b, tol=35):
        d = abs((a - b + 180) % 360 - 180)
        return abs(d - 180) < tol

    if len(big) == 2:
        return "gable" if opposite(azimuths[0], azimuths[1]) else "complex"
    if len(big) in (3, 4) and tilts.std() < 8:
        return "hip"
    if len(big) >= 4 and tilts.std() >= 8:
        return "mansard"
    return "complex"

Returning "complex" rather than guessing is what keeps the classification honest. In any real stock a meaningful share of roofs are genuinely not one of the four textbook forms.

Explanation

Why the flat/pitched split does most of the work

Flat roofs make every downstream calculation the simple case: roof area equals footprint area, azimuth is a free choice, volume from a prism to the eaves is exact, and shadow from a prism is exact. On the reference dataset that covered 63.5% of buildings. Getting the remaining distinctions right matters much less to the totals than getting this one right does.

Why percentiles rather than min and max

A single return on a chimney, an aerial or a bird sets the maximum; a return that slipped through the gutter sets the minimum. The 10th and 90th percentiles are stable against both, at the cost of slightly underestimating the true relief of a very steep roof — which does not change the classification.

Why RANSAC rather than a clustering method

Roof planes are exactly what RANSAC is for: a small number of dominant linear models in data with structured outliers. Clustering on surface normals also works and needs normals estimated first, which on sparse airborne data is noisier than fitting the planes directly.

Why to compare against an existing model

Roof type codes in national LoD2 products follow a published convention, so a comparison is a confusion matrix rather than an opinion. The reference tile's distribution — 1,238 flat, 411 of the next commonest type, then a long tail across six more — is also a useful prior: if your classifier produces 40% "complex", something is wrong with the segmentation parameters rather than with the roofs.

Triage of four roof segmentation parameters — footprint inset, RANSAC threshold at both extremes, and minimum point count — with what each breaks.
Two of the four failures produce plausible planes that are not roof planes.

Edge cases or notes

  • Point density decides what is possible. Below about 4 points/m², plane segmentation is unreliable.
  • Parapets look like relief. A flat roof with a 1 m parapet fails a naive relief test.
  • Dormers add small planes. Filter by point count before classifying.
  • Buildings with several parts have several roofs; segment per part.
  • Solar panels are a plane at a different tilt from the roof beneath.
  • Vegetation classification is imperfect. Check for points above the plausible roof.
  • Terraces share ridges. Clipping to the footprint splits one ridge across several buildings.
  • Report the unclassifiable count. It is a quality measure, not an embarrassment.

FAQ

How do I tell a flat roof from a pitched one?

The 10th-to-90th percentile height range of first returns inside an inset footprint. Below about 0.8 m is flat.

Why inset the footprint?

Wall returns, guttering and the ground just outside contaminate every statistic. Buffering inwards by half a metre to a metre removes them.

How do I find the roof planes?

RANSAC plane fitting, removing inliers each round. A residual threshold of 0.10–0.20 m suits typical airborne LiDAR.

How do I distinguish a gable from a hip?

Two dominant planes with opposite azimuths is a gable; three or four planes of similar tilt is a hip. Two tilts per side is a mansard.

What point density do I need?

At least about 4 points per square metre for plane segmentation; the flat/pitched split works on less.

What if the roof is not one of the standard forms?

Label it "complex". Forcing an exotic roof into a category is worse than admitting the classifier could not decide.