"""Spatial joins in GeoPandas: sjoin, predicates, sjoin_nearest, clip and overlay on Greater Manchester.

13,000+ NaPTAN bus stops and OpenStreetMap schools joined to ONS wards and 2021 LSOAs, with point-in-polygon
counts, predicate behaviour on shared boundaries, nearest-neighbour joins with max_distance, a 400 m
walking-catchment overlay, and a timing comparison of a naive loop against the spatial index.

Part 2 of the TechnicalGIS "Geospatial Python essentials" series:
https://www.technicalgis.com/tutorials/spatial-join-geopandas/
"""
import json
import os
import sys
import time
from importlib.metadata import version
from pathlib import Path

for var in ("PROJ_LIB", "PROJ_DATA", "GDAL_DATA"):   # stray installs (e.g. PostGIS) break the wheels' own copies
    os.environ.pop(var, None)

import numpy as np
import pandas as pd
import geopandas as gpd
import requests
import shapely
import matplotlib
matplotlib.use("Agg")   # render to files only
import matplotlib.pyplot as plt

plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 140)
DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)

ARCGIS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
WARDS = ARCGIS + "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
LADS = ARCGIS + "Local_Authority_Districts_May_2024_Boundaries_UK_BFC/FeatureServer/0/query"
LSOAS = ARCGIS + "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BFC_V10/FeatureServer/0/query"
NOMIS = "https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv"   # Census 2021 TS001, usual residents
NAPTAN = "https://naptan.api.dft.gov.uk/v1/access-nodes"                  # ATCO area 180 = Greater Manchester
OVERPASS = ["https://overpass-api.de/api/interpreter", "https://overpass.private.coffee/api/interpreter"]
GM_BOROUGHS = ["Bolton", "Bury", "Manchester", "Oldham", "Rochdale", "Salford", "Stockport", "Tameside",
               "Trafford", "Wigan"]
IN_GM = "(" + ",".join(f"'{b}'" for b in GM_BOROUGHS) + ")"


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


def ons_query(url, where, fields, out_file):
    """Page through an ONS ArcGIS feature service, asking for British National Grid coordinates."""
    if not out_file.exists():
        parts, offset = [], 0
        while True:
            r = requests.get(url, params={"where": where, "outFields": fields, "outSR": 27700, "f": "geojson",
                                          "resultOffset": offset, "resultRecordCount": 500}, timeout=180)
            r.raise_for_status()
            page = gpd.GeoDataFrame.from_features(r.json()["features"], crs=27700)
            parts.append(page)
            offset += len(page)
            if len(page) < 500:
                break
        pd.concat(parts, ignore_index=True).to_file(out_file, driver="GPKG")
    return gpd.read_file(out_file)


def best_of(fn, repeats=3):
    best = np.inf
    for _ in range(repeats):
        t0 = time.perf_counter()
        res = fn()
        best = min(best, time.perf_counter() - t0)
    return best, res


section("versions")
print(f"python {sys.version.split()[0]} | geopandas {version('geopandas')} | shapely {shapely.__version__} "
      f"(GEOS {shapely.geos_version_string}) | pyogrio {version('pyogrio')} | pandas {pd.__version__}")

# ---------------------------------------------------------------- Step 1: data
section("Step 1: polygons (ONS) and points (NaPTAN, OpenStreetMap), all in EPSG:27700")
boroughs = ons_query(LADS, f"LAD24NM IN {IN_GM}", "LAD24CD,LAD24NM", DATA / "gm_lad24.gpkg")
wards = ons_query(WARDS, f"LAD24NM IN {IN_GM}", "WD24CD,WD24NM,LAD24NM", DATA / "gm_wards_dec2024.gpkg")
lsoa = ons_query(LSOAS, " OR ".join(f"LSOA21NM LIKE '{b} %'" for b in GM_BOROUGHS), "LSOA21CD,LSOA21NM",
                 DATA / "gm_lsoa21.gpkg")
wards["area_km2"] = wards.area / 1e6
census_file = DATA / "census2021_ts001_lsoa.csv"
if not census_file.exists():
    r = requests.get(NOMIS, params={"date": "latest", "geography": "TYPE151", "measures": 20100,
                                    "c2021_restype_3": 0, "select": "geography_code,obs_value"}, timeout=300)
    r.raise_for_status()
    census_file.write_bytes(r.content)
