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 the LoD2 tile of The Hague used here, 63.5% of the solids have a ridge height equal to their eaves height — flat — 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, pandas as pd, geopandas as gpd, laspy

las = laspy.read("tile.laz")
keep = (las.return_number == 1) & (las.classification == 6)          # AHN: 6 = building
pts = gpd.GeoDataFrame({"z": np.asarray(las.z[keep])},
                       geometry=gpd.points_from_xy(las.x[keep], las.y[keep]), crs=b.crs)
inner = b.assign(geometry=b.buffer(-0.75))                            # keep clear of walls and gutters
inner = inner[~inner.is_empty]
z = gpd.sjoin(pts, inner[["id", "geometry"]], predicate="within").groupby("id")["z"]
relief = z.quantile(0.9) - z.quantile(0.1)

b["relief_m"] = b["id"].map(relief)
b["roof_class"] = np.where(b.relief_m < 0.8, "flat", "pitched")
k = b[b.relief_m.notna()]
lod2 = pd.Series(np.where(k.height_ridge_m == k.height_eaves_m, "flat", "pitched"),  # the model's own heights
                 index=k.index, name="LoD2")
print(pd.crosstab(lod2, k.roof_class.rename("relief")))
print(f"agree for {(lod2 == k.roof_class).mean():.1%} of {len(k)} solids with roof points ({len(b) - len(k)} without)")

relief   flat  pitched
LoD2                  
flat      755      121
pitched    91      584
agree for 86.3% of 1551 solids with roof points (439 without)

Percentiles rather than min and max, because a chimney, an aerial, a parapet or an overhanging branch sets the extremes. Checked against the city model's own eaves and ridge heights, one number agrees on flat versus pitched for 86.3% of the 1,551 solids that have roof returns: 121 flat roofs read as pitched and 91 pitched ones as flat. The other 439 solids are parts too small or too narrow to keep any return once the footprint is inset by 0.75 m.

Clip inset, filter returns, measure relief, fit planes

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: without the 0.75 m inset, flat-versus-pitched agreement with the city model falls from 87.3% to 64.6% (Example 1).

2. Keep the right returns

First returns only, and only the building class if the survey has one. An overhanging tree adds a second surface above the roof and turns a gable into a nonsense. Check what the classes are before filtering: AHN4 has no vegetation classes at all — trees are in class 1, unclassified — so excluding the ASPRS vegetation codes 3–5 removes nothing. Keeping class 6 (building) is the filter that works; in this dense city centre, with few trees over roofs, it moved agreement only from 86.2% to 87.3%.

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 — the CityGML 2.0 code list, in which this tile's solids are 1000 flat (1,238 of 1,990), 1030 gable (411), 1130 combination (138) and 1010 shed (110) — and mapping your classes onto them makes the comparison concrete, as Example 3 does.

Relief and RANSAC planes per building

The data used on this page

The LoD2 tile of The Hague's 3D city model published as a CityJSON sample dataset (DenHaag_01.city.json), whose solids carry a roof-type code, and the AHN4 point cloud over it — sub-tile 30DN2_19 from TU Delft's GeoTiles (CC0), cut to the model's extent and saved as tile.laz. b has one row per solid: the ground-surface footprint, the city's roof code, and its eaves and ridge heights above ground.

import json
import os

import geopandas as gpd
import numpy as np
import requests
from shapely.geometry import Polygon
from shapely.ops import unary_union

if not os.path.exists("model.city.json"):
    reply = requests.get("https://3d.bk.tudelft.nl/opendata/cityjson/3dcities/v2.0/DenHaag_01.city.json", timeout=300)
    reply.raise_for_status()
    open("model.city.json", "wb").write(reply.content)
d = json.load(open("model.city.json"))
V = np.array(d["vertices"]) * d["transform"]["scale"] + d["transform"]["translate"]

rows_b = []
for oid, o in d["CityObjects"].items():
    a = o.get("attributes", {})
    for g in o.get("geometry", []):
        if g["type"] != "Solid" or "RelativeEavesHeight" not in a:
            continue
        sem = g["semantics"]
        ground = [Polygon(V[surface[0]][:, :2])
                  for shell, values in zip(g["boundaries"], sem["values"])
                  for surface, v in zip(shell, values)
                  if v is not None and sem["surfaces"][v]["type"] == "GroundSurface"]
        rows_b.append({"id": oid, "roof_type": a.get("roofType"), "height_eaves_m": a["RelativeEavesHeight"],
                       "height_ridge_m": a["RelativeRidgeHeight"], "geometry": unary_union(ground)})
b = gpd.GeoDataFrame(rows_b, crs=28992)

