Extruded buildings float above or sink into the ground

Problem statement

The extrusion worked, the heights are right, and the buildings hover ten metres above the terrain — or the ground floor of every one of them is buried. On a slope, half the block is in the air and the other half is underground.

There are four distinct causes and they need different fixes: the base height was never set, the base came from the wrong vertical datum, the base is a single value for a sloping footprint, or the model and the terrain are in different height systems entirely. All four produce a picture that is obviously wrong and gives no hint which.

Quick answer

Diagnose by the size and sign of the offset:

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

with rasterio.open("dtm.tif") as src:
    terrain = []
    for geom in b.geometry:
        # a ring just outside the walls: a lidar DTM has no ground under a roof
        arr, _ = mask(src, [geom.buffer(3).difference(geom)], crop=True, filled=False)
        v = arr.compressed()
        terrain.append(np.median(v) if v.size else np.nan)

b["terrain_z"] = terrain
b["offset_m"] = b["base_z"] - b["terrain_z"]
print(b["offset_m"].describe().round(2))

count    1935.00
mean       -0.09
std         1.09
min        -4.54
25%        -0.19
50%        -0.09
75%        -0.03
max        13.95
Name: offset_m, dtype: float64

Base set, vertical datum, terrain under the footprint

1. Check whether a base was set at all

trimesh.creation.extrude_polygon produces a prism from z = 0 to z = height. If nothing translated it, every building sits at the vertical origin — sea level for a projected CRS with orthometric heights, or the ellipsoid for a geographic one.

import trimesh

meshes = [trimesh.creation.extrude_polygon(g, h) for g, h in zip(b.geometry, b.height_m)]
print([round(m.bounds[0][2], 2) for m in meshes[:10]])   # all 0.0 means no base was applied

[np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0)]

2. Check the vertical datum

A constant offset of tens of metres is a geoid separation. A model whose z values are heights above a national datum, drawn against a terrain in ellipsoidal height, sits that far underground.

import pyproj
from pyproj import Transformer

t = Transformer.from_crs("EPSG:7415", "EPSG:4979", always_xy=True)   # RD + NAP → ETRS89 ellipsoidal
print("no geoid grid installed:", round(t.transform(78_642, 457_940, 0.0)[2], 2))

pyproj.network.set_network_enabled(True)                              # fetch grids from cdn.proj.org
t = Transformer.from_crs("EPSG:7415", "EPSG:4979", always_xy=True)
print("with the NLGEO2018 grid:", round(t.transform(78_642, 457_940, 0.0)[2], 2))

egm = Transformer.from_crs("EPSG:9518", "EPSG:4979", always_xy=True)  # WGS84 + EGM2008 height
places = {"The Hague": (4.30, 52.08), "London": (-0.13, 51.50), "Edinburgh": (-3.19, 55.95),
          "Paris": (2.35, 48.86), "Madrid": (-3.70, 40.40), "Rome": (12.50, 41.90),
          "Athens": (23.70, 37.98), "Stockholm": (18.07, 59.33), "Helsinki": (24.94, 60.17),
          "south of Sri Lanka": (80.0, 5.0)}
print({k: round(egm.transform(x, y, 0.0)[2], 1) for k, (x, y) in places.items()})

no geoid grid installed: 0.0
with the NLGEO2018 grid: 43.4
{'The Hague': 43.6, 'London': 46.1, 'Edinburgh': 53.0, 'Paris': 44.6, 'Madrid': 51.7, 'Rome': 48.4, 'Athens': 38.5, 'Stockholm': 23.2, 'Helsinki': 17.8, 'south of Sri Lanka': -105.7}

Without the grid, PROJ falls back to a transformation that ignores the vertical part and returns the height unchanged — 0 m in, 0 m out, no warning. The page's first version printed exactly this line and would have shown 0.0. With the grid, the separation here is 43.40 m. Across Europe (EGM2008, a global geoid, so to a few decimetres) it is 46 m in London and 53 m in Edinburgh, 38–52 m round the Mediterranean, but only 18–23 m in Stockholm and Helsinki — and −106 m south of Sri Lanka.

