Source code: Two-step floating catchment area (2SFCA) in Python: spatial access to GP practices in Cornwall

"""Spatial accessibility to GP practices in Cornwall with the two-step floating catchment area
method (2SFCA) and its enhanced, distance-decay version (E2SFCA), on OSMnx drive times.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/spatial-accessibility-2sfca-python/
Part 1 of "Healthcare accessibility analysis in Python". Part 2 (location-allocation-python) reuses
the cached road graph and demand points written to data/ by this script.

Data: NHS England "Patients Registered at a GP Practice" (October 2026) and "General Practice
Workforce" (August 2026) (OGL v3); postcodes.io (ONS Postcode Directory, OGL v3); ONS Census 2021
TS001 via Nomis and ONS Open Geography boundaries and centroids (OGL v3); OpenStreetMap via OSMnx (ODbL).
"""
import json
import sys
import time
import warnings
import zipfile
from importlib.metadata import version
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import requests
import osmnx as ox
import osmium
import shapely
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import dijkstra
from access import Access, weights
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm, BoundaryNorm
from matplotlib.lines import Line2D

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9,
                     "figure.facecolor": "white", "axes.facecolor": "white"})

DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
CRS = "EPSG:27700"          # British National Grid, metres
LAD = "E06000052"           # Cornwall (the Isles of Scilly, E06000053, have no road link and are left out)
BUFFER_M = 25_000           # supply AND demand are taken from 25 km beyond the boundary (edge effects)
SNAP_KMH = 20               # speed for the straight-line hop from a point to its nearest road node
UA = {"User-Agent": "TechnicalGIS-tutorial (https://www.technicalgis.com)"}

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
LAD_URL = f"{ONS}/Local_Authority_Districts_May_2026_Boundaries_UK_BGC/FeatureServer/0/query"
OA_PWC_URL = f"{ONS}/OA_December_2021_EW_PWC_V4/FeatureServer/0/query"
OA_URL = f"{ONS}/Output_Areas_2021_EW_BGC_V2/FeatureServer/0/query"
RUC_URL = f"{ONS}/LSOA_2021_EW_BSC_V4_RUC/FeatureServer/0/query"
TS001 = "https://www.nomisweb.co.uk/output/census/2021/census2021-ts001.zip"
GP_MAP = "https://files.digital.nhs.uk/8E/6B4899/gp-reg-pat-prac-map.zip"          # practice -> postcode, ICB
GP_LIST = "https://files.digital.nhs.uk/FB/7851F6/gp-reg-pat-prac-all.zip"         # registered patients
GP_WORKFORCE = "https://files.digital.nhs.uk/49/2ABF5A/GPWPracticeCSV.082026.zip"  # GP FTE by practice
GEOFABRIK = [f"https://download.geofabrik.de/europe/united-kingdom/england/{c}-latest.osm.pbf"
             for c in ("cornwall", "devon")]
NOT_DRIVE = {"abandoned", "bridleway", "bus_guideway", "construction", "corridor", "cycleway", "elevator",
             "escalator", "footway", "no", "path", "pedestrian", "planned", "platform", "proposed", "raceway",
             "razed", "rest_area", "service", "services", "steps", "track"}
NOT_SERVICE = {"alley", "driveway", "emergency_access", "parking", "parking_aisle", "private"}


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


def download(url, path):
    if not path.exists():
        r = requests.get(url, timeout=600, headers=UA)
        r.raise_for_status()
        path.write_bytes(r.content)
    return path


def arcgis(url, path, **params):
    """Query an ArcGIS FeatureServer layer page by page and cache the result as GeoJSON."""
    if not path.exists():
        feats, offset = [], 0
        while True:
            q = {"outFields": "*", "outSR": 27700, "f": "geojson", "resultOffset": offset,
                 "resultRecordCount": 2000, "where": "1=1", **params}
            r = requests.post(url, data=q, timeout=300, headers=UA)
            r.raise_for_status()
            page = r.json()["features"]
            feats += page
            if len(page) < 2000:
                break
            offset += len(page)
        path.write_text(json.dumps({"type": "FeatureCollection", "features": feats}))
    return gpd.read_file(path).set_crs(CRS, allow_override=True)


def envelope(x0, y0, x1, y1):
    return {"geometry": f"{x0},{y0},{x1},{y1}", "geometryType": "esriGeometryEnvelope", "inSR": 27700,
            "spatialRel": "esriSpatialRelIntersects"}


def drivable(tags):
    """The same rules as OSMnx's network_type="drive" filter, applied to a way's tags."""
    hw = tags.get("highway")
    return (hw is not None and hw not in NOT_DRIVE and tags.get("area") != "yes"
            and tags.get("motor_vehicle") != "no" and tags.get("motorcar") != "no"
            and tags.get("access") != "private" and tags.get("service") not in NOT_SERVICE)