# AHN4 point cloud: sheet 30DN2, sub-tile 19, from TU Delft's GeoTiles (CC0), cut to the model tile
if not os.path.exists("tile.laz"):
    import laspy
    if not os.path.exists("30DN2_19.LAZ"):
        with requests.get("https://geotiles.citg.tudelft.nl/AHN4_T/30DN2_19.LAZ", stream=True, timeout=600) as r:
            r.raise_for_status()
            with open("30DN2_19.LAZ", "wb") as f:
                for chunk in r.iter_content(1 << 20):
                    f.write(chunk)
    las = laspy.read("30DN2_19.LAZ")
    x0, y0, x1, y1 = b.total_bounds
    inside = (las.x >= x0 - 5) & (las.x <= x1 + 5) & (las.y >= y0 - 5) & (las.y <= y1 + 5)
    out = laspy.LasData(las.header)
    out.points = las.points[inside]
    out.write("tile.laz")
print(len(b), "solids;", b.roof_type.value_counts().to_dict())
print(f"ridge height equal to eaves height: {(b.height_ridge_m == b.height_eaves_m).mean():.1%}")

1990 solids; {'1000': 1238, '1030': 411, '1130': 138, '1010': 110, '1040': 40, '1120': 34, '1070': 13, '1060': 5, '9999': 1}
ridge height equal to eaves height: 63.5%

Example 1 — clip, filter and compute relief

import numpy as np, pandas as pd, geopandas as gpd, laspy