census = pd.read_csv(census_file).rename(columns={"GEOGRAPHY_CODE": "LSOA21CD", "OBS_VALUE": "residents"})
lsoa = lsoa.merge(census, on="LSOA21CD", how="left")

stops_csv = DATA / "naptan_gm.csv"
if not stops_csv.exists():
    r = requests.get(NAPTAN, params={"dataFormat": "csv", "atcoAreaCodes": "180"}, timeout=300)
    r.raise_for_status()
    stops_csv.write_bytes(r.content)
naptan = pd.read_csv(stops_csv, low_memory=False)
naptan = naptan[(naptan.Status == "active") & (naptan.StopType == "BCT") & naptan.Easting.notna()]
stops = gpd.GeoDataFrame(naptan[["ATCOCode", "CommonName", "LocalityName"]].reset_index(drop=True),
                         geometry=gpd.points_from_xy(naptan.Easting, naptan.Northing), crs=27700)

schools_json = DATA / "osm_schools_gm.json"
if not schools_json.exists():
    q = ('[out:json][timeout:180];area["name"="Greater Manchester"]["boundary"]->.a;'
         'nwr["amenity"="school"](area.a);out center tags;')
    for attempt in range(6):   # public Overpass servers are busy: retry, alternating two instances
        r = requests.post(OVERPASS[attempt % 2], data={"data": q}, timeout=300,
                          headers={"User-Agent": "TechnicalGIS-tutorial/1.0 (https://www.technicalgis.com)"})
        if r.ok and r.json().get("elements"):
            break
        time.sleep(20)
    r.raise_for_status()
    assert r.json()["elements"], "Overpass returned no schools"
    schools_json.write_bytes(r.content)
els = json.loads(schools_json.read_text(encoding="utf-8"))["elements"]
rows = [{"osm_id": f"{e['type']}/{e['id']}", "name": e.get("tags", {}).get("name", ""),
         "lon": e.get("lon", e.get("center", {}).get("lon")), "lat": e.get("lat", e.get("center", {}).get("lat"))}
        for e in els]
schools = pd.DataFrame(rows).dropna(subset=["lon", "lat"])
schools = gpd.GeoDataFrame(schools, geometry=gpd.points_from_xy(schools.lon, schools.lat), crs=4326).to_crs(27700)

for name, g in [("boroughs", boroughs), ("wards", wards), ("LSOAs 2021", lsoa), ("bus stops", stops),
                ("schools (OSM)", schools)]:
    print(f"{name:14s} {len(g):6,d} rows  {g.geom_type.iloc[0]:12s} CRS {g.crs.to_string()}")
assert wards.crs == stops.crs == schools.crs == lsoa.crs, "reproject first: sjoin needs one CRS"
print(f"Census 2021 residents in the LSOAs: {lsoa.residents.sum():,.0f}")

# ---------------------------------------------------------------- Step 2: point in polygon
section("Step 2: point in polygon - which ward is each bus stop in?")
joined = gpd.sjoin(stops, wards[["WD24CD", "WD24NM", "LAD24NM", "geometry"]], how="left", predicate="within")
print(joined[["ATCOCode", "CommonName", "WD24NM", "LAD24NM", "index_right"]].head(4).to_string())
print(f"rows in: {len(stops):,}  rows out: {len(joined):,}  stops matched to no ward: {joined.WD24CD.isna().sum()}")
outside = joined[joined.WD24CD.isna()]
if len(outside):
    print("unmatched stops are in:", outside.LocalityName.value_counts().head(5).to_dict())

counts = joined.groupby("WD24CD").size().rename("n_stops")
wards = wards.merge(counts, on="WD24CD", how="left").fillna({"n_stops": 0})
wards["stops_per_km2"] = wards.n_stops / wards.area_km2
print(f"wards: {len(wards)}, stops per ward: median {wards.n_stops.median():.0f}, "
      f"min {wards.n_stops.min():.0f}, max {wards.n_stops.max():.0f}")
