--- title: "3D-print a geological ground model: one watertight STL per stratum" description: "Turn kriged stratum surfaces into a multi-material 3D print: one watertight STL per geological unit, checked to fit exactly, plus a colour 3MF." url: https://www.technicalgis.com/tutorials/3d-print-geological-ground-model/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-11-15 updated: 2026-11-15 language: en-GB --- # 3D-print a geological ground model: one watertight STL per stratum > Turn kriged stratum surfaces into a multi-material 3D print: one watertight STL per geological unit, checked to fit exactly, plus a colour 3MF. - Level: Advanced - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/3d-print-geological-ground-model/) Series · Part 7 of 7 ## AGS data in GIS 1. [Read and validate AGS4 ground investigation data with Python](https://www.technicalgis.com/tutorials/read-validate-ags4-python/) 2. [Boreholes on a map: AGS4 locations in British National Grid with QGIS](https://www.technicalgis.com/tutorials/ags4-boreholes-map-qgis/) 3. [3D borehole logs from AGS4 data in PyVista and QGIS](https://www.technicalgis.com/tutorials/3d-borehole-logs-pyvista-qgis/) 4. [A geological cross-section from AGS4 boreholes and LIDAR in Python and QGIS](https://www.technicalgis.com/tutorials/geological-cross-section-ags4-python/) 5. [A 3D ground model from AGS4 boreholes: kriged stratum surfaces in Python](https://www.technicalgis.com/tutorials/geological-ground-model-surfaces-python/) 6. [Visualising SPT N-values from AGS4 data in Python](https://www.technicalgis.com/tutorials/spt-n-values-ags-python/) 7. 3D-print a geological ground model: one watertight STL per stratum ## Key takeaways 1. Build one solid per unit from the stratum surfaces: a heightfield solid from the plinth up to each surface, then subtract the one below. Subtract first and trim to the footprint second, or the shared side walls leave needle-thin triangles on the outside. 2. Decide what is too thin to print in millimetres, not metres. At 1:9,450 with three times vertical exaggeration, 0.6 mm is 1.89 m of strata: about half the superficial deposits' area is handed to the unit below, and the thick head and landslide deposits stay. 3. Prove the print before you write it: every part watertight and consistently oriented, the parts adding up to the ground solid to within 0.0000%, and no volume shared between any two parts. 4. Render the previews from the STL files read back from disk, so the pictures show exactly what the slicer will get. ## Before you start - Part 5 of this series run first: this script reads its output/model/ surfaces - Python 3.10+ with numpy, rasterio, shapely, scipy, manifold3d, pyvista and matplotlib - A multi-material printer (Bambu Lab AMS, Prusa MMU or similar) for the colour print, or filament changes on a single extruder **A ground model you can hold in your hand explains a site faster than any section.** In [Part 5](https://www.technicalgis.com/tutorials/geological-ground-model-surfaces-python/) we kriged the A417 boreholes into a stack of stratum surfaces under the LIDAR ground, and tested them. Now we print them: a 20 cm block of the Cotswold escarpment at Crickley Hill, with the superficial deposits, Great Oolite, Inferior Oolite, Bridport Sand and Lias each printed as a separate part in its own colour. The approach is the one from the [3D-printed maps with Python](https://www.technicalgis.com/tutorials/multi-colour-3d-print-terrain/) series: build one watertight solid per colour with mesh booleans, check every part, and write STLs and a 3MF. What is new is that the colours are not drawn on the surface; they are the geology, and they go all the way through. Cut the print in half and you see the section. The steps are: 1. read the stratum surfaces in the print window; 2. trace the footprint of the model; 3. convert to print heights, and decide what is too thin to print; 4. build one solid per unit; 5. check every part before writing; 6. write one STL per unit and a colour 3MF; 7. render the previews from the STL files. The code is [`print_ground_model.py`](https://www.technicalgis.com/code/3d-print-geological-ground-model/). It reads `output/model/` from Part 5 and runs in about seven seconds. ![Rendered 3D print of the Crickley Hill ground model, 198 by 145 by 71 mm, looking north-east up the escarpment. Brown superficial deposits cover the western slope, orange Inferior Oolite forms the scarp and plateau, pale yellow Great Oolite caps the south-east, and the side walls show a band of olive Bridport Sand over blue-grey Lias mudstones.](https://www.technicalgis.com/_astro/01-printable-ground-model.CoJmyzOQ_Z1dYfNR.webp) _The print, rendered from the five STL files. The side walls show the layers the boreholes found under the escarpment: Inferior Oolite over Bridport Sand over Lias, dipping gently to the south-east._ ## Step 1: the stratum surfaces in the print window Part 5 wrote each surface as a GeoTIFF on the same 10 m grid, plus a `model.json` listing the units from the top down. The top of each unit is the base of the one above, and the bottom unit runs down to a flat model base at 80 m OD. The script checks that before it reads anything, so a change to Part 5 that broke the stack would stop here, not in the slicer. `print_ground_model.py` ```python meta = json.loads((args.model / "model.json").read_text(encoding="utf-8")) units, base_m = meta["units"], float(meta["model_base_mOD"]) for upper, lower in zip(units, units[1:]): assert upper["base"] == lower["top"], f"{upper['code']} base and {lower['code']} top must be the same surface" assert units[-1]["base"] is None, "the bottom unit runs down to the flat model base" grids = [] for u in units: with rasterio.open(args.model / f"{u['top']}.tif") as src: win = from_bounds(*WINDOW, transform=src.transform).round_offsets().round_lengths() win = win.intersection(Window(0, 0, src.width, src.height)) grids.append(src.read(1, window=win).astype("float64")) res, wt = src.res[0], src.window_transform(win) S = np.stack(grids) # (surface, row, col), top surface first valid = np.isfinite(S).all(axis=0) ``` The window, 1.9 by 1.4 km, takes in the scarp face, the plateau and the vale below, where the boreholes are densest. Part 5 only modelled ground within 300 m of a logged hole, so the model covers 89.6% of the window. Part 5 already made sure the surfaces never cross, but a print is unforgiving, so the script applies the rule again (`np.minimum.accumulate` down the stack) and counts how many cells it had to change. None did: `Output` ```text === Step 1: read the stratum surfaces in the print window === window 392400-394300 E, 215000-216400 N: 190 x 140 cells of 10 m | model domain covers 89.6% of it 5 units, model base 80 m OD | cells needing re-clamping 0, below the model base (raised to it) 0 top of SUP Superficial deposits 134.0 to 295.4 m OD top of GOG Great Oolite Group 127.8 to 294.3 m OD top of IOG Inferior Oolite Group 127.8 to 293.1 m OD top of BDS Bridport Sand Formation 127.8 to 238.7 m OD top of LIAS Lias mudstones (Whitby Mudstone and older) 127.8 to 212.6 m OD ``` Where a unit has been eroded away, its top and base are the same surface. That is why so many units share a minimum of 127.8 m OD: down in the vale they all have zero thickness on top of the Lias. ## Step 2: the footprint The print follows the edge of the model, not the rectangle of the window. Each mesh vertex sits at a cell centre, so a grid square only holds real data if all four of its corners are valid. The script turns those squares into a polygon, pulls the edge in by one cell and simplifies it by no more than that. This gives clean straight walls instead of a pixel staircase, without ever reaching outside the modelled ground: `print_ground_model.py` ```python scale = LONG_MM / ((max(rows, cols) - 1) * res) # mm on the print per metre on the ground cell = res * scale quad = valid[:-1, :-1] & valid[1:, :-1] & valid[:-1, 1:] & valid[1:, 1:] quad_to_mm = rasterio.Affine(cell, 0, 0, 0, -cell, (rows - 1) * cell) foot = unary_union([shape(g) for g, _ in features.shapes(quad.astype("uint8"), mask=quad, transform=quad_to_mm)]) foot = foot.buffer(-cell, join_style="mitre").simplify(cell) pieces = sorted(getattr(foot, "geoms", [foot]), key=lambda g: g.area, reverse=True) foot = orient(pieces[0], 1.0) # outer ring anticlockwise, holes clockwise ``` `Output` ```text === Step 2: the footprint of the model domain === horizontal scale 1:9,450 | one 10 m cell = 1.06 mm on the print footprint 2.28 km2 = 255.3 cm2 on the print, 0 hole(s) | 0 smaller piece(s) left out (0 m2) ``` A 10 m cell is 1.06 mm on the print, about two and a half nozzle widths, so a 0.4 mm nozzle can show every cell of the grid. ## Step 3: print heights, and what is too thin to print First the scale. The print is 200 mm on its long side, so 198 by 145 mm fits a common 220 mm build plate with room for a brim, and a 10 m cell is 1.06 mm: two and a half nozzle widths. Heights become millimetres above the print bed: metres above the model base, times the scale, times a vertical exaggeration of three, on top of a 3 mm plinth. Three times is the exaggeration of the 3D views in Parts 3 and 5, so the print matches the pictures. It makes the Bridport Sand a clear 9 mm band and keeps more of the superficial deposits, while the scarp still reads as a slope rather than a cliff. Then comes the most important decision in the script. A layer thinner than about three print layers will not print as a recognisable band: the slicer either drops it or prints a speckle of one colour over another. So any unit thinner than `MIN_MM = 0.6` mm in a cell is handed to the unit below. The script raises that unit’s base up to its top, so the next real unit down fills the space: `print_ground_model.py` ```python Z = (S - base_m) * scale * VE + BASE_MM for k, u in enumerate(units[:-1]): t = Z[k] - Z[k + 1] thin = (t > 1e-6) & (t < MIN_MM) old = Z[k + 1].copy() for j in range(k + 1, len(units)): Z[j] = np.where(thin & (Z[j] >= old - 1e-9), Z[k], Z[j]) ``` The inner loop matters where units have been eroded. If the Great Oolite is absent, its base is the same surface as the top of the Inferior Oolite, so both are raised together and the Inferior Oolite takes the space. `Output` ```text === Step 3: print heights, vertical exaggeration x3, nothing thinner than 0.6 mm === 1 m of ground = 0.317 mm on the print, so 0.6 mm = 1.89 m of strata SUP present under 96.0% of the footprint | thinner than 0.6 mm under 48.5% (1,094,700 m2) -> absorbed by the unit below GOG present under 26.1% of the footprint | thinner than 0.6 mm under 1.4% (32,400 m2) -> absorbed by the unit below IOG present under 75.6% of the footprint | thinner than 0.6 mm under 0.4% (8,200 m2) -> absorbed by the unit below BDS present under 80.9% of the footprint | thinner than 0.6 mm under 0.6% (12,800 m2) -> absorbed by the unit below ``` At this scale 0.6 mm is 1.89 m of ground. The bedrock units are almost all thicker than that, and lose at most 1.4% of the footprint. The superficial deposits are different: they are present almost everywhere, but over much of the plateau they are a metre or so of topsoil and weathered stone (the median logged thickness is 1.0 m). They are thinner than 1.89 m under 48.5% of the footprint, so about half their area prints in the colour of the rock beneath and half survives. What stays is the thick head and landslide debris on the scarp slope and the deeper patches on the plateau, which is what an engineer would want to see. For comparison, a 160 mm print at twice exaggeration would have lost 62% of their area. ![Two panels. Top: a west-to-east section through the five STL files at true print scale, about 200 mm long and up to 64 mm high, showing a thin brown superficial layer over blue-grey Lias in the west, and an olive Bridport Sand band under a thick orange Inferior Oolite in the east, with a small pale yellow wedge of Great Oolite at the eastern end. Bottom: histograms of print thickness for four units, with a dashed red line at 0.6 mm; the superficial deposits have a tall spike below the line, the Bridport Sand peaks at about 8.7 mm and the Inferior Oolite spreads from 17 to 22 mm.](https://www.technicalgis.com/_astro/03-print-layers-section.C5Iuphhk_Z1a4tGB.webp) _(a) A section cut through the STL files themselves, 85 mm from the front edge, at true print scale. (b) How thick each unit prints where it is present, before the 0.6 mm rule. Only the superficial deposits have much of their area below the line._ ## Step 4: one solid per unit Each unit is built from two heightfield solids, using the same `heightfield_solid` helper as the [multi-colour terrain print](https://www.technicalgis.com/tutorials/multi-colour-3d-print-terrain/). `cum[k]` runs from the bottom of the plinth up to surface `k`, so unit `k` is `cum[k]` minus `cum[k + 1]`. The footprint is extruded into a prism and intersected with each part: `print_ground_model.py` ```python prism = mf.Manifold.extrude(mf.CrossSection(contours), float(Z.max()) + 2).translate((0, 0, -1)) cum = [to_manifold(*heightfield_solid(z, cell)) for z in Z] # plinth bottom up to each surface raw = [(cum[k] - cum[k + 1]) ^ prism for k in range(len(units) - 1)] + [cum[-1] ^ prism] whole = cum[0] ^ prism solids = {} for u, m in zip(units, raw): solids[u["code"]], n = drop_slivers(m) ``` The order of the operations matters. If you trim each solid to the footprint first and subtract afterwards, every solid has the same side walls, and subtracting two coincident walls leaves needle-thin triangles on the outside of the print. Subtracting first and trimming last cuts all the side walls in one clean operation. Where two surfaces touch, as they do wherever a unit thins to zero, the booleans can leave zero-volume slivers. `drop_slivers` splits each result into connected pieces and keeps only those over 0.05 mm³: `Output` ```text === Step 4: cumulative solids and booleans === whole model 197.9 x 145.0 x 71.4 mm, including the 3 mm plinth SUP 4 separate piece(s) GOG 1 separate piece(s) IOG 1 separate piece(s) BDS 1 separate piece(s) LIAS 1 separate piece(s) ``` The superficial deposits are in four pieces: the large spread on the scarp slope, a strip along the north-eastern edge of the plateau and two small patches on it. All four are real deposits thicker than 1.89 m. ![The five STL parts rendered pulled apart vertically by 14 mm each. From the top: a brown superficial-deposits sheet over the western slope, a pale yellow Great Oolite cap in the south-east, an orange Inferior Oolite slab, an olive Bridport Sand slab and a blue-grey Lias block with the plinth. Each part's underside matches the top of the part below.](https://www.technicalgis.com/_astro/02-exploded-view.BN8TTWEy_1KB5fN.webp) _The five STL files pulled apart by 14 mm. Each is a separate solid whose underside is exactly the top of the part below, so they fit together with no gaps._ ## Step 5: check every part before writing Nothing is written until every part passes. `check_mesh` counts edges: in a watertight mesh every edge is shared by exactly two triangles, and in a consistently oriented one each edge is used once in each direction. It also works out the volume from the triangles, which must be positive (a negative volume means the mesh is inside out): `print_ground_model.py` ```python def check_mesh(verts, tris): """Every edge shared by exactly two triangles, used once in each direction, and a positive volume.""" 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 ``` Then the parts are checked against each other. Their volumes must add up to the whole ground solid, and no two parts may share any volume. If either test fails, the script stops: `print_ground_model.py` ```python total, whole = sum(m[2] for m in meshes.values()), whole.volume() overlap = sum((solids[a] ^ solids[b]).volume() for i, a in enumerate(codes) for b in codes[i + 1:]) if abs(total - whole) > 0.001 * whole or overlap > 0.001 * whole: raise RuntimeError("the parts do not add up to the ground solid") ``` `Output` ```text === Step 5: check every part before writing === sum of parts 1450.952 cm3 | ground solid 1450.952 cm3 | difference +0.0000% | pairwise overlap 0.0000 mm3 ``` The parts fill the model exactly, with nothing missing and nothing printed twice. ## Step 6: one STL per unit and a colour 3MF Each part is written as a binary STL with the structured-dtype writer from [DEM to STL](https://www.technicalgis.com/tutorials/dem-to-stl-python/), numbered from the top down so they sort in order in the slicer. All five also go into one 3MF, with the unit names and the same colours as Part 5’s figures: `Output` ```text === Step 6: write one STL per unit and a colour 3MF === file unit triangles cm3 PLA g watertight 01_SUP.stl Superficial deposits 46,676 35.62 44.2 True 02_GOG.stl Great Oolite Group 24,120 23.63 29.3 True 03_IOG.stl Inferior Oolite Group 72,058 305.67 379.0 True 04_BDS.stl Bridport Sand Formation 76,698 172.61 214.0 True 05_LIAS.stl Lias mudstones (Whitby Mudstone and older) 93,058 913.43 1132.6 True total 312,610 1450.95 1799.2 (solid PLA; at 15% infill a slicer will use much less) ground_model.3mf 2.9 MB | 5 parts: Superficial deposits, Great Oolite Group, Inferior Oolite Group, Bridport Sand Formation, Lias mudstones (Whitby Mudstone and older) ``` The volumes are for solid parts. With two or three walls and 15% infill, the Lias block uses a fraction of its 1.1 kg, and the thin upper parts are mostly wall anyway. If you want to save more, raise the model base in Part 5: everything below the lowest surface, 127.8 m OD, is plain Lias plinth. ## Step 7: previews from the STL files Every picture in this post is rendered from the STL files read back from disk with `pv.read`, not from the arrays in memory. If the files were wrong, the pictures would be too. The section in figure 3 is cut through the STL meshes with `manifold3d`’s `slice`. `Output` ```text === Step 7: previews rendered from the STL files === read back 5 STLs: 197.9 x 145.0 x 71.4 mm figures written: 01-printable-ground-model.png, 02-exploded-view.png, 03-print-layers-section.png, 04-slicer-style-top.png finished in 6.9 s ``` ![The print seen from directly above, north up. Brown superficial deposits cover the west, orange Inferior Oolite forms a broad band from north to south through the middle, and pale yellow Great Oolite covers the south-east, crossed by a straight cutting that exposes the orange Inferior Oolite beneath.](https://www.technicalgis.com/_astro/04-slicer-style-top.Cs3H0C9R_1skjTr.webp) _Seen from above, the print is a geological map. The Great Oolite on the plateau, the Inferior Oolite on the scarp and the head and landslide deposits on the slope below all match the pattern the boreholes found in Parts 3 and 5._ ## Printing tips - **Loading the parts:** drag all five STLs into PrusaSlicer, Bambu Studio or OrcaSlicer together and accept the offer to load them as **one object with multiple parts**, so they stay aligned. Or open `ground_model.3mf`, which keeps the names and positions. Then give each part its filament. - **Colours:** follow the geological map if you can: a brown or tan for superficial deposits, cream for the Great Oolite, orange for the Inferior Oolite, olive or khaki for the Bridport Sand and a slate blue-grey for the Lias. What matters most is that neighbouring units contrast. - **Layer height:** 0.2 mm suits the 0.6 mm rule. If you print at 0.12 mm, you could lower `MIN_MM` to about 0.4 mm. - **Purge:** the Lias block is large, so use “flush into infill” to absorb most of the purge from the colour changes. - **Single extruder:** the units are not flat bands, so a filament change at one height will not reproduce them. Print the Lias block and the upper parts separately and glue them, or print in one colour and paint the side walls. ## Frequently asked questions ### How do I 3D print a geological model? Build one surface per geological contact, for example by kriging borehole picks as in Part 5. Turn each surface into a heightfield solid from a flat base, subtract the solid below to get each unit, check that every part is watertight and that they add up to the whole, and write one STL per unit. Load the STLs into your slicer as parts of one object and give each unit its own filament. ### Why do the units go all the way through, instead of just colouring the surface? Because the side walls are the most useful part of a geological print: they are cross-sections. Solid units also make the print robust, and you can cut it with a saw along any line to show a section, as figure 3 does on screen. ### What vertical exaggeration should I use for a geological print? Enough that the thinnest unit you care about is at least three or four print layers thick, without turning slopes into cliffs. Three times suits this site: the Bridport Sand, mostly about 28 m thick, prints at about 8.7 mm (figure 3b), half the superficial deposits’ area survives, and it matches the 3D views in Parts 3 and 5. More would keep more superficial deposits but would make the scarp look far steeper than it is. ### Can I print the model at a larger size? Yes. Change `LONG_MM`. At 250 mm the scale is about 1:7,600, a 10 m cell is 1.3 mm and, at three times exaggeration, 0.6 mm is 1.5 m of strata. For anything bigger than your print bed, the [tiling tutorial](https://www.technicalgis.com/tutorials/tile-large-3d-print-terrain/) cuts a model into jigsaw tiles, and the same cut works on each unit. ### Why does the script check the parts when manifold3d already makes valid meshes? Because the check is cheap and the print is expensive. The edge and volume checks catch an inside-out part or a gap at once, and the sum and overlap checks catch a mistake in the stack that would leave two parts claiming the same space. ## Next steps - Back to [Part 5: the ground model](https://www.technicalgis.com/tutorials/geological-ground-model-surfaces-python/), to change the surfaces or extend the model before printing. - [Part 6: SPT N-values](https://www.technicalgis.com/tutorials/spt-n-values-ags-python/) adds the test results, showing which of the printed units are soft and which refuse. - [Multi-colour 3D-printed terrain](https://www.technicalgis.com/tutorials/multi-colour-3d-print-terrain/) explains the one-solid-per-colour approach in more detail, including thin colour caps that save filament. - [From DEM to 3D-printable STL](https://www.technicalgis.com/tutorials/dem-to-stl-python/) starts the 3D-printed maps series from the beginning, with the heightfield solid and the STL writer used here. - [Tiling a large 3D-printed terrain model](https://www.technicalgis.com/tutorials/tile-large-3d-print-terrain/) for a model of the whole A417 route. ## Data and licences - **A417 Missing Link ground investigation** (BGS deposited data items 173984 and 173975; client Highways England, now National Highways), downloaded through the BGS AGS borehole index and export APIs in Part 5. Contains British Geological Survey materials © UKRI 2026, licensed under the [Open Government Licence v3.0](https://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/). - **Environment Agency LIDAR Composite DTM 1 m**, the top surface of the model, from the Defra Data Services Platform WCS. © Environment Agency copyright and/or database right, licensed under the Open Government Licence v3.0. Every figure carries the credit. - Code: [`print_ground_model.py`](https://www.technicalgis.com/code/3d-print-geological-ground-model/), tested with Python 3.12.0, numpy 2.5.4, rasterio 1.5.2, shapely 2.2.0, manifold3d 3.5.5 and pyvista 0.49.2. It runs entirely offline on Part 5’s output. ## Key terms - **Heightfield solid**: A closed mesh made from a grid of heights: a top surface, a flat bottom and walls around the edge. Each unit here is the difference between two of them. - **Watertight**: A mesh with no holes: every edge is shared by exactly two triangles. Slicers need watertight parts to know what is inside. - **3MF**: A zipped 3D printing format that keeps several named parts together, in position and with colours, in one file. - **Vertical exaggeration**: Stretching heights relative to horizontal distances so that relief and thin layers show at small scales. Here heights are tripled. ## Data for this tutorial - [BGS AGS borehole data](https://www.technicalgis.com/data-sources/united-kingdom/#src-bgs-ags-borehole-data) · Donated AGS ground investigation data, Open Government Licence - [GIS data sources for the United Kingdom](https://www.technicalgis.com/data-sources/united-kingdom/) · Including Environment Agency LIDAR, Open Government Licence Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). [← Part 6 Visualising SPT N-values from AGS4 data in Python](https://www.technicalgis.com/tutorials/spt-n-values-ags-python/) --- Source: https://www.technicalgis.com/tutorials/3d-print-geological-ground-model/