def roof_points(las_path, buildings, classes=(6,), exclude=(), inset=0.75, min_points=30):
    las = laspy.read(las_path)
    keep = np.asarray(las.return_number == 1)
    if classes:
        keep &= np.isin(las.classification, classes)
    if exclude:
        keep &= ~np.isin(las.classification, exclude)
    # laspy's scaled z is a view, not an array: GeoDataFrame reads it as a scalar without asarray
    pts = gpd.GeoDataFrame({"z": np.asarray(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

las = laspy.read("tile.laz")
classes, counts = np.unique(np.asarray(las.classification), return_counts=True)
print("AHN4 classes in the tile:", {int(c): f"{n:,}" for c, n in zip(classes, counts)})

lod2_flat = b.set_index("id").eval("height_ridge_m == height_eaves_m")
for label, kw in {"first returns, not class 3-5 (first version)": dict(classes=(), exclude=(3, 4, 5)),
                  "first returns, class 6 (building)": dict(),
                  "class 6, no inset": dict(inset=0)}.items():
    stats, points = roof_points("tile.laz", b, **kw)
    flat = stats.relief_m < 0.8
    agree = (flat == lod2_flat[stats.index]).mean()
    print(f"{label:45} {len(stats):,} buildings with 30+ returns; flat {flat.mean():.1%}; "
          f"agrees with LoD2 {agree:.1%}")

AHN4 classes in the tile: {1: '2,013,216', 2: '8,415,418', 6: '7,454,532'}
first returns, not class 3-5 (first version)  1,392 buildings with 30+ returns; flat 49.8%; agrees with LoD2 86.2%
first returns, class 6 (building)             1,375 buildings with 30+ returns; flat 50.5%; agrees with LoD2 87.3%
class 6, no inset                             1,880 buildings with 30+ returns; flat 29.6%; agrees with LoD2 64.6%

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. On AHN4's 27 roof returns per m², 0.15 m is the setting Example 4 found best.

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"

import time

stats, points = roof_points("tile.laz", b)
# CityGML 2.0 roof type codes as the model uses them; 1070 (pavilion) is four planes to a point
CODE = {"1000": "flat", "1010": "shed", "1030": "gable", "1040": "hip", "1070": "hip",
        "1060": "mansard", "1120": "complex", "1130": "complex"}
codes = b.set_index("id")["roof_type"]

start, rows = time.perf_counter(), []
for bid, g in points.groupby("id"):
    if bid not in stats.index:
        continue
    xyz = np.c_[g.geometry.x - g.geometry.x.mean(), g.geometry.y - g.geometry.y.mean(), g.z]
    relief = stats.relief_m[bid]
    planes = segment_planes(xyz) if relief >= 0.8 else []
    rows.append({"id": bid, "lod2": CODE.get(codes[bid], "other"), "lidar": classify_roof(relief, planes)})
r = pd.DataFrame(rows)
print(f"{len(r):,} buildings in {time.perf_counter() - start:.0f} s")
print(pd.crosstab(r.lod2, r.lidar, margins=True))
print(f"same class as the LoD2 model: {(r.lod2 == r.lidar).mean():.1%}; "
      f"pitched roofs only: {(r.lod2 == r.lidar)[r.lod2 != 'flat'].mean():.1%}")

1,375 buildings in 27 s
lidar    complex  flat  gable  hip  mansard  shed   All
lod2                                                   
complex       16    18     19   75       19     2   149
flat          19   651     13    8        0    31   722
gable         45    20    270   35       16    13   399
hip            3    15      6   23        2     1    50
mansard        1     1      0    2        1     0     5
other          0     1      0    0        0     0     1
shed           5    28      0    2        0    14    49
All           89   734    308  145       38    61  1375
same class as the LoD2 model: 70.9%; pitched roofs only: 49.6%

Against the city model the classifier gets 70.9% of roofs into the same class, but that is mostly the flat ones; on pitched roofs it is 49.6%. Gables are the strong case (270 of 399). The rules for the rest are too coarse for real roofs: 75 of the 149 roofs the model calls combinations come out as "hip", because three or four similar planes is what a combination roof often is, and 28 of the 49 shed roofs come out flat. 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.

Example 4 — how sensitive the parameters are

The parameters below were measured, not assumed — on the buildings the LoD2 model calls gable roofs, with the relief test passed:

rng = np.random.default_rng(0)
area = b.set_index("id").buffer(-0.75).area
density = stats.n / area[stats.index]
print(f"roof returns per m²: median {density.median():.1f}")

gables = [bid for bid in stats.index if codes[bid] == "1030" and stats.relief_m[bid] >= 0.8]
groups = {bid: g for bid, g in points.groupby("id") if bid in set(gables)}

def gable_rate(threshold=0.15, keep=1.0):
    hits, n_planes = 0, []
    for bid, g in groups.items():
        g = g.sample(frac=keep, random_state=0) if keep < 1 else g
        xyz = np.c_[g.geometry.x - g.geometry.x.mean(), g.geometry.y - g.geometry.y.mean(), g.z]
        planes = segment_planes(xyz, residual_threshold=threshold)
        n_planes.append(len([p for p in planes if p["n_points"] >= 25 and p["tilt_deg"] > 5]))
        hits += classify_roof(stats.relief_m[bid], planes) == "gable"
    return hits / len(groups), np.median(n_planes)

for t in (0.05, 0.15, 0.40):
    rate, planes = gable_rate(threshold=t)
    print(f"threshold {t:.2f} m: gable {rate:.1%}, median {planes:.0f} planes")
for pts_m2 in (8, 4, 2):
    rate, planes = gable_rate(keep=pts_m2 / density.median())
    print(f"thinned to ~{pts_m2} returns/m²: gable {rate:.1%}, median {planes:.0f} planes")

roof returns per m²: median 27.0
threshold 0.05 m: gable 28.1%, median 3 planes
threshold 0.15 m: gable 70.9%, median 2 planes
threshold 0.40 m: gable 63.5%, median 2 planes
thinned to ~8 returns/m²: gable 79.3%, median 2 planes
thinned to ~4 returns/m²: gable 66.4%, median 2 planes
thinned to ~2 returns/m²: gable 26.0%, median 1 planes

At 0.05 m the planes split — a median of three per gable roof — and gables fall to 28.1%; at 0.40 m they fall only to 63.5%. Density is less simple than a threshold. Thinning the cloud to 8 returns per m² raised the gable rate to 79.3%, because the minimum of 25 inliers then filters out dormers and chimneys that it keeps at full density; at 4 per m² the rate is 66.4%, and at 2 it collapses to 26.0%. The inlier minimum is a point count, so it means a different area at every density — scale it with the density rather than fixing it.

The flat-or-pitched split does most of the work

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 tile's distribution — 1,238 flat, 411 gable, then a long tail across seven more codes — is also a useful prior: if your classifier produces 40% "complex", something is wrong with the segmentation parameters rather than with the roofs.

Point density, parapets and solar panels

• Point density decides what is possible. Thinned to 2 returns/m², gables were found for 26.0% of roofs against 66.4% at 4 (Example 4).
• 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.

Internal links

• Where building heights come from, and how wrong they are — the other thing LiDAR gives you
• Level of detail explained for 3D city models — what this gets you towards
• How to read a CityJSON city model in Python — the model to validate against
• How to estimate rooftop solar potential in Python — the main consumer of roof planes
• How to read LAS and LAZ files in Python with laspy — loading the point cloud
• How to filter a point cloud by class and return — the return and class filtering
• How to build a canopy height model in Python — separating vegetation properly
• LiDAR surfaces explained: DSM, DTM and CHM — what the returns represent

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. On AHN4 a residual threshold of 0.15 m worked best of 0.05, 0.15 and 0.40 m.

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?

On AHN4 roofs thinned for the test, gables were found for 66.4% at 4 returns per m² and 26.0% at 2 — plan on at least 4. The flat/pitched split needs only enough points for two percentiles.

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.