"""Hot spot analysis (Getis-Ord Gi*) and spatial autocorrelation (Moran's I, LISA) in Python
with PySAL, on 2025 STATS19 road collisions in Birmingham.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/hot-spot-analysis-getis-ord-python/
Data: DfT STATS19 road safety data (OGL v3), ONS boundaries (OGL v3), OpenStreetMap via Geofabrik (ODbL).
"""
import warnings
from importlib.metadata import version
from pathlib import Path
import sys

import numpy as np
import pandas as pd
import geopandas as gpd
import requests
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
from shapely.geometry import Polygon
import libpysal
from libpysal.weights import Queen, DistanceBand
import esda

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

DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
SEED = 42
PERMS = 9999
np.random.seed(SEED)   # esda.Moran has no seed argument; this makes the global tests repeatable
LAD = "E08000025"   # Birmingham
CRS = "EPSG:27700"  # British National Grid, metres

STATS19 = "https://data.dft.gov.uk/road-accidents-safety-data/dft-road-casualty-statistics-collision-2025.csv"
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
LAD_URL = f"{ONS}/Local_Authority_Districts_May_2026_Boundaries_UK_BGC/FeatureServer/0/query"
LSOA_URL = f"{ONS}/Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query"
OSM_PBF = "https://download.geofabrik.de/europe/united-kingdom/england/west-midlands-latest.osm.pbf"


def download(url, path, **params):
    if not path.exists():
        r = requests.get(url, params=params or None, timeout=600, headers={"User-Agent": "TechnicalGIS-tutorial"})
        r.raise_for_status()
        path.write_bytes(r.content)
    return path


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


section("versions")
print(f"python {sys.version.split()[0]} | geopandas {version('geopandas')} | libpysal {version('libpysal')} | "
      f"esda {version('esda')} | numpy {version('numpy')} | pandas {version('pandas')}")

# ---------------------------------------------------------------- Step 1: load the data
section("STATS19 collisions")
csv = download(STATS19, DATA / "dft-road-casualty-statistics-collision-2025.csv")
df = pd.read_csv(csv, low_memory=False)
print(f"{len(df):,} collisions in Great Britain in 2025 "
      f"({df.date.str[3:5].nunique()} months, {df.date.min()} .. {df.date.max()})")
bham = df[df.local_authority_ons_district == LAD].dropna(subset=["location_easting_osgr"])
pts = gpd.GeoDataFrame(bham, geometry=gpd.points_from_xy(bham.location_easting_osgr, bham.location_northing_osgr),
                       crs=CRS)
sev = pts.collision_severity.map({1: "fatal", 2: "serious", 3: "slight"}).value_counts()
print(f"Birmingham: {len(pts):,} collisions | fatal {sev.get('fatal', 0)}, serious {sev.get('serious', 0)}, "
      f"slight {sev.get('slight', 0)}")

boundary = gpd.read_file(download(LAD_URL, DATA / "birmingham_lad.geojson", where=f"LAD26CD='{LAD}'",
                                  outFields="LAD26CD,LAD26NM", outSR=27700, f="geojson"))
boundary = boundary.set_crs(CRS, allow_override=True)
city = boundary.geometry.iloc[0]
print(f"boundary: {boundary.LAD26NM.iloc[0]}, {city.area / 1e6:.1f} km2")


# ---------------------------------------------------------------- Step 2: hexagon grid and counts
def hex_grid(poly, width):
    """Flat-topped hexagons `width` metres across the flats, kept if their centre is inside `poly`."""
    side = width / np.sqrt(3)
    dx, dy = 1.5 * side, width
    x0, y0, x1, y1 = poly.bounds
    angles = np.radians(np.arange(0, 360, 60))
    cells = []
    for i, cx in enumerate(np.arange(x0 - side, x1 + side, dx)):
        for cy in np.arange(y0 - width + (dy / 2 if i % 2 else 0), y1 + width, dy):
            # round to the millimetre so neighbouring hexagons share exactly the same vertices
            cells.append(Polygon(zip(np.round(cx + side * np.cos(angles), 3), np.round(cy + side * np.sin(angles), 3))))
    g = gpd.GeoDataFrame(geometry=cells, crs=CRS)
    return g[g.centroid.within(poly)].reset_index(drop=True)


