Source code: Emboss a GPX route on 3D-printed terrain with Python: raised, engraved or two-colour inlay

"""Emboss a GPX route on a 3D-printable terrain model: the Pennine Way from Edale to Crowden.

Reads a GPX track, drapes it over Environment Agency LIDAR terrain, and writes three printable
solids: the route raised 0.6 mm above the terrain, the route engraved 0.8 mm into it, and a
two-part inlay (terrain plus route) for a multi-material printer. Every mesh is checked for
watertightness and orientation before it is written as a binary STL, and the inlay is also
written as one 3MF file with two parts.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/gpx-route-3d-print-terrain/

Data: Environment Agency LIDAR Composite DTM 1 m via the Defra WCS (Open Government Licence v3);
the Pennine Way (Edale to Crowden) route relation from OpenStreetMap
(© OpenStreetMap contributors, Open Database Licence).
"""
import sys
import time
import zipfile
import xml.etree.ElementTree as ET
from importlib.metadata import version
from pathlib import Path

import numpy as np
import requests
import rasterio
from pyproj import Transformer
from scipy import ndimage
from shapely.geometry import LineString, box
from shapely.geometry.polygon import orient
from shapely.ops import linemerge
import manifold3d as mf
import pyvista as pv
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource

plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9,
                     "figure.facecolor": "white", "axes.facecolor": "white"})
pv.OFF_SCREEN = True

DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
UA = {"User-Agent": "TechnicalGIS-tutorial (https://www.technicalgis.com)"}

RELATION = 4080403                                      # OSM: "Pennine Way (Edale to Crowden)"
AOI = (403_000, 384_000, 415_000, 400_000)              # 12 x 16 km, British National Grid
RES = 25.0                                              # metres per cell in the model
FETCH = 5.0                                             # metres per cell requested from the WCS
TILE = 4000                                             # metres per WCS request
WCS = ("https://environment.data.gov.uk/spatialdata/lidar-composite-digital-terrain-model-dtm-1m/wcs"
       "?service=WCS&version=2.0.1&request=GetCoverage&format=image/tiff"
       "&CoverageId=13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m")

MODEL_LONG_MM = 170.0       # printed length of the longer (north-south) side
Z_EXAGGERATION = 2.0
BASE_MM = 4.0
ROUTE_WIDTH_MM = 1.2        # three 0.4 mm extrusion lines
RAISE_MM = 0.6              # route height above the terrain (three 0.2 mm layers)
DEPTH_MM = 0.8              # groove depth below the terrain

RED, SAND = "#d53e4f", "#e8dcc0"


def section(title):
    print(f"\n=== {title} ===")


def get(url, path, tries=6):
    """Download once and cache. The Defra WCS rate-limits bursts (HTTP 429), so back off and retry."""
    if not path.exists():
        for attempt in range(tries):
            try:
                r = requests.get(url, timeout=300, headers=UA)
                r.raise_for_status()
                if not r.headers.get("Content-Type", "").startswith("image/"):
                    raise requests.RequestException(f"not an image: {r.text[:200]}")
                break
            except requests.RequestException as err:
                if attempt == tries - 1:
                    raise
                wait = getattr(err.response, "headers", {}).get("Retry-After")
                time.sleep(float(wait) if wait and wait.isdigit() else 20 * (attempt + 1))
        path.write_bytes(r.content)
        time.sleep(1)
    return path


# ------------------------------------------------------------------ mesh helpers
def grid_triangles(r, c):
    """Two anticlockwise (from above) triangles per grid cell, as in the DEM-to-STL tutorial."""
    idx = np.arange(r * c).reshape(r, c)
    a, b = idx[:-1, :-1].ravel(), idx[:-1, 1:].ravel()
    d, e = idx[1:, :-1].ravel(), idx[1:, 1:].ravel()
    return np.concatenate([np.column_stack([a, d, b]), np.column_stack([b, d, e])])


def heightfield_solid(Z, cell_mm):
    """Closed solid: the surface Z (mm), a flat base at z = 0 and four side walls."""
    rows, cols = Z.shape
    X, Y = np.meshgrid(np.arange(cols) * cell_mm, (rows - 1 - np.arange(rows)) * cell_mm)
    top = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
    bot = top.copy()
    bot[:, 2] = 0.0
    n = rows * cols
    tri_top = grid_triangles(rows, cols)
    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])])
    return np.vstack([top, bot]), np.vstack([tri_top, tri_top[:, ::-1] + n, walls])