cols = ["WD24NM", "LAD24NM", "n_stops", "area_km2", "stops_per_km2"]
print("densest wards:\n" + wards.nlargest(3, "stops_per_km2")[cols].round(1).to_string(index=False))
print("sparsest wards:\n" + wards.nsmallest(3, "stops_per_km2")[cols].round(1).to_string(index=False))
by_borough = joined.groupby("LAD24NM").size()
print("stops per borough:", by_borough.sort_values(ascending=False).to_dict())

# ---------------------------------------------------------------- Step 3: predicates
section("Step 3: predicates - points vs polygons, and polygons vs polygons")
gw = wards[["WD24CD", "geometry"]]
for pred in ["within", "intersects"]:
    print(f"stops {pred:10s} wards: {len(gpd.sjoin(stops, gw, predicate=pred)):,} pairs")
print(f"wards contains   stops: {len(gpd.sjoin(gw, stops, predicate='contains')):,} pairs (same pairs, other way round)")
print(f"valid predicates: {sorted(stops.sindex.valid_query_predicates - {None})}")

lsoa_pts = lsoa.copy()
lsoa_pts["geometry"] = lsoa.representative_point()
pred_pairs = {}
for pred in ["intersects", "overlaps", "touches", "within", "covered_by"]:
    pred_pairs[pred] = gpd.sjoin(lsoa[["LSOA21CD", "geometry"]], gw, predicate=pred)
pred_pairs["rep. point within"] = gpd.sjoin(lsoa_pts[["LSOA21CD", "geometry"]], gw, predicate="within")
for pred, j in pred_pairs.items():
    per = j.groupby("LSOA21CD").size()
    print(f"LSOAs {pred:18s} wards: {len(j):5,d} pairs, {per.size:5,d} LSOAs matched, "
          f"{(per > 1).sum():5,d} matched to 2+ wards")
inter = gpd.overlay(lsoa[["LSOA21CD", "geometry"]], gw, how="intersection", keep_geom_type=True)
inter["a"] = inter.area
max_share = (inter.a / inter.groupby("LSOA21CD").a.transform("sum")).groupby(inter.LSOA21CD).max()
print(f"overlay check: {(max_share < 0.99).sum():,} of {lsoa.LSOA21CD.nunique():,} LSOAs have under 99% of their "
      f"area in one ward (they genuinely straddle a ward boundary); "
      f"{(max_share >= 0.99).sum():,} sit in one ward apart from slivers")

# ---------------------------------------------------------------- Step 4: how= and counts per LSOA
section("Step 4: how='left' / 'inner' / 'right', and stops per 1,000 residents by LSOA")
gl = lsoa[["LSOA21CD", "LSOA21NM", "residents", "geometry"]]
for how in ["inner", "left", "right"]:
    j = gpd.sjoin(stops, gl, how=how, predicate="within")
    print(f"how={how:6s}: {len(j):6,d} rows, geometry from the {'right' if how == 'right' else 'left'} frame "
          f"({j.geom_type.iloc[0]})")
per_lsoa = gpd.sjoin(stops, gl, predicate="within").groupby("LSOA21CD").size().rename("n_stops")
lsoa = lsoa.merge(per_lsoa, on="LSOA21CD", how="left").fillna({"n_stops": 0})
lsoa["stops_per_1000"] = lsoa.n_stops / lsoa.residents * 1000
print(f"LSOAs with no bus stop at all: {(lsoa.n_stops == 0).sum():,} of {len(lsoa):,} "
      f"(home to {lsoa.loc[lsoa.n_stops == 0, 'residents'].sum():,.0f} residents)")
print(f"stops per 1,000 residents: median {lsoa.stops_per_1000.median():.1f}, "
      f"Greater Manchester overall {lsoa.n_stops.sum() / lsoa.residents.sum() * 1000:.1f}")

# ---------------------------------------------------------------- Step 5: sjoin_nearest
section("Step 5: sjoin_nearest - the nearest bus stop to every school")
schools = gpd.sjoin(schools, boroughs[["LAD24NM", "geometry"]], predicate="within").drop(columns="index_right")
print(f"OSM schools inside the ten boroughs: {len(schools):,}")
s_cols = ["ATCOCode", "CommonName", "geometry"]
t, near = best_of(lambda: gpd.sjoin_nearest(schools, stops[s_cols], how="left", distance_col="dist_m"), 1)
print(f"sjoin_nearest, no limit: {len(near):,} rows for {len(schools):,} schools in {t * 1000:.0f} ms "
      f"(ties add rows: {near.index.duplicated().sum()})")
