Source code: Boreholes on a map: AGS4 locations in British National Grid with QGIS

"""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 AGS
File Utilities API exports those boreholes as AGS4 files, one ZIP per request. A request URL
can 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 start
from a clean, merged set of tables with load_clean().

Data: A417 Missing Link ground investigation (BGS deposited data item 173984; client Highways
England, contractor Geotechnical Engineering Ltd), from the BGS AGS download service under the
Open Government Licence v3.0. Contains British Geological Survey materials © UKRI 2026.
"""
import io
import time
import zipfile
from pathlib import Path

import pandas as pd
import requests
from python_ags4 import AGS4

ITEM = 173984                       # BGS deposited data item: HE551505 A417 Missing Link GI
INDEX = "https://ogcapi.bgs.ac.uk/collections/agsboreholeindex/items"
EXPORT = "https://agsapi.bgs.ac.uk/v7/ags_export/"
BATCH = 40
UA = {"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 out

"""Boreholes on a map: AGS4 LOCA to points in British National Grid, checked against LIDAR.

Part 2 of the TechnicalGIS series "AGS data in GIS". Turns the LOCA group of the A417 Missing
Link ground investigation into a GeoDataFrame in British National Grid (EPSG:27700), shows the
two coordinate mistakes that send boreholes to the wrong place, summarises each hole's geology
from GEOL (depth to rock, the unit at rockhead, the deepest unit reached), checks every ground
level against Environment Agency LIDAR, draws two maps and writes a GeoPackage with a
ready-made HTML popup for QGIS (styled by qgis_borehole_map.py).

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/ags4-boreholes-map-qgis/

Data: A417 Missing Link ground investigation (BGS deposited data item 173984; client Highways
England, contractor Geotechnical Engineering Ltd), from the BGS AGS download service under the
Open 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 time
from html import escape
from importlib.metadata import version
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import requests
import rasterio
from rasterio.transform import from_origin
from pyproj import Transformer
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
from matplotlib.lines import Line2D

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)             # 4 x 4 km around the holes, British National Grid
RES = 2.0                                               # metres per cell requested from the WCS
TILE = 1000                                             # metres per WCS request (500 x 500 cells)
WCS = ("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")

# One colour per kind of hole, shared with the QGIS style in qgis_borehole_map.py
TYPE_COLOURS = {"RC": "#1f78b4", "TP": "#e6550d", "CPT": "#31a354", "RO": "#756bb1", "CP": "#d62728"}


def section(title):
    print(f"\n=== {title} ===")


def get(url, path, tries=6):
    """Download once and cache. The Defra WCS rate-limits bursts (HTTP 429), so back off and retry."""
    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)
        time.sleep(1)
    return path


print(f"geopandas {gpd.__version__} | pyproj {version('pyproj')} | python-ags4 {version('python-ags4')}")

# ------------------------------------------------------------------ Step 1: LOCA to points
section("Step 1: LOCA to points in British National Grid")
ags = A.load_clean(DATA)
loca = ags["LOCA"]
holes = gpd.GeoDataFrame(loca, geometry=gpd.points_from_xy(loca.LOCA_NATE, loca.LOCA_NATN), crs="EPSG:27700")
print(f"{len(holes)} holes | CRS {holes.crs.name} ({holes.crs.to_string()})")
print(f"extent E {holes.total_bounds[0]:,.0f} to {holes.total_bounds[2]:,.0f}, "
      f"N {holes.total_bounds[1]:,.0f} to {holes.total_bounds[3]:,.0f}")
centre = holes.dissolve().centroid.to_crs("EPSG:4326").iloc[0]
print(f"centre in WGS84: {centre.y:.5f} N, {abs(centre.x):.5f} W (Crickley Hill and Birdlip)")

# ------------------------------------------------------------------ Step 2: two classic mistakes
section("Step 2: what goes wrong with the wrong CRS")
first = loca.iloc[0]
to_wgs = Transformer.from_crs("EPSG:27700", "EPSG:4326", always_xy=True)
lon, lat = to_wgs.transform(first.LOCA_NATN, first.LOCA_NATE)       # easting and northing swapped
print(f"{first.LOCA_ID} with easting and northing swapped lands at {lat:.3f} N, {lon:.3f} W: "
      f"in the sea off north Wales, {np.hypot(first.LOCA_NATE - first.LOCA_NATN, first.LOCA_NATN - first.LOCA_NATE) / 1000:,.0f} km away")
wrong = gpd.GeoDataFrame(geometry=gpd.points_from_xy(loca.LOCA_NATE, loca.LOCA_NATN), crs="EPSG:4326")
print(f"tagged as EPSG:4326 instead, x runs to {wrong.total_bounds[2]:,.0f} 'degrees'; "
      f"valid longitudes stop at 180, so nothing can be drawn")

# ------------------------------------------------------------------ Step 3: geology per hole
section("Step 3: summarise GEOL for each hole")
geol = ags["GEOL"].copy()
geol["rank"] = geol.UNIT.map(A.UNITS["rank"])
rock = geol[geol["rank"] > 0].groupby("LOCA_ID")
summary = pd.DataFrame({
    "rockhead_depth": rock.GEOL_TOP.min(),
    "rockhead_unit": rock.UNIT.first(),
    "deepest_unit": rock.UNIT.last(),
})
abbr = ags["ABBR"].query("ABBR_HDNG == 'LOCA_TYPE'").set_index("ABBR_CODE").ABBR_DESC
water = ags["WSTG"].groupby("LOCA_ID").WSTG_DPTH.min().rename("first_water_strike")
counts = pd.DataFrame({"n_samples": ags["SAMP"].groupby("LOCA_ID").size(),
                       "n_spt": ags["ISPT"].groupby("LOCA_ID").size()})


def log_html(g):
    """The hole's log as a small HTML table, one row per unit, for the QGIS popup."""
    runs = g.groupby((g.UNIT != g.UNIT.shift()).cumsum()).agg(
        unit=("UNIT", "first"), top=("GEOL_TOP", "min"), base=("GEOL_BASE", "max"))
    rows = "".join(f'<tr><td style="background:{A.UNITS.colour[r.unit]};width:10px"></td>'
                   f"<td>{r.top:.2f}&ndash;{r.base:.2f} m</td><td>{escape(A.UNITS.name[r.unit])}</td></tr>"
                   for r in runs.itertuples())
    return f'<table style="font-size:11px;border-collapse:collapse">{rows}</table>'


holes = holes.join(summary, on="LOCA_ID").join(water, on="LOCA_ID").join(counts, on="LOCA_ID")
holes[["n_samples", "n_spt"]] = holes[["n_samples", "n_spt"]].fillna(0).astype(int)
holes["type_name"] = holes.LOCA_TYPE.map(abbr)
holes["rockhead_mOD"] = (holes.LOCA_GL - holes.rockhead_depth).round(2)
holes["rockhead_name"] = holes.rockhead_unit.map(A.UNITS.name)
holes["log_html"] = holes.LOCA_ID.map(geol.groupby("LOCA_ID").apply(log_html, include_groups=False))
print(holes.type_name.value_counts().to_string())
print(f"\n{holes.rockhead_depth.notna().sum()} holes reach bedrock; depth to rock "
      f"{holes.rockhead_depth.min():.2f} to {holes.rockhead_depth.max():.2f} m (median {holes.rockhead_depth.median():.2f} m)")
print("unit at rockhead:", holes.rockhead_unit.value_counts().to_dict())
print(holes.loc[holes.LOCA_ID.isin(["DSRC415", "RC516", "CP217"]),
                ["LOCA_ID", "type_name", "LOCA_GL", "LOCA_FDEP", "rockhead_depth", "rockhead_mOD", "rockhead_name",
                 "first_water_strike", "n_samples"]].to_string(index=False))

# ------------------------------------------------------------------ Step 4: ground level against LIDAR
section("Step 4: check LOCA_GL against Environment Agency LIDAR")
X0, Y0, X1, Y1 = AOI
dtm = np.full((int((Y1 - Y0) / RES), int((X1 - X0) / RES)), np.nan, dtype="float32")
for x in range(X0, X1, TILE):
    for y in range(Y0, Y1, TILE):
        f = get(f"{WCS}&subset=E({x},{x + TILE})&subset=N({y},{y + TILE})&SCALEFACTOR={1 / RES}",
                DATA / f"ea_dtm_{x}_{y}_{TILE}m_tile_{RES:.0f}m.tif")
        with rasterio.open(f) as src:
            a = src.read(1).astype("float32")
            a[a < -1000] = np.nan
            r, c = int(round((Y1 - src.bounds.top) / RES)), int(round((src.bounds.left - X0) / RES))
            dtm[r:r + a.shape[0], c:c + a.shape[1]] = a
transform = from_origin(X0, Y1, RES, RES)
print(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 | "
      f"no-data cells: {np.isnan(dtm).sum():,}")
profile = dict(driver="GTiff", width=dtm.shape[1], height=dtm.shape[0], count=1, dtype="float32",
               crs="EPSG:27700", transform=transform, nodata=np.nan, compress="deflate")
with rasterio.open(OUT / "a417_dtm_2m.tif", "w", **profile) as dst:
    dst.write(dtm, 1)
ls = LightSource(azdeg=315, altdeg=45)
shade = ls.hillshade(np.nan_to_num(dtm, nan=np.nanmean(dtm)), vert_exag=2, dx=RES, dy=RES)
with rasterio.open(OUT / "a417_hillshade.tif", "w", **{**profile, "dtype": "uint8", "nodata": None}) as dst:
    dst.write((shade * 255).astype("uint8"), 1)

col = ((holes.LOCA_NATE - X0) / RES).astype(int)
row = ((Y1 - holes.LOCA_NATN) / RES).astype(int)
holes["lidar_gl"] = dtm[row, col].round(2)
holes["gl_minus_lidar"] = (holes.LOCA_GL - holes.lidar_gl).round(2)
d = holes.gl_minus_lidar
print(f"LOCA_GL minus LIDAR: median {d.median():+.2f} m | middle half {d.quantile(.25):+.2f} to {d.quantile(.75):+.2f} m | "
      f"within 0.5 m: {(d.abs() <= 0.5).sum()} of {d.notna().sum()} | within 1 m: {(d.abs() <= 1).sum()}")
odd = holes.loc[d.abs() > 1, ["LOCA_ID", "type_name", "LOCA_CLST", "LOCA_GL", "lidar_gl", "gl_minus_lidar"]]
print("more than 1 m apart:")
print(odd.sort_values("gl_minus_lidar").to_string(index=False))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3.8), gridspec_kw={"width_ratios": [1.1, 1]})
for t, c in TYPE_COLOURS.items():
    s = holes[holes.LOCA_TYPE == t]
    ax1.scatter(s.lidar_gl, s.LOCA_GL, s=12, color=c, label=f"{t} {abbr.get(t, '')}", edgecolor="white", lw=0.3)