3. Sample the terrain under each footprint, not at its centroid

A centroid can fall in a courtyard, on the far side of a river, or outside the raster. Masking by the polygon is more robust — but a DTM made from lidar is empty under the roof, so mask a ring just outside the walls. On this tile the cells inside the footprints leave 211 buildings without a value even with a one-cell buffer (Example 1); a 3 m ring leaves 55.

4. Choose the statistic deliberately

• Minimum — the building never floats; the uphill side is buried.
• Median — split the difference; both sides are out by half the slope.
• Per-vertex — a sloping base, which needs a mesh rather than a prism.

For most work the median is right. For rendering, the minimum looks better because a floating building is more obviously wrong than a buried one.

5. Handle the sloping case properly if it matters

On steep ground, a flat base over a 20 m footprint at 5% is a metre out at each end. The honest fix is to extrude from a sloping base — build the prism with per-vertex base heights sampled from the DTM.

6. Check the terrain and the model share a CRS

Including the vertical part. EPSG:7415 is not EPSG:28992: the first carries NAP heights, the second is horizontal only, and reprojecting between them silently drops or keeps z depending on the tool.

7. In a web map, remember the base is relative to the terrain

With a terrain source enabled, fill-extrusion-base is measured from the terrain surface, so 0 is correct. Copying a base from the flat case buries every building.

Robust terrain sampling and datum 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 base_z, the height of the solid's ground surface in the model (metres above NAP) — the base the city itself gave each building. 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"], "base_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.base_z.min():.2f}–{b.base_z.max():.2f} m NAP, "
      f"median {b.base_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 — robust terrain sampling with a report

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

def terrain_under(buildings, dtm_path, stat="median", min_cells=4):
    stats, problems = [], []
    with rasterio.open(dtm_path) as src:
        if buildings.crs != src.crs:
            raise ValueError(f"CRS mismatch: buildings {buildings.crs}, dtm {src.crs}")
        for i, geom in enumerate(buildings.geometry):
            try:
                arr, _ = mask(src, [geom], crop=True, filled=False)
            except ValueError:
                problems.append((i, "outside the raster"))
                stats.append(np.nan)
                continue
            v = arr.compressed()
            if v.size < min_cells:
                # fall back to a small buffer so slivers still get a value
                arr, _ = mask(src, [geom.buffer(src.res[0])], crop=True, filled=False)
                v = arr.compressed()
            if v.size == 0:
                problems.append((i, "all nodata"))
                stats.append(np.nan)
            else:
                stats.append(float(np.median(v) if stat == "median" else v.min()))
    return np.array(stats), problems

base, problems = terrain_under(b, "dtm.tif")
print(f"{np.isnan(base).sum()} buildings without a terrain value; "
      f"{len(problems)} reported problems")
d = base - b["base_z"]
print(f"against the model's own base: median {np.nanmedian(d):+.2f} m, "
      f"within 0.5 m for {(np.abs(d) <= 0.5).sum()} of {np.isfinite(d).sum()}")

211 buildings without a terrain value; 211 reported problems
against the model's own base: median +0.14 m, within 0.5 m for 1551 of 1779

Sampling inside the footprint, with a one-cell buffer as the fallback, leaves 211 of the 1,990 buildings with no terrain value: the cells a lidar DTM has near a wall are the few it kept at the edge of the hole. Where it does find a value, it agrees with the model's own base to half a metre for 1,551 of 1,779. Without the buffer it is 545 (the extrusion guide counts them); the 3 m ring in the quick answer leaves 55.

Example 2 — detect a datum offset rather than guessing

import numpy as np

def diagnose(base_z, terrain_z, label):
    base_z, terrain_z = np.asarray(base_z, float), np.asarray(terrain_z, float)
    diff = base_z - terrain_z
    med = np.nanmedian(diff)
    q1, q3 = np.nanpercentile(diff, [25, 75])
    far = np.abs(diff - med) > 2
    print(f"{label}: offset median {med:.2f} m, std {np.nanstd(diff):.2f}, IQR {q3 - q1:.2f}, "
          f"{far.sum()} more than 2 m from the median")
    if np.all(base_z == 0):
        print("  every base is 0 — the base was never set")
    elif abs(med) > 20 and q3 - q1 < 2:
        print("  constant large offset — a vertical datum mismatch, not a sampling problem")
    elif q3 - q1 > 1 and abs(med) < 1:
        print("  offset varies and averages to zero — terrain slope under a flat base")
    elif far.any():
        print(f"  consistent apart from {far.sum()} buildings — look at those one by one")

