"""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"))))