Estimating Roof Pitch and Aspect from Normals

TL;DR: Run filters.normal with always_up: true on class-6 points, cluster each building’s points by normal direction to separate roof planes, fit a plane to each group with a least-squares SVD, and report pitch as degrees(arccos(|nz|)) and aspect as degrees(atan2(nx, ny)) mod 360 — a south-facing 35° plane comes out as pitch 35, aspect 180.

# Context and Motivation

This guide is part of Building Extraction from LiDAR. Solar potential studies, insurance roof surveys and LOD2 city models all need the same two numbers per roof face: how steep it is and which way it faces. A footprint with a single height cannot answer either. The surface normal at each roof point already encodes both, and PDAL computes it from the same kind of neighbourhood analysis used for planarity: the normal is the eigenvector of the smallest eigenvalue.

Per-point normals are noisy, though, and a roof has several faces. The useful output is per plane, which means grouping points into faces first and then estimating one normal per face from all of its points.

# Prerequisites and Assumptions

• Class-6 roof points, from the building extraction workflow or a vendor delivery.
• Building IDs per point, either ClusterID from segmentation or a footprint polygon burned in with filters.overlay.
• PDAL 2.4+ with filters.normal; Python with NumPy, pandas and scikit-learn (for DBSCAN on normal vectors).
• A projected CRS whose Y axis points to grid north. Aspect is measured relative to grid north; the difference from true north (grid convergence) is usually under two degrees but matters for precise solar work.

# Step-by-Step Implementation

# Step 1 — Compute normals facing upward

always_up: true flips any normal with a negative Z component, so every roof normal points to the sky and aspects are not reversed at random.

{ "type": "filters.normal", "knn": 12, "always_up": true, "where": "Classification == 6" }

# Step 2 — Group each building’s points by normal direction

Within one building, points on the same face share a normal. DBSCAN on the three normal components with a small eps (0.08, roughly 5 degrees) separates faces; add scaled X and Y if parallel faces on opposite sides of a courtyard must stay apart.

# Step 3 — Fit a plane per face

The normal averaged over a face is decent, but an SVD plane fit over the face’s coordinates is better, because it uses every point at once rather than averaging noisy local estimates.

# Step 4 — Convert the normal to pitch and aspect

Pitch is the angle between the normal and vertical. Aspect is the compass bearing of the normal’s horizontal component, measured clockwise from north.

# Step 5 — Filter and report

Drop faces with fewer than about 30 points or an area under 4 m², and flag faces with pitch under 5° as flat, where aspect is meaningless.

# Complete Working Example

