"""Tile a large 3D-printed terrain model: the central Lake District as six interlocking jigsaw tiles. A 45 x 30 km Copernicus GLO-30 DEM is scaled once for the whole model (so seams match), the plan is cut into a 3 x 2 grid of 150 mm tiles with jigsaw knobs and a printing clearance, and every tile is proved watertight and non-overlapping. A second set of tiles uses straight edges with holes for 1.75 mm filament dowels instead. TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/tile-large-3d-print-terrain/ Data: Copernicus DEM GLO-30 (© DLR e.V. 2010-2014 and © Airbus Defence and Space GmbH 2014-2018, provided under COPERNICUS by the European Union and ESA). """ import sys from importlib.metadata import version from pathlib import Path import numpy as np import rasterio from rasterio.transform import from_origin from rasterio.warp import reproject, Resampling from shapely.geometry import Point, box from shapely.geometry.polygon import orient from shapely.ops import unary_union 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 OUT = Path("output"); OUT.mkdir(exist_ok=True) DEM_URLS = [f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N54_00_W00{w}_00_DEM/" f"Copernicus_DSM_COG_10_N54_00_W00{w}_00_DEM.tif" for w in (4, 3)] AOI = (310_000, 495_000, 355_000, 525_000) # 45 x 30 km of the central Lake District, BNG RES = 40.0 NX, NY = 3, 2 # tiles across and up TILE_MM = 150.0 # printed size of each square tile Z_EXAGGERATION = 2.0 BASE_MM = 5.0 CLEARANCE_MM = 0.2 # total gap between mating edges KNOB_NECK_MM, KNOB_R_MM, KNOB_CENTRE_MM = 6.0, 4.0, 6.0 PIN_HOLE_D_MM, PIN_DEPTH_MM = 1.9, 10.0 # for 1.75 mm filament used as dowels X0, Y0, X1, Y1 = AOI COLS, ROWS = int((X1 - X0) / RES) + 1, int((Y1 - Y0) / RES) + 1 # vertices on cell centres DST = from_origin(X0 - RES / 2, Y1 + RES / 2, RES, RES) def section(title): print(f"\n=== {title} ===") # ------------------------------------------------------------------ mesh helpers (as in Part 1) def grid_triangles(r, c): 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, x0=0.0, y0=0.0): rows, cols = Z.shape X, Y = np.meshgrid(x0 + np.arange(cols) * cell_mm, y0 + (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(verts, tris): 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 return bool((und == 2).all()), bool((drc == 1).all()), 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", " 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 prism(poly, height): polys = getattr(poly, "geoms", [poly]) contours = [] for p in polys: p = orient(p, 1.0) contours += [np.array(p.exterior.coords)[:-1]] + [np.array(h.coords)[:-1] for h in p.interiors] return mf.Manifold.extrude(mf.CrossSection(contours), height + 10).translate((0, 0, -5)) 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: one DEM from two source tiles section("Step 1: mosaic two GLO-30 tiles onto a 40 m grid") dem = np.full((ROWS, COLS), np.nan, dtype="float32") with rasterio.Env(GDAL_HTTP_MAX_RETRY="4", GDAL_HTTP_RETRY_DELAY="5"): for url in DEM_URLS: with rasterio.open(url) as src: reproject(rasterio.band(src, 1), dem, dst_transform=DST, dst_crs="EPSG:27700", resampling=Resampling.bilinear, dst_nodata=np.nan, init_dest_nodata=False) print(f"{url.rsplit('/', 1)[1]:<44s} NaNs left: {int(np.isnan(dem).sum()):,}") print(f"grid {COLS} x {ROWS} vertices at {RES:.0f} m | {np.nanmin(dem):.1f} - {np.nanmax(dem):.1f} m") # ------------------------------------------------------------------ Step 2: one scale for every tile section("Step 2: one horizontal and vertical scale for the whole model") scale = NX * TILE_MM / (X1 - X0) # 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"scale 1:{1000 / scale:,.0f} | cell {cell_mm:.2f} mm | whole model {W_mm:.0f} x {D_mm:.0f} x {Z.max():.1f} mm") print(f"{NX} x {NY} tiles of {TILE_MM:.0f} x {TILE_MM:.0f} mm = {TILE_MM / scale / 1000:.0f} x " f"{TILE_MM / scale / 1000:.0f} km each") assert abs(D_mm - NY * TILE_MM) < 1e-6 # ------------------------------------------------------------------ Step 3: jigsaw outlines section("Step 3: interlocking outlines with a printing clearance") def knob(x, y, direction): """A jigsaw knob rooted on a tile edge at (x, y), pointing 'E' or 'N'.""" dx, dy = (1, 0) if direction == "E" else (0, 1) cx, cy = x + dx * KNOB_CENTRE_MM, y + dy * KNOB_CENTRE_MM half = KNOB_NECK_MM / 2 if direction == "E": neck = box(x - 0.01, y - half, cx, y + half) else: neck = box(x - half, y - 0.01, x + half, cy) return neck.union(Point(cx, cy).buffer(KNOB_R_MM, quad_segs=16)) squares = {(i, j): box(i * TILE_MM, j * TILE_MM, (i + 1) * TILE_MM, (j + 1) * TILE_MM) for i in range(NX) for j in range(NY)} exact = dict(squares) for (i, j) in squares: for frac in (1 / 3, 2 / 3): if i < NX - 1: # knobs from this tile into its east neighbour k = knob((i + 1) * TILE_MM, (j + frac) * TILE_MM, "E") exact[(i, j)] = exact[(i, j)].union(k) exact[(i + 1, j)] = exact[(i + 1, j)].difference(k) if j < NY - 1: # and into its north neighbour k = knob((i + frac) * TILE_MM, (j + 1) * TILE_MM, "N") exact[(i, j)] = exact[(i, j)].union(k) exact[(i, j + 1)] = exact[(i, j + 1)].difference(k) whole = box(0, 0, W_mm, D_mm) pair_overlap = sum(exact[a].intersection(exact[b]).area for a in exact for b in exact if a < b) print(f"exact partition: sum of tile areas {sum(p.area for p in exact.values()):,.3f} mm2 | " f"whole model {whole.area:,.3f} mm2 | overlap between tiles {pair_overlap:.6f} mm2") # shrink every outline by half the clearance, so that each joint has the full clearance between parts outlines = {k: p.buffer(-CLEARANCE_MM / 2, join_style="mitre") for k, p in exact.items()} gaps = [outlines[a].distance(outlines[b]) for a in outlines for b in outlines if a < b and exact[a].touches(exact[b])] print(f"clearance between neighbouring tiles: min {min(gaps):.3f} mm, max {max(gaps):.3f} mm " f"over {len(gaps)} joints") name = {(i, j): f"{'ABC'[i]}{j + 1}" for i in range(NX) for j in range(NY)} # ------------------------------------------------------------------ Step 4: cut the tiles section("Step 4: cut six tiles and prove they fit") def tile_solid(poly): """Heightfield over the outline's bounding box (plus a cell), cut to the outline.""" minx, miny, maxx, maxy = poly.bounds c0, c1 = max(int(minx / cell_mm) - 1, 0), min(int(np.ceil(maxx / cell_mm)) + 2, COLS) r0 = max(ROWS - 1 - int(np.ceil(maxy / cell_mm)) - 1, 0) r1 = min(ROWS - int(miny / cell_mm) + 1, ROWS) hf = to_manifold(*heightfield_solid(Z[r0:r1, c0:c1], cell_mm, c0 * cell_mm, (ROWS - r1) * cell_mm)) return hf ^ prism(poly, Z.max()) tiles = {k: tile_solid(p) for k, p in outlines.items()} whole_terrain = to_manifold(*heightfield_solid(Z, cell_mm)) for k, m in tiles.items(): v, t = from_manifold(m) wt, ori, vol = check_mesh(v, t) assert wt and ori and vol > 0 lo, hi = v.min(axis=0), v.max(axis=0) print(f"tile {name[k]}: {len(t):>7,} triangles | watertight {wt} | oriented {ori} | {vol / 1000:6.1f} cm3 | " f"footprint {hi[0] - lo[0]:.1f} x {hi[1] - lo[1]:.1f} mm, height {hi[2]:.1f} mm") write_stl(OUT / f"lake_district_jigsaw_{name[k]}.stl", v, t, f"TechnicalGIS Lake District tile {name[k]} - Copernicus GLO-30") keys = list(tiles) overlap = sum((tiles[a] ^ tiles[b]).volume() for i, a in enumerate(keys) for b in keys[i + 1:]) lost = whole_terrain.volume() - sum(m.volume() for m in tiles.values()) print(f"pairwise overlap between tiles: {overlap:.4f} mm3") print(f"whole model {whole_terrain.volume() / 1000:.2f} cm3 | tiles {sum(m.volume() for m in tiles.values()) / 1000:.2f} cm3 | " f"removed by the clearance {lost / 1000:.2f} cm3") # ------------------------------------------------------------------ Step 5: the dowel-pin alternative section("Step 5: straight edges with holes for filament dowels") def pin_hole(x, y, axis): """A horizontal hole, PIN_DEPTH_MM either side of a tile edge, half way up the base.""" cyl = mf.Manifold.cylinder(2 * PIN_DEPTH_MM, PIN_HOLE_D_MM / 2, circular_segments=24, center=True) cyl = cyl.rotate((0, 90, 0)) if axis == "x" else cyl.rotate((90, 0, 0)) return cyl.translate((x, y, BASE_MM / 2)) holes, centres = [], [] for (i, j) in squares: for frac in (1 / 3, 2 / 3): if i < NX - 1: holes.append(pin_hole((i + 1) * TILE_MM, (j + frac) * TILE_MM, "x")) centres.append(((i + 1) * TILE_MM, (j + frac) * TILE_MM)) if j < NY - 1: holes.append(pin_hole((i + frac) * TILE_MM, (j + 1) * TILE_MM, "y")) centres.append(((i + frac) * TILE_MM, (j + 1) * TILE_MM)) all_holes = mf.Manifold.batch_boolean(holes, mf.OpType.Add) pinned = {k: tile_solid(sq.buffer(-CLEARANCE_MM / 2, join_style="mitre")) - all_holes for k, sq in squares.items()} for k, m in pinned.items(): v, t = from_manifold(m) wt, ori, vol = check_mesh(v, t) assert wt and ori and vol > 0 write_stl(OUT / f"lake_district_pinned_{name[k]}.stl", v, t, f"TechnicalGIS Lake District pinned tile {name[k]} - Copernicus GLO-30") print(f"{len(holes)} dowels, each {2 * PIN_DEPTH_MM - CLEARANCE_MM:.1f} mm long; every pinned tile watertight " f"and oriented: True") print("pin holes per tile: " + ", ".join( f"{name[k]} {sum(sq.boundary.distance(Point(c)) < 1e-6 for c in centres)}" for k, sq in squares.items())) wall = (BASE_MM - PIN_HOLE_D_MM) / 2 print(f"plastic above and below each hole: {wall:.2f} mm") # ------------------------------------------------------------------ Figures section("figures") km = 1 / scale / 1000 def to_km(x, y): return X0 / 1000 + np.asarray(x) * km, Y0 / 1000 + np.asarray(y) * km ext = (X0 / 1000 - RES / 2000, X1 / 1000 + RES / 2000, Y0 / 1000 - RES / 2000, Y1 / 1000 + RES / 2000) hs = LightSource(azdeg=315, altdeg=45).hillshade(dem, vert_exag=1.5, dx=RES, dy=RES) fig, ax = plt.subplots(figsize=(10, 6.8)) im = ax.imshow(dem, cmap="Spectral_r", extent=ext) ax.imshow(hs, cmap="gray", extent=ext, alpha=0.35) for k, p in exact.items(): xs, ys = to_km(*p.exterior.xy) ax.plot(xs, ys, color="black", lw=1.0) cx, cy = to_km(*np.array(squares[k].centroid.coords[0])) ax.text(cx, cy, name[k], ha="center", va="center", fontsize=16, weight="bold", color="white", bbox=dict(boxstyle="round", fc="black", alpha=0.45, lw=0)) for label, e, n in [("Scafell Pike", 321_540, 507_210), ("Helvellyn", 334_210, 515_140)]: ax.plot(e / 1000, n / 1000, "k^", ms=6) ax.annotate(label, (e / 1000, n / 1000), xytext=(5, 5), textcoords="offset points", fontsize=8) plt.colorbar(im, ax=ax, shrink=0.7, label="Elevation (m)") ax.set_xlabel("Easting (km, British National Grid)") ax.set_ylabel("Northing (km)") ax.set_title(f"Copernicus GLO-30, 40 m: {NX} x {NY} tiles of {TILE_MM / scale / 1000:.0f} km " f"({TILE_MM:.0f} mm printed)") fig.savefig(OUT / "01-tile-layout.png"); plt.close(fig) # a close-up of one joint, showing the clearance fig, axes = plt.subplots(1, 2, figsize=(10, 4.6)) jx, jy = TILE_MM, TILE_MM / 3 for ax, (half_w, title) in zip(axes, [(14, "One jigsaw knob (tile A1 into B1)"), (1.2, "Zoomed to the 0.2 mm clearance")]): for k, col in [((0, 0), "#3288bd"), ((1, 0), "#f46d43")]: p = outlines[k] for g in getattr(p, "geoms", [p]): ax.fill(*g.exterior.xy, color=col, alpha=0.35, lw=0) ax.plot(*g.exterior.xy, color=col, lw=1.2) ax.plot(*exact[k].exterior.xy, color="black", lw=0.6, ls=":") cx = jx + (KNOB_CENTRE_MM + KNOB_R_MM if half_w < 2 else KNOB_CENTRE_MM) ax.set_xlim(cx - half_w, cx + half_w) ax.set_ylim(jy - half_w * 0.8, jy + half_w * 0.8) ax.set_aspect("equal") ax.set_title(title) ax.set_xlabel("x on the print (mm)") axes[0].set_ylabel("y on the print (mm)") axes[0].legend(handles=[plt.Line2D([], [], color="#3288bd", lw=2, label="tile A1"), plt.Line2D([], [], color="#f46d43", lw=2, label="tile B1"), plt.Line2D([], [], color="black", lw=0.8, ls=":", label="exact cut line")], fontsize=8, loc="lower left") fig.savefig(OUT / "02-joint-clearance.png"); plt.close(fig) def save_cropped(p, png, pad=20): img = p.screenshot(return_img=True) ink = np.where((img < 250).any(axis=2)) r0, r1 = max(ink[0].min() - pad, 0), min(ink[0].max() + pad, img.shape[0]) c0, c1 = max(ink[1].min() - pad, 0), min(ink[1].max() + pad, img.shape[1]) plt.imsave(OUT / png, img[r0:r1, c0:c1]) p.close() def look(p, centre, dist, elev_deg, azim_deg): el, az = np.radians(elev_deg), np.radians(azim_deg) c = np.asarray(centre, dtype=float) 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 clim = (float(Z.min()), float(Z.max())) p = pv.Plotter(off_screen=True, window_size=(2000, 1300)) p.set_background("white") for k in tiles: mesh = pv.read(OUT / f"lake_district_jigsaw_{name[k]}.stl").translate((k[0] * 8.0, k[1] * 8.0, 0)) p.add_mesh(mesh, scalars=mesh.points[:, 2], cmap="Spectral_r", clim=clim, show_scalar_bar=False, ambient=0.2, diffuse=0.8) p.enable_lightkit() look(p, (W_mm / 2 + 8, D_mm / 2 + 4, 5), 1.75 * W_mm, 38, -70) save_cropped(p, "03-stl-preview-tiles.png") fig_imgs = [] p = pv.Plotter(off_screen=True, window_size=(1800, 900), shape=(1, 2), border=False) p.subplot(0, 0) p.set_background("white") for k, off, colour in [((0, 0), -4.0, "#d9c8a9"), ((1, 0), 4.0, "#9ecae1")]: mesh = pv.read(OUT / f"lake_district_jigsaw_{name[k]}.stl") p.add_mesh(mesh.clip_box((TILE_MM - 30, TILE_MM + 30, 25, 75, -1, 100), invert=False).translate((off, 0, 0)), color=colour, ambient=0.25, diffuse=0.75) p.add_text("Jigsaw knob, tiles A1 and B1 pulled 8 mm apart", font_size=12, color="#222222") p.enable_lightkit() look(p, (TILE_MM, 50, 6), 125, 58, -95) p.subplot(0, 1) p.set_background("white") mesh = pv.read(OUT / f"lake_district_pinned_{name[(1, 0)]}.stl") p.add_mesh(mesh.clip_box((TILE_MM - 1, TILE_MM + 20, 82, 118, -1, 100), invert=False), color="#9ecae1", ambient=0.25, diffuse=0.75) p.add_text("Dowel hole, west edge of pinned tile B1", font_size=12, color="#222222") p.enable_lightkit() look(p, (TILE_MM, 100, 8), 80, 16, 205) save_cropped(p, "04-joint-closeups.png") print("wrote 01-tile-layout, 02-joint-clearance, 03-stl-preview-tiles, 04-joint-closeups")