How to extrude building footprints into 3D in Python

Problem statement

Extrusion is the cheapest 3D there is: a polygon, a base height, a top height, and a closed prism. It is what turns a footprint layer into something you can compute shadows, volumes and visibility from, and it is the step every 2.5D workflow starts with.

The arithmetic is trivial. What is not trivial is everything around it — which height to extrude to, what the base is when the ground slopes, what to do with holes and multipart geometries, and whether the result is a watertight solid or a bag of triangles that looks right and measures nothing.

This guide extrudes a real footprint layer correctly, checks the result, and writes it out in the formats the next step needs.

Quick answer

import numpy as np, trimesh
from shapely.geometry import Polygon

def extrude(polygon: Polygon, base: float, top: float) -> trimesh.Trimesh:
    """A watertight prism from a 2D polygon, with holes handled."""
    mesh = trimesh.creation.extrude_polygon(polygon, height=top - base)
    mesh.apply_translation([0, 0, base])
    return mesh

meshes = [extrude(geom, base, base + h)
          for geom, base, h in zip(b.geometry, b.ground_z, b.height_m)]
scene = trimesh.util.concatenate(meshes)
print(f"{sum(m.is_watertight for m in meshes)} of {len(meshes)} watertight; "
      f"{sum(m.volume for m in meshes):,.0f} m³ to the eaves")

1990 of 1990 watertight; 627,257 m³ to the eaves

trimesh.creation.extrude_polygon triangulates the polygon including its interior rings, caps both ends and welds the walls, so the result is watertight by construction. Hand-rolled extrusions usually are not.

Projected CRS, terrain base, chosen top, multipolygons

1. Work in a projected CRS

Extrusion mixes horizontal and vertical units. In degrees, a 12 m extrusion of a polygon whose coordinates are around 4.3 and 52.1 — this tile, in WGS84 — produces a building twelve degrees tall. Reproject first.

2. Get the base height from the terrain, not from zero

A building extruded from z = 0 sits at sea level — in this part of The Hague that is a median 5.96 m below the street. Sample the DTM and use a robust statistic, the median, or the minimum if you want the building to sit on its lowest ground contact.

Sample it around the footprint, not inside it. A DTM made from lidar has no ground under a roof: AHN leaves those cells empty, and inside these footprints it has a median of 3 cells where the area would hold 98. For 545 of the 1,990 solids there is not one, and a pipeline that skips empty samples silently drops a quarter of the buildings.

import numpy as np, rasterio
from rasterio.mask import mask

def ground_height(dtm_path, geom, ring_m=3, stat="median"):
    """Terrain in a ring just outside the walls: a lidar DTM is empty under the roof."""
    with rasterio.open(dtm_path) as src:
        arr, _ = mask(src, [geom.buffer(ring_m).difference(geom)], crop=True, filled=False)
    values = arr.compressed()
    if not values.size:
        return np.nan
    return float(np.median(values) if stat == "median" else values.min())

inside, spread = [], []
with rasterio.open("dtm.tif") as src:
    for geom in b.geometry:
        arr, _ = mask(src, [geom], crop=True, filled=False)
        inside.append(arr.compressed().size)
        ring, _ = mask(src, [geom.buffer(3).difference(geom)], crop=True, filled=False)
        spread.append(np.ptp(ring.compressed()) if ring.compressed().size else np.nan)
inside, spread = np.array(inside), np.array(spread)
print(f"DTM cells inside a footprint: median {np.median(inside):.0f}, "
      f"where the area would hold {np.median(b.area / 0.25):.0f}; none at all for {(inside == 0).sum()} of {len(b)}")
print(f"terrain in the 3 m ring spans a median {np.nanmedian(spread):.2f} m, "
      f"{np.nanpercentile(spread, 95):.2f} m at the 95th percentile")
print(f"extruded from z = 0, {(b.height_m < b.ground_z).mean():.0%} of the buildings would have their eaves below the ground")

DTM cells inside a footprint: median 3, where the area would hold 98; none at all for 545 of 1990
terrain in the 3 m ring spans a median 0.39 m, 1.87 m at the 95th percentile
extruded from z = 0, 55% of the buildings would have their eaves below the ground

3. Decide which top height you are extruding to

Eaves gives a flat-roofed block that under-represents volume; ridge over-represents it. On this tile two thirds of the solids have flat roofs, where the two are the same, and on the pitched ones the ridge is a median 3.15 m above the eaves (where building heights come from has the breakdown). Example 1 measures what that does to the volume. If the choice matters, say which you used.

