"""The same cross-section in QGIS: a print layout with a map and an elevation profile along section A-I. Part 4 of the TechnicalGIS series "AGS data in GIS". Run cross_section.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_section.py It loads output/a417_section.gpkg and the 2 m DTM, sets their elevation properties so QGIS's elevation profile tools can draw them (the DTM as a ground line, each GEOL interval as a point extruded upwards by its thickness), checks QGIS's own ground profile against ground_profile.csv, builds a print layout and saves output/a417_section.qgz. The layout is exported to output/_capture/06-qgis-section-layout.png. Open the project in QGIS Desktop and choose View > Elevation Profile, then Capture Curve from Feature on the section line, to explore the same section interactively. Tested with QGIS 3.40.12-Bratislava. Data: A417 Missing Link ground investigation (BGS deposited data items 173984 and 173975), Open Government Licence v3.0, contains British Geological Survey materials © UKRI 2026. Environment Agency LIDAR Composite DTM (Open Government Licence v3). """ import csv from pathlib import Path from qgis.core import ( Qgis, QgsApplication, QgsCategorizedSymbolRenderer, QgsCoordinateReferenceSystem, QgsFillSymbol, QgsHillshadeRenderer, QgsLayout, QgsLayoutExporter, QgsLayoutItemElevationProfile, QgsLayoutItemLabel, QgsLayoutItemMap, QgsLayoutPoint, QgsLayoutSize, QgsLineSymbol, QgsMapLayerElevationProperties, QgsMarkerSymbol, QgsPalLayerSettings, QgsPoint, QgsPrintLayout, QgsProfileExporter, QgsProfileRequest, QgsProject, QgsProperty, QgsRasterLayer, QgsRectangle, QgsRendererCategory, QgsSingleSymbolRenderer, QgsSymbolLayer, QgsTextFormat, QgsUnitTypes, QgsVectorLayer, QgsVectorLayerSimpleLabeling, ) from qgis.PyQt.QtGui import QColor, QFont HERE = Path(__file__).parent OUT = HERE / "output" CREDIT = ("Contains British Geological Survey materials © UKRI (OGL v3). Environment Agency LIDAR Composite DTM " "© Environment Agency and database right (OGL v3).") HALF = 60 # corridor half-width, as in cross_section.py Z_RANGE = (95, 305) def main() -> None: qgs = QgsApplication([], False) qgs.initQgis() print("QGIS", Qgis.version()) project = QgsProject.instance() project.setCrs(QgsCoordinateReferenceSystem("EPSG:27700")) gpkg = OUT / "a417_section.gpkg" layer = {name: QgsVectorLayer(f"{gpkg}|layername={name}", title, "ogr") for name, title in [ ("section_line", "Section line A-I"), ("corridor", f"Corridor ±{HALF} m"), ("bends", "Bends"), ("holes", "Projected holes"), ("intervals", "GEOL intervals"), ("contacts", "Picked contacts")]} dtm = QgsRasterLayer(str(OUT / "a417_section_dtm_2m.tif"), "LIDAR DTM 2 m") for name, lyr in layer.items(): print(f" {name:13s} valid={lyr.isValid()} features={lyr.featureCount()}") print(f" dtm valid={dtm.isValid()} {dtm.width()} x {dtm.height()}") # --- 2D styles for the map dtm.setRenderer(QgsHillshadeRenderer(dtm.dataProvider(), 1, 315, 45)) layer["corridor"].setRenderer(QgsSingleSymbolRenderer(QgsFillSymbol.createSimple( {"color": "255,212,0,70", "outline_style": "no"}))) layer["section_line"].setRenderer(QgsSingleSymbolRenderer(QgsLineSymbol.createSimple( {"color": "#c62828", "width": "0.6"}))) layer["holes"].setRenderer(QgsSingleSymbolRenderer(QgsMarkerSymbol.createSimple( {"name": "circle", "color": "#111111", "outline_color": "white", "size": "1.6"}))) layer["bends"].setRenderer(QgsSingleSymbolRenderer(QgsMarkerSymbol.createSimple( {"name": "circle", "color": "0,0,0,0", "outline_style": "no", "size": "0.1"}))) label = QgsPalLayerSettings() label.fieldName = "letter" fmt = QgsTextFormat() fmt.setFont(QFont("Arial", 9, QFont.Bold)) fmt.setSize(9) fmt.setColor(QColor("#c62828")) fmt.buffer().setEnabled(True) fmt.buffer().setSize(0.8) label.setFormat(fmt) label.dist = 1.5 layer["bends"].setLabeling(QgsVectorLayerSimpleLabeling(label)) layer["bends"].setLabelsEnabled(True) # --- Elevation properties: what the elevation profile draws, and how dtm.elevationProperties().setEnabled(True) dtm.elevationProperties().profileLineSymbol().setColor(QColor("#3b3024")) dtm.elevationProperties().profileLineSymbol().setWidth(0.5) iv = layer["intervals"] units = {} for f in iv.getFeatures(): units.setdefault(f["unit"], (f["unit_name"], f["colour"])) iv.setRenderer(QgsCategorizedSymbolRenderer("unit", [ QgsRendererCategory(u, QgsMarkerSymbol.createSimple({"color": c, "size": "1.2"}), n) for u, (n, c) in units.items()])) ep = iv.elevationProperties() ep.setClamping(Qgis.AltitudeClamping.Absolute) # z is the base of the interval in m OD ep.setExtrusionEnabled(True) # ... and the line rises by the interval's thickness ep.dataDefinedProperties().setProperty(QgsMapLayerElevationProperties.Property.ExtrusionHeight, QgsProperty.fromField("thickness")) ep.setRespectLayerSymbology(False) strip = QgsLineSymbol.createSimple({"width": "2.2", "capstyle": "flat"}) strip.symbolLayer(0).setDataDefinedProperty(QgsSymbolLayer.Property.StrokeColor, QgsProperty.fromField("colour")) ep.setProfileLineSymbol(strip) cp = layer["contacts"].elevationProperties() cp.setClamping(Qgis.AltitudeClamping.Absolute) cp.setRespectLayerSymbology(False) cp.setProfileMarkerSymbol(QgsMarkerSymbol.createSimple( {"name": "circle", "color": "white", "outline_color": "#111111", "size": "1.8"})) project.addMapLayers([*layer.values(), dtm], False) tree = project.layerTreeRoot() for lyr in (layer["bends"], layer["contacts"], layer["holes"], layer["section_line"], layer["intervals"], layer["corridor"], dtm): tree.addLayer(lyr) tree.findLayer(layer["intervals"].id()).setItemVisibilityChecked(False) tree.findLayer(layer["contacts"].id()).setItemVisibilityChecked(False) # --- QGIS's own ground profile along the line, compared with the one cross_section.py sampled curve = next(layer["section_line"].getFeatures()).geometry().constGet().clone() # OH413 sits on the very end of the line; carry the curve 30 m on so QGIS keeps it (chainage is unchanged) (x1, y1), (x2, y2) = [(curve.xAt(i), curve.yAt(i)) for i in (curve.numPoints() - 2, curve.numPoints() - 1)] leg = ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5 curve.addVertex(QgsPoint(x2 + 30 * (x2 - x1) / leg, y2 + 30 * (y2 - y1) / leg)) request = QgsProfileRequest(curve) request.setCrs(project.crs()) request.setTransformContext(project.transformContext()) request.setTolerance(1) exporter = QgsProfileExporter([dtm], request, Qgis.ProfileExportType.DistanceVsElevationTable) exporter.run() table = exporter.toLayers()[0] qgis_profile = sorted((f["distance"], f["elevation"]) for f in table.getFeatures()) with open(OUT / "ground_profile.csv", newline="") as fh: rows = [(float(r["chainage"]), float(r["ground_mOD"])) for r in csv.DictReader(fh)] py_ch, py_z = zip(*rows) def interp(x): i = min(range(len(py_ch)), key=lambda k: abs(py_ch[k] - x)) return py_z[i] diffs = [z - interp(d) for d, z in qgis_profile if d <= py_ch[-1]] diffs.sort(key=abs) on_line = [z for d, z in qgis_profile if d <= py_ch[-1]] print(f"QGIS ground profile: {len(on_line)} points, {min(on_line):.1f} to {max(on_line):.1f} m OD " f"| against ground_profile.csv: median |difference| {abs(diffs[len(diffs) // 2]):.2f} m, " f"95% within {abs(diffs[int(len(diffs) * 0.95)]):.2f} m") # --- Print layout: the map above, the elevation profile below layout = QgsPrintLayout(project) layout.initializeDefaults() layout.setName("Section A-I") page = layout.pageCollection().page(0) page.setPageSize(QgsLayoutSize(297, 210, QgsUnitTypes.LayoutMillimeters)) project.layoutManager().addLayout(layout) def text(s, x, y, w, h, size, bold=False, colour="#222222"): item = QgsLayoutItemLabel(layout) tf = QgsTextFormat() tf.setFont(QFont("Arial", size, QFont.Bold if bold else QFont.Normal)) tf.setSize(size) tf.setColor(QColor(colour)) item.setTextFormat(tf) item.setText(s) item.attemptMove(QgsLayoutPoint(x, y, QgsUnitTypes.LayoutMillimeters)) item.attemptResize(QgsLayoutSize(w, h, QgsUnitTypes.LayoutMillimeters)) layout.addLayoutItem(item) return item text("Section A-I up the Cotswold escarpment at Crickley Hill, drawn by QGIS", 10, 6, 277, 8, 13, bold=True) ext = layer["corridor"].extent() ext.grow(160) mp = QgsLayoutItemMap(layout) mp.setLayers([layer["bends"], layer["holes"], layer["section_line"], layer["corridor"], dtm]) mp.attemptMove(QgsLayoutPoint(10, 16, QgsUnitTypes.LayoutMillimeters)) mp.attemptResize(QgsLayoutSize(120, 72, QgsUnitTypes.LayoutMillimeters)) mp.setExtent(ext) mp.setFrameEnabled(True) layout.addLayoutItem(mp) text("The map frame shows the 2 m LIDAR hillshade, the section line and the ±60 m corridor. " "The elevation profile item below uses the same line as its profile curve with a 60 m tolerance, so it " "draws every GEOL interval in the corridor at its true elevation: each interval is a point at its base, " "extruded upwards by its thickness and coloured from the 'colour' field. White dots are the contacts " "picked by cross_section.py. Vertical exaggeration is about ×5.5.", 140, 18, 147, 50, 9) prof = QgsLayoutItemElevationProfile(layout) prof.setCrs(project.crs()) prof.setProfileCurve(curve.clone()) prof.setTolerance(HALF) prof.setLayers([dtm, layer["intervals"], layer["contacts"]]) # later layers draw on top prof.attemptMove(QgsLayoutPoint(10, 92, QgsUnitTypes.LayoutMillimeters)) prof.attemptResize(QgsLayoutSize(277, 102, QgsUnitTypes.LayoutMillimeters)) plot = prof.plot() plot.setXMinimum(-20) plot.setXMaximum(curve.length() + 20) plot.setYMinimum(Z_RANGE[0]) plot.setYMaximum(Z_RANGE[1]) plot.xAxis().setGridIntervalMajor(500) plot.xAxis().setGridIntervalMinor(100) plot.xAxis().setLabelInterval(500) plot.yAxis().setGridIntervalMajor(50) plot.yAxis().setGridIntervalMinor(10) plot.yAxis().setLabelInterval(50) layout.addLayoutItem(prof) text("distance along the section (m)", 120, 194, 80, 5, 8) text("m OD", 2, 92, 10, 5, 8) text(CREDIT, 10, 200, 277, 6, 7, colour="#444444") project.write(str(OUT / "a417_section.qgz")) print(f"saved output/a417_section.qgz: {len(units)} units in the profile ({', '.join(units)})") cap = OUT / "_capture" cap.mkdir(exist_ok=True) settings = QgsLayoutExporter.ImageExportSettings() settings.dpi = 200 settings.generateWorldFile = False result = QgsLayoutExporter(layout).exportToImage(str(cap / "06-qgis-section-layout.png"), settings) print("layout exported to output/_capture/06-qgis-section-layout.png" if result == QgsLayoutExporter.Success else f"layout export failed ({result})") qgs.exitQgis() if __name__ == "__main__": main()