def count_points(units, points):
    joined = gpd.sjoin(points[["geometry"]], units[["geometry"]], predicate="within")
    return joined.groupby("index_right").size().reindex(units.index, fill_value=0).to_numpy()


section("hexagon grid")
hexes = hex_grid(city, 500)
hexes["n"] = count_points(hexes, pts)
print(f"{len(hexes):,} hexagons 500 m across ({hexes.area.iloc[0] / 1e6:.3f} km2 each) cover "
      f"{hexes.area.sum() / 1e6:.1f} km2")
print(f"{hexes.n.sum():,} of {len(pts):,} collisions fall in the grid; "
      f"{(hexes.n == 0).mean():.1%} of hexagons have none; max {hexes.n.max()} in one hexagon")
print(f"mean {hexes.n.mean():.2f}, variance {hexes.n.var():.2f} (variance/mean {hexes.n.var() / hexes.n.mean():.1f})")

# OpenStreetMap roads for context (and later as an exposure measure)
pbf = download(OSM_PBF, DATA / "west-midlands-latest.osm.pbf")
ROAD_TYPES = ("motorway", "trunk", "primary", "secondary", "tertiary", "unclassified", "residential",
              "living_street", "motorway_link", "trunk_link", "primary_link", "secondary_link", "tertiary_link")
wgs = boundary.to_crs(4326).total_bounds
roads = gpd.read_file(pbf, layer="lines", columns=["highway"], bbox=tuple(wgs), engine="pyogrio",
                      where="highway IN ({})".format(",".join(f"'{t}'" for t in ROAD_TYPES))).to_crs(CRS)
roads = roads[roads.intersects(city)]
major = roads[roads.highway.isin(["motorway", "trunk", "primary"])]


centre = gpd.GeoSeries.from_xy([-1.9026], [52.4797], crs=4326).to_crs(CRS).iloc[0]   # Victoria Square


def basemap(ax, title=None, label_centre=False):
    boundary.boundary.plot(ax=ax, color="#333333", lw=0.7)
    major.plot(ax=ax, color="#555555", lw=0.35, alpha=0.6)
    ax.set_axis_off()
    if title:
        ax.set_title(title, fontsize=10)
    if label_centre:
        ax.plot(centre.x, centre.y, marker="o", mfc="none", mec="#111111", ms=9, mew=1.4)
        ax.annotate("city centre", (centre.x, centre.y), xytext=(-60, 30), textcoords="offset points", fontsize=8,
                    arrowprops={"arrowstyle": "-", "color": "#111111", "lw": 0.8})


def legend_below(ax, handles, ncol=3):
    ax.legend(handles=handles, loc="upper center", bbox_to_anchor=(0.5, 0.0), ncol=ncol, frameon=False, fontsize=8)


fig, axes = plt.subplots(1, 2, figsize=(12, 5.9), constrained_layout=True)
basemap(axes[0], f"{len(pts):,} injury collisions, 2025")
for code, colour, size, label in [(3, "#3288bd", 2, "slight"), (2, "#fdae61", 5, "serious"), (1, "#d53e4f", 14, "fatal")]:
    s = pts[pts.collision_severity == code]
    s.plot(ax=axes[0], color=colour, markersize=size, label=f"{label} ({len(s):,})", alpha=0.8, linewidth=0)
legend_below(axes[0], axes[0].get_legend_handles_labels()[0])
for h in axes[0].get_legend().legend_handles:
    h.set_sizes([30])
hexes.plot(ax=axes[1], column="n", cmap="YlOrRd", legend=True, edgecolor="white", linewidth=0.1,
           legend_kwds={"label": "collisions per 500 m hexagon", "shrink": 0.7})