lim = [holes.lidar_gl.min() - 5, holes.lidar_gl.max() + 5]
ax1.plot(lim, lim, color="#555555", lw=0.8)
ax1.set(xlim=lim, ylim=lim, xlabel="LIDAR DTM at the hole (m OD)", ylabel="LOCA_GL in the AGS file (m OD)",
        title="Ground level: AGS against LIDAR")
ax1.legend(fontsize=7, frameon=False)
ax2.hist(d.clip(-3, 3), bins=np.arange(-3, 3.01, 0.2), color="#3f7f9f")
ax2.axvline(0, color="#555555", lw=0.8)
ax2.set(xlabel="LOCA_GL minus LIDAR (m, clipped at ±3)", ylabel="holes",
        title=f"{(d.abs() <= 0.5).sum()} of {d.notna().sum()} holes agree within 0.5 m")
for ax in (ax1, ax2):
    ax.spines[["top", "right"]].set_visible(False)
fig.savefig(OUT / "01-ground-level-check.png")
plt.close(fig)

# ------------------------------------------------------------------ Step 5: maps
section("Step 5: maps")
extent = (X0, X1, Y0, Y1)


def basemap(ax):
    ax.imshow(shade, cmap="gray", extent=extent, vmin=0, vmax=1.15)
    ax.imshow(dtm, cmap="terrain", extent=extent, alpha=0.25, vmin=np.nanmin(dtm) - 40, vmax=np.nanmax(dtm) + 40)
    ax.set_xlim(X0, X1); ax.set_ylim(Y0, Y1)
    ax.set_xticks(range(X0, X1 + 1, 1000), [f"{v / 1000:.0f}" for v in range(X0, X1 + 1, 1000)], fontsize=7)
    ax.set_yticks(range(Y0, Y1 + 1, 1000), [f"{v / 1000:.0f}" for v in range(Y0, Y1 + 1, 1000)], fontsize=7)
    ax.set_xlabel("easting (km, EPSG:27700)", fontsize=8); ax.set_ylabel("northing (km)", fontsize=8)
    ax.text(0.01, 0.01, "Contains Environment Agency information © Environment Agency and database right; "
            "BGS © UKRI (OGL v3)", transform=ax.transAxes, fontsize=5.5, color="#333333",
            bbox=dict(facecolor="white", alpha=0.7, lw=0, pad=1.5))


