"""Borehole logs in the QGIS 3D map view, on LIDAR terrain cut away in front of the escarpment holes. Part 3 of the TechnicalGIS series "AGS data in GIS". Run borehole_logs_3d.py first, then run this with the QGIS Python launcher, for example on Windows: "C:\\Program Files\\QGIS 3.40.12\\bin\\python-qgis-ltr.bat" qgis_3d_logs.py It writes, beside the outputs of borehole_logs_3d.py: - a417_hillshade_5m.tif hillshade of the LIDAR DTM - a417_terrain_texture.tif the hillshade over OpenStreetMap, rendered once at 1 m as the terrain's texture - a417_dtm_cutaway.tif the DTM with the ground in front of a stepped cut line set to no data - a417_logs_3d.gpkg new layers 'columns' (one prism per GEOL interval in the holes near the cut) and 'cut_face' (the vertical face of the cut) - a417_logs_3d.qgz the project, with the cutaway DTM as its 3D terrain Open the project in QGIS Desktop, choose View > 3D Map Views > New 3D Map View, then run setup_3d_view() in the Python console (see the bottom of this file) to apply the view settings. Tested with QGIS 3.40.12-Bratislava. Data: A417 Missing Link ground investigation (BGS deposited data item 173984), Open Government Licence v3.0, contains British Geological Survey materials © UKRI 2026. Environment Agency LIDAR Composite DTM (Open Government Licence v3). Basemap © OpenStreetMap contributors (ODbL). """ from pathlib import Path import numpy as np from osgeo import gdal from qgis.core import ( Qgis, QgsAbstract3DSymbol, QgsApplication, QgsCategorizedSymbolRenderer, QgsCoordinateReferenceSystem, QgsDatumTransform, QgsFeature, QgsField, QgsFields, QgsGeometry, QgsLineSymbol, QgsMapRendererSequentialJob, QgsMapSettings, QgsMapThemeCollection, QgsPointXY, QgsProject, QgsProperty, QgsRasterDemTerrainProvider, QgsRasterLayer, QgsRectangle, QgsRendererCategory, QgsVector3D, QgsVectorFileWriter, QgsVectorLayer, ) from qgis._3d import ( QgsDirectionalLightSettings, QgsPhongMaterialSettings, QgsPolygon3DSymbol, QgsRuleBased3DRenderer, QgsVectorLayer3DRenderer, ) from qgis.PyQt.QtCore import QMetaType, QSize from qgis.PyQt.QtGui import QColor, QImage, QPainter gdal.UseExceptions() OUT = Path(__file__).parent / "output" # The cut follows terrain-tile edges, which QGIS meshes exactly at every level of detail, so the cut edge # stays straight. QGIS tiles the 3D extent from its north-west corner in squares of 2,200 m, halving at each # level, so y = 215700 and 215975 and x = 393075 are tile edges from the third level down. The cut runs # along the foot of the escarpment (CP102, CP202, CP104, CP105), then steps up to the crest (DSRC109, RC507, # DSRC302, OH405); the ground south of it is cut away. SECTION = [(391700.0, 215700.0), (393075.0, 215700.0), (393075.0, 215975.0), (393900.0, 215975.0)] SECTION_WIDTH = 60.0 # m; holes this close to the cut are projected onto its face CUT_LIP = 10.0 # m; two DTM cells of ground are kept beyond the cut, so the tile-edge vertices have data FACE_BASE = 100.0 # m OD, below the deepest log on the section COLUMN_RADIUS = 8.0 # m; the holes are 0.1-0.2 m wide, but must show across 2 km COLUMN_OFFSET = 6.0 # m in front of the face, so each log stands mostly proud of it EXAGGERATION = 2.0 # vertical exaggeration of the terrain, the logs and the face alike VIEW_EXTENT = (391700.0, 215300.0, 393900.0, 216800.0) SKY = "#a9c8e6" # plain background; the capture swaps it for a gradient FACE_COLOUR = "#8a7660" THEME = "Terrain texture" OSM_URL = "type=xyz&url=https://tile.openstreetmap.org/{z}/{x}/{y}.png&zmin=0&zmax=19" def south_of_cut(x, y, lip=0.0): """True where the ground is cut away; lip moves the cut that far into the cut-away side.""" (_, y_foot), (x_step, _), (_, y_crest), _ = SECTION return ((x < x_step + lip) & (y < y_foot - lip)) | ((x >= x_step + lip) & (y < y_crest - lip)) def section_line() -> QgsGeometry: return QgsGeometry.fromPolylineXY([QgsPointXY(x, y) for x, y in SECTION]) def phong(colour: str) -> QgsPhongMaterialSettings: """Mostly ambient, with a little diffuse light for shape, so the prisms keep the legend's colours.""" m = QgsPhongMaterialSettings() c = QColor(colour) m.setDiffuse(c) m.setAmbient(c) m.setSpecular(QColor("black")) m.setAmbientCoefficient(0.8) m.setDiffuseCoefficient(0.35) m.setSpecularCoefficient(0.0) return m def flat(colour: str) -> QgsPhongMaterialSettings: """An unlit material: ambient only, so every face of the shape has the same colour.""" m = QgsPhongMaterialSettings() m.setAmbient(QColor(colour)) m.setDiffuse(QColor("black")) m.setSpecular(QColor("black")) m.setAmbientCoefficient(1.0) m.setDiffuseCoefficient(0.0) m.setSpecularCoefficient(0.0) return m def write_layer(layer: QgsVectorLayer, name: str) -> QgsVectorLayer: opts = QgsVectorFileWriter.SaveVectorOptions() opts.layerName = name opts.actionOnExistingFile = QgsVectorFileWriter.ActionOnExistingFile.CreateOrOverwriteLayer gpkg = OUT / "a417_logs_3d.gpkg" err = QgsVectorFileWriter.writeAsVectorFormatV3(layer, str(gpkg), QgsProject.instance().transformContext(), opts) assert err[0] == QgsVectorFileWriter.WriterError.NoError, err return QgsVectorLayer(f"{gpkg}|layername={name}", name, "ogr") def make_rasters(dtm_path: Path) -> tuple[Path, Path]: hillshade = OUT / "a417_hillshade_5m.tif" gdal.DEMProcessing(str(hillshade), str(dtm_path), "hillshade", azimuth=315, altitude=40, computeEdges=True) ds = gdal.Open(str(dtm_path)) z = ds.GetRasterBand(1).ReadAsArray().astype("float32") gx, dx, _, gy, _, dy = ds.GetGeoTransform() xs = gx + dx * (np.arange(ds.RasterXSize) + 0.5) ys = gy + dy * (np.arange(ds.RasterYSize) + 0.5) xx, yy = np.meshgrid(xs, ys) z[south_of_cut(xx, yy, lip=CUT_LIP) | ~np.isfinite(z)] = -9999 cutaway = OUT / "a417_dtm_cutaway.tif" out = gdal.GetDriverByName("GTiff").CreateCopy(str(cutaway), ds) out.GetRasterBand(1).SetNoDataValue(-9999) out.GetRasterBand(1).WriteArray(z) out = None return hillshade, cutaway def make_columns(intervals: QgsVectorLayer) -> QgsVectorLayer: """One prism footprint per interval for the holes near the section, projected onto the cut face. base_mOD is where each prism starts and thickness how far it is extruded; offset_m records how far the hole was moved to reach the section line. """ fields = QgsFields() for name, t in [("LOCA_ID", QMetaType.Type.QString), ("unit", QMetaType.Type.QString), ("unit_name", QMetaType.Type.QString), ("colour", QMetaType.Type.QString), ("base_mOD", QMetaType.Type.Double), ("thickness", QMetaType.Type.Double), ("offset_m", QMetaType.Type.Double)]: fields.append(QgsField(name, t)) mem = QgsVectorLayer("Polygon?crs=EPSG:27700", "columns", "memory") mem.dataProvider().addAttributes(fields) mem.updateFields() line = section_line() feats = [] for f in intervals.getFeatures(): top = f.geometry().constGet().pointN(0) collar = QgsGeometry.fromPointXY(QgsPointXY(top.x(), top.y())) offset = collar.distance(line) if offset > SECTION_WIDTH: continue chain = line.lineLocatePoint(collar) on_line = line.interpolate(chain).asPoint() bearing = line.interpolateAngle(chain) # radians clockwise from north # right-hand normal of the line, i.e. towards the viewer on the cut-away side x = on_line.x() + COLUMN_OFFSET * np.cos(bearing) y = on_line.y() - COLUMN_OFFSET * np.sin(bearing) out = QgsFeature(fields) out.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(x, y)).buffer(COLUMN_RADIUS, 8)) out.setAttributes([f["LOCA_ID"], f["unit"], f["unit_name"], f["colour"], f["base_mOD"], round(f["top_mOD"] - f["base_mOD"], 2), round(offset, 1)]) feats.append(out) mem.dataProvider().addFeatures(feats) return write_layer(mem, "columns") def make_cut_face(dtm: QgsRasterLayer) -> QgsVectorLayer: """A vertical face just in front of the section line, from FACE_BASE up to the ground, in 5 m strips.""" mem = QgsVectorLayer("Polygon?crs=EPSG:27700&field=ground_mOD:double&field=height:double", "cut_face", "memory") _, (x_step, _), (_, y_crest), _ = SECTION # The crest face carries on 30 m behind the step, to back the corner where the terrain mesh thins out corner = QgsGeometry.fromPolylineXY([QgsPointXY(x_step - 30, y_crest), QgsPointXY(x_step, y_crest)]) feats = [] for line in (section_line(), corner): for s in np.arange(0, line.length(), 5.0): s1 = min(s + 5.0, line.length()) a, b = line.interpolate(s).asPoint(), line.interpolate(s1).asPoint() mid = line.interpolate((s + s1) / 2).asPoint() length = float(np.hypot(b.x() - a.x(), b.y() - a.y())) ux, uy = (b.y() - a.y()) / length, -(b.x() - a.x()) / length # unit normal towards the viewer # The 3D terrain meshes the DTM coarsely, so its edge can sit a little back from the cut; the face # reaches up to the highest ground in the 40 m behind the cut, which closes any gap z, ok = dtm.dataProvider().sample(mid, 1) if not ok: continue behind = [dtm.dataProvider().sample(QgsPointXY(mid.x() - ux * d, mid.y() - uy * d), 1) for d in range(5, 45, 5)] top = max([z] + [zb for zb, okb in behind if okb]) near, far = 0.2, 0.7 # m in front of the cut, so the face hides the DTM's ragged lip ring = [QgsPointXY(a.x() + ux * near, a.y() + uy * near), QgsPointXY(b.x() + ux * near, b.y() + uy * near), QgsPointXY(b.x() + ux * far, b.y() + uy * far), QgsPointXY(a.x() + ux * far, a.y() + uy * far)] f = QgsFeature(mem.fields()) f.setGeometry(QgsGeometry.fromPolygonXY([ring])) f.setAttributes([round(z, 2), round(top - FACE_BASE, 2)]) feats.append(f) mem.dataProvider().addFeatures(feats) return write_layer(mem, "cut_face") def make_texture(project: QgsProject, hillshade: QgsRasterLayer, osm: QgsRasterLayer) -> Path: """Render the hillshade over OpenStreetMap once, at 1 m per pixel, as the terrain's texture.""" x0, y0, x1, y1 = VIEW_EXTENT w, h = int(x1 - x0), int(y1 - y0) ms = QgsMapSettings() ms.setDestinationCrs(QgsCoordinateReferenceSystem("EPSG:27700")) ms.setTransformContext(project.transformContext()) ms.setExtent(QgsRectangle(*VIEW_EXTENT)) ms.setOutputSize(QSize(w, h)) ms.setBackgroundColor(QColor("white")) ms.setLayers([hillshade, osm]) job = QgsMapRendererSequentialJob(ms) job.start() job.waitForFinished() img = job.renderedImage().convertToFormat(QImage.Format.Format_RGB32) ptr = img.constBits() ptr.setsize(img.sizeInBytes()) bgra = np.frombuffer(ptr, np.uint8).reshape(h, img.bytesPerLine() // 4, 4)[:, :w] rgb = bgra[:, :, 2::-1].copy() # The lip of ground kept beyond the cut is painted the colour of the face xx, yy = np.meshgrid(x0 + np.arange(w) + 0.5, y1 - np.arange(h) - 0.5) rgb[south_of_cut(xx, yy)] = QColor(FACE_COLOUR).getRgb()[:3] path = OUT / "a417_terrain_texture.tif" ds = gdal.GetDriverByName("GTiff").Create(str(path), w, h, 3, gdal.GDT_Byte, ["COMPRESS=DEFLATE"]) ds.SetGeoTransform((x0, 1.0, 0.0, y1, 0.0, -1.0)) ds.SetProjection(QgsCoordinateReferenceSystem("EPSG:27700").toWkt()) for band in range(3): ds.GetRasterBand(band + 1).WriteArray(rgb[:, :, band]) ds = None return path def extruded(colour: str) -> QgsPolygon3DSymbol: """A polygon extruded from base_mOD (absolute, m OD) by its thickness field.""" sym = QgsPolygon3DSymbol() sym.setAltitudeClamping(Qgis.AltitudeClamping.Absolute) props = sym.dataDefinedProperties() props.setProperty(QgsAbstract3DSymbol.Property.Height, QgsProperty.fromExpression(f'"base_mOD" * {EXAGGERATION}')) props.setProperty(QgsAbstract3DSymbol.Property.ExtrusionHeight, QgsProperty.fromExpression(f'"thickness" * {EXAGGERATION}')) sym.setDataDefinedProperties(props) sym.setMaterialSettings(phong(colour)) return sym def main() -> None: qgs = QgsApplication([], False) qgs.initQgis() print("QGIS", Qgis.version()) project = QgsProject.instance() project.setCrs(QgsCoordinateReferenceSystem("EPSG:27700")) gpkg = OUT / "a417_logs_3d.gpkg" intervals = QgsVectorLayer(f"{gpkg}|layername=intervals", "GEOL intervals (3D lines)", "ogr") dtm = QgsRasterLayer(str(OUT / "a417_dtm_5m.tif"), "LIDAR DTM 5 m") print(f"intervals valid={intervals.isValid()} ({intervals.featureCount()}) | dtm valid={dtm.isValid()}") hs_path, cut_path = make_rasters(OUT / "a417_dtm_5m.tif") hillshade = QgsRasterLayer(str(hs_path), "LIDAR hillshade") hillshade.setBlendMode(QPainter.CompositionMode.CompositionMode_Multiply) hillshade.renderer().setOpacity(0.7) cutaway = QgsRasterLayer(str(cut_path), "LIDAR DTM, cut away") osm = QgsRasterLayer(OSM_URL, "OpenStreetMap", "wms") # Store the most accurate Web Mercator -> BNG operation, so QGIS doesn't ask which one to use for OSM bng, merc = QgsCoordinateReferenceSystem("EPSG:27700"), QgsCoordinateReferenceSystem("EPSG:3857") op = min((o for o in QgsDatumTransform.operations(merc, bng) if o.isAvailable and o.accuracy > 0), key=lambda o: o.accuracy) ctx = project.transformContext() ctx.addCoordinateOperation(merc, bng, op.proj) project.setTransformContext(ctx) # The 3D view drapes a local copy of the hillshade over OpenStreetMap, rendered once here texture = QgsRasterLayer(str(make_texture(project, hillshade, osm)), "Terrain texture (OSM and hillshade)") columns = make_columns(intervals) face = make_cut_face(dtm) print(f"hillshade valid={hillshade.isValid()} | cutaway valid={cutaway.isValid()} | " f"osm valid={osm.isValid()} | texture valid={texture.isValid()}") holes = sorted({f["LOCA_ID"] for f in columns.getFeatures()}) print(f"columns: {columns.featureCount()} prisms in {len(holes)} holes within {SECTION_WIDTH:g} m of the " f"section, radius {COLUMN_RADIUS:g} m | cut face: {face.featureCount()} strips from " f"{FACE_BASE:g} m OD to ground") print("holes on the section:", ", ".join(holes)) # Units in the order they appear, with the colours written by borehole_logs_3d.py units = {} for f in intervals.getFeatures(): units.setdefault(f["unit"], (f["unit_name"], f["colour"])) # 2D style for the map canvas: the same colours, categorised by unit cats = [QgsRendererCategory(u, QgsLineSymbol.createSimple({"color": c, "width": "1.2"}), n) for u, (n, c) in units.items()] intervals.setRenderer(QgsCategorizedSymbolRenderer("unit", cats)) # 3D style: one rule per unit, each an extruded prism in the unit's colour root = QgsRuleBased3DRenderer.Rule(None) for u, (name, colour) in units.items(): root.appendChild(QgsRuleBased3DRenderer.Rule(extruded(colour), f"\"unit\" = '{u}'", name)) columns.setRenderer3D(QgsRuleBased3DRenderer(root)) face.setRenderer3D(QgsVectorLayer3DRenderer(extruded_face())) # The cutaway DTM is the project's terrain, so every new 3D view starts on it terrain = QgsRasterDemTerrainProvider() terrain.setLayer(cutaway) project.elevationProperties().setTerrainProvider(terrain) layers = [columns, face, intervals, texture, hillshade, osm, cutaway, dtm] project.addMapLayers(layers, False) tree = project.layerTreeRoot() for layer in layers: tree.addLayer(layer) for layer in (intervals, texture, cutaway, dtm): tree.findLayer(layer.id()).setItemVisibilityChecked(False) # Only the pre-rendered texture is draped on the terrain theme = QgsMapThemeCollection.MapThemeRecord() theme.setLayerRecords([QgsMapThemeCollection.MapThemeLayerRecord(texture)]) project.mapThemeCollection().insert(THEME, theme) project.write(str(OUT / "a417_logs_3d.qgz")) print(f"saved output/a417_logs_3d.qgz: {len(units)} 3D rules ({', '.join(units)}), terrain {cutaway.name()}") qgs.exitQgis() def extruded_face() -> QgsPolygon3DSymbol: sym = QgsPolygon3DSymbol() sym.setAltitudeClamping(Qgis.AltitudeClamping.Absolute) sym.setOffset(FACE_BASE * EXAGGERATION) props = sym.dataDefinedProperties() # The extra 2 m tucks the top of the face under the edge of the ground props.setProperty(QgsAbstract3DSymbol.Property.ExtrusionHeight, QgsProperty.fromExpression(f'"height" * {EXAGGERATION} + 2')) sym.setDataDefinedProperties(props) sym.setMaterialSettings(flat(FACE_COLOUR)) return sym def setup_3d_view(settings, project=None) -> list[tuple[str, str]]: """Apply the view settings to an open 3D map view, e.g. iface.mapCanvases3D()[0].mapSettings(). Returns (name, colour) for the units in the logs in front of the cut, in stratigraphic order, for a legend. """ project = project or QgsProject.instance() layer = {lyr.name(): lyr for lyr in project.mapLayers().values()} settings.setExtent(QgsRectangle(*VIEW_EXTENT)) settings.setLayers([layer["columns"], layer["cut_face"]]) settings.setTerrainMapTheme(THEME) settings.setTerrainVerticalScale(EXAGGERATION) # the prisms and the face are scaled to match settings.setMapTileResolution(1024) settings.setMaxTerrainScreenError(0.5) settings.setMaxTerrainGroundError(0.5) settings.setTerrainShadingEnabled(False) # the hillshade already shades the ground settings.setBackgroundColor(QColor(SKY)) settings.setFieldOfView(35) key, fill = QgsDirectionalLightSettings(), QgsDirectionalLightSettings() key.setDirection(QgsVector3D(0.6, -1.0, -0.8)) key.setIntensity(0.9) fill.setDirection(QgsVector3D(-0.7, -0.4, 0.5)) fill.setIntensity(0.45) settings.setLightSources([key, fill]) order = ["TOPS", "MGR", "HEAD", "SLIP", "WHL", "HMB", "GOG", "FE", "SALS", "ASLS", "BLPL", "IOG", "BDS", "WHM", "MRB", "DYS", "CHAM", "NR"] seen = {} for f in layer["columns"].getFeatures(): seen[f["unit"]] = (f["unit_name"], f["colour"]) return [seen[u] for u in order if u in seen] if __name__ == "__main__": main()