Source code: 3D borehole logs from AGS4 data in PyVista and QGIS
Runnable Python scripts, requirements and console output for the tutorial “3D borehole logs from AGS4 data in PyVista and QGIS”.
Last updated
These are the exact scripts used to produce every figure and number in 3D borehole logs from AGS4 data in PyVista and QGIS. Download them, install the requirements into a fresh virtual environment and run them from this folder.
- a417_ags.py
- borehole_logs_3d.py
- qgis_3d_logs.py
- requirements.txt
- console.txt: the console output recorded when we ran it
python -m venv .venv.venv\Scripts\activate # macOS/Linux: source .venv/bin/activatepip install -r requirements.txtpython a417_ags.pya417_ags.py
"""Download, merge and repair the A417 Missing Link ground investigation AGS4 data.
Shared by every part of the TechnicalGIS series "AGS data in GIS"; each part keeps an identical copy.
The British Geological Survey (BGS) holds AGS files donated by the ground investigation industry.Its OGC API lists every borehole it can export (collection "agsboreholeindex"), and its AGSFile Utilities API exports those boreholes as AGS4 files, one ZIP per request. A request URLcan only carry so many borehole IDs, so we ask for them in batches of 40.
repair() applies the fixes worked through step by step in Part 1, so Parts 2 and 3 can startfrom a clean, merged set of tables with load_clean().
Data: A417 Missing Link ground investigation (BGS deposited data item 173984; client HighwaysEngland, contractor Geotechnical Engineering Ltd), from the BGS AGS download service under theOpen Government Licence v3.0. Contains British Geological Survey materials © UKRI 2026."""import ioimport timeimport zipfilefrom pathlib import Path
import pandas as pdimport requestsfrom python_ags4 import AGS4
ITEM = 173984 # BGS deposited data item: HE551505 A417 Missing Link GIINDEX = "https://ogcapi.bgs.ac.uk/collections/agsboreholeindex/items"EXPORT = "https://agsapi.bgs.ac.uk/v7/ags_export/"BATCH = 40UA = {"User-Agent": "TechnicalGIS-tutorial (https://www.technicalgis.com)"}
# Geological units met in the A417 logs, youngest at the top. Bedrock codes are BGS Lexicon codes# (GEOL_BGS); rank gives the expected order down a borehole, and an undivided group spans the ranks of its# formations (rank to last). Superficial deposits have rank 0.UNITS = pd.DataFrame([ # code, name, group, rank, colour ("TOPS", "Topsoil", "Superficial", 0, "#5b4636"), ("MGR", "Made ground", "Superficial", 0, "#9e9e9e"), ("HEAD", "Head", "Superficial", 0, "#d9b77e"), ("SLIP", "Landslide deposits", "Superficial", 0, "#c98c5a"), ("WHL", "White Limestone Formation", "Great Oolite Group", 1, "#f7f1bd"), ("HMB", "Hampen Formation", "Great Oolite Group", 2, "#e4d79c"), ("GOG", "Great Oolite Group (undivided)", "Great Oolite Group", 1, "#eee3a6"), ("FE", "Fuller's Earth Formation", "Great Oolite Group", 3, "#9fb08f"), ("SALS", "Salperton Limestone Formation", "Inferior Oolite Group", 4, "#f4cd72"), ("ASLS", "Aston Limestone Formation", "Inferior Oolite Group", 5, "#e9b157"), ("BLPL", "Birdlip Limestone Formation", "Inferior Oolite Group", 6, "#d98f3e"), ("IOG", "Inferior Oolite Group (undivided)", "Inferior Oolite Group", 4, "#e6a855"), ("BDS", "Bridport Sand Formation", "Lias Group", 7, "#d2c48c"), ("WHM", "Whitby Mudstone Formation", "Lias Group", 8, "#7d93ab"), ("MRB", "Marlstone Rock Formation", "Lias Group", 9, "#b0703f"), ("DYS", "Dyrham Formation", "Lias Group", 10, "#a3b8cc"), ("CHAM", "Charmouth Mudstone Formation", "Lias Group", 11, "#55698a"), ("NR", "No recovery or not classified", "Other", 0, "#f2f2f2"),], columns=["code", "name", "group", "rank", "colour"]).set_index("code")UNITS["last"] = UNITS["rank"]UNITS.loc[["GOG", "IOG"], "last"] = [3, 6]
# Where GEOL_BGS is blank, the project's own code (GEOL_GEOL) still tells us the unit.FROM_GEOL = {"TOP": "TOPS", "MG": "MGR", "HDD": "HEAD", "GOG": "GOG", "IOG": "IOG", "SLIP": "SLIP"}
def borehole_index(item=ITEM): """Every borehole in one BGS deposited data item, from the BGS OGC API (GeoJSON).""" r = requests.get(INDEX, params={"f": "json", "limit": 1000, "item_index_id": item}, headers=UA, timeout=120) r.raise_for_status() rows = [f["properties"] for f in r.json()["features"]] return pd.DataFrame(rows).sort_values("bgs_loca_id").reset_index(drop=True)
def _get(url, params, tries=5): for attempt in range(tries): try: r = requests.get(url, params=params, headers=UA, timeout=300) r.raise_for_status() return r except requests.RequestException: if attempt == tries - 1: raise time.sleep(10 * (attempt + 1))
def download(data_dir, item=ITEM): """Download the item's boreholes as AGS4 files, once. Returns the .ags paths in order.""" raw = Path(data_dir) / f"bgs_{item}" raw.mkdir(parents=True, exist_ok=True) if not list(raw.glob("*.ags")): ids = borehole_index(item).bgs_loca_id.tolist() for i in range(0, len(ids), BATCH): r = _get(EXPORT, {"bgs_loca_id": ";".join(ids[i:i + BATCH])}) with zipfile.ZipFile(io.BytesIO(r.content)) as z: for name in z.namelist(): (raw / f"batch{i // BATCH + 1:02d}_{Path(name).name}").write_bytes(z.read(name)) time.sleep(1) return sorted(raw.glob("*.ags"))
def merge(paths): """Read several AGS4 files from one project and stack them, group by group.
Each table keeps its UNIT and TYPE rows (from the first file) followed by every DATA row; rows repeated in every file (PROJ, TRAN, ABBR, DICT and so on) are kept once. """ tables, headings = {}, {} for p in paths: t, h = AGS4.AGS4_to_dataframe(str(p)) for group, df in t.items(): if group not in tables: tables[group], headings[group] = df, list(h[group]) continue headings[group] += [c for c in h[group] if c not in headings[group]] tables[group] = pd.concat([tables[group], df[df.HEADING == "DATA"]], ignore_index=True) for group, df in tables.items(): df = df.reindex(columns=headings[group]).fillna("") tables[group] = df.drop_duplicates(ignore_index=True) return tables, headings
def repair(tables, headings): """Apply the Part 1 fixes to merged A417 tables, in place, and return them.""" # Rule 1: a degree sign double-encoded as "Β°" (Greek capital beta + degree); the logs write "deg". for df in tables.values(): df.replace("Β°", "deg", regex=True, inplace=True)
# Rule 10a: three DICT headings defined twice; keep the later, archive descriptions. d = tables["DICT"] tables["DICT"] = d[~d.duplicated(["HEADING", "DICT_TYPE", "DICT_GRP", "DICT_HDNG"], keep="last") | (d.HEADING != "DATA")].reset_index(drop=True)
# Rules 7, 9, 19a, 19b: rename the BGS identifiers to valid user-defined headings and define them. new_dict = [] for group, old, new, desc in [("PROJ", "BGS_PROJ_ID", "PROJ_BGSI", "BGS export project identifier"), ("LOCA", "BGS_LOCA_ID", "LOCA_BGSI", "BGS borehole index identifier")]: tables[group] = tables[group].rename(columns={old: new}) tables[group].loc[tables[group].HEADING == "TYPE", new] = "X" headings[group] = [new if c == old else c for c in headings[group]] new_dict.append({"HEADING": "DATA", "DICT_TYPE": "HEADING", "DICT_GRP": group, "DICT_HDNG": new, "DICT_STAT": "OTHER", "DICT_DTYP": "X", "DICT_DESC": desc})
# The *_TEXT columns hold the depths and diameters exactly as originally written ("0.50"), where the # plain columns were reformatted ("0.5"). Put the original text back, restore the standard data type # and drop the extra column. std, _ = AGS4.AGS4_to_dataframe(str(Path(AGS4.__file__).parent / "Standard_dictionary_v4_0_3.ags")) dtype = std["DICT"].query("DICT_TYPE == 'HEADING'").set_index(["DICT_GRP", "DICT_HDNG"]).DICT_DTYP for group in list(tables): df = tables[group] for col in [c for c in headings[group] if c.endswith("_TEXT")]: base = col[:-5] data = df.HEADING == "DATA" same = pd.to_numeric(df.loc[data, base]).round(6).equals(pd.to_numeric(df.loc[data, col]).round(6)) assert same, f"{group}.{col} disagrees with {base}" df.loc[data, base] = df.loc[data, col] df.loc[df.HEADING == "TYPE", base] = dtype[(group, base)] tables[group] = df = df.drop(columns=col) headings[group].remove(col)
# Rule 10c: in the child groups LOCA_ID was exported with TYPE "X", so it no longer matches the # parent key in LOCA (TYPE "ID"). Restore the standard type. for group, df in tables.items(): if group != "LOCA" and "LOCA_ID" in df: df.loc[df.HEADING == "TYPE", "LOCA_ID"] = "ID"
tables["DICT"] = pd.concat([tables["DICT"], pd.DataFrame(new_dict)], ignore_index=True).fillna("")
# Rule 7 (only checkable once the headings above are defined): headings must follow the order of # the standard dictionary, with the file's own DICT headings after them in DICT order. own = tables["DICT"].query("HEADING == 'DATA' and DICT_TYPE == 'HEADING'") defined = pd.concat([std["DICT"].query("DICT_TYPE == 'HEADING'"), own]) for group in tables: listed = defined[defined.DICT_GRP == group].DICT_HDNG.drop_duplicates().tolist() order = {h: i for i, h in enumerate(listed)} headings[group] = ["HEADING"] + sorted(headings[group][1:], key=lambda h: order.get(h, len(order))) tables[group] = tables[group][headings[group]]
tables["TRAN"].loc[tables["TRAN"].HEADING == "DATA", "TRAN_DESC"] += \ "; three BGS exports merged and repaired (TechnicalGIS AGS data in GIS, Part 1)" return tables, headings
def write(tables, headings, ags_path, metadata): """Write a clean AGS4 file with the FILE folder AGS Rule 20 expects beside it.""" ags_path = Path(ags_path) ags_path.parent.mkdir(parents=True, exist_ok=True) AGS4.dataframe_to_AGS4(tables, headings, str(ags_path), warnings=False) fset = ags_path.parent / "FILE" / "BGSFileSet01" fset.mkdir(parents=True, exist_ok=True) (fset / "BGS_download_metadata.txt").write_bytes(Path(metadata).read_bytes()) return ags_path
def numeric(tables): """DATA rows only, with numbers as numbers. python-ags4 reads the TYPE row to decide; it leaves "U" (variable decimal places) as text, so convert those columns too.""" out = {} for g, df in tables.items(): types = df[df.HEADING == "TYPE"].iloc[0] n = AGS4.convert_to_numeric(df) for col in types.index[types == "U"]: n[col] = pd.to_numeric(n[col], errors="coerce") out[g] = n return out
def unit_code(geol): """One unit code per GEOL row: the BGS Lexicon code, else the project code, else NR.""" code = geol.GEOL_BGS.where(geol.GEOL_BGS != "", geol.GEOL_GEOL.map(FROM_GEOL)) return code.fillna("NR")
def in_order(codes): """True if a top-to-bottom list of bedrock unit codes never steps back to a younger unit.""" reached = 0 for code in codes: if UNITS.at[code, "last"] < reached: return False reached = max(reached, UNITS.at[code, "rank"]) return True
def load_clean(data_dir): """Download (once), merge and repair the A417 data. Returns numeric DataFrames by group.""" tables, headings = repair(*merge(download(data_dir))) out = numeric(tables) out["GEOL"] = out["GEOL"].sort_values(["LOCA_ID", "GEOL_TOP"], ignore_index=True) out["GEOL"]["UNIT"] = unit_code(out["GEOL"]) return outborehole_logs_3d.py
"""3D borehole logs from AGS4: strata cylinders coloured by GEOL, in PyVista and QGIS.
Part 3 of the TechnicalGIS series "AGS data in GIS". Turns every GEOL interval of the A417 MissingLink ground investigation into a short cylinder at the hole's British National Grid position and itstrue elevation (ground level from LOCA_GL minus depth), coloured by formation, and draws them underEnvironment Agency LIDAR terrain. Then it checks whether neighbouring logs tell the same story at thebase of the Inferior Oolite, finds a block that does not, and shows what its logs record. Finally itwrites LineStringZ intervals to a GeoPackage for the QGIS 3D map view (qgis_3d_logs.py).
TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/3d-borehole-logs-pyvista-qgis/
Data: A417 Missing Link ground investigation (BGS deposited data item 173984; client HighwaysEngland, contractor Geotechnical Engineering Ltd), from the BGS AGS download service under theOpen Government Licence v3.0. Contains British Geological Survey materials © UKRI 2026.Environment Agency LIDAR Composite DTM 1 m via the Defra WCS (Open Government Licence v3)."""import itertoolsimport timefrom importlib.metadata import versionfrom pathlib import Path
import numpy as npimport pandas as pdimport geopandas as gpdimport requestsimport rasteriofrom shapely.geometry import LineString, Pointimport pyvista as pvimport matplotlibmatplotlib.use("Agg")import matplotlib.pyplot as pltfrom matplotlib.colors import to_rgbfrom matplotlib.lines import Line2Dfrom matplotlib.patches import Patch
import a417_ags as A
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9, "figure.facecolor": "white", "axes.facecolor": "white"})pd.set_option("display.width", 120)
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)"}
AOI = (391_600, 212_800, 395_600, 216_800) # same 4 x 4 km window as Part 2RES = 5.0 # 5 m terrain is plenty under 112 holesVE = 3.0 # vertical exaggeration for the 3D viewsCLUSTER = (393_320, 216_060, 250) # centre and radius (m) of the close-up at the top of Crickley HillWCS = ("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")CREDIT = ("BGS AGS data © UKRI; LIDAR © Environment Agency (both OGL v3). " f"Vertical exaggeration x{VE:.0f}.")
def section(title): print(f"\n=== {title} ===")
def get(url, path, tries=6): """Download once and cache, backing off when the Defra WCS is busy.""" 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) return path
print(f"pyvista {pv.__version__} | vtk {version('vtk')} | geopandas {gpd.__version__} | " f"python-ags4 {version('python-ags4')}")
# ------------------------------------------------------------------ Step 1: intervals with elevationssection("Step 1: GEOL intervals with elevations")ags = A.load_clean(DATA)loca = ags["LOCA"].set_index("LOCA_ID")geol = ags["GEOL"].copy()geol["UNIT"] = geol.UNIT.fillna("NR")geol = geol.join(loca[["LOCA_NATE", "LOCA_NATN", "LOCA_GL"]], on="LOCA_ID")geol["top_mOD"] = geol.LOCA_GL - geol.GEOL_TOPgeol["base_mOD"] = geol.LOCA_GL - geol.GEOL_BASEgeol["thick"] = geol.GEOL_BASE - geol.GEOL_TOPgeol = geol[geol.thick > 0].reset_index(drop=True)print(f"{len(geol)} intervals in {geol.LOCA_ID.nunique()} holes | {geol.thick.sum():,.0f} m of logged ground | " f"elevations {geol.base_mOD.min():.1f} to {geol.top_mOD.max():.1f} m OD")used = [u for u in A.UNITS.index if u in set(geol.UNIT)]print(geol.groupby("UNIT").agg(intervals=("LOCA_ID", "size"), holes=("LOCA_ID", "nunique"), metres=("thick", "sum")) .reindex(used).round(1).to_string())
# ------------------------------------------------------------------ Step 2: terrainsection("Step 2: LIDAR terrain")X0, Y0, X1, Y1 = AOIf = get(f"{WCS}&subset=E({X0},{X1})&subset=N({Y0},{Y1})&SCALEFACTOR={1 / RES}", DATA / f"ea_dtm_{X0}_{Y0}_{X1 - X0}m_{RES:.0f}m.tif")with rasterio.open(f) as src: dtm = src.read(1).astype("float32") dtm[dtm < -1000] = np.nan left, top, dtm_transform = src.bounds.left, src.bounds.top, src.transformprint(f"DTM {dtm.shape[1]} x {dtm.shape[0]} cells of {RES:.0f} m, {np.nanmin(dtm):.1f} to {np.nanmax(dtm):.1f} m OD")
# Local coordinates keep the numbers small for the graphics card: metres east and north of the AOI corner.xs = left - X0 + RES / 2 + np.arange(dtm.shape[1]) * RESys = top - Y0 - RES / 2 - np.arange(dtm.shape[0]) * RESgx, gy = np.meshgrid(xs, ys)z = np.nan_to_num(dtm, nan=np.nanmin(dtm))terrain = pv.StructuredGrid(gx, gy, z * VE)terrain["elevation"] = z.ravel(order="F")
# ------------------------------------------------------------------ Step 3: one cylinder per intervalsection("Step 3: build the strata cylinders")
def cylinders(intervals, radius): """One cylinder per GEOL interval, with an RGB colour per face so the whole set is a single mesh.""" pieces = [] for r in intervals.itertuples(): c = pv.Cylinder(center=(r.LOCA_NATE - X0, r.LOCA_NATN - Y0, (r.top_mOD + r.base_mOD) / 2 * VE), direction=(0, 0, 1), radius=radius, height=r.thick * VE, resolution=24) rgb = (np.array(to_rgb(A.UNITS.colour[r.UNIT])) * 255).astype(np.uint8) c.cell_data["rgb"] = np.tile(rgb, (c.n_cells, 1)) pieces.append(c) return pv.merge(pieces)
logs = cylinders(geol, radius=22)print(f"{len(geol)} cylinders merged into one mesh: {logs.n_points:,} points, {logs.n_cells:,} faces")
def compose(img, units, path, title, note=""): """Put a PyVista screenshot beside a proper legend (PyVista's own legend tints the text).""" ink = np.where((img[..., :3] < 245).any(axis=2)) # trim the empty white margins img = img[max(ink[0].min() - 10, 0):ink[0].max() + 10, max(ink[1].min() - 10, 0):ink[1].max() + 10] h, w = img.shape[:2] fig = plt.figure(figsize=(11, 11 * h / w * 0.86)) ax = fig.add_axes([0, 0.035, 0.8, 0.905]) ax.imshow(img) ax.set_axis_off() fig.text(0.01, 0.975, title, fontsize=11, weight="bold", va="top") if note: fig.text(0.01, 0.005, note, fontsize=7, color="#444444") handles = [Patch(facecolor=A.UNITS.colour[u], edgecolor="#555555", lw=0.4, label=A.UNITS.name[u]) for u in units] fig.legend(handles=handles, loc="center left", bbox_to_anchor=(0.79, 0.5), fontsize=7.5, frameon=False, title="GEOL unit", title_fontsize=8, alignment="left") fig.savefig(path, dpi=150) plt.close(fig)
# ------------------------------------------------------------------ Step 4: the overviewsection("Step 4: render the overview")pl = pv.Plotter(off_screen=True, window_size=(1800, 1150))pl.set_background("white")pl.add_mesh(terrain, scalars="elevation", cmap="gist_earth", clim=[60, 330], opacity=0.4, smooth_shading=True, specular=0.0, ambient=0.3, show_scalar_bar=False)pl.add_mesh(terrain.contour(np.arange(90, 310, 10) * VE, scalars=terrain.points[:, 2]), color="#6f6a5f", line_width=1, opacity=0.5)pl.add_mesh(logs, scalars="rgb", rgb=True, smooth_shading=True, ambient=0.35)for lid in ["CP217", "RC516", "DSRC301"]: p = loca.loc[lid] pl.add_point_labels([(p.LOCA_NATE - X0, p.LOCA_NATN - Y0, (p.LOCA_GL + 30) * VE)], [lid], font_size=20, point_size=1, shape_opacity=0.8, always_visible=True, text_color="black", shape_color="white")pl.camera_position = [(-900, -1200, 1950 * VE / 3 + 1500), (2200, 2700, 200 * VE), (0, 0, 1)]pl.camera.zoom(1.5)img = pl.screenshot(return_img=True)pl.close()compose(img, used, OUT / "01-3d-borehole-logs.png", "A417 Missing Link: 112 borehole logs under LIDAR terrain, looking north-east up Crickley Hill", CREDIT + " 10 m ground contours.")print("rendered 01-3d-borehole-logs.png")
# ------------------------------------------------------------------ Step 5: do neighbouring logs agree?section("Step 5: do neighbouring logs agree? The base of the Inferior Oolite")rank = geol.UNIT.map(A.UNITS["rank"])contacts = []for lid, g in geol[rank > 0].groupby("LOCA_ID", sort=False): bds = g[g.UNIT == "BDS"] if bds.empty: continue b = bds.iloc[0] above = g[(g.UNIT.isin(["SALS", "ASLS", "BLPL", "IOG"])) & (g.GEOL_TOP < b.GEOL_TOP)] if above.empty: continue contacts.append(dict(LOCA_ID=lid, E=b.LOCA_NATE, N=b.LOCA_NATN, GL=b.LOCA_GL, depth=b.GEOL_TOP, z=b.top_mOD, above=above.iloc[-1].GEOL_DESC, below=b.GEOL_DESC))contacts = pd.DataFrame(contacts).set_index("LOCA_ID")print(f"{len(contacts)} holes log Inferior Oolite limestone resting on Bridport Sand; the contact lies at " f"{contacts.z.min():.1f} to {contacts.z.max():.1f} m OD")lith = contacts.below.str.extract(r"\b(SILT|CLAY|MUDSTONE|SILTSTONE|SAND|SANDSTONE|LIMESTONE)\b", expand=False)print("first Bridport Sand material below the contact:", lith.value_counts().to_dict())
pairs = []for a, b in itertools.combinations(contacts.itertuples(), 2): s = np.hypot(a.E - b.E, a.N - b.N) if s <= 100: pairs.append(dict(a=a.Index, b=b.Index, metres_apart=round(s), dz=round(a.z - b.z, 1)))pairs = pd.DataFrame(pairs).sort_values("metres_apart")pairs["agree_3m"] = pairs.dz.abs() <= 3print(f"\n{len(pairs)} pairs of holes within 100 m of each other; {pairs.agree_3m.sum()} agree within 3 m")print(pairs.to_string(index=False))
# A hole is "offset" if its contact sits more than 8 m from the median of its neighbours within 100 m.off = {}for lid, c in contacts.iterrows(): nb = contacts[(np.hypot(contacts.E - c.E, contacts.N - c.N) <= 100) & (contacts.index != lid)] if len(nb) >= 2: off[lid] = c.z - nb.z.median()off = pd.Series(off)offset = off[off.abs() > 8].index.tolist()print("\ncontact minus the median of its neighbours within 100 m (holes with 2+ neighbours):")print(off.round(1).sort_values().to_string())rest = pairs[~pairs.a.isin(offset) & ~pairs.b.isin(offset)]print(f"offset by more than 8 m: {', '.join(offset)}. Without them, {len(rest)} pairs agree to within " f"{rest.dz.abs().max():.1f} m (median {rest.dz.abs().median():.1f} m)")
# What do the offset holes record that their neighbours do not?steep = geol.GEOL_DESC.str.extract(r"bedded \((\d+)\s*-\s*(\d+)\s*deg\)", expand=True).astype(float)geol["bedding_dip"] = steep[1]geol["slickensides"] = geol.GEOL_DESC.str.contains("slickensid", case=False)cx, cy, cr = CLUSTERin_cluster = contacts[np.hypot(contacts.E - cx, contacts.N - cy) <= cr].index.tolist()ev = geol[geol.LOCA_ID.isin(in_cluster) & (geol.bedding_dip.notna() | geol.slickensides)]print(f"\n{len(in_cluster)} holes reach the contact within {cr} m of the close-up centre; their logs record:")for r in ev.itertuples(): what = f"bedding dipping up to {r.bedding_dip:.0f} deg" if pd.notna(r.bedding_dip) else "slickensides" print(f" {r.LOCA_ID:8s} {r.GEOL_TOP:6.2f}-{r.GEOL_BASE:6.2f} m {A.UNITS.name[r.UNIT]:28s} {what}")
# Strip logs of the cluster at true elevation, west to east.cl = contacts.loc[in_cluster].sort_values("E")fig, ax = plt.subplots(figsize=(10, 5.2))for i, (lid, c) in enumerate(cl.iterrows()): for r in geol[geol.LOCA_ID == lid].itertuples(): ax.bar(i, r.thick, bottom=r.base_mOD, width=0.62, color=A.UNITS.colour[r.UNIT], edgecolor="none") ax.plot([i - 0.42, i + 0.42], [c.z, c.z], color="black", lw=1.6) ax.text(i + 0.36, c.z + 0.8, f"{c.z:.1f}", fontsize=7, ha="left", va="bottom") e = ev[ev.LOCA_ID == lid] for r in e[e.bedding_dip.notna()].itertuples(): ax.plot(i, (r.top_mOD + r.base_mOD) / 2, marker="$\u2215$", ms=11, color="#b2182b") for r in e[e.slickensides].itertuples(): ax.plot(i, (r.top_mOD + r.base_mOD) / 2, marker="x", ms=7, mew=2, color="#b2182b")names = [f"{lid}\n{'offset' if lid in offset else ''}" for lid in cl.index]ax.set_xticks(range(len(cl)), names, fontsize=7, rotation=0)for t, lid in zip(ax.get_xticklabels(), cl.index): t.set_color("#b2182b" if lid in offset else "black")ax.set_ylim(170, 245)ax.set_xlim(-0.8, len(cl) - 0.2)ax.set_ylabel("elevation (m OD)")ax.set_title(f"Top of Crickley Hill: {len(cl)} holes within {cr} m, west to east. Black bars mark the base of the " "Inferior Oolite (m OD)", fontsize=9)ax.spines[["top", "right"]].set_visible(False)handles = [Patch(facecolor=A.UNITS.colour[u], label=A.UNITS.name[u]) for u in used if u in set(geol.loc[geol.LOCA_ID.isin(cl.index), "UNIT"])]handles += [Line2D([], [], color="black", lw=1.6, label="base of Inferior Oolite"), Line2D([], [], ls="", marker="$\u2215$", ms=9, color="#b2182b", label="log records bedding dipping 30-50°"), Line2D([], [], ls="", marker="x", ms=6, mew=1.6, color="#b2182b", label="log records slickensides")]ax.legend(handles=handles, fontsize=6.5, loc="center left", bbox_to_anchor=(1.0, 0.5), frameon=False)fig.savefig(OUT / "02-contact-check-strip-logs.png")plt.close(fig)
# The same holes in 3D, with the contact drawn as a black disc and the ground as 2 m contours.box = (cx - X0 - 330, cx - X0 + 330, cy - Y0 - 330, cy - Y0 + 330, -1e6, 1e6)near = loca[np.hypot(loca.LOCA_NATE - cx, loca.LOCA_NATN - cy) <= 330].indexclose_logs = cylinders(geol[geol.LOCA_ID.isin(near)], radius=5)local = terrain.clip_box(box, invert=False)pl = pv.Plotter(off_screen=True, window_size=(1800, 1150))pl.set_background("white")pl.add_mesh(local.contour(np.arange(200, 300, 2) * VE, scalars=local.points[:, 2]), color="#8f8a7e", line_width=1.2)pl.add_mesh(close_logs, scalars="rgb", rgb=True, smooth_shading=True, ambient=0.25)for lid, c in contacts.loc[contacts.index.isin(near)].iterrows(): pl.add_mesh(pv.Disc(center=(c.E - X0, c.N - Y0, c.z * VE), inner=0, outer=11, normal=(0, 0, 1), c_res=32), color="#b2182b" if lid in offset else "black")labels = contacts.loc[contacts.index.isin(near)]pl.add_point_labels(np.c_[labels.E - X0, labels.N - Y0, (labels.GL + 6) * VE], [f"{lid} {c.z:.0f}" for lid, c in labels.iterrows()], font_size=15, point_size=1, shape_opacity=0.75, always_visible=True, text_color="black", shape_color="white")pl.camera_position = [(cx - X0 - 150, cy - Y0 - 900, 360 * VE), (cx - X0, cy - Y0, 205 * VE), (0, 0, 1)]img = pl.screenshot(return_img=True)pl.close()compose(img, [u for u in used if u in set(geol.loc[geol.LOCA_ID.isin(near), "UNIT"])], OUT / "03-contact-check-3d.png", "Top of Crickley Hill, looking north: discs mark the base of the Inferior Oolite (red = offset holes)", CREDIT + " 2 m ground contours. Labels give the contact in m OD.")print("rendered 03-contact-check-3d.png")
# ------------------------------------------------------------------ Step 6: LineStringZ for QGIS 3Dsection("Step 6: write LineStringZ intervals for QGIS")lines = gpd.GeoDataFrame( {"LOCA_ID": geol.LOCA_ID, "unit": geol.UNIT, "unit_name": geol.UNIT.map(A.UNITS.name), "colour": geol.UNIT.map(A.UNITS.colour), "top_depth": geol.GEOL_TOP, "base_depth": geol.GEOL_BASE, "top_mOD": geol.top_mOD.round(2), "base_mOD": geol.base_mOD.round(2), "GEOL_DESC": geol.GEOL_DESC}, geometry=[LineString([(r.LOCA_NATE, r.LOCA_NATN, r.top_mOD), (r.LOCA_NATE, r.LOCA_NATN, r.base_mOD)]) for r in geol.itertuples()], crs="EPSG:27700")collars = gpd.GeoDataFrame(loca[["LOCA_TYPE", "LOCA_GL", "LOCA_FDEP"]].reset_index(), geometry=[Point(r.LOCA_NATE, r.LOCA_NATN, r.LOCA_GL) for r in loca.itertuples()], crs="EPSG:27700")contact_pts = gpd.GeoDataFrame(contacts[["depth", "z"]].assign(offset=contacts.index.isin(offset)).reset_index(), geometry=[Point(c.E, c.N, c.z) for c in contacts.itertuples()], crs="EPSG:27700")gpkg = OUT / "a417_logs_3d.gpkg"gpkg.unlink(missing_ok=True)lines.to_file(gpkg, layer="intervals", driver="GPKG")collars.to_file(gpkg, layer="collars", driver="GPKG")contact_pts.to_file(gpkg, layer="base_inferior_oolite", driver="GPKG")check = gpd.read_file(gpkg, layer="intervals")print(f"wrote {gpkg.as_posix()}: 'intervals' {len(check)} {check.geom_type.iloc[0]}s with z: {check.has_z.all()}, " f"'collars' {len(collars)} points, 'base_inferior_oolite' {len(contact_pts)} points")with rasterio.open(OUT / "a417_dtm_5m.tif", "w", driver="GTiff", width=dtm.shape[1], height=dtm.shape[0], count=1, dtype="float32", crs="EPSG:27700", transform=dtm_transform, nodata=np.nan, compress="deflate") as dst: dst.write(dtm, 1)print("terrain for QGIS 3D: a417_dtm_5m.tif")print("figures written:", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))qgis_3d_logs.py
"""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 withthe QGIS Python launcher, for example on Windows: "C:\\Program Files\\QGIS 3.40.12\\bin\\python-qgis-ltr.bat" qgis_3d_logs.pyIt 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 terrainOpen the project in QGIS Desktop, choose View > 3D Map Views > New 3D Map View, then runsetup_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 Licencev3.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 npfrom osgeo import gdalfrom 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, QSizefrom 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 faceCUT_LIP = 10.0 # m; two DTM cells of ground are kept beyond the cut, so the tile-edge vertices have dataFACE_BASE = 100.0 # m OD, below the deepest log on the sectionCOLUMN_RADIUS = 8.0 # m; the holes are 0.1-0.2 m wide, but must show across 2 kmCOLUMN_OFFSET = 6.0 # m in front of the face, so each log stands mostly proud of itEXAGGERATION = 2.0 # vertical exaggeration of the terrain, the logs and the face alikeVIEW_EXTENT = (391700.0, 215300.0, 393900.0, 216800.0)SKY = "#a9c8e6" # plain background; the capture swaps it for a gradientFACE_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()