near = near[~near.index.duplicated()]
d = near.dist_m
print(f"distance to the nearest stop: median {d.median():.0f} m, 90th pct {d.quantile(0.9):.0f} m, max {d.max():.0f} m")
far = near.nlargest(3, "dist_m")[["name", "LAD24NM", "CommonName", "dist_m"]]
print("furthest:\n" + far.round(0).to_string(index=False))
lim = gpd.sjoin_nearest(schools, stops[s_cols], how="left", max_distance=400, distance_col="dist_m")
lim = lim[~lim.index.duplicated()]
print(f"max_distance=400: {lim.ATCOCode.isna().sum()} schools have no stop within 400 m "
      f"({lim.ATCOCode.isna().mean() * 100:.1f}%), returned with NaN instead of a far-away match")
t_lim, _ = best_of(lambda: gpd.sjoin_nearest(schools, stops[s_cols], max_distance=400), 1)
print(f"time with max_distance: {t_lim * 1000:.0f} ms")
k400 = gpd.sjoin(schools, gpd.GeoDataFrame(geometry=stops.buffer(400), crs=27700), predicate="intersects")
print(f"stops within 400 m of each school (buffer + sjoin): median {k400.groupby(level=0).size().median():.0f}")

# ---------------------------------------------------------------- Step 6: clip and overlay
section("Step 6: clip and overlay - areas more than 400 m from a bus stop")
manc = boroughs[boroughs.LAD24NM == "Manchester"]
t_clip, clipped = best_of(lambda: gpd.clip(stops, manc), 1)
print(f"gpd.clip(stops, Manchester): {len(clipped):,} stops in {t_clip * 1000:.0f} ms; "
      f"sjoin gives {len(gpd.sjoin(stops, manc[['geometry']], predicate='intersects')):,}")
catch = gpd.GeoDataFrame(geometry=[stops.buffer(400).union_all()], crs=27700)
t_ov, beyond = best_of(lambda: gpd.overlay(lsoa[["LSOA21CD", "LSOA21NM", "residents", "geometry"]], catch,
                                           how="difference"), 1)
beyond["beyond_km2"] = beyond.area / 1e6
lsoa = lsoa.merge(beyond[["LSOA21CD", "beyond_km2"]], on="LSOA21CD", how="left").fillna({"beyond_km2": 0})
lsoa["beyond_pct"] = lsoa.beyond_km2 / (lsoa.area / 1e6) * 100
gm_area = boroughs.area.sum() / 1e6
print(f"overlay(how='difference') in {t_ov:.1f} s: {lsoa.beyond_km2.sum():,.1f} km2 of {gm_area:,.1f} km2 "
      f"({lsoa.beyond_km2.sum() / gm_area * 100:.1f}%) is more than 400 m in a straight line from a bus stop")
print(f"LSOAs with more than half their area beyond 400 m: {(lsoa.beyond_pct > 50).sum()} "
      f"(home to {lsoa.loc[lsoa.beyond_pct > 50, 'residents'].sum():,.0f} residents)")

# ---------------------------------------------------------------- Step 7: timing
section("Step 7: timing - naive loop vs vectorised vs spatial index")
ward_codes = wards.WD24CD.to_numpy()
pts = list(stops.geometry)
# fresh copies: earlier sjoins "prepared" the ward polygons in place, which would flatter the loop
ward_geoms = list(shapely.from_wkb(shapely.to_wkb(wards.geometry.values)))


def naive_loop(points, polys):
    out, tests = [], 0
    for p in points:                     # every point against every ward until a hit, in pure Python
        code = None
        for g, c in zip(polys, ward_codes):
            tests += 1
            if g.contains(p):
                code = c
                break
        out.append(code)
    return out, tests