"""Per-face roof pitch and aspect from PDAL normals."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pandas as pd
import pdal
from sklearn.cluster import DBSCAN


def roof_points(src: Path) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.range", "limits": "Classification[6:6]"},
        {"type": "filters.normal", "knn": 12, "always_up": True},
        {"type": "filters.cluster", "tolerance": 1.0, "min_points": 50, "is3d": True},
    ]}))
    p.execute()
    return p.arrays[0]


def fit_plane(xyz: np.ndarray) -> np.ndarray:
    centred = xyz - xyz.mean(axis=0)
    _, _, vt = np.linalg.svd(centred, full_matrices=False)
    n = vt[2]
    return n if n[2] >= 0 else -n


def pitch_aspect(n: np.ndarray) -> tuple[float, float]:
    pitch = float(np.degrees(np.arccos(np.clip(abs(n[2]), 0.0, 1.0))))
    aspect = float(np.degrees(np.arctan2(n[0], n[1])) % 360.0)
    return pitch, aspect


def faces(points: np.ndarray, eps: float = 0.08, min_pts: int = 30) -> pd.DataFrame:
    rows = []
    for bid in np.unique(points["ClusterID"]):
        if bid == 0:
            continue
        b = points[points["ClusterID"] == bid]
        normals = np.column_stack([b["NormalX"], b["NormalY"], b["NormalZ"]])
        labels = DBSCAN(eps=eps, min_samples=10).fit_predict(normals)
        for face in np.unique(labels[labels >= 0]):
            sel = b[labels == face]
            if len(sel) < min_pts:
                continue
            xyz = np.column_stack([sel["X"], sel["Y"], sel["Z"]])
            n = fit_plane(xyz)
            pitch, aspect = pitch_aspect(n)
            # Plan area of the face from its points, corrected for slope.
            plan = len(sel) / max(len(b), 1) * np.ptp(b["X"]) * np.ptp(b["Y"])
            rows.append({"building": int(bid), "face": int(face), "points": len(sel),
                         "pitch_deg": round(pitch, 1),
                         "aspect_deg": round(aspect, 0) if pitch >= 5 else np.nan,
                         "flat": pitch < 5,
                         "area_m2": round(plan / max(np.cos(np.radians(pitch)), 0.2), 1)})
    return pd.DataFrame(rows)


if __name__ == "__main__":
    table = faces(roof_points(Path("tile_5840_2710_bldg.laz")))
    print(table.head(12).to_string(index=False))

Example output for a hip-roofed house and a flat-roofed garage:

 building  face  points  pitch_deg  aspect_deg   flat  area_m2
       14     0     612       31.8       178.0  False     41.2
       14     1     588       32.4       358.0  False     39.6
       14     2     204       33.1        88.0  False     13.8
       14     3     197       32.6       269.0  False     13.1
       15     0     341        1.9         NaN   True     22.4

# Key Parameter Table

# Verification

• Symmetric roofs agree. Opposite faces of a gable or hip roof should have pitches within a degree or two and aspects 180° apart. Large asymmetry on a roof you know is symmetric points to a face that absorbed ridge or eave points.
• Plane residuals. The RMS distance of face points to the fitted plane should be a few centimetres. Residuals above 0.15 m mean two faces were merged.
• Known roofs. Check two or three roofs against drawings or a site visit; typical residential pitches of 25 to 45 degrees give an immediate sanity check.

def plane_rms(xyz: np.ndarray, n: np.ndarray) -> float:
    d = (xyz - xyz.mean(axis=0)) @ n
    return float(np.sqrt(np.mean(d ** 2)))

# Gotchas and Edge Cases

Aspect flips on near-flat roofs. Tiny tilts from drainage falls produce arbitrary aspects. That is why aspect is reported only above 5 degrees of pitch.

Dormers join the main face. A dormer with the same aspect but different pitch has a slightly different normal; with a generous eps it merges into the main face. Lower eps, or cluster on pitch and aspect angles instead of raw normal components.

Grid north is not true north. Aspect is measured against the projected CRS’s Y axis. In UTM, grid convergence reaches a couple of degrees toward zone edges. For solar yield it rarely matters; for anything precise, correct with the convergence angle from pyproj.

Walls contaminate faces. Wall returns along eaves have near-horizontal normals and pitches near 90 degrees. Filter points with NormalZ < 0.3 before clustering; no roof you care about is steeper than about 72 degrees.

# Frequently Asked Questions

Why fit a plane instead of averaging point normals?

Point normals come from small neighbourhoods and are noisy, especially at edges. A plane fitted to all of a face’s points uses far more information, gives a more accurate orientation and also provides residuals that tell you whether the face was segmented correctly.

What does always_up do in filters.normal?

The sign of an eigenvector is arbitrary, so half the normals could point into the roof. always_up flips any normal with a negative vertical component, making all roof normals point upward and aspects consistent.

How accurate is LiDAR roof pitch?

On well-sampled faces, within one or two degrees of design pitch. Small faces, low density and dormers degrade it; report the point count per face so users can judge.

Can I compute solar potential from these values?

Pitch, aspect and area are the main geometric inputs to solar yield models. Shading from trees and neighbouring buildings also matters, which requires a DSM-based horizon analysis on top of these per-face values.

# Related

• Building Extraction from LiDAR — classifying the roof points used here
• Detecting Planar Roofs with Covariance Features — the same neighbourhood analysis applied to classification
• Extracting Building Footprints from LiDAR — the outlines these faces sit inside
• Generating Slope and Aspect Rasters with gdaldem — the raster equivalent for terrain
• DBSCAN Segmentation with filters.dbscan — the clustering method used on normals