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:
arr, _ = mask(src, [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))
| offset | cause |
|---|---|
| equal to the base_z, terrain near 0 | base never set โ buildings at z = 0 |
| a constant 30โ60 m | vertical datum mismatch (geoid separation) |
| varies with slope, mean near 0 | one base per footprint on sloping ground |
| a constant few metres | different DTM or a different reference epoch |
Step-by-step solution
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.
print([round(m.bounds[0][2], 2) for m in meshes[:10]]) # all 0.0 means no base was applied
2. Check the vertical datum
A constant offset of tens of metres is a geoid separation. The Netherlands is about 43 m; Britain about 45โ50 m; parts of the Mediterranean and the Baltic differ by similar amounts. A model whose z values are heights above a national datum, drawn against a terrain in ellipsoidal height, sits that far underground.
from pyproj import Transformer
t = Transformer.from_crs("EPSG:7415", "EPSG:4979", always_xy=True) # NAP โ ellipsoidal
print(t.transform(78_642, 457_940, 0.0)[2])
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 and taking a statistic of the cells inside is both more robust and more honest.
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.
Code examples
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")
Example 2 โ detect a datum offset rather than guessing
import numpy as np
diff = b["model_base_z"] - b["terrain_z"]
print(f"offset: median {np.nanmedian(diff):.2f} m, "
f"std {np.nanstd(diff):.2f} m, range {np.nanmin(diff):.1f} to {np.nanmax(diff):.1f}")
if abs(np.nanmedian(diff)) > 20 and np.nanstd(diff) < 2:
print("constant large offset โ this is a vertical datum mismatch, not a sampling problem")
elif np.nanstd(diff) > 1 and abs(np.nanmedian(diff)) < 1:
print("offset varies and averages to zero โ this is terrain slope under a flat base")
A small standard deviation with a large median is the signature of a datum problem; the reverse is a slope problem. Printing both saves an hour of reading code.
Example 3 โ a sloping base
import numpy as np, trimesh, rasterio
from shapely.geometry import Polygon
def extrude_on_terrain(poly: Polygon, height, dtm_src):
"""Base follows the terrain under each footprint vertex."""
coords = np.array(poly.exterior.coords[:-1])
base_z = np.array([v[0] for v in dtm_src.sample(coords)], dtype=float)
n = len(coords)
verts = np.vstack([np.c_[coords, base_z],
np.c_[coords, base_z + height]])
faces = []
for i in range(n):
j = (i + 1) % n
faces += [[i, j, j + n], [i, j + n, i + n]]
cap = trimesh.creation.triangulate_polygon(poly)[1]
faces += [[a, b, c] for a, b, c in cap] # bottom
faces += [[a + n, c + n, b + n] for a, b, c in cap] # top, reversed
return trimesh.Trimesh(vertices=verts, faces=np.array(faces), process=True)
The reversed winding on the top cap is what keeps the solid watertight; without it the normals disagree and the volume is wrong.
Explanation
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 most global terrain products use heights above the WGS84 ellipsoid. The two differ by the geoid separation โ roughly 43 m in the Netherlands, 45โ55 m over much of western Europe, and negative in parts of the Indian Ocean. 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. Masking by the polygon and taking a statistic over the cells inside is barely more code and cannot do either.
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.
Edge cases or notes
- 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 some DTMs mask buildings as nodata. Mask by the polygon and take a statistic of the cells inside.
How do I handle a steep site?
Extrude from a sloping base, sampling the DTM at each footprint vertex, and reverse the winding on the top cap so the solid stays 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.