4. Handle multipart footprints and holes

A MultiPolygon is several buildings or several parts of one; extrude each part and keep the relationship. Interior rings are courtyards and must stay holes — extrude_polygon handles them, a naive wall-building loop does not. (This tile has none, which is exactly how a courtyard bug survives testing; the first version of Example 2 dropped them from its caps.)

5. Check the result is watertight

bad = [i for i, m in enumerate(meshes) if not m.is_watertight]
print(f"{len(bad)} of {len(meshes)} meshes are not watertight")

0 of 1990 meshes are not watertight

Anything not watertight has no volume and will fail a boolean operation later.

6. Decide how to write it out

• GeoPackage with PolygonZ — keeps it in the GIS world, loses the solid.
• CityJSON — a real solid with semantics, readable by 3D tools.
• glTF / OBJ / PLY — for viewers and renderers.
• A height column on the footprints — for web maps that extrude client-side.

The last is usually the right answer for delivery: the browser extrudes, and the payload stays 2D.

7. Keep the attributes

The mesh loses everything except geometry. Carry an id array alongside so the extruded model can be joined back to the footprint table.

Extrusion with terrain and checks

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), and the AHN4 terrain model under it at 0.5 m from PDOK (CC0). b has one row per solid: its ground-surface footprint, the eaves and ridge heights above ground the city recorded, and ground_z, the height of the solid's ground surface in the model (metres above NAP). height_m is the eaves height.

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"]
        ground_z = [V[surface[0]][:, 2].min()
                    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"], "ground_z": min(ground_z),
                       "geometry": unary_union(ground)})
b = gpd.GeoDataFrame(rows_b, crs=28992)
b["height_m"] = b["height_eaves_m"]

# AHN4 terrain, 0.5 m, from PDOK's WCS: the tile plus a margin
if not os.path.exists("dtm.tif"):
    reply = requests.get("https://service.pdok.nl/rws/ahn/wcs/v1_0", timeout=300, params=[
        ("SERVICE", "WCS"), ("VERSION", "2.0.1"), ("REQUEST", "GetCoverage"), ("COVERAGEID", "dtm_05m"),
        ("SUBSET", "x(78200,79100)"), ("SUBSET", "y(457550,458330)"), ("FORMAT", "image/tiff")])
    reply.raise_for_status()
    open("dtm.tif", "wb").write(reply.content)
print(f"{len(b)} solids; ground in the model {b.ground_z.min():.2f}–{b.ground_z.max():.2f} m NAP, "
      f"median {b.ground_z.median():.2f}; median eaves {b.height_m.median():.2f} m above it")

1990 solids; ground in the model 3.00–22.22 m NAP, median 5.96; median eaves 5.61 m above it

Example 1 — the full pipeline, with terrain and checks

import geopandas as gpd, numpy as np, pandas as pd, rasterio, trimesh
from rasterio.mask import mask

def extrude_layer(buildings, dtm_path, height_col="height_m", crs="EPSG:28992", ring_m=3):
    b = buildings.to_crs(crs).explode(index_parts=False).reset_index(drop=True)

    with rasterio.open(dtm_path) as src:
        bases = []
        for geom in b.geometry:
            # the ring outside the walls: AHN has no ground under a roof
            arr, _ = mask(src, [geom.buffer(ring_m).difference(geom)], crop=True, filled=False)
            vals = arr.compressed()
            bases.append(float(np.median(vals)) if vals.size else np.nan)
    b["ground_z"] = bases

    meshes, ids = [], []
    for row in b.itertuples():
        if not np.isfinite(row.ground_z) or not np.isfinite(getattr(row, height_col)):
            continue
        m = trimesh.creation.extrude_polygon(row.geometry, getattr(row, height_col))
        m.apply_translation([0, 0, row.ground_z])
        meshes.append(m)
        ids.append(row.Index)

    b["extruded"] = b.index.isin(ids)
    return meshes, ids, b

model_z = b.set_index("id")["ground_z"]
meshes, ids, t = extrude_layer(b.drop(columns="ground_z"), "dtm.tif")
print(f"extruded {len(meshes):,} of {len(t):,}; no terrain in the ring for {(~t.extruded).sum()}; "
      f"not watertight: {sum(not m.is_watertight for m in meshes)}")
eaves = sum(m.volume for m in meshes)
ridge = sum(trimesh.creation.extrude_polygon(g, h).volume
            for g, h in zip(t.geometry[ids], t.height_ridge_m[ids]))