def build_drive_graph(area, folder):
    """Drive network for `area` (BNG polygon) from Geofabrik extracts, via an .osm file read by OSMnx."""
    poly = gpd.GeoSeries([area], crs=CRS).to_crs(4326).iloc[0]
    osm_file = folder / "drive_network.osm"
    keep_nodes, ways = set(), {}
    for url in GEOFABRIK:
        pbf = download(url, folder / url.rsplit("/", 1)[1])
        for w in osmium.FileProcessor(str(pbf), osmium.osm.NODE | osmium.osm.WAY).with_locations():
            if not w.is_way() or w.id in ways or not drivable(w.tags):
                continue
            xy = np.array([(n.lon, n.lat) for n in w.nodes])
            if shapely.contains_xy(poly, xy[:, 0], xy[:, 1]).any():   # keep ways that touch the study area
                ways[w.id] = ([n.ref for n in w.nodes], dict(w.tags))
                keep_nodes.update(n.ref for n in w.nodes)
    written = set()
    with osmium.SimpleWriter(str(osm_file), overwrite=True) as writer:
        for url in GEOFABRIK:
            for n in osmium.FileProcessor(str(folder / url.rsplit("/", 1)[1]), osmium.osm.NODE):
                if n.id in keep_nodes and n.id not in written:
                    writer.add_node(n)
                    written.add(n.id)
        for wid, (refs, tags) in sorted(ways.items()):
            writer.add_way(osmium.osm.mutable.Way(id=wid, version=1, nodes=refs, tags=tags))
    print(f"wrote {len(ways):,} drivable ways and {len(written):,} nodes to {osm_file.name}")
    G = ox.graph_from_xml(osm_file, simplify=True)
    G = ox.routing.add_edge_speeds(G)           # km/h from maxspeed, else the mean for that road type
    G = ox.routing.add_edge_travel_times(G)     # seconds
    G = ox.truncate.largest_component(G, strongly=True)
    return ox.project_graph(G, to_crs=CRS)


def read_zip_csv(path, name_part):
    with zipfile.ZipFile(path) as z:
        name = next(n for n in z.namelist() if name_part in n and n.endswith(".csv"))
        return pd.read_csv(z.open(name), low_memory=False, encoding="utf-8-sig", encoding_errors="replace")


section("versions")
print(f"python {sys.version.split()[0]} | osmnx {version('osmnx')} | networkx {version('networkx')} | "
      f"access {version('access')} | geopandas {version('geopandas')} | scipy {version('scipy')} | "
      f"numpy {version('numpy')} | pandas {version('pandas')}")

# ---------------------------------------------------------------- Step 1: study area and demand
section("study area")
lads = arcgis(LAD_URL, DATA / "lads.geojson", **envelope(130000, 0, 300000, 140000))
cornwall = lads[lads.LAD26CD == LAD].geometry.iloc[0]
study = cornwall.buffer(BUFFER_M)
print(f"Cornwall: {cornwall.area / 1e6:,.0f} km2 | with a {BUFFER_M / 1000:.0f} km buffer: {study.area / 1e6:,.0f} km2 "
      f"(much of it sea)")

pop = read_zip_csv(download(TS001, DATA / "census2021-ts001.zip"), "-oa.csv")
pop = pop.rename(columns={"geography code": "OA21CD", "Residence type: Total; measures: Value": "pop"})
pwc = arcgis(OA_PWC_URL, DATA / "oa_pwc.geojson", **envelope(*study.bounds))
pwc = pwc[pwc.within(study)].merge(pop[["OA21CD", "pop"]], on="OA21CD").reset_index(drop=True)
oas = arcgis(OA_URL, DATA / "cornwall_oas.geojson", where="LSOA21NM LIKE 'Cornwall %'")
pwc["in_cornwall"] = pwc.OA21CD.isin(oas.OA21CD)
print(f"{len(pwc):,} output areas in the buffered study area, {pwc['pop'].sum():,} residents (Census 2021)")
print(f"  in Cornwall: {pwc.in_cornwall.sum():,} OAs, {pwc.loc[pwc.in_cornwall, 'pop'].sum():,} residents; "
      f"in the buffer (Devon): {(~pwc.in_cornwall).sum():,} OAs, {pwc.loc[~pwc.in_cornwall, 'pop'].sum():,} residents")

# ---------------------------------------------------------------- Step 2: supply (GP practices)
section("GP practices")
pmap = read_zip_csv(download(GP_MAP, DATA / "gp-reg-pat-prac-map.zip"), "map")
plist = read_zip_csv(download(GP_LIST, DATA / "gp-reg-pat-prac-all.zip"), "all")
gpw = read_zip_csv(download(GP_WORKFORCE, DATA / "GPWPracticeCSV.082026.zip"), "Detailed")
print(f"England: {len(pmap):,} practices in the October 2026 list; {len(gpw):,} in the August 2026 workforce file")