def to_manifold(verts, tris):
    m = mf.Manifold(mf.Mesh(vert_properties=np.ascontiguousarray(verts, dtype=np.float32),
                            tri_verts=np.ascontiguousarray(tris, dtype=np.uint32)))
    assert m.status() == mf.Error.NoError, m.status()
    return m


def from_manifold(m):
    mesh = m.to_mesh()
    return np.asarray(mesh.vert_properties)[:, :3].astype("float64"), np.asarray(mesh.tri_verts).astype("int64")


def check_mesh(name, verts, tris):
    """Watertight: every undirected edge used twice. Oriented: every directed edge used once."""
    e = np.concatenate([tris[:, [0, 1]], tris[:, [1, 2]], tris[:, [2, 0]]])
    n = len(verts)
    und = np.unique(np.sort(e, axis=1) @ np.array([n, 1]), return_counts=True)[1]
    drc = np.unique(e @ np.array([n, 1]), return_counts=True)[1]
    t = verts[tris]
    vol = np.einsum("ij,ij->i", t[:, 0], np.cross(t[:, 1], t[:, 2])).sum() / 6
    ok = bool((und == 2).all() and (drc == 1).all() and vol > 0)
    print(f"{name:<12s} {len(tris):>9,} triangles | watertight {bool((und == 2).all())} | "
          f"oriented {bool((drc == 1).all())} | volume {vol / 1000:7.2f} cm3")
    return ok, vol


def write_stl(path, verts, tris, header):
    t = verts[tris]
    nrm = np.cross(t[:, 1] - t[:, 0], t[:, 2] - t[:, 0])
    ln = np.linalg.norm(nrm, axis=1, keepdims=True)
    rec = np.zeros(len(tris), dtype=[("normal", "<f4", 3), ("v", "<f4", (3, 3)), ("attr", "<u2")])
    rec["normal"] = np.divide(nrm, ln, out=np.zeros_like(nrm), where=ln > 0)
    rec["v"] = t
    with open(path, "wb") as f:
        f.write(header.encode("ascii")[:80].ljust(80, b" "))
        f.write(np.uint32(len(rec)).tobytes())
        f.write(rec.tobytes())


def write_3mf(path, parts):
    """Minimal 3MF: one mesh object per part, each with its own display colour."""
    ns = "http://schemas.microsoft.com/3dmanufacturing/core/2015/02"
    xml = [f'<?xml version="1.0" encoding="UTF-8"?>\n<model unit="millimeter" xml:lang="en-GB" xmlns="{ns}">',
           "<resources>", '<basematerials id="1">']
    xml += [f'<base name="{name}" displaycolor="{colour}FF"/>' for name, _, _, colour in parts]
    xml.append("</basematerials>")
    for i, (name, v, t, _) in enumerate(parts):
        xml.append(f'<object id="{i + 2}" name="{name}" type="model" pid="1" pindex="{i}"><mesh><vertices>')
        xml.append("".join(f'<vertex x="{x:.4f}" y="{y:.4f}" z="{z:.4f}"/>' for x, y, z in v))
        xml.append("</vertices><triangles>")
        xml.append("".join(f'<triangle v1="{a}" v2="{b}" v3="{c}"/>' for a, b, c in t))
        xml.append("</triangles></mesh></object>")
    xml += ["</resources>", "<build>"] + [f'<item objectid="{i + 2}"/>' for i in range(len(parts))]
    xml += ["</build>", "</model>"]
    with zipfile.ZipFile(path, "w", zipfile.ZIP_DEFLATED) as z:
        z.writestr("[Content_Types].xml",
                   '<?xml version="1.0" encoding="UTF-8"?><Types xmlns="http://schemas.openxmlformats.org/'
                   'package/2006/content-types"><Default Extension="rels" ContentType="application/'
                   'vnd.openxmlformats-package.relationships+xml"/><Default Extension="model" '
                   'ContentType="application/vnd.ms-package.3dmanufacturing-3dmodel+xml"/></Types>')
        z.writestr("_rels/.rels",
                   '<?xml version="1.0" encoding="UTF-8"?><Relationships xmlns="http://schemas.openxmlformats.org/'
                   'package/2006/relationships"><Relationship Target="/3D/3dmodel.model" Id="rel0" '
                   'Type="http://schemas.microsoft.com/3dmanufacturing/2013/01/3dmodel"/></Relationships>')
        z.writestr("3D/3dmodel.model", "\n".join(xml))