print(f"volume to the eaves {eaves:,.0f} m³, to the ridge {ridge:,.0f} m³ ({ridge / eaves - 1:+.1%})")

# the DTM base against the ground surface the model itself records
diff = t["ground_z"] - t["id"].map(model_z)
print(f"DTM base minus model ground: median {diff.median():+.2f} m, "
      f"within 0.5 m for {(diff.abs() <= 0.5).mean():.1%}; "
      f"{(diff < -2).sum()} solids sit more than 2 m above the terrain in the model")

# are those standing on another solid? its eaves or ridge would be at their ground height
raised = t[diff < -2]
on_roof = 0
for r in raised.itertuples():
    below = b[(b.id != r.id) & b.intersects(r.geometry.buffer(-0.2))]
    tops = pd.concat([below.ground_z + below.height_eaves_m, below.ground_z + below.height_ridge_m])
    on_roof += bool(((tops - model_z[r.id]).abs() < 1).any())
print(f"{on_roof} of {len(raised)} stand on the roof of the solid beneath them")

extruded 1,935 of 1,990; no terrain in the ring for 55; not watertight: 0
volume to the eaves 622,219 m³, to the ridge 737,951 m³ (+18.6%)
DTM base minus model ground: median +0.09 m, within 0.5 m for 86.8%; 18 solids sit more than 2 m above the terrain in the model
16 of 18 stand on the roof of the solid beneath them

Sampled around the walls, 1,935 of the 1,990 solids get a base; the 55 without one are small parts wedged between other buildings with no open ground within 3 m. Against the ground height the model itself records, the DTM base is a median 0.09 m higher and within half a metre for 86.8% of them. The 18 that sit more than 2 m above the terrain in the model are not errors: 16 of them stand exactly on the eaves or ridge of the solid beneath, so their ground surface is another building's roof. Extruded from the DTM, they would be driven down through the building beneath; a part on a roof needs its base from the model, not the terrain. To the ridge instead of the eaves, the same buildings hold 18.6% more volume.

Example 2 — extrude to a PolygonZ layer instead of a mesh

from shapely.geometry import Polygon, MultiPolygon
import geopandas as gpd

def prism_surfaces(poly: Polygon, base: float, top: float):
    """Ground, roof and wall faces as 3D polygons — a GIS-friendly 'solid'."""
    def ring3(coords, z):
        return [(x, y, z) for x, y in coords]

    holes = lambda z: [ring3(r.coords, z) for r in poly.interiors]   # courtyards stay open
    faces = [Polygon(ring3(poly.exterior.coords, base), holes(base)),
             Polygon(ring3(poly.exterior.coords, top), holes(top))]
    for ring in [poly.exterior, *poly.interiors]:
        pts = list(ring.coords)
        for a, c in zip(pts, pts[1:]):
            faces.append(Polygon([(a[0], a[1], base), (c[0], c[1], base),
                                  (c[0], c[1], top), (a[0], a[1], top)]))
    return faces

rows = []
for row in t[t.extruded].itertuples():
    for f in prism_surfaces(row.geometry, row.ground_z, row.ground_z + row.height_m):
        rows.append({"id": row.id, "geometry": f})
surfaces = gpd.GeoDataFrame(rows, crs=t.crs)
surfaces.to_file("buildings_3d.gpkg", layer="surfaces", driver="GPKG")
print(f"{len(surfaces):,} faces for {surfaces.id.nunique():,} buildings; "
      f"{sum(len(g.interiors) for g in t.geometry)} footprints have a courtyard")

13,860 faces for 1,935 buildings; 0 footprints have a courtyard

Useful when the consumer is a GIS rather than a 3D engine. Note it is a collection of faces, not a solid: nothing enforces that they close.

Example 3 — the web-map answer: do not extrude at all

import os

web = t.loc[t.extruded, ["id", "height_m", "geometry"]].to_crs(4326)
web.to_file("buildings.geojson", driver="GeoJSON", COORDINATE_PRECISION=6)

scene = trimesh.util.concatenate(meshes)
scene.export("buildings.glb")
import gzip
size = lambda f: f"{os.path.getsize(f) / 1e6:.2f} MB, {len(gzip.compress(open(f, 'rb').read())) / 1e3:.0f} kB gzipped"
print(f"GeoJSON {size('buildings.geojson')}, "
      f"{sum(len(g.exterior.coords) + sum(len(r.coords) for r in g.interiors) for g in web.geometry):,} vertices")