prac = pmap[pmap.PRACTICE_POSTCODE.str.match(r"^(TR|PL|EX)\d")][
    ["PRACTICE_CODE", "PRACTICE_NAME", "PRACTICE_POSTCODE", "ICB_NAME"]].copy()
prac = prac.merge(plist[["CODE", "NUMBER_OF_PATIENTS"]].rename(columns={"CODE": "PRACTICE_CODE"}), how="left")
gpw["TOTAL_GP_FTE"] = pd.to_numeric(gpw.TOTAL_GP_FTE, errors="coerce")
prac = prac.merge(gpw[["PRAC_CODE", "TOTAL_GP_FTE"]].rename(columns={"PRAC_CODE": "PRACTICE_CODE"}), how="left")

geo_cache = DATA / "postcodes.json"
coords = json.loads(geo_cache.read_text()) if geo_cache.exists() else {}
todo = [p for p in prac.PRACTICE_POSTCODE.unique() if p not in coords]
for i in range(0, len(todo), 100):   # postcodes.io takes up to 100 postcodes per request
    r = requests.post("https://api.postcodes.io/postcodes", json={"postcodes": todo[i:i + 100]}, headers=UA, timeout=60)
    r.raise_for_status()
    for res in r.json()["result"]:
        rr = res["result"]
        coords[res["query"]] = [rr["eastings"], rr["northings"]] if rr else None
geo_cache.write_text(json.dumps(coords))
prac["xy"] = prac.PRACTICE_POSTCODE.map(coords)
print(f"{len(prac)} practices with TR, PL or EX postcodes; {prac.xy.isna().sum()} postcodes not found")
prac = prac.dropna(subset=["xy"])
prac = gpd.GeoDataFrame(prac, geometry=gpd.points_from_xy(prac.xy.str[0], prac.xy.str[1]), crs=CRS)
prac = prac[prac.within(study)].drop(columns="xy").reset_index(drop=True)
prac["in_cornwall"] = prac.within(cornwall)
print(f"{len(prac)} practices in the buffered study area: {prac.in_cornwall.sum()} in Cornwall, "
      f"{(~prac.in_cornwall).sum()} in Devon")
missing = prac.TOTAL_GP_FTE.isna() | (prac.TOTAL_GP_FTE <= 0)
print(f"practices with no GP FTE in the workforce file: {missing.sum()}"
      + (f" ({', '.join(prac.loc[missing, 'PRACTICE_NAME'].str.title())})" if missing.any() else ""))
# Fill any gaps from the list size at the study-area ratio of patients per GP FTE
ratio = prac.loc[~missing, "NUMBER_OF_PATIENTS"].sum() / prac.loc[~missing, "TOTAL_GP_FTE"].sum()
prac.loc[missing, "TOTAL_GP_FTE"] = prac.loc[missing, "NUMBER_OF_PATIENTS"] / ratio
c = prac[prac.in_cornwall]
print(f"Cornwall: {c.NUMBER_OF_PATIENTS.sum():,} registered patients, {c.TOTAL_GP_FTE.sum():,.1f} GP FTE "
      f"({c.TOTAL_GP_FTE.sum() / c.NUMBER_OF_PATIENTS.sum() * 1000:.3f} per 1,000 registered patients); "
      f"{ratio:,.0f} patients per GP FTE across the study area")
print(f"Cornwall practice size: median {c.NUMBER_OF_PATIENTS.median():,.0f} patients, "
      f"median {c.TOTAL_GP_FTE.median():.1f} GP FTE")

# ---------------------------------------------------------------- Step 3: road network and travel times
section("road network (OpenStreetMap -> OSMnx)")
graph_file = DATA / "cornwall_drive.graphml"
t0 = time.time()
if not graph_file.exists():
    G = build_drive_graph(study, DATA)
    ox.save_graphml(G, graph_file)
G = ox.load_graphml(graph_file)
print(f"{G.number_of_nodes():,} nodes, {G.number_of_edges():,} edges (largest strongly connected component); "
      f"{time.time() - t0:.0f} s")
speeds = pd.Series([d["speed_kph"] for *_, d in G.edges(data=True)])
print(f"edge speeds: median {speeds.median():.0f} km/h; {(speeds >= 96).mean():.1%} of edges at 96 km/h or more")

nodes = list(G.nodes)
pos = {n: i for i, n in enumerate(nodes)}
edges = pd.DataFrame([(pos[a], pos[b], d["travel_time"]) for a, b, d in G.edges(data=True)], columns=["u", "v", "t"])
edges = edges.groupby(["u", "v"], as_index=False).t.min()   # keep the fastest of any parallel edges
W = csr_matrix((edges.t, (edges.u, edges.v)), shape=(len(nodes), len(nodes)))