def read_gpx(path):
    """Return an (n, 2) array of lon/lat from the track points (or route points) of any GPX 1.0/1.1 file."""
    root = ET.parse(path).getroot()
    ns = root.tag[1:].split("}")[0] if root.tag.startswith("{") else ""
    q = (lambda t: f"{{{ns}}}{t}") if ns else (lambda t: t)
    pts = list(root.iter(q("trkpt"))) or list(root.iter(q("rtept")))
    return np.array([[float(p.get("lon")), float(p.get("lat"))] for p in pts])


section("versions")
print(f"python {sys.version.split()[0]} | numpy {np.__version__} | rasterio {rasterio.__version__} | "
      f"shapely {version('shapely')} | manifold3d {version('manifold3d')} | pyvista {pv.__version__}")

# ------------------------------------------------------------------ Step 1: a GPX route
section("Step 1: the route as GPX")
osm = requests.get(f"https://api.openstreetmap.org/api/0.6/relation/{RELATION}/full.json",
                   timeout=180, headers=UA).json()["elements"]
nodes = {e["id"]: (e["lon"], e["lat"]) for e in osm if e["type"] == "node"}
ways = {e["id"]: e["nodes"] for e in osm if e["type"] == "way"}
rel = next(e for e in osm if e["type"] == "relation")
line = linemerge([LineString([nodes[n] for n in ways[m["ref"]]]) for m in rel["members"] if m["type"] == "way"])
print(f"OSM relation {RELATION}: '{rel['tags']['name']}', {len(ways)} ways -> {line.geom_type}")
coords = np.array(line.coords)
if coords[0, 1] > coords[-1, 1]:                        # walk it south to north: Edale first
    coords = coords[::-1]
gpx = ET.Element("gpx", version="1.1", creator="TechnicalGIS", xmlns="http://www.topografix.com/GPX/1/1")
trk = ET.SubElement(gpx, "trk")
ET.SubElement(trk, "name").text = rel["tags"]["name"]
seg = ET.SubElement(trk, "trkseg")
for lon, lat in coords:
    ET.SubElement(seg, "trkpt", lat=f"{lat:.7f}", lon=f"{lon:.7f}")
gpx_path = DATA / "pennine_way_edale_crowden.gpx"
ET.ElementTree(gpx).write(gpx_path, xml_declaration=True, encoding="UTF-8")

lonlat = read_gpx(gpx_path)                             # from here on, any GPX file works
to_bng = Transformer.from_crs("EPSG:4326", "EPSG:27700", always_xy=True)
E, N = to_bng.transform(lonlat[:, 0], lonlat[:, 1])
route = LineString(np.column_stack([E, N]))
print(f"read {gpx_path.name}: {len(lonlat):,} track points")
print(f"start {E[0]:,.0f} E {N[0]:,.0f} N | end {E[-1]:,.0f} E {N[-1]:,.0f} N | length {route.length / 1000:.2f} km")
print(f"inside the model area: {box(*AOI).contains(route)}")

# ------------------------------------------------------------------ Step 2: terrain
section("Step 2: Environment Agency LIDAR DTM, fetched at 5 m and averaged to 25 m")
X0, Y0, X1, Y1 = AOI
k = int(RES / FETCH)                                    # 5 x 5 blocks of 5 m cells per 25 m cell
fine = np.full((int((Y1 - Y0) / FETCH), int((X1 - X0) / FETCH)), np.nan)
n_tiles = 0
for x in range(X0, X1, TILE):
    for y in range(Y0, Y1, TILE):
        f = get(f"{WCS}&subset=E({x},{x + TILE})&subset=N({y},{y + TILE})&SCALEFACTOR={1 / FETCH}",
                DATA / f"ea_dtm_{x}_{y}_{TILE}m_tile_{FETCH:.0f}m.tif")
        with rasterio.open(f) as src:
            a = src.read(1).astype("float64")
            a[a < -1000] = np.nan
            r, c = int(round((Y1 - src.bounds.top) / FETCH)), int(round((src.bounds.left - X0) / FETCH))
            fine[r:r + a.shape[0], c:c + a.shape[1]] = a
        n_tiles += 1
