Source code: From DEM to 3D-printable STL with Python and NumPy

"""From DEM to 3D-printable STL with Python - TechnicalGIS tutorial.

Streams a Copernicus GLO-30 DEM tile for the Edale valley (Peak District),
reprojects it to British National Grid, scales it to a printable size and
writes a watertight binary STL using nothing but NumPy.
Run with:  python dem_to_stl.py
"""
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import rioxarray  # noqa: F401  (registers the .rio accessor)
from matplotlib.colors import LightSource
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from rasterio.enums import Resampling

HERE = Path(__file__).parent
OUT = HERE / "output"
OUT.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})

DEM_URL = ("https://copernicus-dem-30m.s3.amazonaws.com/"
           "Copernicus_DSM_COG_10_N53_00_W002_00_DEM/Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif")
AOI_BNG = (406_000, 380_000, 418_000, 389_000)   # 12 x 9 km: Kinder plateau, Edale, Mam Tor
MODEL_WIDTH_MM = 150.0       # printed size of the longest side
Z_EXAGGERATION = 2.0         # vertical exaggeration makes relief readable at this scale
BASE_MM = 4.0                # solid plinth below the lowest point


def section(title: str) -> None:
    print(f"\n=== {title} ===")


# --- Step 1: get the DEM in a projected CRS ---------------------------------------
dem = rioxarray.open_rasterio(DEM_URL, masked=True).squeeze("band", drop=True)
dem = dem.rio.clip_box(-1.95, 53.30, -1.70, 53.45).load()            # generous lon/lat window
dem = dem.rio.reproject("EPSG:27700", resolution=30, resampling=Resampling.bilinear)
dem = dem.rio.clip_box(*AOI_BNG)
z = dem.values.astype("float64")
section("DEM")
print("shape (rows, cols):", z.shape, "| resolution:", dem.rio.resolution())
print("NaNs:", int(np.isnan(z).sum()))
print(f"elevation range: {np.nanmin(z):.1f} - {np.nanmax(z):.1f} m")

# --- Step 2: scale real-world metres to millimetres -------------------------------
rows, cols = z.shape
cell_m = abs(dem.rio.resolution()[0])
width_m, height_m = (cols - 1) * cell_m, (rows - 1) * cell_m
scale = MODEL_WIDTH_MM / max(width_m, height_m)        # mm per metre
section("model scale")
print(f"ground extent: {width_m / 1000:.2f} x {height_m / 1000:.2f} km")
print(f"horizontal scale: 1:{1 / (scale / 1000):,.0f}")
print(f"cell size on the print: {cell_m * scale:.3f} mm")

x = np.arange(cols) * cell_m * scale
y = (rows - 1 - np.arange(rows)) * cell_m * scale    # row 0 is north -> largest y
X, Y = np.meshgrid(x, y)
Z = BASE_MM + (z - np.nanmin(z)) * scale * Z_EXAGGERATION
print(f"model size: {x.max():.1f} x {y.max():.1f} x {Z.max():.1f} mm (W x D x H)")


# --- Step 3: triangulate a closed (watertight) solid -------------------------------
def grid_triangles(r: int, c: int) -> np.ndarray:
    """Vertex indices for two counter-clockwise triangles per grid cell (viewed from above)."""
    idx = np.arange(r * c).reshape(r, c)
    a, b = idx[:-1, :-1].ravel(), idx[:-1, 1:].ravel()   # top-left, top-right
    d, e = idx[1:, :-1].ravel(), idx[1:, 1:].ravel()     # bottom-left, bottom-right
    # y decreases with row, so (a, d, b) and (b, d, e) are anticlockwise from above
    return np.concatenate([np.column_stack([a, d, b]), np.column_stack([b, d, e])])


top_v = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
bot_v = top_v.copy()
bot_v[:, 2] = 0.0
verts = np.vstack([top_v, bot_v])
n = rows * cols
tri_top = grid_triangles(rows, cols)
tri_bottom = tri_top[:, ::-1] + n                         # reversed winding -> faces down