def snap(points):
    """Nearest graph node for each point, the hop to it in minutes and its length in metres."""
    ids, dist = ox.distance.nearest_nodes(G, points.geometry.x, points.geometry.y, return_dist=True)
    dist = np.asarray(dist)
    return np.array([pos[i] for i in ids]), dist / 1000 / SNAP_KMH * 60, dist


d_node, d_hop, d_dist = snap(pwc)
far = d_dist > 2000                                   # e.g. Lundy: no road link to the mainland
if far.any():
    print(f"dropping {far.sum()} OA(s) more than 2 km from the network: {', '.join(pwc.OA21CD[far])} "
          f"({pwc['pop'][far].sum()} residents, in Cornwall: {pwc.in_cornwall[far].sum()})")
    pwc = pwc[~far].reset_index(drop=True)
    d_node, d_hop, d_dist = snap(pwc)
s_node, s_hop, s_dist = snap(prac)
print(f"snap distance: OA centroids median {np.median(d_dist):.0f} m, 95th percentile {np.percentile(d_dist, 95):.0f} m, "
      f"max {d_dist.max():.0f} m; practices median {np.median(s_dist):.0f} m, max {s_dist.max():.0f} m")

t0 = time.time()
# Dijkstra from each practice on the REVERSED graph gives the time from every node to that practice
to_prac = dijkstra(W.T.tocsr(), directed=True, indices=s_node) / 60        # practices x nodes, minutes
T = to_prac[:, d_node].T + d_hop[:, None] + s_hop[None, :]                 # OAs x practices, minutes
print(f"travel-time matrix: {T.shape[0]:,} OAs x {T.shape[1]} practices in {time.time() - t0:.1f} s; "
      f"unreachable pairs: {np.isinf(T).sum()}")
pd.DataFrame(T, index=pwc.OA21CD, columns=prac.PRACTICE_CODE).to_parquet(OUT / "oa_practice_minutes.parquet")
pwc.assign(node=[nodes[i] for i in d_node], hop_min=d_hop).to_file(DATA / "demand_points.gpkg")
prac.assign(node=[nodes[i] for i in s_node], hop_min=s_hop).to_file(DATA / "practices.gpkg")

# ---------------------------------------------------------------- Step 4: nearest practice
section("nearest practice")
inC = pwc.in_cornwall.to_numpy()
popC = pwc.loc[inC, "pop"].to_numpy()


def wstats(x, w=popC):
    order = np.argsort(x)
    cum = np.cumsum(w[order]) / w.sum()
    return np.average(x, weights=w), x[order][np.searchsorted(cum, 0.5)]


pwc["nearest_min"] = T.min(axis=1)
pwc["nearest_code"] = prac.PRACTICE_CODE.to_numpy()[T.argmin(axis=1)]
near = pwc.loc[inC, "nearest_min"].to_numpy()
m, med = wstats(near)
print(f"Cornwall residents' drive time to the nearest GP practice: mean {m:.1f} min, median {med:.1f} min, "
      f"max {near.max():.1f} min")
for t in (10, 15, 20):
    print(f"  more than {t} min: {popC[near > t].sum():,} residents ({popC[near > t].sum() / popC.sum():.1%})")
cross = ~pwc.loc[inC, "nearest_code"].isin(prac.loc[prac.in_cornwall, "PRACTICE_CODE"])
print(f"Cornwall OAs whose nearest practice is in Devon: {cross.sum()} ({popC[cross.to_numpy()].sum():,} residents)")

# ---------------------------------------------------------------- Step 5: 2SFCA and E2SFCA with access
section("2SFCA and E2SFCA (access package)")
cost = (pd.DataFrame(T, index=pwc.OA21CD, columns=prac.PRACTICE_CODE).stack().rename("minutes")
        .rename_axis(["origin", "dest"]).reset_index())
cost = cost[cost.minutes <= 60]
print(f"cost table: {len(cost):,} OA-practice pairs within 60 min")
STEP = {10: 1.0, 20: 0.68, 30: 0.22}        # Luo and Qi (2009) zone weights
SIGMA = 15                                  # Gaussian decay: weight 0.61 at 15 min, 0.14 at 30 min
A = Access(demand_df=pwc[["OA21CD", "pop"]], demand_index="OA21CD", demand_value="pop",
           supply_df=prac[["PRACTICE_CODE", "TOTAL_GP_FTE"]], supply_index="PRACTICE_CODE",
           supply_value="TOTAL_GP_FTE", cost_df=cost, cost_origin="origin", cost_dest="dest", cost_name="minutes")