print(f"{n_tiles} WCS tiles -> {fine.shape[1]} x {fine.shape[0]} cells of {FETCH:.0f} m | "
      f"no-data: {np.isnan(fine).sum():,}")
rows, cols = fine.shape[0] // k, fine.shape[1] // k
with np.errstate(all="ignore"):
    dem = np.nanmean(fine.reshape(rows, k, cols, k), axis=(1, 3))   # area mean, not point samples
gaps = np.isnan(dem)
print(f"block mean {k} x {k} -> {cols} x {rows} cells of {RES:.0f} m | empty cells filled: {gaps.sum():,}")
dem = dem[tuple(ndimage.distance_transform_edt(gaps, return_distances=False, return_indices=True))]
print(f"elevation {dem.min():.1f} - {dem.max():.1f} m")

# ------------------------------------------------------------------ Step 3: drape and profile
section("Step 3: drape the route on the DEM")
d = np.arange(0, route.length, RES)
pts = np.array([route.interpolate(s).coords[0] for s in d])
rc = np.vstack([(Y1 - pts[:, 1]) / RES - 0.5, (pts[:, 0] - X0) / RES - 0.5])   # fractional row, col
elev = ndimage.map_coordinates(dem, rc, order=1)
climb = np.diff(elev)
top_i = int(np.argmax(elev))
print(f"sampled every {RES:.0f} m: {len(d):,} points")
print(f"lowest {elev.min():.0f} m | highest {elev.max():.0f} m at {d[top_i] / 1000:.1f} km "
      f"({pts[top_i, 0]:,.0f} E {pts[top_i, 1]:,.0f} N)")
print(f"total ascent {climb[climb > 0].sum():,.0f} m | total descent {-climb[climb < 0].sum():,.0f} m")

# ------------------------------------------------------------------ Step 4: model scale
section("Step 4: model scale")
scale = MODEL_LONG_MM / ((max(rows, cols) - 1) * RES)          # mm per metre
cell_mm = RES * scale
Z = BASE_MM + (dem - dem.min()) * scale * Z_EXAGGERATION
W_mm, D_mm = (cols - 1) * cell_mm, (rows - 1) * cell_mm
print(f"horizontal scale 1:{1000 / scale:,.0f} | cell on the print {cell_mm:.3f} mm")
print(f"model {W_mm:.1f} x {D_mm:.1f} x {Z.max():.1f} mm | route {ROUTE_WIDTH_MM} mm wide = "
      f"{ROUTE_WIDTH_MM / scale:.0f} m on the ground")


def to_mm(xy):
    return np.column_stack([(xy[:, 0] - (X0 + RES / 2)) * scale, (xy[:, 1] - (Y0 + RES / 2)) * scale])


route_mm = LineString(to_mm(np.array(route.coords))).simplify(0.02)
band = route_mm.buffer(ROUTE_WIDTH_MM / 2, quad_segs=6).intersection(box(0.5, 0.5, W_mm - 0.5, D_mm - 0.5))
band = orient(band, 1.0)
contours = [np.array(band.exterior.coords)[:-1]] + [np.array(h.coords)[:-1] for h in band.interiors]
print(f"route line: {len(route.coords):,} -> {len(route_mm.coords):,} vertices after a 0.02 mm simplify")
print(f"route outline: {band.geom_type} with {sum(len(c) for c in contours):,} vertices, "
      f"{len(band.interiors)} holes, area {band.area:.1f} mm2")

# ------------------------------------------------------------------ Step 5: booleans
section("Step 5: raised, engraved and inlay solids")
t0 = time.perf_counter()
terrain = to_manifold(*heightfield_solid(Z, cell_mm))
terrain_up = to_manifold(*heightfield_solid(Z + RAISE_MM, cell_mm))
terrain_down = to_manifold(*heightfield_solid(Z - DEPTH_MM, cell_mm))
prism = mf.Manifold.extrude(mf.CrossSection(contours), Z.max() + 10).translate((0, 0, -5))
ribbon = prism ^ terrain_up                           # the route column, cut off 0.6 mm above the surface
raised = terrain + ribbon                             # single colour, route stands proud
engraved = terrain - (prism - terrain_down)           # single colour, 0.8 mm groove
inlay = ribbon - terrain_down                         # second body that fills the groove
meshes = {}
for name, m in {"terrain": terrain, "raised": raised, "engraved": engraved, "route inlay": inlay}.items():
    v, t = from_manifold(m)
    ok, vol = check_mesh(name, v, t)
    assert ok, name
    meshes[name] = (v, t, vol)
