Generating Terrain Meshes from DEM Rasters for Digital Twins

This page turns a DEM raster — a GeoTIFF of elevations — into a watertight triangle mesh a digital twin can render, using rasterio to read the grid, numpy to build vertices, and trimesh to triangulate, simplify, and export glTF or OBJ. The path is: read the raster and its nodata mask, sample every valid cell into an XYZ vertex on an explicit CRS (here EPSG:32618 horizontal with vertical EPSG:5703), stitch the regular grid into triangles, optionally decimate to a TIN to cut triangle count, and write a .glb or .obj ready to drape imagery over.

You hit this the moment a terrain surface has to be geometry rather than a raster: a Cesium or game-engine twin needs a mesh to occlude buildings, cast shadows, and let the camera skim the ground, and analysis like cut-and-fill needs a continuous triangulated surface. The DEM already exists — produced by the digital elevation model workflows that rasterize, void-fill, and hydro-flatten a point cloud — so this guide starts from a clean GeoTIFF and ends at a mesh, and the two things that break it are mishandled nodata (which pulls spikes down to the fill value) and a CRS left implicit (which mis-scales or mis-places the surface).

DEM raster to triangulated terrain mesh A regular grid of DEM cells with elevation values becomes a lattice of vertices, each carrying its easting, northing, and height; adjacent vertices are joined into two triangles per cell to form a continuous terrain mesh. DEM cells (GeoTIFF) verticesE, N, H triangulated mesh EPSG:32618 + 5703
Each valid DEM cell becomes a vertex carrying easting, northing, and orthometric height; adjacent vertices are stitched into two triangles per cell.

Prerequisites

  • Python 3.10+ with rasterio>=1.3 (bundles GDAL 3.6+), numpy>=1.24, trimesh>=4.0, and pyproj>=3.6. Install via conda install -c conda-forge rasterio trimesh pyproj numpy so GDAL and PROJ are one matched build.
  • A single-band elevation GeoTIFF with a populated nodata value and a declared CRS. The examples use EPSG:32618 (UTM 18N, metres) for the horizontal grid and orthometric heights on EPSG:5703 (NAVD88), i.e. the compound EPSG:32618+5703. Do not mesh a geographic raster (EPSG:4326) — degree cells are not square metres and the mesh comes out sheared.
  • A DEM that is already void-filled. Interior nodata holes become gaps or spikes in the mesh; fill them at the raster stage, not here.
  • For very large rasters, enough RAM to hold the elevation array plus roughly three times its vertex-array footprint, or the windowed approach in step 5.

Step-by-Step

1. Read the DEM and its nodata mask

Open the GeoTIFF, read band 1, and capture the affine transform, the CRS, and the nodata value together — the transform maps pixel row/column to easting/northing and is what georeferences every vertex.

python
import rasterio
import numpy as np

with rasterio.open("dtm_18N_filled.tif") as src:
    elev = src.read(1).astype(np.float64)
    transform = src.transform          # affine: (col,row) -> (E,N) in EPSG:32618
    nodata = src.nodata
    crs = src.crs

print("CRS:", crs, "| nodata:", nodata, "| shape:", elev.shape)
assert crs is not None and crs.is_projected, "DEM must be a projected metric CRS (EPSG:32618)"
valid = np.isfinite(elev) & (elev != nodata)
print("valid cells:", int(valid.sum()), "of", elev.size)

2. Sample the grid into XYZ vertices in the raster CRS

Build one vertex per cell. rasterio.transform.xy converts row/column indices to eastings and northings; the height is the elevation value. Keep a lookup from (row, col) to vertex index so the triangulation in step 3 can reference neighbours.

python
from rasterio.transform import xy

rows, cols = np.nonzero(valid)
eastings, northings = xy(transform, rows, cols, offset="center")
heights = elev[rows, cols]

vertices = np.column_stack([eastings, northings, heights]).astype(np.float64)  # EPSG:32618+5703

# Dense index grid: -1 where nodata, else the vertex row in `vertices`.
vidx = np.full(elev.shape, -1, dtype=np.int64)
vidx[rows, cols] = np.arange(len(vertices))
print("vertices:", len(vertices))

3. Triangulate the regular grid into two triangles per cell

Every 2×2 block of valid vertices forms a quad that splits into two triangles. Vectorize over the grid: take the top-left corner of each cell, require all four corners valid, and emit the two triangles with consistent (counter-clockwise, upward) winding so normals point to the sky.

python
tl = vidx[:-1, :-1]     # top-left of each cell
tr = vidx[:-1, 1:]
bl = vidx[1:, :-1]
br = vidx[1:, 1:]

full = (tl >= 0) & (tr >= 0) & (bl >= 0) & (br >= 0)   # all four corners present
tl, tr, bl, br = tl[full], tr[full], bl[full], br[full]

faces = np.vstack([
    np.column_stack([tl, bl, br]),   # lower triangle
    np.column_stack([tl, br, tr]),   # upper triangle
])
print("faces:", len(faces))

The index arithmetic is mechanical, but the diagonal is not. Splitting every cell the same way is the standard approach and it is anisotropic: a ridge or a kerb that happens to run against the chosen diagonal is spanned by two triangles that average across it, losing up to half the feature’s height, while the same feature aligned with the diagonal survives intact. For terrain at metre resolution this rarely matters; for a 0.25 m grid carrying kerbs, ditch banks, and crowned roads it produces a systematic bias whose direction is fixed by your code. If the features you care about have a dominant orientation, choose the diagonal that matches it — or split each cell into four triangles around a centre vertex interpolated from the four corners, which costs twice the geometry and removes the directional bias entirely.