t_naive, (naive_res, n_tests) = best_of(lambda: naive_loop(pts, ward_geoms), 1)
prepared = list(shapely.from_wkb(shapely.to_wkb(wards.geometry.values)))
shapely.prepare(prepared)
t_prep, (prep_res, _) = best_of(lambda: naive_loop(pts, prepared), 1)


def vectorised_no_index():
    res = np.full(len(stops), None, dtype=object)
    arr = stops.geometry.values
    for g, c in zip(prepared, ward_codes):   # one vectorised test per ward, but every point every time
        res[shapely.contains(g, arr)] = c
    return res


t_vec, vec_res = best_of(vectorised_no_index)


def with_sjoin():
    return gpd.sjoin(stops, wards[["WD24CD", "geometry"]], how="left", predicate="within")


def tree_on_wards():
    tree = shapely.STRtree(wards.geometry.values)
    return tree.query(stops.geometry.values, predicate="within")


def tree_on_stops():          # what sjoin(predicate="within") does internally: index the points, ask each ward
    tree = shapely.STRtree(stops.geometry.values)
    return tree.query(wards.geometry.values, predicate="contains")


t_sj, sj_res = best_of(with_sjoin, 5)
t_tw, _ = best_of(tree_on_wards, 5)
t_ts, _ = best_of(tree_on_stops, 5)
sj_codes = sj_res.WD24CD.to_numpy()
same = all((pd.Series(r).fillna("-").to_numpy() == pd.Series(sj_codes).fillna("-").to_numpy()).all()
           for r in (naive_res, prep_res, vec_res))
print(f"{len(stops):,} stops x {len(wards)} wards; all methods give identical ward codes: {same}")
print(f"the naive loop ran {n_tests:,} contains tests ({n_tests / len(stops):.0f} per stop on average)")
timings = {"naive Python loop": t_naive, "naive loop, prepared polygons": t_prep,
           "shapely.contains per ward (prepared)": t_vec, "STRtree on wards, query stops": t_tw,
           "STRtree on stops, query wards": t_ts, "gpd.sjoin(predicate='within')": t_sj}
for k, v in timings.items():
    print(f"  {k:38s} {v:8.3f} s   {t_naive / v:6.1f}x the loop's speed")
t_build, _ = best_of(lambda: shapely.STRtree(stops.geometry.values))
print(f"building an STRtree for {len(stops):,} stops: {t_build * 1000:.1f} ms")
wards[["WD24CD", "WD24NM", "LAD24NM", "n_stops", "area_km2", "stops_per_km2", "geometry"]].to_file(
    OUT / "gm_wards_bus_stops.gpkg", layer="wards", driver="GPKG")

# ---------------------------------------------------------------- figures
section("figures")
fig, ax = plt.subplots(figsize=(8, 6.4))
bins = np.quantile(wards.stops_per_km2, np.linspace(0, 1, 7))   # sextile class breaks
norm = matplotlib.colors.BoundaryNorm(bins, 256)
wards.plot("stops_per_km2", ax=ax, cmap="Spectral_r", norm=norm, edgecolor="white", linewidth=0.25, legend=True,
           legend_kwds={"label": "bus stops per km² (sextile classes)", "shrink": 0.6,
                        "ticks": bins, "format": "%.1f"})
boroughs.boundary.plot(ax=ax, color="#222222", lw=0.7)
ax.set_axis_off()
ax.set_title(f"{len(stops):,} NaPTAN bus stops joined to {len(wards)} wards with gpd.sjoin(predicate='within')")
fig.savefig(OUT / "01-bus-stops-per-ward.png")
plt.close(fig)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained", gridspec_kw={"width_ratios": [1, 1.1]})
names = list(pred_pairs)
vals = [len(pred_pairs[k]) for k in names]
bars = axes[0].barh(names, vals, color=["#d53e4f", "#fdae61", "#fee08b", "#66c2a5", "#3288bd", "#5e4fa2"])
axes[0].axvline(len(lsoa), color="#333333", ls="--", lw=0.8)
axes[0].annotate(f"{len(lsoa):,} LSOAs", (len(lsoa), -0.6), ha="center", va="bottom", fontsize=8)
for b, v in zip(bars, vals):
    axes[0].annotate(f"{v:,}", (v, b.get_y() + b.get_height() / 2), xytext=(4, 0), textcoords="offset points",
                     va="center", fontsize=8)