A.two_stage_fca(name="2sfca", max_cost=30)
A.enhanced_two_stage_fca(name="e2sfca", max_cost=30, weight_fn=weights.step_fn(STEP))
A.enhanced_two_stage_fca(name="g2sfca", max_cost=60, weight_fn=weights.gaussian(SIGMA))
acc = A.access_df.reindex(pwc.OA21CD).fillna(0) * 1000                    # GP FTE per 1,000 residents
for name in ("2sfca", "e2sfca", "g2sfca"):
    pwc[name] = acc[f"{name}_TOTAL_GP_FTE"].to_numpy()


def fca(T, demand, supply, weight):
    """Vectorised 2SFCA: weight is a function of the travel-time matrix (OAs x practices)."""
    Wt = weight(T)
    ratio = supply / (Wt.T @ demand)          # step 1: supply-to-demand ratio for each practice
    return Wt @ ratio * 1000                  # step 2: sum the ratios each OA can reach


pop_all, fte = pwc["pop"].to_numpy(float), prac.TOTAL_GP_FTE.to_numpy(float)
binary = lambda t, d=30: (t <= d).astype(float)
step = lambda t: np.select([t <= 10, t <= 20, t <= 30], [1.0, 0.68, 0.22], 0.0)
gauss = lambda t: np.where(t <= 60, np.exp(-t ** 2 / (2 * SIGMA ** 2)), 0.0)
for name, fn in [("2sfca", binary), ("e2sfca", step), ("g2sfca", gauss)]:
    mine = fca(T, pop_all, fte, fn)
    print(f"{name}: numpy vs access max abs difference {np.abs(mine - pwc[name]).max():.2e}")

total = prac.loc[prac.in_cornwall, "TOTAL_GP_FTE"].sum() / popC.sum() * 1000
print(f"Cornwall overall: {total:.3f} GP FTE per 1,000 residents (Cornwall practices / Cornwall residents)")
for name in ("2sfca", "e2sfca", "g2sfca"):
    x = pwc.loc[inC, name].to_numpy()
    m, med = wstats(x)
    q10, q90 = np.percentile(x, [10, 90])
    print(f"{name:7s} mean {m:.3f} | median {med:.3f} | OA 10th-90th pct {q10:.3f}-{q90:.3f} | "
          f"min {x.min():.3f} max {x.max():.3f} | CV {x.std() / x.mean():.2f}")
r = pd.Series(near).corr(pd.Series(pwc.loc[inC, "e2sfca"].to_numpy()), method="spearman")
print(f"Spearman rank correlation, nearest-practice time vs E2SFCA (Cornwall OAs): {r:.2f}")

# ---------------------------------------------------------------- Step 6: catchment threshold
section("catchment threshold sensitivity (2SFCA)")
pop_all, fte = pwc["pop"].to_numpy(float), prac.TOTAL_GP_FTE.to_numpy(float)
base30 = pwc.loc[inC, "2sfca"].to_numpy()
thresholds = [10, 15, 20, 30, 45, 60]
by_threshold = {}
for d in thresholds:
    x = fca(T, pop_all, fte, lambda t, d=d: binary(t, d))[inC]
    by_threshold[d] = x
    zero = popC[x == 0].sum()
    rho = pd.Series(x).corr(pd.Series(base30), method="spearman")
    print(f"{d:2d} min: OA 10th-90th pct {np.percentile(x, 10):.3f}-{np.percentile(x, 90):.3f} | "
          f"CV {x.std() / x.mean():.2f} | residents with no practice in reach {zero:,} | Spearman vs 30 min {rho:.2f}")

# ---------------------------------------------------------------- Step 7: edge effects
section("edge effects at the Cornwall-Devon border")
devon = lads[~lads.LAD26CD.isin([LAD, "E06000053"])].union_all()
border = cornwall.boundary.intersection(devon.buffer(200))
pwc["border_km"] = pwc.distance(border) / 1000
pC = prac.in_cornwall.to_numpy()
variants = {
    "Cornwall only": fca(T[inC][:, pC], popC.astype(float), fte[pC], step),
    "+ Devon practices only": fca(T[inC], popC.astype(float), fte, step),
    "+ Devon practices and residents": pwc.loc[inC, "e2sfca"].to_numpy(),
}
bands = pd.cut(pwc.loc[inC, "border_km"], [0, 5, 10, 20, 40, 200], right=False,
               labels=["0-5", "5-10", "10-20", "20-40", "40+"])
edge = pd.DataFrame({k: v for k, v in variants.items()}).assign(pop=popC, band=bands.values)
edge_tab = edge.groupby("band", observed=True).apply(
    lambda g: pd.Series({k: np.average(g[k], weights=g["pop"]) for k in variants} | {"residents": g["pop"].sum()}))