Two triangles per cell, and why the diagonal direction matters A regular DEM grid is triangulated by splitting every cell along one diagonal, giving two triangles per cell with vertex indices derived from row times column count plus column. When a ridge crest runs along the chosen diagonal the crest survives; when the diagonal cuts across the crest, the two triangles average it away and the ridge loses height. 0 1 2 3 v = i·ncols + j (v, v+1, v+ncols) and (v+1, v+ncols+1, v+ncols) terrain ridge crest a diagonal along the crest keeps its height the other diagonal averages it away — up to 0.4 m lost
The index arithmetic is the easy half. The diagonal is a modelling choice: a fixed direction quietly flattens every crest that runs the other way.

4. Build the mesh and (optionally) decimate to a TIN

Assemble a trimesh object, shift to a local origin so the large UTM eastings do not overflow float32 on export, and quadric-decimate to a TIN. Terrain is smooth over most of its area, so a target of 20–40% of the faces typically holds sub-decimetre vertical error while shrinking the payload.

python
import trimesh

origin = vertices.mean(axis=0)                 # local origin in EPSG:32618+5703
mesh = trimesh.Trimesh(vertices=vertices - origin, faces=faces, process=True)
trimesh.repair.fix_normals(mesh)

print("pre-decimation:", len(mesh.faces), "faces, watertight:", mesh.is_watertight)

target = max(5000, len(mesh.faces) // 4)       # ~4x reduction to a TIN
tin = mesh.simplify_quadric_decimation(target)
print("post-decimation:", len(tin.faces), "faces")
np.save("terrain_origin_utm.npy", origin)      # re-add at placement time
Regular grid sampling versus a TIN over the same profile A terrain profile sampled on a regular grid places a vertex under every cell, so most vertices land on flat ground where they describe nothing. The same profile as a triangulated irregular network keeps vertices only at breaks of slope, reproducing the shape with a fraction of the geometry. regular grid — a vertex under every cell eight of these fourteen sit on flat ground and describe nothing TIN — vertices only at breaks of slope eight vertices, the same crest and toe, and a stated maximum vertical error
Decimation to a TIN is not a quality trade so much as a redistribution: the vertices move to where the surface actually changes.

5. Reproject or export with an explicit CRS

For an OBJ handed to CAD, keep the metric EPSG:32618 vertices. For a web twin, reproject the local-origin vertices back to full UTM, then export glTF; the saved origin plus EPSG:32618+5703 is what a downstream tiler bakes into a 3D Tiles root transform. Never let an exporter guess the frame.

python
from pyproj import Transformer

world = tin.vertices + origin                  # back to EPSG:32618+5703 metres

# OBJ for CAD: metric UTM, Z-up, CRS recorded in a sidecar (.obj carries none itself).
trimesh.Trimesh(world, tin.faces).export("terrain_utm18n.obj")
with open("terrain_utm18n.prj", "w") as f:
    f.write("EPSG:32618+5703")

# glTF for the web twin: keep local metres, tag the source CRS in extras.
glb = trimesh.Trimesh(tin.vertices, tin.faces)
glb.metadata["source_crs"] = "EPSG:32618+5703"
glb.export("terrain.glb")
print("exported terrain_utm18n.obj and terrain.glb")

Expected Output & Verification

After step 5 you have terrain.glb, terrain_utm18n.obj, and the saved origin. Verify the mesh re-imports, has upward normals, and covers the raster footprint:

python
import trimesh, numpy as np

m = trimesh.load("terrain_utm18n.obj", force="mesh")
assert m.is_winding_consistent, "winding broke on export"
assert (m.face_normals[:, 2] > 0).mean() > 0.95, "normals not predominantly up — fix winding"

world = m.vertices
print(f"extent E {world[:,0].ptp():.1f} m  N {world[:,1].ptp():.1f} m  "
      f"H {world[:,2].min():.1f}–{world[:,2].max():.1f} m")

Expected console output for a 1 km² 0.5 m DTM tile decimated to a TIN:

text
CRS: EPSG:32618 | nodata: -9999.0 | shape: (2000, 2000)
valid cells: 3987214 of 4000000
vertices: 3987214
faces: 7960000
pre-decimation: 7960000 faces, watertight: False
post-decimation: 1990000 faces
extent E 999.5 m  N 999.5 m  H 4.2–61.8 m

The extent should match the raster’s ground coverage and the height range should match the DEM’s min/max — if heights are off by ~30 m, the vertical datum was dropped somewhere; confirm EPSG:5703 travelled through. A terrain mesh is rarely fully watertight at the tile edge, which is expected; stitch adjacent tiles at delivery rather than forcing a closed volume here.

One habit is worth adopting before any of these errors can occur: write the raster’s CRS, cell size, nodata value and vertex count into the mesh sidecar at the moment of export, not afterwards. Every failure below is easier to diagnose when the mesh states what it was built from.

Common Errors

Terrain has deep spikes pulling down to a flat floor. The nodata sentinel (-9999) was meshed as a real elevation because the mask in step 1 missed it — often because the raster stores nodata as NaN, not the header value, or vice versa. Test both: np.isfinite(elev) & (elev != nodata), and confirm src.nodata is populated before trusting it.

ValueError: could not broadcast or a sheared, stretched mesh. The DEM was in a geographic CRS (EPSG:4326), so rasterio.transform.xy returned degrees and the vertices span fractions of a unit while heights span tens of metres. Reproject the raster to a projected metric CRS such as EPSG:32618 before meshing; assert crs.is_projected as in step 1.

Model renders far from the origin, jitters, or the exporter warns about precision. Raw UTM eastings (hundreds of thousands of metres) were written into a float32 glTF buffer. Subtract the local origin before export (step 4) and re-add it only when placing the tile, exactly as the saved terrain_origin_utm.npy allows.

Back to Digital Elevation Model Workflows.