fig, ax = plt.subplots(figsize=(7.5, 7.5))
basemap(ax)
for t, c in TYPE_COLOURS.items():
    s = holes[holes.LOCA_TYPE == t]
    ax.scatter(s.geometry.x, s.geometry.y, s=22, color=c, edgecolor="white", lw=0.6,
               label=f"{abbr.get(t, t)} ({len(s)})", zorder=3)
placed = []  # label the deepest holes first, skipping any within 180 m of a label already placed
for _, h in holes[holes.LOCA_FDEP >= 40].sort_values("LOCA_FDEP", ascending=False).iterrows():
    x, y = h.geometry.x, h.geometry.y
    if any(np.hypot(x - px, y - py) < 180 for px, py in placed):
        continue
    placed.append((x, y))
    ax.annotate(h.LOCA_ID, (x, y), xytext=(4, 3), textcoords="offset points", fontsize=6,
                bbox=dict(facecolor="white", alpha=0.75, lw=0, pad=0.8), zorder=4)
ax.legend(loc="upper right", fontsize=7, framealpha=0.9)
ax.set_title("A417 Missing Link: 112 exploratory holes (labels on the deepest, 40 m or more)")
fig.savefig(OUT / "02-borehole-map.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(7.5, 7.5))
basemap(ax)
none = holes[holes.rockhead_unit.isna()]
ax.scatter(none.geometry.x, none.geometry.y, s=14, facecolor="white", edgecolor="#555555", lw=0.6, zorder=3)
order = [u for u in A.UNITS.index if u in set(holes.rockhead_unit.dropna())]
for u in order:
    s = holes[holes.rockhead_unit == u]
    ax.scatter(s.geometry.x, s.geometry.y, s=30, color=A.UNITS.colour[u], edgecolor="#222222", lw=0.6, zorder=3)
handles = [Line2D([], [], marker="o", ls="", markersize=6, markerfacecolor=A.UNITS.colour[u], markeredgecolor="#222222",
                  label=f"{A.UNITS.name[u]} ({(holes.rockhead_unit == u).sum()})") for u in order]
handles.append(Line2D([], [], marker="o", ls="", markersize=5, markerfacecolor="white", markeredgecolor="#555555",
                      label=f"bedrock not reached ({len(none)})"))
ax.legend(handles=handles, loc="upper right", fontsize=6.5, framealpha=0.9, title="first bedrock unit", title_fontsize=7)
ax.set_title("The formation at rockhead in each hole")
fig.savefig(OUT / "03-rockhead-unit-map.png")
plt.close(fig)

# ------------------------------------------------------------------ Step 6: GeoPackage for QGIS
section("Step 6: write a GeoPackage for QGIS")
keep = ["LOCA_ID", "LOCA_TYPE", "type_name", "LOCA_CLST", "LOCA_GL", "LOCA_FDEP", "lidar_gl", "gl_minus_lidar",
        "rockhead_depth", "rockhead_mOD", "rockhead_unit", "rockhead_name", "deepest_unit", "first_water_strike",
        "n_samples", "n_spt", "LOCA_BGSI", "LOCA_REM", "log_html", "geometry"]
gpkg = OUT / "a417_boreholes.gpkg"
gpkg.unlink(missing_ok=True)
holes[keep].to_file(gpkg, layer="boreholes", driver="GPKG")
check = gpd.read_file(gpkg, layer="boreholes")
print(f"wrote {gpkg.as_posix()}: layer 'boreholes', {len(check)} points, {len(check.columns) - 1} fields, "
      f"CRS {check.crs.to_epsg()}")
print("rasters for the QGIS basemap: a417_dtm_2m.tif, a417_hillshade.tif")
print("figures written:", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))