print("population-weighted mean E2SFCA (GP FTE per 1,000) by distance from the Devon border (km):")
print(edge_tab.round(3).to_string())
full = variants["+ Devon practices and residents"]
for k in ("Cornwall only", "+ Devon practices only"):
    diff = (variants[k] - full) / full * 100
    near_b = bands.isin(["0-5", "5-10"]).to_numpy()
    print(f"{k}: vs full buffer, OAs within 10 km of the border: median {np.median(diff[near_b]):+.1f}%, "
          f"range {diff[near_b].min():+.1f}% to {diff[near_b].max():+.1f}%; beyond 40 km max |diff| "
          f"{np.abs(diff[(edge.band == '40+').to_numpy()]).max():.1f}%")

# ---------------------------------------------------------------- Step 8: urban and rural
section("urban and rural areas")
ruc = arcgis(RUC_URL, DATA / "cornwall_lsoa_ruc.geojson", where="LSOA21NM LIKE 'Cornwall %'",
             outFields="LSOA21CD,RUC21NM,Urban_rura")
pwc = pwc.merge(oas[["OA21CD", "LSOA21CD"]], on="OA21CD", how="left").merge(
    ruc[["LSOA21CD", "Urban_rura", "RUC21NM"]], on="LSOA21CD", how="left")
c = pwc[pwc.in_cornwall]
for ur, g in c.groupby("Urban_rura"):
    print(f"{ur}: {g['pop'].sum():,} residents ({g['pop'].sum() / c['pop'].sum():.0%}) | nearest practice "
          f"{np.average(g.nearest_min, weights=g['pop']):.1f} min | 2SFCA {np.average(g['2sfca'], weights=g['pop']):.3f} | "
          f"E2SFCA {np.average(g.e2sfca, weights=g['pop']):.3f} GP FTE per 1,000")
low = c.e2sfca < np.percentile(c.e2sfca, 10)
print(f"lowest-access decile of OAs (E2SFCA < {np.percentile(c.e2sfca, 10):.3f}): {c.loc[low, 'pop'].sum():,} residents, "
      f"{(c.loc[low, 'Urban_rura'] == 'Urban').mean():.0%} urban, mean nearest-practice time "
      f"{np.average(c.loc[low, 'nearest_min'], weights=c.loc[low, 'pop']):.1f} min")
pwc.drop(columns="geometry").to_csv(DATA / "oa_access.csv", index=False)

# ---------------------------------------------------------------- Figures
section("figures")
cmap_oa = oas.merge(pwc.drop(columns="geometry"), on="OA21CD")
cmap_oa["density"] = cmap_oa["pop"] / (cmap_oa.area / 1e6)
devon_lads = lads[~lads.LAD26CD.isin([LAD, "E06000053"])]
x0, y0, x1, y1 = cornwall.bounds
XLIM, YLIM = (x0 - 3000, x1 + 22000), (y0 - 3000, y1 + 3000)
pc, pd_ = prac[prac.in_cornwall], prac[~prac.in_cornwall]


def frame(ax, title, xlim=XLIM, ylim=YLIM):
    devon_lads.boundary.plot(ax=ax, color="#9aa0a0", linewidth=0.5)
    gpd.GeoSeries([cornwall]).boundary.plot(ax=ax, color="#333333", linewidth=0.6)
    ax.set_xlim(xlim); ax.set_ylim(ylim); ax.set_axis_off(); ax.set_title(title, fontsize=10)
    sx, sy = xlim[0] + 3000, ylim[0] + 3000
    ax.plot([sx, sx + 10000], [sy, sy], color="black", lw=2)
    ax.text(sx + 5000, sy + 2500, "10 km", ha="center", fontsize=8)


def fte_size(f):
    return 6 + f * 4


# Figure 1: study area, demand and supply
fig, ax = plt.subplots(figsize=(9, 6.2))
cmap_oa.plot(ax=ax, column="density", cmap="YlGnBu", norm=matplotlib.colors.LogNorm(10, 10000), linewidth=0,
             legend=True, legend_kwds={"label": "Residents per km² (Census 2021, log scale)", "shrink": 0.6})
pwc[~pwc.in_cornwall].plot(ax=ax, color="#7a7a7a", markersize=1.5, label="Devon OA centroids (buffer demand)")
gpd.GeoSeries([study], crs=CRS).boundary.plot(ax=ax, color="#555555", linestyle="--", linewidth=0.8)
pc.plot(ax=ax, color="#d95f02", edgecolor="black", linewidth=0.4, markersize=fte_size(pc.TOTAL_GP_FTE), zorder=4)
pd_.plot(ax=ax, facecolor="none", edgecolor="#d95f02", linewidth=1, markersize=fte_size(pd_.TOTAL_GP_FTE), zorder=4)
frame(ax, f"Cornwall: {pwc.in_cornwall.sum():,} output areas and {pC.sum()} GP practices, "
          f"plus {(~pC).sum()} practices in the 25 km Devon buffer")