diagnose(b["base_z"], b["terrain_z"], "model vs AHN")
diagnose(np.zeros(len(b)), b["terrain_z"], "base never set")
# the same AHN terrain, had it been delivered in ellipsoidal heights
x, y = b.centroid.x.values, b.centroid.y.values
ell = np.array(t.transform(x, y, b["terrain_z"].values)[2])
diagnose(b["base_z"], ell, "terrain in ellipsoidal heights")

model vs AHN: offset median -0.09 m, std 1.09, IQR 0.15, 43 more than 2 m from the median
  consistent apart from 43 buildings — look at those one by one
base never set: offset median -6.10 m, std 2.55, IQR 2.64, 709 more than 2 m from the median
  every base is 0 — the base was never set
terrain in ellipsoidal heights: offset median -43.49 m, std 1.09, IQR 0.15, 43 more than 2 m from the median
  constant large offset — a vertical datum mismatch, not a sampling problem

The first version of this check keyed on the standard deviation, and on the real model it called a perfectly consistent layer "terrain slope under a flat base": 18 rooftop parts, which stand on other buildings 13 m above the ground, pushed the standard deviation to 1.09 m while the interquartile range is 0.15 m. It also said nothing at all about bases that were never set. The IQR separates the cases, and the stragglers are listed rather than averaged in. Converting the same terrain to ellipsoidal heights with the NLGEO2018 grid moves the median to −43.49 m and leaves the IQR untouched — the signature of a datum problem.

Example 3 — a sloping base

The first version of this example sampled the DTM at each footprint vertex and built the prism by hand. On AHN that cannot work, before the geometry is even reached:

with rasterio.open("dtm.tif") as src:
    pts = np.vstack([np.asarray(g.exterior.coords) for g in b.geometry])
    vals = np.array([v[0] for v in src.sample(pts)])
    print(f"{(vals == src.nodata).mean():.0%} of {len(pts):,} footprint vertices sample a nodata cell "
          f"({src.nodata:.3g})")

69% of 12,240 footprint vertices sample a nodata cell (3.4e+38)

A vertex sits on the wall, and the wall is where the DTM's hole under the building begins. The nodata value then becomes a base height of 3.4 × 10³⁸ m. The hand-built prism had a second bug: it took the cap triangles from triangulate_polygon and used their indices as if they numbered the footprint's own vertices, which they do not — on this tile it raised IndexError: index 8 is out of bounds for axis 0 with size 8. Filling the DTM under the buildings first, and bending a prism that extrude_polygon has already made watertight, avoids both:

import numpy as np, trimesh, rasterio
from rasterio.fill import fillnodata
from rasterio.transform import rowcol
from shapely.geometry import Polygon

with rasterio.open("dtm.tif") as src:
    dtm = src.read(1, masked=True)
    ground = fillnodata(dtm.filled(0), mask=~dtm.mask, max_search_distance=100)   # interpolate under roofs
    ground_transform = src.transform

def extrude_on_terrain(poly: Polygon, height, ground, transform):
    """A prism whose base follows the terrain at each vertex; the topology is extrude_polygon's."""
    m = trimesh.creation.extrude_polygon(poly, height)
    v = m.vertices.copy()
    rows, cols = rowcol(transform, v[:, 0], v[:, 1])
    v[:, 2] += ground[rows, cols]              # bottom at the ground, top at ground + height
    return trimesh.Trimesh(v, m.faces, process=False)