print(f"[booleans and checks: {time.perf_counter() - t0:.1f} s]")
overlap = (engraved ^ inlay).volume()
print(f"engraved + inlay = {(meshes['engraved'][2] + meshes['route inlay'][2]) / 1000:.3f} cm3 | "
      f"raised = {meshes['raised'][2] / 1000:.3f} cm3 | overlap between the two parts {overlap:.4f} mm3")
print(f"raised - terrain = {(meshes['raised'][2] - meshes['terrain'][2]):.1f} mm3 "
      f"(route area {band.area:.1f} mm2 x {RAISE_MM} mm = {band.area * RAISE_MM:.1f} mm3)")
print(f"route inlay: {len(inlay.decompose())} piece(s)")

# ------------------------------------------------------------------ Step 6: files
section("Step 6: write STL and 3MF")
files = {"raised": "pennine_way_raised.stl", "engraved": "pennine_way_engraved.stl",
         "route inlay": "pennine_way_route_inlay.stl"}
for name, fn in files.items():
    v, t, _ = meshes[name]
    write_stl(OUT / fn, v, t, f"TechnicalGIS Pennine Way {name} - EA LIDAR DTM, (c) OSM contributors")
    size = (OUT / fn).stat().st_size
    with open(OUT / fn, "rb") as fh:
        fh.seek(80)
        n = int(np.frombuffer(fh.read(4), "<u4")[0])
    print(f"{fn:<29s} {size / 1e6:5.1f} MB | {n:,} facets | 84 + 50 x facets = file size: {size == 84 + 50 * n}")
write_3mf(OUT / "pennine_way_inlay.3mf", [("terrain", *meshes["engraved"][:2], SAND),
                                          ("route", *meshes["route inlay"][:2], RED)])
with zipfile.ZipFile(OUT / "pennine_way_inlay.3mf") as z:
    model = ET.fromstring(z.read("3D/3dmodel.model"))
objs = model.findall(".//{http://schemas.microsoft.com/3dmanufacturing/core/2015/02}object")
print(f"pennine_way_inlay.3mf         {(OUT / 'pennine_way_inlay.3mf').stat().st_size / 1e6:5.1f} MB | "
      f"parts: {', '.join(o.get('name') for o in objs)}")

# ------------------------------------------------------------------ Figures
section("figures")
hs = LightSource(azdeg=315, altdeg=45).hillshade(dem, vert_exag=2, dx=RES, dy=RES)
fig, ax = plt.subplots(figsize=(6.2, 7.6))
ext = (X0 / 1000, X1 / 1000, Y0 / 1000, Y1 / 1000)
im = ax.imshow(dem, cmap="Spectral_r", extent=ext, alpha=0.85)
ax.imshow(hs, cmap="gray", extent=ext, alpha=0.35)
ax.plot(E / 1000, N / 1000, color=RED, lw=2.2, label="Pennine Way (OpenStreetMap)")
ax.plot(E / 1000, N / 1000, color="white", lw=0.6)
for (x, y), lab, ha in [((E[0], N[0]), "Edale", "right"), ((E[-1], N[-1]), "Crowden", "right"),
                         ((pts[top_i, 0], pts[top_i, 1]), f"highest point {elev.max():.0f} m", "left")]:
    ax.plot(x / 1000, y / 1000, "o", ms=6, mfc="white", mec="black")
    ax.annotate(lab, (x / 1000, y / 1000), xytext=(-8 if ha == "right" else 8, 4), textcoords="offset points",
                ha=ha, fontsize=8.5, bbox=dict(boxstyle="round,pad=0.2", fc="white", ec="none", alpha=0.8))
ax.set_xlabel("Easting (km, British National Grid)")
ax.set_ylabel("Northing (km)")
ax.legend(loc="lower left", fontsize=8)
plt.colorbar(im, ax=ax, shrink=0.6, label="Elevation (m)")
ax.set_title("Pennine Way, Edale to Crowden, on the EA LIDAR DTM (25 m)")
fig.savefig(OUT / "01-route-map.png"); plt.close(fig)