handles = [Line2D([], [], marker="o", ls="", markerfacecolor="#d95f02", markeredgecolor="black",
                  markersize=np.sqrt(fte_size(f)), label=f"{f} GP FTE") for f in (5, 10, 20)]
handles += [Line2D([], [], marker="o", ls="", markerfacecolor="none", markeredgecolor="#d95f02", markersize=6,
                   label="Practice in Devon"),
            Line2D([], [], marker="o", ls="", color="#7a7a7a", markersize=2, label="Devon OA centroid"),
            Line2D([], [], color="#555555", ls="--", label="25 km buffer")]
ax.legend(handles=handles, loc="upper left", fontsize=8, frameon=True)
ax.text(251000, 43500, "Plymouth", ha="center", fontsize=8, color="#333333")
fig.savefig(OUT / "01-study-area.png"); plt.close(fig)

# Figure 2: nearest practice drive time
fig, axes = plt.subplots(1, 2, figsize=(11, 4.6), gridspec_kw={"width_ratios": [1.7, 1]})
bins = [0, 2, 4, 6, 8, 10, 15]
cmap_oa.plot(ax=axes[0], column="nearest_min", cmap="YlOrRd", norm=BoundaryNorm(bins, 256), linewidth=0,
             legend=True, legend_kwds={"label": "Drive time to nearest practice (min)", "shrink": 0.7,
                                       "ticks": bins})
pc.plot(ax=axes[0], color="black", markersize=4, zorder=4)
pd_.plot(ax=axes[0], facecolor="none", edgecolor="black", markersize=6, linewidth=0.6, zorder=4)
frame(axes[0], "Free-flow drive time to the nearest GP practice")
order = np.argsort(near)
axes[1].plot(near[order], np.cumsum(popC[order]) / popC.sum() * 100, color="#1543b0", lw=2)
for t in (5, 10):
    share = popC[near <= t].sum() / popC.sum() * 100
    axes[1].plot([t], [share], "o", color="#d95f02")
    axes[1].annotate(f"{share:.1f}% within {t} min", (t, share), xytext=(t + 0.5, share - 8), fontsize=8)
axes[1].set_xlabel("Drive time to nearest practice (min)"); axes[1].set_ylabel("Cornwall residents (cumulative %)")
axes[1].set_ylim(0, 101); axes[1].grid(alpha=0.3); axes[1].set_title("Share of residents within a drive time", fontsize=10)
fig.savefig(OUT / "02-nearest-practice-time.png"); plt.close(fig)

# Figure 3: 2SFCA and E2SFCA, diverging around the Cornwall-wide ratio
fig, axes = plt.subplots(1, 2, figsize=(12, 4.9))
norm = TwoSlopeNorm(vmin=-1.3, vcenter=0, vmax=1.3)
for ax, (col, title) in zip(axes, [("2sfca", "2SFCA, 30-minute catchment"),
                                    ("e2sfca", "E2SFCA, 10/20/30-minute zones (1, 0.68, 0.22)")]):
    vals = np.log2(cmap_oa[col].clip(lower=0.05) / total)
    cmap_oa.assign(v=vals).plot(ax=ax, column="v", cmap="PuOr", norm=norm, linewidth=0)
    pc.plot(ax=ax, color="black", markersize=3, zorder=4)
    pd_.plot(ax=ax, facecolor="none", edgecolor="black", markersize=5, linewidth=0.5, zorder=4)
    frame(ax, title)
sm = plt.cm.ScalarMappable(norm=norm, cmap="PuOr")
ticks = np.array([0.3, 0.45, total, 1.0, 1.4])
cb = fig.colorbar(sm, ax=axes, shrink=0.75, ticks=np.log2(ticks / total), pad=0.01)
cb.ax.set_yticklabels([f"{t:.2f}" + (" (Cornwall)" if t == total else "") for t in ticks])
cb.set_label("GP FTE per 1,000 residents (log scale)")
fig.savefig(OUT / "03-2sfca-e2sfca.png"); plt.close(fig)

# Figure 4: decay functions and catchment threshold
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
t = np.linspace(0, 60, 601)
axes[0].plot(t, binary(t), color="#1543b0", lw=2, label="2SFCA: binary, 30 min")
axes[0].plot(t, step(t), color="#d95f02", lw=2, label="E2SFCA: step 1 / 0.68 / 0.22")
axes[0].plot(t, gauss(t), color="#1b9e77", lw=2, ls="--", label=f"Gaussian, σ = {SIGMA} min, cut at 60")
axes[0].set_xlabel("Drive time (min)"); axes[0].set_ylabel("Weight"); axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3); axes[0].set_title("Distance-decay functions", fontsize=10)
bp = axes[1].boxplot([by_threshold[d] for d in thresholds], positions=range(len(thresholds)), whis=(5, 95),
                     showfliers=False, patch_artist=True, widths=0.55)