sloped = [extrude_on_terrain(g, h, ground, ground_transform) for g, h in zip(b.geometry, b.height_m)]
tilt = np.array([np.ptp(m.vertices[m.vertices[:, 2] < m.vertices[:, 2].mean(), 2]) for m in sloped])
print(f"{sum(m.is_watertight for m in sloped)} of {len(sloped)} watertight; "
      f"base varies across the footprint by a median {np.median(tilt):.2f} m, "
      f"95th percentile {np.percentile(tilt, 95):.2f} m, max {tilt.max():.2f} m")

1990 of 1990 watertight; base varies across the footprint by a median 0.17 m, 95th percentile 1.15 m, max 4.53 m

Moving vertices vertically never changes which vertices a face joins, so a solid that was watertight stays watertight — there is no winding to get wrong.

Geoid separation and centroid sampling

Why the geoid separation is the commonest large offset

National height datums are defined relative to mean sea level, which follows the geoid. GNSS, and lidar or photogrammetry processed straight from it, gives heights above the ellipsoid; global DEMs such as Copernicus GLO-30 are referenced to a geoid (EGM2008) instead. The two differ by the geoid separation — 43.4 m at The Hague, 38–53 m across most of western and southern Europe, 18–23 m round the northern Baltic, and −106 m south of Sri Lanka (step 2). Nothing in a file forces the two to be labelled, so the mismatch is invisible until buildings sink.

Why the centroid is a bad sampling point

A building's centroid can be outside its own footprint when the shape is concave, which for a U-shaped block means sampling the courtyard. It can also be a nodata cell in a DTM that masks buildings out — and in AHN that is every cell under a roof. A statistic over a ring just outside the walls is barely more code and avoids both.

Why a flat base is a decision and not a bug

Buildings have level floors; the ground does not. A real building on a slope has a plinth, a stepped foundation or a sunken side, and none of that is in a footprint layer. A flat base is therefore a reasonable model — but which flat base matters, and the minimum, median and mean answer different questions.

Why web maps invert the rule

Renderers with terrain support place extrusions on the terrain surface, so the base is an offset from the ground rather than an absolute height. That is the opposite convention from a desktop 3D scene, and code moved between the two without a change buries or floats everything.

Basements, DTM versus DSM and bridges

• Basements need a negative base. base = terrain − depth.
• DTM versus DSM. Sampling a DSM under a building gives the roof.
• Some DTMs mask buildings as nodata. Buffer outwards to find ground.
• Bridges and buildings over water have no terrain beneath them.
• Check the CRS on both layers, including the vertical component.
• Terrain exaggeration in viewers is display only.
• Report the buildings you could not place. NaN is better than zero.
• Record the statistic used. Minimum and median give different volumes.

Internal links

• How to extrude building footprints into 3D in Python — where the base is set
• Vertical datums explained — the constant-offset case
• 3D spatial data models explained: 2.5D, meshes and solids — why z needs its own CRS
• How to publish a 3D building layer to a web map — the inverted base convention
• How to extract raster values at points with rasterio — sampling the DTM
• LiDAR surfaces explained: DSM, DTM and CHM — sampling the right surface
• How to do a datum shift transformation in Python — converting between height systems
• Raster and vector data do not align — the horizontal version of the same problem

FAQ

Why are my extruded buildings floating above the terrain?

Most often no base height was applied, so the prisms start at z = 0, or the model and the terrain use different vertical datums.

What is a 40 m constant offset?

A geoid separation. National height datums and ellipsoidal heights differ by tens of metres — about 43 m in the Netherlands.

Should I use the minimum or the median terrain height?

The median splits the error on a slope; the minimum guarantees the building does not float and buries the uphill side. Pick one and record it.

Why does sampling the centroid fail?

A concave footprint's centroid can lie outside it, and lidar DTMs such as AHN leave every cell under a roof as nodata. Take a statistic over a ring just outside the walls.

How do I handle a steep site?

Extrude from a sloping base: fill the DTM under the buildings, extrude the footprint as usual, then lift each vertex by the terrain height beneath it. Moving vertices vertically keeps the solid watertight.

Why do my web map buildings sink into a hill?

With a terrain source enabled, fill-extrusion-base is relative to the terrain. A base copied from the flat case is subtracted twice.