fig, ax = plt.subplots(figsize=(8, 3.2))
ax.fill_between(d / 1000, elev, elev.min() - 30, color="#abdda4", alpha=0.6)
ax.plot(d / 1000, elev, color="#3288bd", lw=1.2)
ax.set_xlim(0, d[-1] / 1000); ax.set_ylim(elev.min() - 30, elev.max() + 40)
ax.set_xlabel("Distance from Edale (km)"); ax.set_ylabel("Elevation (m)")
sec = ax.secondary_yaxis("right", functions=(lambda e: BASE_MM + (e - dem.min()) * scale * Z_EXAGGERATION,
                                               lambda z: (z - BASE_MM) / (scale * Z_EXAGGERATION) + dem.min()))
sec.set_ylabel("Height on the print (mm)")
ax.annotate(f"{elev.max():.0f} m", (d[top_i] / 1000, elev.max()), xytext=(0, 6), textcoords="offset points",
            ha="center", fontsize=8)
ax.grid(alpha=0.3)
ax.set_title("Elevation profile sampled from the DTM every 25 m")
fig.savefig(OUT / "02-elevation-profile.png"); plt.close(fig)

# vertical sections through the real meshes: move the y = y_cut plane to the origin, rotate it flat, slice
y_cut = float(np.median(to_mm(np.array(route.coords))[:, 1]))
x_hit = band.intersection(LineString([(0, y_cut), (W_mm, y_cut)])).centroid.x
zc = float(Z[int(round((D_mm - y_cut) / cell_mm)), int(round(x_hit / cell_mm))])
fig, axes = plt.subplots(1, 3, figsize=(9.5, 3.0), sharey=True)
for ax, (title, parts) in zip(axes, [(f"Raised {RAISE_MM} mm", [(raised, SAND)]),
                                    (f"Engraved {DEPTH_MM} mm", [(engraved, SAND)]),
                                    ("Two-part inlay", [(engraved, SAND), (inlay, RED)])]):
    for m, colour in parts:
        sl = m.translate((0, -y_cut, 0)).rotate((90, 0, 0)).slice(0.0)
        for poly in sl.to_polygons():
            p = np.vstack([poly, poly[:1]])
            ax.fill(p[:, 0], -p[:, 1], color=colour, ec="#555555", lw=0.5)   # the rotation maps z to -y
    ax.set_xlim(x_hit - 4, x_hit + 4)
    ax.set_title(title)
    ax.set_xlabel("x on the print (mm)")
    ax.set_aspect("equal")
axes[0].set_ylabel("z (mm)")
axes[0].set_ylim(zc - 3, zc + 2)
fig.suptitle(f"West-east sections through the meshes where the route crosses y = {y_cut:.1f} mm", y=1.02)
fig.savefig(OUT / "03-route-sections.png"); plt.close(fig)
print(f"sections at y = {y_cut:.1f} mm, route centre x = {x_hit:.1f} mm, terrain z = {zc:.1f} mm")


def render(parts, png, zoom=1.0, focal=None, elev_deg=35, azim_deg=-70, size=(1600, 1150)):
    p = pv.Plotter(off_screen=True, window_size=size)
    p.set_background("white")
    for fn, colour in parts:
        p.add_mesh(pv.read(OUT / fn), color=colour, smooth_shading=False, ambient=0.25, diffuse=0.75,
                   specular=0.1)
    p.enable_lightkit()
    c = np.array(focal if focal is not None else (W_mm / 2, D_mm / 2, 0))
    dist = 1.25 * max(W_mm, D_mm) / zoom
    el, az = np.radians(elev_deg), np.radians(azim_deg)
    p.camera_position = [tuple(c + dist * np.array([np.cos(el) * np.cos(az), np.cos(el) * np.sin(az), np.sin(el)])),
                         tuple(c), (0, 0, 1)]
    p.camera.view_angle = 30
    p.screenshot(OUT / png)
    p.close()


render([(files["engraved"], SAND), (files["route inlay"], RED)], "04-stl-preview-inlay.png")
fx, fy = to_mm(np.array([[pts[top_i, 0], pts[top_i, 1]]]))[0]
render([(files["raised"], SAND)], "05-stl-preview-raised-closeup.png", zoom=4.5,
       focal=(fx, fy, float(Z.max()) - 2), elev_deg=40, azim_deg=-60)
print("wrote 01-route-map, 02-elevation-profile, 03-route-sections, 04-stl-preview-inlay, "
      "05-stl-preview-raised-closeup")