basemap(axes[1], "Counts per hexagon", label_centre=True)
fig.savefig(OUT / "01-collisions-hexagons.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 3: spatial weights and global Moran's I
section("spatial weights")
w_queen = Queen.from_dataframe(hexes, use_index=False)
print(f"queen contiguity: {w_queen.n} units, mean {w_queen.mean_neighbors:.2f} neighbours "
      f"(min {w_queen.min_neighbors}, max {w_queen.max_neighbors}), islands {len(w_queen.islands)}")
cent = hexes.copy(); cent["geometry"] = hexes.centroid
w_d1 = DistanceBand.from_dataframe(cent, threshold=1001, binary=True, silence_warnings=True, use_index=False)
w_d2 = DistanceBand.from_dataframe(cent, threshold=2001, binary=True, silence_warnings=True, use_index=False)
for name, w in [("distance band 1 km", w_d1), ("distance band 2 km", w_d2)]:
    print(f"{name}: mean {w.mean_neighbors:.1f} neighbours (min {w.min_neighbors}, max {w.max_neighbors})")

section("global Moran's I (counts, queen, row-standardised)")
y = hexes.n.to_numpy().astype(float)
w_queen.transform = "r"
mi = esda.Moran(y, w_queen, permutations=PERMS)
print(f"I = {mi.I:.3f}  E[I] = {mi.EI:.4f}  z(normal) = {mi.z_norm:.1f}  "
      f"z(permutation) = {mi.z_sim:.1f}  pseudo p = {mi.p_sim:.4f}")
print(f"largest of {PERMS} permuted I values: {mi.sim.max():.3f}")
for name, w in [("distance band 1 km", w_d1), ("distance band 2 km", w_d2)]:
    w.transform = "r"
    m = esda.Moran(y, w, permutations=PERMS)
    print(f"{name}: I = {m.I:.3f}, z(permutation) = {m.z_sim:.1f}, pseudo p = {m.p_sim:.4f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), constrained_layout=True)
z = (y - y.mean()) / y.std()
lag = libpysal.weights.lag_spatial(w_queen, z)
axes[0].scatter(z, lag, s=6, color="#3288bd", alpha=0.5, linewidth=0)
b = np.polyfit(z, lag, 1)
xx = np.linspace(z.min(), z.max(), 2)
axes[0].plot(xx, b[0] * xx + b[1], color="#d53e4f", lw=1.5, label=f"slope = Moran's I = {mi.I:.3f}")
axes[0].axhline(0, color="#999999", lw=0.6); axes[0].axvline(0, color="#999999", lw=0.6)
for (tx, ty, t) in [(0.97, 0.97, "HH"), (0.03, 0.97, "LH"), (0.03, 0.03, "LL"), (0.97, 0.03, "HL")]:
    axes[0].text(tx, ty, t, transform=axes[0].transAxes, ha="right" if tx > 0.5 else "left",
                 va="top" if ty > 0.5 else "bottom", fontsize=10, color="#666666")
axes[0].set_xlabel("collisions per hexagon (z-score)")
axes[0].set_ylabel("spatial lag: mean z of neighbours")
axes[0].set_title("Moran scatterplot", fontsize=10)
axes[0].legend(frameon=False, loc="upper center")
axes[1].hist(mi.sim, bins=60, color="#99d594", edgecolor="none")
axes[1].axvline(mi.I, color="#d53e4f", lw=2)
axes[1].annotate(f"observed I = {mi.I:.3f}", (mi.I, 0), xytext=(-8, 150), textcoords="offset points",
                 ha="right", color="#d53e4f")
axes[1].set_xlabel(f"Moran's I under {PERMS:,} random permutations")
axes[1].set_ylabel("count")
axes[1].set_title("Permutation reference distribution", fontsize=10)
fig.savefig(OUT / "02-moran-scatterplot.png")
plt.close(fig)


# ---------------------------------------------------------------- Step 4: Getis-Ord Gi* with FDR
def gi_star(values, w):
    """Gi* z-scores and two-sided pseudo p-values: double the one-tailed ("directed") permutation p."""
    g = esda.G_Local(values, w, transform="B", permutations=PERMS, star=True, seed=SEED, keep_simulations=False,
                     alternative="directed")
    return g.Zs, np.minimum(1.0, 2 * g.p_sim)


def classify(zs, p, alpha_levels=(0.01, 0.05, 0.10)):
    """ArcGIS-style bins: +/-3 = 99 %, +/-2 = 95 %, +/-1 = 90 % confidence, 0 = not significant."""
    out = np.zeros(len(zs), dtype=int)
    for level, a in zip((3, 2, 1), alpha_levels):
        out[(out == 0) & (p <= a)] = level
    return out * np.sign(zs).astype(int)


def fdr_bins(zs, p, alpha=0.05):
    thr = esda.fdr(p, alpha)
    return np.where(p <= thr, np.sign(zs).astype(int) * 3, 0), thr


section("Getis-Ord Gi* (queen + self, binary weights)")
w_queen.transform = "b"
zs, p = gi_star(y, w_queen)
hexes["gi_z"], hexes["gi_p"] = zs, p
raw = classify(zs, p)
fdr, thr = fdr_bins(zs, p)
print(f"two-sided pseudo p-values: min {p.min():.4f} (floor 2/(1+{PERMS}) = {2 / (1 + PERMS):.4f})")
for lvl, name in [(3, "99%"), (2, "95%"), (1, "90%")]:
    print(f"hot {name}: {(raw == lvl).sum():4d}   cold {name}: {(raw == -lvl).sum():4d}")
print(f"p <= 0.05 without correction: {(p <= 0.05).sum()} hexagons "
      f"(expected false positives under the null: {0.05 * len(p):.0f})")
print(f"FDR (Benjamini-Hochberg, alpha 0.05) p threshold = {thr:.4f}: "
      f"{(fdr > 0).sum()} hot, {(fdr < 0).sum()} cold hexagons")
print(f"Bonferroni threshold would be {0.05 / len(p):.6f}, below the smallest attainable p, "
      f"so it would find nothing with {PERMS} permutations")
g2 = esda.G_Local(y, w_queen, transform="B", permutations=PERMS, star=True, seed=SEED, keep_simulations=False,
                  alternative="two-sided")
print(f"esda alternative='two-sided' instead: smallest p {g2.p_sim.min():.4f}, FDR threshold "
      f"{esda.fdr(g2.p_sim, 0.05):.4f}, {(g2.p_sim <= esda.fdr(g2.p_sim, 0.05)).sum()} significant hexagons")
hot = hexes[fdr > 0]
print(f"FDR hot spots hold {hot.n.sum():,} collisions ({hot.n.sum() / hexes.n.sum():.1%}) "
      f"in {hot.area.sum() / 1e6:.1f} km2 ({len(hot) / len(hexes):.1%} of the city)")
hexes["gi_fdr"] = fdr

GI_COLOURS = {-3: "#3288bd", -2: "#66c2a5", -1: "#abdda4", 0: "#eeeeee", 1: "#fee08b", 2: "#fdae61", 3: "#d53e4f"}
GI_LABELS = {3: "hot 99%", 2: "hot 95%", 1: "hot 90%", 0: "not significant", -1: "cold 90%", -2: "cold 95%",
             -3: "cold 99%"}


def plot_classes(ax, gdf, classes, colours, labels, title, label_centre=False, extra=()):
    gdf.plot(ax=ax, color=[colours[c] for c in classes], edgecolor="white", linewidth=0.1)
    basemap(ax, title, label_centre)
    present = [k for k in labels if (np.asarray(classes) == k).any()]
    handles = [Patch(color=colours[k], label=f"{labels[k]} ({(np.asarray(classes) == k).sum()})") for k in present]
    legend_below(ax, handles + list(extra), ncol=min(4, len(handles) + len(extra)))


fig, axes = plt.subplots(1, 2, figsize=(12, 5.6), constrained_layout=True)
plot_classes(axes[0], hexes, raw, GI_COLOURS, GI_LABELS, "Gi* at 90/95/99%, no correction")
plot_classes(axes[1], hexes, fdr, GI_COLOURS, {3: "hot spot (FDR 5%)", 0: "not significant",
                                                 -3: "cold spot (FDR 5%)"}, "Gi* after FDR correction", True)
fig.savefig(OUT / "03-gi-star-fdr.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 5: LISA (local Moran's I)
section("local Moran's I (LISA, queen, row-standardised)")
w_queen.transform = "r"
lisa = esda.Moran_Local(y, w_queen, permutations=PERMS, seed=SEED, keep_simulations=False, alternative="directed")
lp = np.minimum(1.0, 2 * lisa.p_sim)
lthr = esda.fdr(lp, 0.05)
quad = np.where(lp <= lthr, lisa.q, 0)
QUAD = {1: "High-High", 2: "Low-High", 3: "Low-Low", 4: "High-Low", 0: "not significant"}
print(f"FDR threshold {lthr:.4f}; " + ", ".join(f"{QUAD[k]} {(quad == k).sum()}" for k in (1, 2, 3, 4)))
print(f"uncorrected p <= 0.05: " + ", ".join(f"{QUAD[k]} {((lp <= 0.05) & (lisa.q == k)).sum()}" for k in (1, 2, 3, 4)))
both = (quad == 1) & (fdr > 0)
print(f"High-High LISA and Gi* hot spots agree on {both.sum()} hexagons; "
      f"Gi*-only {((fdr > 0) & (quad != 1)).sum()}, LISA-only {((quad == 1) & (fdr <= 0)).sum()}")
LISA_COLOURS = {1: "#d53e4f", 2: "#abdda4", 3: "#3288bd", 4: "#fdae61", 0: "#eeeeee"}
fig, axes = plt.subplots(1, 2, figsize=(12, 5.9), constrained_layout=True)
plot_classes(axes[0], hexes, np.where(lp <= 0.05, lisa.q, 0), LISA_COLOURS, QUAD, "LISA, p <= 0.05 uncorrected")
plot_classes(axes[1], hexes, quad, LISA_COLOURS, QUAD, "LISA after FDR correction", True)
fig.savefig(OUT / "04-lisa-clusters.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 6: the weights matter
section("Gi* hot spots under different weights (FDR 5%)")
results = {}
for name, w in [("queen (6 neighbours)", w_queen), ("distance band 1 km (18)", w_d1), ("distance band 2 km (~60)", w_d2)]:
    w.transform = "b"
    zz, pp = gi_star(y, w)
    cls, t = fdr_bins(zz, pp)
    results[name] = cls
    print(f"{name:26s} hot {(cls > 0).sum():4d}  cold {(cls < 0).sum():4d}  "
          f"hot area {(cls > 0).sum() * hexes.area.iloc[0] / 1e6:5.1f} km2  FDR p threshold {t:.4f}")
names = list(results)
for a, b2 in [(0, 1), (0, 2), (1, 2)]:
    A, B = results[names[a]] > 0, results[names[b2]] > 0
    print(f"hot-spot overlap (Jaccard) {names[a].split(' (')[0]} vs {names[b2].split(' (')[0]}: {(A & B).sum() / (A | B).sum():.2f}")
fig, axes = plt.subplots(1, 3, figsize=(15, 5.2), constrained_layout=True)
for ax, (name, cls) in zip(axes, results.items()):
    plot_classes(ax, hexes, cls, GI_COLOURS, {3: "hot", 0: "n.s.", -3: "cold"}, name)
fig.savefig(OUT / "05-weights-sensitivity.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 7: MAUP
section("modifiable areal unit problem")
lsoa = gpd.read_file(download(LSOA_URL, DATA / "birmingham_lsoa21.geojson", where="LSOA21NM LIKE 'Birmingham%'",
                              outFields="LSOA21CD,LSOA21NM", outSR=27700, f="geojson"))
lsoa = lsoa.set_crs(CRS, allow_override=True)
lsoa["n"] = count_points(lsoa, pts)
units = {"hexagons 250 m": hex_grid(city, 250), "hexagons 500 m": hexes[["geometry", "n"]],
         "hexagons 1 km": hex_grid(city, 1000), f"LSOAs ({len(lsoa)})": lsoa}
maup = {}
for name, u in units.items():
    if "n" not in u:
        u["n"] = count_points(u, pts)
    w = Queen.from_dataframe(u, use_index=False, silence_warnings=True)
    yy = u.n.to_numpy().astype(float)
    w.transform = "r"
    m = esda.Moran(yy, w, permutations=PERMS)
    w.transform = "b"
    zz, pp = gi_star(yy, w)
    cls, _ = fdr_bins(zz, pp)
    maup[name] = (u, cls)
    hot_area = u.area[cls > 0].sum() / 1e6
    print(f"{name:16s} units {len(u):5d}  mean count {yy.mean():5.2f}  Moran's I {m.I:.3f} (z {m.z_sim:5.1f})  "
          f"FDR hot {(cls > 0).sum():4d} units = {hot_area:5.1f} km2 holding {yy[cls > 0].sum() / yy.sum():.1%} of collisions")
fig, axes = plt.subplots(1, 4, figsize=(17, 4.8), constrained_layout=True)
for ax, (name, (u, cls)) in zip(axes, maup.items()):
    plot_classes(ax, u, cls, GI_COLOURS, {3: "hot", 0: "n.s.", -3: "cold"}, name)
fig.savefig(OUT / "06-maup.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 8: exposure (collisions per km of road)
section("normalising by exposure: road length")
pieces = gpd.overlay(roads[["geometry"]], hexes[["geometry"]].reset_index(names="hex"), how="intersection",
                     keep_geom_type=True)
hexes["road_km"] = pieces.assign(km=pieces.length / 1000).groupby("hex").km.sum().reindex(hexes.index, fill_value=0)
print(f"{len(roads):,} OSM road ways; {hexes.road_km.sum():,.0f} km of road inside the grid "
      f"(dual carriageways are mapped as two ways and count twice)")
print(f"road length per hexagon: median {hexes.road_km.median():.1f} km, "
      f"{(hexes.road_km < 1).sum()} hexagons under 1 km")
r = np.corrcoef(hexes.n, hexes.road_km)[0, 1]
print(f"correlation between collisions and road length per hexagon: r = {r:.2f}")
sub = hexes[hexes.road_km >= 1].reset_index(drop=True)
ws = Queen.from_dataframe(sub, use_index=False, silence_warnings=True)
keep = [i for i in range(ws.n) if i not in set(ws.islands)]
sub = sub.loc[keep].reset_index(drop=True)
ws = Queen.from_dataframe(sub, use_index=False, silence_warnings=True)
e, bkm = sub.n.to_numpy().astype(float), sub.road_km.to_numpy()
print(f"{len(sub)} hexagons with >= 1 km of road; city-wide rate {e.sum() / bkm.sum():.3f} collisions per road-km")
mr = esda.Moran_Rate(e, bkm, ws, adjusted=True, transformation="r", permutations=PERMS)
print(f"Moran's I of the empirical-Bayes rate: {mr.I:.3f} (z {mr.z_sim:.1f}, pseudo p {mr.p_sim:.4f})")
rate = esda.Moran_Local_Rate(e, bkm, ws, adjusted=True, transformation="r", permutations=PERMS, seed=SEED,
                             keep_simulations=False, alternative="directed")
rp = np.minimum(1.0, 2 * rate.p_sim)
rq = np.where(rp <= esda.fdr(rp, 0.05), rate.q, 0)
print("EB-rate LISA (FDR 5%): " + ", ".join(f"{QUAD[k]} {(rq == k).sum()}" for k in (1, 2, 3, 4)))
ws.transform = "b"
cz, cp = gi_star(e, ws)
cnt_hot = fdr_bins(cz, cp)[0] > 0
rate_hot = rq == 1
print(f"count Gi* hot spots on this subset: {cnt_hot.sum()}; EB-rate High-High: {rate_hot.sum()}; "
      f"in both: {(cnt_hot & rate_hot).sum()}")
sub["rate"] = e / bkm
top = sub.sort_values("n", ascending=False).head(10)
print(f"10 busiest hexagons: median {top.road_km.median():.1f} road-km vs {sub.road_km.median():.1f} median of all {len(sub)}")
fig, axes = plt.subplots(1, 2, figsize=(12, 5.9), constrained_layout=True)
excluded = Patch(facecolor="white", edgecolor="#999999", hatch="////", label="< 1 km of road (excluded)")
for ax in axes:
    hexes[hexes.road_km < 1].plot(ax=ax, color="white", edgecolor="#bbbbbb", linewidth=0.1, hatch="////")
plot_classes(axes[0], sub, np.where(cnt_hot, 3, 0), GI_COLOURS, {3: "hot spot", 0: "n.s."},
             "Gi* on raw counts (FDR 5%)", True, [excluded])
plot_classes(axes[1], sub, rq, LISA_COLOURS, QUAD, "LISA on collisions per road-km (EB, FDR 5%)", False, [excluded])
fig.savefig(OUT / "07-exposure-rate.png")
plt.close(fig)

hexes.to_file(OUT / "birmingham_hexagons_gi.gpkg", driver="GPKG")
print("\nwrote figures 01-07 and output/birmingham_hexagons_gi.gpkg")