axes[0].invert_yaxis()
axes[0].set_xlabel("LSOA-ward pairs returned by sjoin")
axes[0].set_title("Same layers, different predicates")
multi = pred_pairs["intersects"].groupby("LSOA21CD").size()
ex = multi[multi == multi.max()].index[0]
ex_geom = lsoa[lsoa.LSOA21CD == ex]
hit = wards[wards.WD24CD.isin(pred_pairs["intersects"].loc[pred_pairs["intersects"].LSOA21CD == ex, "WD24CD"])]
ax = axes[1]
hit.plot(ax=ax, column="WD24NM", cmap="Spectral", edgecolor="white", alpha=0.75, legend=True,
         legend_kwds={"loc": "upper left", "bbox_to_anchor": (1, 1), "frameon": False, "fontsize": 7})
ex_geom.boundary.plot(ax=ax, color="black", lw=1.6)
lsoa_pts[lsoa_pts.LSOA21CD == ex].plot(ax=ax, color="black", marker="*", markersize=60)
ax.set_axis_off()
ax.set_title(f"{ex_geom.LSOA21NM.iloc[0]}: 'intersects' matches {len(hit)} wards;\n"
             "the representative point (star) picks one")
fig.savefig(OUT / "02-predicates.png")
plt.close(fig)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), layout="constrained", gridspec_kw={"width_ratios": [1.3, 1]})
ax = axes[0]
boroughs.plot(ax=ax, color="#eeeeee", edgecolor="#888888", lw=0.6)
srt = near.sort_values("dist_m")
sc = ax.scatter(srt.geometry.x, srt.geometry.y, c=srt.dist_m.clip(upper=800), cmap="Spectral_r", s=7, vmin=0, vmax=800)
fig.colorbar(sc, ax=ax, shrink=0.7, label="distance to nearest bus stop (m, capped at 800)")
ax.set_axis_off()
ax.set_title(f"sjoin_nearest: {len(near):,} OSM schools to {len(stops):,} bus stops")
axes[1].hist(d, bins=np.arange(0, d.max() + 25, 25), color="#3288bd")
axes[1].axvline(400, color="#d53e4f", ls="--")
axes[1].annotate(f"max_distance=400\n{lim.ATCOCode.isna().sum()} schools unmatched", (400, axes[1].get_ylim()[1] * 0.8),
                 xytext=(6, 0), textcoords="offset points", color="#d53e4f", fontsize=8)
axes[1].set_xlabel("distance to nearest bus stop (m)")
axes[1].set_ylabel("schools")
fig.savefig(OUT / "03-sjoin-nearest-schools.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(8, 6.4))
lsoa.plot("beyond_pct", ax=ax, cmap="Spectral_r", edgecolor="none", legend=True, vmin=0, vmax=100,
          legend_kwds={"label": "% of LSOA area more than 400 m from a bus stop", "shrink": 0.6})
boroughs.boundary.plot(ax=ax, color="#222222", lw=0.7)
ax.set_axis_off()
ax.set_title("overlay(how='difference'): LSOAs minus the 400 m bus-stop catchment")
fig.savefig(OUT / "04-overlay-400m-catchment.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(8, 3.6), layout="constrained")
ks = list(timings)
bars = ax.barh(ks, [timings[k] for k in ks], color=["#d53e4f", "#f46d43", "#fdae61", "#abdda4", "#66c2a5", "#3288bd"])
for b, k in zip(bars, ks):
    ax.annotate(f"{timings[k]:.3f} s" if timings[k] < 1 else f"{timings[k]:.1f} s", (timings[k], b.get_y() + 0.4),
                xytext=(4, 0), textcoords="offset points", va="center", fontsize=8)
ax.set_xlim(0, max(timings.values()) * 1.18)
ax.invert_yaxis()
ax.set_xlabel("seconds (best of repeated runs)")
ax.set_title(f"Point in polygon: {len(stops):,} bus stops x {len(wards)} wards")
fig.savefig(OUT / "05-timing.png")
plt.close(fig)
print("saved:", ", ".join(sorted(p.name for p in OUT.glob("*.png"))))