# Walls: walk the outline clockwise and join each top edge to the bottom.
idx = np.arange(n).reshape(rows, cols)
ring = np.concatenate([idx[0, :], idx[1:, -1], idx[-1, -2::-1], idx[-2:0:-1, 0]])
p, q = ring, np.roll(ring, -1)
walls = np.concatenate([np.column_stack([p, q, p + n]), np.column_stack([q, q + n, p + n])])
triangles = np.vstack([tri_top, tri_bottom, walls])
section("mesh")
print(f"vertices: {len(verts):,}  triangles: {len(triangles):,} "
      f"(surface {len(tri_top):,}, base {len(tri_bottom):,}, walls {len(walls):,})")

# Watertight check: every undirected edge must be shared by exactly two triangles,
# and every directed edge must appear exactly once (consistent orientation).
edges = np.concatenate([triangles[:, [0, 1]], triangles[:, [1, 2]], triangles[:, [2, 0]]])
_, und_counts = np.unique(np.sort(edges, axis=1), axis=0, return_counts=True)
_, dir_counts = np.unique(edges, axis=0, return_counts=True)
print("watertight:", bool((und_counts == 2).all()), "| consistently oriented:", bool((dir_counts == 1).all()))

tri_xyz = verts[triangles]                                 # (n_tri, 3 vertices, xyz)
normals = np.cross(tri_xyz[:, 1] - tri_xyz[:, 0], tri_xyz[:, 2] - tri_xyz[:, 0])
volume_mm3 = np.einsum("ij,ij->i", tri_xyz[:, 0], np.cross(tri_xyz[:, 1], tri_xyz[:, 2])).sum() / 6
print(f"enclosed volume: {volume_mm3 / 1000:.1f} cm³ (positive = outward-facing normals)")

# --- Step 4: write a binary STL with a NumPy structured array -----------------------
lengths = np.linalg.norm(normals, axis=1, keepdims=True)
unit = np.divide(normals, lengths, out=np.zeros_like(normals), where=lengths > 0)
record = np.dtype([("normal", "<f4", 3), ("v", "<f4", (3, 3)), ("attr", "<u2")])
data = np.zeros(len(triangles), dtype=record)
data["normal"], data["v"] = unit, tri_xyz
stl_path = OUT / "edale_kinder_mam_tor.stl"
with open(stl_path, "wb") as f:
    f.write(b"TechnicalGIS Edale DEM - Copernicus GLO-30".ljust(80, b" "))
    f.write(np.uint32(len(data)).tobytes())
    f.write(data.tobytes())
section("STL")
print(f"wrote {stl_path.name}: {stl_path.stat().st_size / 1e6:.1f} MB, {len(data):,} facets")

# Read it back to prove the file is valid (and to draw the preview from the file itself).
raw = np.fromfile(stl_path, dtype=np.uint8)
count = int(np.frombuffer(raw[80:84].tobytes(), "<u4")[0])
back = np.frombuffer(raw[84:].tobytes(), dtype=record)
assert count == len(back) == len(triangles)
print("read back:", count, "facets | bounding box (mm):",
      np.round(back["v"].reshape(-1, 3).min(0), 1), "->", np.round(back["v"].reshape(-1, 3).max(0), 1))

# --- Step 5: previews ------------------------------------------------------------------
ls = LightSource(azdeg=315, altdeg=45)
fig, ax = plt.subplots(figsize=(7.5, 5.8))
rgb = ls.shade(z, cmap=plt.cm.gist_earth, blend_mode="soft", vert_exag=2, dx=cell_m, dy=cell_m,
               vmin=np.nanmin(z) - 150, vmax=np.nanmax(z) + 50)