for b in bp["boxes"]:
    b.set_facecolor("#c6dbef"); b.set_edgecolor("#1543b0")
zero10 = popC[by_threshold[10] == 0].sum() / popC.sum()
axes[1].set_xlabel(f"Catchment threshold (at 10 min, {zero10:.1%} of residents have no practice in reach: access 0, not shown)",
                   fontsize=8)
axes[1].axvspan(3.5, 5.5, color="#eeeeee", zorder=0)
axes[1].text(4.5, 1.25, "catchment wider than\nthe 25 km buffer", ha="center", fontsize=8, color="#555555")
axes[1].set_xticks(range(len(thresholds)), [f"{d} min" for d in thresholds])
axes[1].set_ylabel("2SFCA, GP FTE per 1,000 residents"); axes[1].grid(alpha=0.3, axis="y")
axes[1].set_title("2SFCA by catchment threshold (Cornwall OAs, 5th–95th pct)", fontsize=10)
fig.savefig(OUT / "04-decay-and-threshold.png"); plt.close(fig)

# Figure 5: edge effects
fig, axes = plt.subplots(1, 3, figsize=(13, 4.3), gridspec_kw={"width_ratios": [1, 1, 0.9]})
enorm = TwoSlopeNorm(vmin=-2, vcenter=0, vmax=2)
cw = cmap_oa[cmap_oa.in_cornwall].copy()
idx = pd.Series(np.arange(inC.sum()), index=pwc.loc[inC, "OA21CD"])
for ax, k in zip(axes[:2], ("Cornwall only", "+ Devon practices only")):
    lr = np.log2(variants[k] / full)[idx[cw.OA21CD].to_numpy()]
    cw.assign(d=lr).plot(ax=ax, column="d", cmap="PuOr", norm=enorm, linewidth=0)
    pd_.plot(ax=ax, facecolor="none", edgecolor="black", markersize=5, linewidth=0.5, zorder=4)
    frame(ax, f"{k} ÷ full buffer", xlim=(195000, XLIM[1]), ylim=(38000, YLIM[1]))
ratio_ticks = [0.25, 0.5, 0.8, 1, 1.25, 2, 4]
cb = fig.colorbar(plt.cm.ScalarMappable(norm=enorm, cmap="PuOr"), ax=axes[:2], shrink=0.75, pad=0.01, extend="both",
                  ticks=np.log2(ratio_ticks))
cb.ax.set_yticklabels([f"×{t:g}" for t in ratio_ticks])
cb.set_label("E2SFCA ratio to the full-buffer result (log scale)")
for k, colr, ls in zip(variants, ("#d95f02", "#7570b3", "#1b9e77"), ("-", "--", "-")):
    axes[2].plot(edge_tab.index.astype(str), edge_tab[k], marker="o", color=colr, ls=ls, label=k)
axes[2].set_xlabel("Distance from the Devon border (km)"); axes[2].set_ylabel("Mean E2SFCA, GP FTE per 1,000")
axes[2].legend(fontsize=8); axes[2].grid(alpha=0.3); axes[2].set_title("Population-weighted mean by distance band", fontsize=10)
fig.savefig(OUT / "05-edge-effects.png"); plt.close(fig)

# Figure 6: nearest-practice time vs E2SFCA
fig, ax = plt.subplots(figsize=(7.5, 4.8))
handles = []
for ur, colr in (("Urban", "#0072B2"), ("Rural", "#E69F00")):
    g = c[c.Urban_rura == ur]
    ax.scatter(g.nearest_min, g.e2sfca, s=g["pop"] / 40, color=colr, alpha=0.55, edgecolor="none")
    handles.append(Line2D([], [], marker="o", ls="", color=colr, markersize=6,
                          label=f"{ur} LSOA ({g['pop'].sum() / c['pop'].sum():.0%} of residents)"))
ax.axhline(total, color="#555555", ls=":", lw=1)
ax.text(14.8, total + 0.02, f"Cornwall overall {total:.2f}", ha="right", fontsize=8, color="#555555")
ax.set_xlabel("Drive time to nearest practice (min)"); ax.set_ylabel("E2SFCA, GP FTE per 1,000 residents")
ax.set_title(f"Proximity and access disagree (Spearman ρ = {r:.2f}); point size = OA population", fontsize=10)
ax.legend(handles=handles, fontsize=8); ax.grid(alpha=0.3)
fig.savefig(OUT / "06-time-vs-access.png"); plt.close(fig)
print("saved", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))