"""Style the A417 boreholes in QGIS: categorised symbols, labels, an HTML map tip and a saved project.

Part 2 of the TechnicalGIS series "AGS data in GIS". Run boreholes_map.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_borehole_map.py
It writes output/a417_boreholes.qgz (open it in QGIS Desktop) and renders output/04-qgis-borehole-map.png.
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).
"""
from pathlib import Path

from qgis.core import (
    Qgis,
    QgsApplication,
    QgsCategorizedSymbolRenderer,
    QgsCoordinateReferenceSystem,
    QgsMapRendererParallelJob,
    QgsMapSettings,
    QgsMarkerSymbol,
    QgsPalLayerSettings,
    QgsProject,
    QgsRasterLayer,
    QgsRectangle,
    QgsRendererCategory,
    QgsTextBufferSettings,
    QgsTextFormat,
    QgsVectorLayer,
    QgsVectorLayerSimpleLabeling,
)
from qgis.PyQt.QtCore import QPoint, QSize, Qt
from qgis.PyQt.QtGui import QColor, QFont, QPainter

OUT = Path(__file__).parent / "output"
CREDIT = ("Contains Environment Agency information © Environment Agency and database right; "
          "BGS © UKRI (OGL v3)")

# Same colours as boreholes_map.py
TYPE_COLOURS = {"RC": "#1f78b4", "TP": "#e6550d", "CPT": "#31a354", "RO": "#756bb1", "CP": "#d62728"}