print(f"glTF    {size('buildings.glb')}, {len(scene.vertices):,} vertices")

GeoJSON 0.62 MB, 92 kB gzipped, 11,925 vertices
glTF    0.63 MB, 196 kB gzipped, 19,980 vertices

MapLibre's fill-extrusion layer takes its height from a property — "fill-extrusion-height": ["get", "height_m"] — and extrudes on the GPU. The payload stays 2D and the client does the work. On this tile the saving is smaller than the first version of this page claimed: raw, the GeoJSON and the binary glTF are the same size, because GeoJSON spells every coordinate out as text; gzipped for the wire the footprints are half the size of the meshes, and the mesh has 1.7 times the vertices. What the height column keeps that the glTF loses is the attributes.

Base height dominates error on slopes

Why extrude_polygon rather than a wall loop

Building the walls by hand is easy; getting the caps right is not. A polygon with interior rings needs a constrained triangulation to cap it, and the triangulation has to be consistent between the top and bottom faces or the solid leaks. trimesh.creation.extrude_polygon uses a proper triangulator and produces a watertight result including holes.

Why the base height dominates the error on sloping ground

A flat base under a building on a 5% slope across a 20 m footprint is a metre out at one end, which is more than most height estimates are wrong by. The Hague is flat, and even there the terrain in the 3 m ring round a building spans a median 0.39 m, and 1.87 m at the 95th percentile (step 2). Using the median of the terrain around the footprint splits the difference; using the minimum guarantees the building does not float but buries the uphill side.

Why extruding to the ridge overstates volume

A prism to the ridge fills the whole roof void. For a symmetric pitched roof, the prism overstates the building's volume by about half the roof volume — on a 6 m building with a 1 m roof rise, roughly 8%. Extruding to the eaves understates it by the same amount. Neither is wrong; the choice must be recorded.

Why the browser is often the right place to extrude

The same 1,935 buildings are 196 kB of gzipped glTF or 92 kB of gzipped GeoJSON with a height column (Example 3) — and the GeoJSON still carries its ids and heights, which the mesh does not. The GPU extrudes it in a frame. Server-side extrusion is for when you need the geometry itself — volumes, shadows, intersections — not for when you need a picture.

Zero heights, invalid footprints and ids

• Zero or negative heights. Clip to a minimum, or drop with a count.
• Invalid footprints. make_valid before extruding; a bowtie polygon triangulates badly.
• Self-touching rings produce degenerate triangles.
• Very small footprints can extrude to slivers; set a minimum area.
• explode before extruding. A MultiPolygon extrudes as separate prisms.
• Keep an id array. Meshes carry no attributes.
• Watch memory. A hundred thousand prisms is a large scene.
• Terrain and buildings need the same vertical datum.

Internal links

• Where building heights come from, and how wrong they are — the height you extrude to
• Extruded buildings float above or sink into the ground — the base height in practice
• 3D spatial data models explained: 2.5D, meshes and solids — what a prism can and cannot represent
• How to calculate building volumes and floor area in Python — measuring the result
• A building mesh is not watertight and its volume is wrong — when the check fails
• How to export a 3D city model to glTF in Python — writing it out
• How to publish a 3D building layer to a web map — the client-side route
• How to fix invalid geometries in GeoPandas — before extruding

FAQ

How do I extrude a polygon to 3D in Python?

trimesh.creation.extrude_polygon(polygon, height) produces a watertight prism including any interior rings; translate it in z to set the base.

What should the base height be?

The terrain at the footprint, sampled from a DTM — the median of the cells in a ring just outside the walls, because a lidar DTM has no ground under a roof. Extruding from zero puts buildings at sea level, which in central The Hague buries most of them.

Should I extrude to the eaves or the ridge?

Whichever your analysis needs, and record it. On a LoD2 tile of The Hague, two thirds of the roofs are flat and the two agree; over all 1,935 buildings the ridge prisms hold 18.6% more volume than the eaves prisms.

Why is my extruded mesh not watertight?

Usually a hand-built wall loop with no caps, or an invalid footprint. Use a proper triangulating extruder and validate the polygon first.

Do I need to extrude for a web map?

No. Ship footprints with a height property and let MapLibre's fill-extrusion do it on the GPU — the payload is a fraction of the size.

How do I keep the attributes on the mesh?

You cannot; meshes carry geometry only. Keep a parallel array of ids so the model can be joined back to the footprint table.