Extracting Building Footprints from LiDAR

TL;DR: Rasterize class-6 points at roughly the point spacing, label connected roof regions, polygonize each with rasterio.features.shapes, close gaps with a mitred buffer out-and-in, snap each polygon to its dominant edge orientation, and write the result to a GeoPackage with eave and ridge heights taken from point percentiles.

# Context and Motivation

This guide is part of Building Extraction from LiDAR, which covers classifying roofs to class 6. Classification is only half of what a planning department or insurer asks for; the other half is a polygon per building they can load into a GIS, measure and join to an address register. Turning points into polygons is where most footprint projects lose quality — jagged staircase edges, holes where a skylight absorbed the laser, two houses fused through a shared garage — and every one of those defects is visible the moment someone overlays the result on an orthophoto.

The approach here is deliberately raster-first. Alpha shapes and concave hulls work directly on points and look elegant, but they are sensitive to density variations along a roof and slow on thousands of buildings. A binary mask at the point spacing, cleaned with morphology and polygonized, is fast, predictable and easy to regularize afterwards.

# Prerequisites and Assumptions

• A LAZ tile where buildings are class 6, produced by the parent workflow or by a vendor.
• PDAL 2.4+ and its Python bindings; rasterio, Shapely 2.x, GeoPandas and NumPy.
• A projected CRS in metres. The regularization step measures angles and lengths in map units.
• Average point spacing known to within a factor of two. At 10 pts/m² spacing is about 0.32 m; at 25 pts/m², 0.2 m.

# Step-by-Step Implementation

# Step 1 — Read only the roof points

Reading class 6 alone keeps the arrays small and avoids masking later.

p = pdal.Pipeline(json.dumps({"pipeline": [
    {"type": "readers.las", "filename": "tile_5840_2710_bldg.laz"},
    {"type": "filters.range", "limits": "Classification[6:6]"},
    {"type": "filters.hag_nn", "count": 2},
]}))

filters.hag_nn needs ground points to compute heights, so place the range filter after it if your tile has ground in class 2 — the version in the complete example does exactly that.

# Step 2 — Rasterize to a binary mask

Choose a cell about equal to the point spacing, burn a 1 into every cell containing a roof return, and dilate by one cell so that single missing cells do not split a roof.

# Step 3 — Label and polygonize

scipy.ndimage.label assigns an integer to each connected region; rasterio.features.shapes turns each labelled region into a polygon in map coordinates using the raster transform. Using labels rather than the binary mask directly means each polygon carries an ID you can join back to the points.

# Step 4 — Close gaps and drop slivers

Buffer each polygon out by one cell with join_style="mitre" and back in by the same amount. This fills notches narrower than two cells without rounding corners. Drop polygons smaller than your minimum building area.

# Step 5 — Regularize the outline

Find the polygon’s dominant edge direction from the minimum rotated rectangle, rotate so that direction is horizontal, snap near-horizontal and near-vertical edges, then rotate back. Buildings with genuinely non-orthogonal walls should skip this step; flag them by how much area regularization changed.

# Step 6 — Attach heights and write

Take eave height as the 10th percentile and ridge height as the 98th percentile of HeightAboveGround of the points inside each polygon, and write to a GeoPackage layer.

# Complete Working Example

"""Class-6 points to regularized footprint polygons with heights."""
from __future__ import annotations

import json
import logging
from pathlib import Path

import geopandas as gpd
import numpy as np
import pdal
import rasterio.features
from rasterio.transform import from_origin
from scipy import ndimage as ndi
from shapely import affinity
from shapely.geometry import Polygon, shape

log = logging.getLogger("footprints")


def roof_points(src: Path) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.range", "limits": "Classification[6:6]"},
    ]}))
    p.execute()
    return p.arrays[0]


def regularize(poly: Polygon, snap_deg: float = 12.0) -> Polygon:
    mrr = poly.minimum_rotated_rectangle
    x, y = mrr.exterior.coords.xy
    edges = [(x[i + 1] - x[i], y[i + 1] - y[i]) for i in range(2)]
    dx, dy = max(edges, key=lambda e: e[0] ** 2 + e[1] ** 2)
    angle = np.degrees(np.arctan2(dy, dx))
    rot = affinity.rotate(poly, -angle, origin="centroid")
    coords = np.asarray(rot.exterior.coords)
    snapped = [coords[0]]
    for pt in coords[1:]:
        prev = snapped[-1]
        seg = np.degrees(np.arctan2(pt[1] - prev[1], pt[0] - prev[0])) % 180
        if min(seg, 180 - seg) < snap_deg:          # near horizontal
            pt = np.array([pt[0], prev[1]])
        elif abs(seg - 90) < snap_deg:              # near vertical
            pt = np.array([prev[0], pt[1]])
        snapped.append(pt)
    out = Polygon(snapped).buffer(0)
    return affinity.rotate(out, angle, origin=poly.centroid) if out.is_valid else poly