ax.imshow(rgb, extent=AOI_BNG[0::2] + AOI_BNG[1::2])
labels = {"Mam Tor": (412_750, 383_640), "Kinder Scout": (408_600, 387_600), "Edale": (412_300, 385_700),
          "Lose Hill": (415_300, 385_300)}
for name, (lx, ly) in labels.items():
    ax.plot(lx, ly, "^", color="black", ms=5)
    ax.annotate(name, (lx, ly), xytext=(5, 5), textcoords="offset points", fontsize=8,
                bbox=dict(boxstyle="round,pad=0.2", fc="white", alpha=0.75, lw=0))
ax.set_title("Copernicus GLO-30 DEM, Edale valley (British National Grid, 30 m)")
ax.ticklabel_format(style="plain")
fig.savefig(OUT / "01-dem-aoi.png")
plt.close(fig)

# Render a decimated copy of the surface (every 3rd vertex) - plotting ~100k facets
# in matplotlib is slow, and the full-resolution file has already been validated above.
step = 3
sub = grid_triangles(len(range(0, rows, step)), len(range(0, cols, step)))
Xs, Ys, Zs = X[::step, ::step], Y[::step, ::step], Z[::step, ::step]
pts = np.column_stack([Xs.ravel(), Ys.ravel(), Zs.ravel()])
facet_xyz = pts[sub]
facet_n = np.cross(facet_xyz[:, 1] - facet_xyz[:, 0], facet_xyz[:, 2] - facet_xyz[:, 0])
fig = plt.figure(figsize=(9, 6.5))
ax = fig.add_subplot(projection="3d")
height = (facet_xyz[:, :, 2].mean(1) - BASE_MM) / (Z.max() - BASE_MM)
colours = plt.cm.gist_earth(0.25 + 0.6 * height)
colours[:, :3] *= (0.45 + 0.55 * ls.shade_normals(facet_n, fraction=1.0))[:, None]
ax.add_collection3d(Poly3DCollection(facet_xyz, linewidths=0, facecolors=colours))
for wall in (np.column_stack([x, np.full(cols, y.max()), Z[0]]), np.column_stack([x, np.zeros(cols), Z[-1]])):
    ax.add_collection3d(Poly3DCollection([np.vstack([wall, [[wall[-1, 0], wall[-1, 1], 0], [wall[0, 0], wall[0, 1], 0]]])],
                                         facecolors="#c9b79c", linewidths=0))
ax.set_xlim(0, x.max())
ax.set_ylim(0, y.max())
ax.set_zlim(0, Z.max())
ax.set_box_aspect((x.max(), y.max(), Z.max()))  # true printed proportions
ax.set_zticks([0, BASE_MM, round(Z.max())])
for axis in (ax.xaxis, ax.yaxis, ax.zaxis):
    axis.pane.set_visible(False)
ax.view_init(elev=28, azim=-60)
ax.set_xlabel("mm")
ax.set_ylabel("mm")
ax.set_zlabel("mm")
ax.set_title(f"STL preview: {x.max():.0f} x {y.max():.0f} mm, {Z_EXAGGERATION}x vertical exaggeration")
fig.savefig(OUT / "02-stl-preview.png")
plt.close(fig)

# Cross-section through Mam Tor showing the plinth and exaggeration.
r = int(round((AOI_BNG[3] - 383_640) / cell_m))
fig, ax = plt.subplots(figsize=(8, 2.6))
ax.fill_between(x, 0, Z[r], color="#c9b79c")
ax.plot(x, Z[r], color="#7f5539", lw=1)
ax.axhline(BASE_MM, color="#495057", ls="--", lw=0.8)
ax.text(2, BASE_MM + 0.6, f"{BASE_MM:.0f} mm plinth", fontsize=8)
ax.set_aspect("equal")
ax.set_xlabel("mm")
ax.set_ylabel("mm")
ax.set_title("Model cross-section (west → east) through Mam Tor")
fig.savefig(OUT / "03-cross-section.png")
plt.close(fig)
print("\nfigures written to", OUT)