# The map tip: a header from LOCA, a few summary fields, then the hole's log as the HTML table
# that boreholes_map.py stored in the log_html field.
MAP_TIP = """<div style="font-family:sans-serif;font-size:12px;min-width:260px">
<b style="font-size:14px">[% "LOCA_ID" %]</b> &middot; [% "type_name" %]<br>
Ground level [% format_number("LOCA_GL", 2) %] m OD (LIDAR [% format_number("lidar_gl", 2) %])<br>
Final depth [% format_number("LOCA_FDEP", 2) %] m
[% CASE WHEN "rockhead_depth" IS NOT NULL
   THEN '&middot; rock at ' || format_number("rockhead_depth", 2) || ' m (' || "rockhead_name" || ')' END %]<br>
[% CASE WHEN "first_water_strike" IS NOT NULL
   THEN 'First water strike ' || format_number("first_water_strike", 2) || ' m<br>' END %]
[% "n_samples" %] samples, [% "n_spt" %] SPTs
<hr style="margin:4px 0">[% "log_html" %]</div>"""


def main() -> None:
    qgs = QgsApplication([], False)
    qgs.initQgis()
    print("QGIS", Qgis.version())
    project = QgsProject.instance()
    project.setCrs(QgsCoordinateReferenceSystem("EPSG:27700"))

    shade = QgsRasterLayer(str(OUT / "a417_hillshade.tif"), "LIDAR hillshade")
    holes = QgsVectorLayer(f"{OUT / 'a417_boreholes.gpkg'}|layername=boreholes", "A417 exploratory holes", "ogr")
    print(f"hillshade valid={shade.isValid()} | boreholes valid={holes.isValid()}, "
          f"{holes.featureCount()} features, CRS {holes.crs().authid()}")

    # One symbol per LOCA_TYPE, labelled with the description from the ABBR group.
    names = {f["LOCA_TYPE"]: f["type_name"] for f in holes.getFeatures()}
    categories = []
    for code, colour in TYPE_COLOURS.items():
        symbol = QgsMarkerSymbol.createSimple({"name": "circle", "color": colour, "outline_color": "white",
                                               "outline_width": "0.3", "size": "2.6"})
        categories.append(QgsRendererCategory(code, symbol, f"{code}: {names.get(code, code)}"))
    holes.setRenderer(QgsCategorizedSymbolRenderer("LOCA_TYPE", categories))

    # Label every hole with its LOCA_ID, haloed so it reads over the hillshade.
    fmt = QgsTextFormat()
    fmt.setFont(QFont("Arial"))
    fmt.setSize(7)
    buffer = QgsTextBufferSettings()
    buffer.setEnabled(True)
    buffer.setSize(0.8)
    buffer.setColor(QColor("white"))
    fmt.setBuffer(buffer)
    label = QgsPalLayerSettings()
    label.fieldName = "LOCA_ID"
    label.placement = Qgis.LabelPlacement.AroundPoint
    label.setFormat(fmt)
    holes.setLabeling(QgsVectorLayerSimpleLabeling(label))
    holes.setLabelsEnabled(True)

    holes.setDisplayExpression('"LOCA_ID"')
    holes.setMapTipTemplate(MAP_TIP)
    project.addMapLayers([shade, holes], False)
    root = project.layerTreeRoot()
    root.addLayer(holes)                    # first in the tree draws on top
    root.addLayer(shade)
    # Credit the data on every map: as layer attribution, and as QGIS Desktop's copyright decoration
    # (View > Decorations > Copyright Label), which the project remembers.
    shade.serverProperties().setAttribution("© Environment Agency, Open Government Licence v3")
    holes.serverProperties().setAttribution("Contains British Geological Survey materials © UKRI, Open Government Licence v3")
    project.writeEntry("CopyrightLabel", "/Enabled", True)
    project.writeEntry("CopyrightLabel", "/Label", CREDIT)
    project.write(str(OUT / "a417_boreholes.qgz"))
    print("saved output/a417_boreholes.qgz with the style, labels and map tip")

    # Render the busiest part of the site, the top of Crickley Hill, at a scale where labels fit.
    settings = QgsMapSettings()
    settings.setLayers([holes, shade])
    settings.setDestinationCrs(project.crs())
    settings.setExtent(QgsRectangle(392_700, 215_500, 393_900, 216_500))
    settings.setOutputSize(QSize(1500, 1250))
    settings.setOutputDpi(150)
    settings.setBackgroundColor(QColor("white"))
    job = QgsMapRendererParallelJob(settings)
    job.start()
    job.waitForFinished()
    image = job.renderedImage()
    painter = QPainter(image)
    painter.setFont(QFont("Arial", 7))
    box = painter.fontMetrics().boundingRect(CREDIT).adjusted(-6, -4, 6, 4)
    box.moveBottomLeft(QPoint(8, image.height() - 8))
    painter.fillRect(box, QColor(255, 255, 255, 220))
    painter.setPen(QColor("#333333"))
    painter.drawText(box, Qt.AlignCenter, CREDIT)
    painter.end()
    image.save(str(OUT / "04-qgis-borehole-map.png"))
    print(f"rendered 04-qgis-borehole-map.png at 1:{settings.scale():,.0f}; render errors: {len(job.errors())}")
    qgs.exitQgis()


if __name__ == "__main__":
    main()