def footprints(src: Path, dst: Path, crs: str = "EPSG:6347",
               cell: float = 0.5, min_area: float = 20.0) -> gpd.GeoDataFrame:
    pts = roof_points(src)
    x0, y1 = np.floor(pts["X"].min()) - cell, np.ceil(pts["Y"].max()) + cell
    cols = int(np.ceil((pts["X"].max() - x0) / cell)) + 2
    rows = int(np.ceil((y1 - pts["Y"].min()) / cell)) + 2
    mask = np.zeros((rows, cols), dtype=bool)
    mask[((y1 - pts["Y"]) / cell).astype(int), ((pts["X"] - x0) / cell).astype(int)] = True
    mask = ndi.binary_closing(mask, iterations=1)
    labels, n = ndi.label(mask)
    log.info("%d connected roof regions", n)

    transform = from_origin(x0, y1, cell, cell)
    col_of = ((pts["X"] - x0) / cell).astype(int)
    row_of = ((y1 - pts["Y"]) / cell).astype(int)
    point_label = labels[row_of, col_of]

    records = []
    for geom, lab in rasterio.features.shapes(labels.astype(np.int32), mask=labels > 0,
                                              transform=transform):
        poly = shape(geom).buffer(cell, join_style="mitre").buffer(-cell, join_style="mitre")
        if poly.is_empty or poly.area < min_area:
            continue
        reg = regularize(poly)
        hag = pts["HeightAboveGround"][point_label == int(lab)]
        records.append({
            "region": int(lab),
            "area_m2": round(reg.area, 1),
            "reg_change": round(abs(reg.area - poly.area) / poly.area, 3),
            "eave_m": round(float(np.percentile(hag, 10)), 2),
            "ridge_m": round(float(np.percentile(hag, 98)), 2),
            "geometry": reg,
        })
    gdf = gpd.GeoDataFrame(records, crs=crs)
    gdf.to_file(dst, layer="footprints", driver="GPKG")
    log.info("%s: %d footprints", dst.name, len(gdf))
    return gdf


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    footprints(Path("tile_5840_2710_bldg.laz"), Path("tile_5840_2710_footprints.gpkg"))

# Key Parameter Table

# Verification

Check the output in three ways before handing it over.

1. Count against classified regions. The number of polygons should equal the number of labelled regions above the area threshold. A large difference means the buffer step is erasing thin buildings.
2. Regularization change. The reg_change column records how much area regularization moved. Values above 0.05 are worth a look; values above 0.15 usually indicate a curved or angled building that should keep its raw outline.
3. Visual overlay. Load the GeoPackage over an orthophoto or a hillshaded DSM and inspect twenty buildings from different parts of the tile.

gdf = gpd.read_file("tile_5840_2710_footprints.gpkg", layer="footprints")
assert gdf.is_valid.all(), "invalid geometries in output"
assert (gdf.ridge_m >= gdf.eave_m).all(), "ridge below eave"
print(gdf.reg_change.describe())

# Gotchas and Edge Cases

LiDAR footprints are roof outlines, not wall lines. Eaves overhang walls by 0.3 to 1 m. Cadastral footprints are usually wall lines, so a LiDAR layer will look systematically larger. If the client needs wall lines, buffer inward by a project-specific overhang and say so in the metadata.

Shared walls fuse neighbours. Semi-detached and terraced houses are one connected roof. Split them with a parcel layer, or with height discontinuities along the ridge line, rather than by eroding the mask, which damages every other building.

Courtyards become filled. Morphological closing fills holes smaller than the structuring element but leaves larger courtyards. Check poly.interiors if the client needs inner rings preserved, and never close with more than two iterations.

Tile edges cut buildings. A building on a tile boundary yields two partial polygons. Run on buffered tiles and keep a polygon only if its centroid falls inside the tile’s nominal extent.

# Frequently Asked Questions

Should I use alpha shapes instead of a raster mask?

Alpha shapes can follow curved walls more closely, but they are sensitive to the alpha parameter and to density changes across a roof, and slow on large tiles. The raster route is faster and more predictable; use alpha shapes only for the minority of buildings that regularization flags as non-rectilinear.

What raster cell size gives the best footprints?

About the average point spacing. Smaller cells leave empty cells inside roofs that closing must repair; larger cells round corners and merge buildings separated by narrow gaps.

How do I get wall lines instead of roof outlines?

LiDAR sees roofs from above, so the outline is the eave line. Buffer inward by the typical overhang for the building stock, or use wall returns from oblique or mobile scans if they exist.

Can this produce 3D building models?

It produces the footprint and two heights, which is enough for a block model at LOD1. Roof shapes need plane segmentation per building, covered in estimating roof pitch and aspect from normals.

# Related

• Building Extraction from LiDAR — the classification workflow that produces class 6
• Detecting Planar Roofs with Covariance Features — improving the roof classification itself
• Estimating Roof Pitch and Aspect from Normals — per-plane geometry inside each footprint
• Extracting Objects from Segment Labels — the same object-table pattern for any class
• LiDAR Classification and Feature Extraction — the section overview