"""Community-scale CTA rail accessibility and distributional diagnostics.

Sources
-------
1. City of Chicago / CTA 'L' Stops (8pix-ypme), deduplicated by map_id.
2. City of Chicago community-area boundaries (igwz-8jzy), 77 areas.
3. Chicago selected socioeconomic indicators, 2008–2012 (kn9c-c2s2).

Design
------
The script samples a regular geographic grid within official community-area
polygons and assigns each cell the Haversine distance to its nearest unique rail
station. It reports citywide threshold coverage, community distributions, an
explicit comparison of high- and low-hardship quartiles, a descriptive
hardship/access rank correlation, and boundary/grid sensitivity. Results are
area-weighted descriptive exposure measures—not network travel times, current
population-weighted access, or causal effects.

Dependencies: Python standard library + Pillow.
"""
import csv
import json
import math
import statistics
from collections import defaultdict
from pathlib import Path

from PIL import Image, ImageDraw, ImageFont


ROOT = Path(__file__).resolve().parent
OUT = ROOT.parent / "public" / "gis"
OUT.mkdir(parents=True, exist_ok=True)

PAPER = "#f5f4ef"
INK = "#171716"
SOFT = "#6e6b65"
HAIR = "#cfcac0"
COBALT = "#2849a8"
CORAL = "#d56754"
GREEN = "#54736a"
PALE_BLUE = "#dbe3f5"


def font(size, bold=False, serif=False):
    candidates = []
    if serif:
        candidates += [
            "C:/Windows/Fonts/georgiab.ttf" if bold else "C:/Windows/Fonts/georgia.ttf",
            "C:/Windows/Fonts/timesbd.ttf" if bold else "C:/Windows/Fonts/times.ttf",
        ]
    candidates += ["C:/Windows/Fonts/arialbd.ttf" if bold else "C:/Windows/Fonts/arial.ttf"]
    for candidate in candidates:
        if Path(candidate).exists():
            return ImageFont.truetype(candidate, size)
    return ImageFont.load_default()


def percentile(values, q):
    ordered = sorted(values)
    pos = (len(ordered) - 1) * q
    lo, hi = math.floor(pos), math.ceil(pos)
    if lo == hi:
        return ordered[lo]
    return ordered[lo] * (hi - pos) + ordered[hi] * (pos - lo)


def haversine_km(a, b):
    lon1, lat1 = map(math.radians, a)
    lon2, lat2 = map(math.radians, b)
    term = math.sin((lat2 - lat1) / 2) ** 2
    term += math.cos(lat1) * math.cos(lat2) * math.sin((lon2 - lon1) / 2) ** 2
    return 6371 * 2 * math.asin(math.sqrt(term))


def point_in_ring(point, ring):
    x, y = point
    inside = False
    for i, a in enumerate(ring):
        b = ring[i - 1]
        if ((a[1] > y) != (b[1] > y)) and x < (b[0] - a[0]) * (y - a[1]) / (b[1] - a[1]) + a[0]:
            inside = not inside
    return inside


def point_in_polygon(point, polygon):
    return point_in_ring(point, polygon[0]) and not any(point_in_ring(point, hole) for hole in polygon[1:])


def point_in_area(point, area):
    x, y = point
    xmin, ymin, xmax, ymax = area["bbox"]
    if not (xmin <= x <= xmax and ymin <= y <= ymax):
        return False
    return any(point_in_polygon(point, polygon) for polygon in area["polygons"])


def rankdata(values):
    ordered = sorted(enumerate(values), key=lambda item: item[1])
    ranks = [0.0] * len(values)
    i = 0
    while i < len(ordered):
        j = i + 1
        while j < len(ordered) and ordered[j][1] == ordered[i][1]:
            j += 1
        average_rank = (i + 1 + j) / 2
        for k in range(i, j):
            ranks[ordered[k][0]] = average_rank
        i = j
    return ranks


def pearson(x, y):
    mx, my = statistics.mean(x), statistics.mean(y)
    numerator = sum((a - mx) * (b - my) for a, b in zip(x, y))
    denominator = math.sqrt(sum((a - mx) ** 2 for a in x) * sum((b - my) ** 2 for b in y))
    return numerator / denominator


def spearman(x, y):
    return pearson(rankdata(x), rankdata(y))


def convex_hull(points):
    def cross(o, a, b):
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

    points = sorted(set(points))
    lower, upper = [], []
    for point in points:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], point) <= 0:
            lower.pop()
        lower.append(point)
    for point in reversed(points):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], point) <= 0:
            upper.pop()
        upper.append(point)
    return lower[:-1] + upper[:-1]


raw_stops = json.loads((ROOT / "data" / "chicago_cta_stops.json").read_text(encoding="utf-8"))
station_groups = defaultdict(list)
for row in raw_stops:
    location = row.get("location") or {}
    map_id = row.get("map_id")
    if map_id and location.get("latitude") and location.get("longitude"):
        station_groups[map_id].append(
            (float(location["longitude"]), float(location["latitude"]), row.get("station_name", "Unknown"))
        )
stations = {
    map_id: (
        statistics.mean(point[0] for point in points),
        statistics.mean(point[1] for point in points),
        points[0][2],
    )
    for map_id, points in station_groups.items()
}
station_points = [(lon, lat) for lon, lat, _ in stations.values()]

geo = json.loads((ROOT / "data" / "chicago_community_areas.geojson").read_text(encoding="utf-8"))
areas = []
for feature in geo["features"]:
    properties = feature["properties"]
    coordinates = feature["geometry"]["coordinates"]
    polygons = coordinates if feature["geometry"]["type"] == "MultiPolygon" else [coordinates]
    all_points = [point for polygon in polygons for ring in polygon for point in ring]
    xs, ys = [point[0] for point in all_points], [point[1] for point in all_points]
    areas.append(
        {
            "ca": int(properties.get("area_num_1") or properties.get("area_numbe")),
            "name": properties["community"].title(),
            "polygons": polygons,
            "bbox": (min(xs), min(ys), max(xs), max(ys)),
        }
    )

xmin = min(area["bbox"][0] for area in areas)
ymin = min(area["bbox"][1] for area in areas)
xmax = max(area["bbox"][2] for area in areas)
ymax = max(area["bbox"][3] for area in areas)
outside_station_count = sum(
    not any(point_in_area((lon, lat), area) for area in areas)
    for lon, lat, _name in stations.values()
)

socio_raw = json.loads((ROOT / "data" / "chicago_socioeconomic.json").read_text(encoding="utf-8"))
socio = {}
for row in socio_raw:
    try:
        ca = int(row["ca"])
        hardship = float(row["hardship_index"])
    except (KeyError, TypeError, ValueError):
        continue
    socio[ca] = {
        "hardship": hardship,
        "poverty": float(row["percent_households_below_poverty"]),
        "unemployment": float(row["percent_aged_16_unemployed"]),
        "income": float(row["per_capita_income_"]),
    }


def sample_city(nx, ny, keep_points=False):
    by_ca = defaultdict(list)
    records = []
    for iy in range(ny):
        lat = ymax - (iy + 0.5) / ny * (ymax - ymin)
        for ix in range(nx):
            lon = xmin + (ix + 0.5) / nx * (xmax - xmin)
            area = next((candidate for candidate in areas if point_in_area((lon, lat), candidate)), None)
            if area is None:
                continue
            distance = min(haversine_km((lon, lat), station) for station in station_points)
            by_ca[area["ca"]].append(distance)
            if keep_points:
                records.append((lon, lat, area["ca"], distance))
    all_distances = [distance for values in by_ca.values() for distance in values]
    return {"nx": nx, "ny": ny, "by_ca": by_ca, "distances": all_distances, "points": records}


city_sample = sample_city(180, 260, keep_points=True)
coarse_sample = sample_city(110, 160)
fine_sample = sample_city(220, 320)
distances = city_sample["distances"]
thresholds = [0.5, 0.8, 1.0, 1.5, 2.0, 3.0]


def coverage(values, threshold):
    return 100 * sum(value <= threshold for value in values) / len(values)


city_coverage = {f"{threshold:.1f}": coverage(distances, threshold) for threshold in thresholds}

community_metrics = []
area_lookup = {area["ca"]: area for area in areas}
for ca, values in city_sample["by_ca"].items():
    if ca not in socio:
        continue
    community_metrics.append(
        {
            "ca": ca,
            "name": area_lookup[ca]["name"],
            "hardship": socio[ca]["hardship"],
            "poverty": socio[ca]["poverty"],
            "income": socio[ca]["income"],
            "mean_km": statistics.mean(values),
            "median_km": statistics.median(values),
            "p90_km": percentile(values, 0.9),
            "within_1km_pct": coverage(values, 1.0),
            "sample_cells": len(values),
        }
    )

hardship_values = [row["hardship"] for row in community_metrics]
hardship_q1 = percentile(hardship_values, 0.25)
hardship_q3 = percentile(hardship_values, 0.75)
low_ca = {row["ca"] for row in community_metrics if row["hardship"] <= hardship_q1}
high_ca = {row["ca"] for row in community_metrics if row["hardship"] >= hardship_q3}
low_distances = [value for ca in low_ca for value in city_sample["by_ca"][ca]]
high_distances = [value for ca in high_ca for value in city_sample["by_ca"][ca]]
low_coverage = {f"{threshold:.1f}": coverage(low_distances, threshold) for threshold in thresholds}
high_coverage = {f"{threshold:.1f}": coverage(high_distances, threshold) for threshold in thresholds}
low_equal_community_1km = statistics.mean(
    row["within_1km_pct"] for row in community_metrics if row["ca"] in low_ca
)
high_equal_community_1km = statistics.mean(
    row["within_1km_pct"] for row in community_metrics if row["ca"] in high_ca
)

x_hardship = [row["hardship"] for row in community_metrics]
y_access = [row["mean_km"] for row in community_metrics]
rho = spearman(x_hardship, y_access)

mean_x, mean_y = statistics.mean(x_hardship), statistics.mean(y_access)
slope = sum((x - mean_x) * (y - mean_y) for x, y in zip(x_hardship, y_access)) / sum(
    (x - mean_x) ** 2 for x in x_hardship
)
intercept = mean_y - slope * mean_x

station_hull = convex_hull(station_points)
hull_xmin, hull_xmax = min(x for x, _ in station_hull), max(x for x, _ in station_hull)
hull_ymin, hull_ymax = min(y for _, y in station_hull), max(y for _, y in station_hull)
hull_distances = []
for iy in range(100):
    lat = hull_ymax - (iy + 0.5) / 100 * (hull_ymax - hull_ymin)
    for ix in range(150):
        lon = hull_xmin + (ix + 0.5) / 150 * (hull_xmax - hull_xmin)
        if point_in_ring((lon, lat), station_hull):
            hull_distances.append(min(haversine_km((lon, lat), station) for station in station_points))
hull_1km = coverage(hull_distances, 1.0)

grid_sensitivity = {
    "110x160": coverage(coarse_sample["distances"], 1.0),
    "180x260": coverage(city_sample["distances"], 1.0),
    "220x320": coverage(fine_sample["distances"], 1.0),
}
access_gaps = sorted(community_metrics, key=lambda row: row["mean_km"], reverse=True)[:8]

# Reproducible community table.
with (OUT / "chicago_community_metrics.csv").open("w", newline="", encoding="utf-8") as handle:
    writer = csv.DictWriter(
        handle,
        fieldnames=["ca", "name", "hardship", "poverty", "income", "mean_km", "median_km", "p90_km", "within_1km_pct", "sample_cells"],
    )
    writer.writeheader()
    for row in sorted(community_metrics, key=lambda item: item["ca"]):
        writer.writerow({key: round(value, 4) if isinstance(value, float) else value for key, value in row.items()})

# Figure 1 — official-boundary distance surface.
W, H = 1500, 950
image = Image.new("RGB", (W, H), "#080b0a")
draw = ImageDraw.Draw(image)
draw.text((70, 55), "CHICAGO / RAPID-TRANSIT ACCESS", font=font(20, True), fill="#d7ff50")
draw.text((70, 92), "Distance to the nearest CTA rail station", font=font(52, True), fill="#f2f0e9")
draw.text((70, 160), "77 official community areas · 144 unique CTA rail stations", font=font(18), fill="#8e9993")
draw.line((70, 205, 1430, 205), fill="#2c3531", width=2)

left, top, right, bottom = 70, 240, 1005, 885
map_scale = min((right - left) / (xmax - xmin), (bottom - top) / (ymax - ymin))
map_width, map_height = (xmax - xmin) * map_scale, (ymax - ymin) * map_scale
map_left = left + (right - left - map_width) / 2
map_top = top + (bottom - top - map_height) / 2


def screen(point):
    return (
        map_left + (point[0] - xmin) * map_scale,
        map_top + (ymax - point[1]) * map_scale,
    )


def distance_colour(value):
    t = max(0, min(1, value / 3.0))
    low, high = (215, 255, 80), (211, 92, 61)
    return tuple(int(low[i] + t * (high[i] - low[i])) for i in range(3))


cell_w = (xmax - xmin) / city_sample["nx"] * map_scale + 1
cell_h = (ymax - ymin) / city_sample["ny"] * map_scale + 1
for lon, lat, _ca, distance in city_sample["points"]:
    sx, sy = screen((lon, lat))
    draw.rectangle((sx - cell_w / 2, sy - cell_h / 2, sx + cell_w / 2, sy + cell_h / 2), fill=distance_colour(distance))

for area in areas:
    for polygon in area["polygons"]:
        for ring_index, ring in enumerate(polygon):
            colour = "#f2f0e9" if ring_index == 0 else "#8e9993"
            draw.line([screen(point) for point in ring], fill=colour, width=1)
for lon, lat, _name in stations.values():
    sx, sy = screen((lon, lat))
    draw.ellipse((sx - 2.5, sy - 2.5, sx + 2.5, sy + 2.5), fill="#080b0a", outline="#f2f0e9")

# Distance and station legends, including stations just beyond the study boundary.
legend_x, legend_y, legend_w = 95, 820, 210
draw.rectangle((legend_x - 14, legend_y - 20, legend_x + legend_w + 16, legend_y + 48), fill="#080b0a", outline="#2c3531")
for index in range(legend_w):
    draw.line((legend_x + index, legend_y, legend_x + index, legend_y + 11), fill=distance_colour(index / legend_w * 3))
draw.text((legend_x, legend_y + 17), "0", font=font(10), fill="#f2f0e9")
draw.text((legend_x + legend_w / 2 - 11, legend_y + 17), "1.5", font=font(10), fill="#f2f0e9")
draw.text((legend_x + legend_w - 19, legend_y + 17), "3+ km", font=font(10), fill="#f2f0e9")
draw.ellipse((350, 824, 358, 832), fill="#080b0a", outline="#f2f0e9")
draw.text((370, 819), f"144 stations · {outside_station_count} outside-boundary points retained", font=font(11), fill="#8e9993")

city_1km = city_coverage["1.0"]
draw.rounded_rectangle((1065, 260, 1430, 735), radius=18, fill="#111714", outline="#2c3531", width=2)
draw.text((1105, 300), "COMMUNITY-AREA READOUT", font=font(12, True), fill="#d7ff50")
draw.text((1105, 350), f"{city_1km:.0f}%", font=font(62, True), fill="#f2f0e9")
draw.text((1105, 420), "sampled land within 1 km", font=font(16), fill="#8e9993")
draw.text((1105, 474), f"Median cell      {statistics.median(distances):.2f} km", font=font(17, True), fill="#f2f0e9")
draw.text((1105, 516), f"90th percentile  {percentile(distances, .9):.2f} km", font=font(17, True), fill="#f2f0e9")
draw.text((1105, 570), f"Low hardship Q1  {low_coverage['1.0']:.1f}%", font=font(17, True), fill="#d7ff50")
draw.text((1105, 612), f"High hardship Q4 {high_coverage['1.0']:.1f}%", font=font(17, True), fill="#d7ff50")
draw.text((1105, 664), f"Spearman ρ       {rho:.2f}", font=font(17, True), fill="#f2f0e9")
draw.text((1105, 700), "Area-weighted · descriptive", font=font(14), fill="#8e9993")
draw.text((70, 912), "METHOD  ·  official polygons  /  Haversine distance grid  /  hardship-quartile comparison", font=font(13, True), fill="#8e9993")
image.save(OUT / "chicago_transit_access.png", optimize=True)

# Figure 2 — accessibility threshold curves.
W, H = 1600, 920
image = Image.new("RGB", (W, H), PAPER)
draw = ImageDraw.Draw(image)
draw.text((78, 58), "FIGURE 02 · THRESHOLD CURVE", font=font(16, True), fill=COBALT)
draw.text((78, 98), "Access conclusions change with the distance threshold", font=font(44, serif=True), fill=INK)
draw.text((80, 166), "Area-weighted share of sampled land within each straight-line distance of a unique CTA rail station.", font=font(18), fill=SOFT)
draw.line((80, 218, 1520, 218), fill=INK, width=2)

plot_left, plot_top, plot_right, plot_bottom = 130, 290, 1460, 745
for percent in range(0, 101, 20):
    y = plot_bottom - percent / 100 * (plot_bottom - plot_top)
    draw.line((plot_left, y, plot_right, y), fill=HAIR, width=1)
    draw.text((75, y - 9), f"{percent}%", font=font(13), fill=SOFT)
for threshold in thresholds:
    x = plot_left + threshold / 3 * (plot_right - plot_left)
    draw.line((x, plot_bottom, x, plot_bottom + 8), fill=INK, width=2)
    draw.text((x - 12, plot_bottom + 18), f"{threshold:g}", font=font(13), fill=SOFT)


def draw_curve(values, colour, width=5):
    points = []
    for threshold in thresholds:
        x = plot_left + threshold / 3 * (plot_right - plot_left)
        y = plot_bottom - values[f"{threshold:.1f}"] / 100 * (plot_bottom - plot_top)
        points.append((x, y))
    draw.line(points, fill=colour, width=width, joint="curve")
    for x, y in points:
        draw.ellipse((x - 7, y - 7, x + 7, y + 7), fill=colour)


draw_curve(city_coverage, INK, 5)
draw_curve(low_coverage, COBALT, 6)
draw_curve(high_coverage, CORAL, 6)
draw.text((plot_left, 785), "Distance to nearest station (km)", font=font(15), fill=SOFT)
draw.line((980, 835, 1020, 835), fill=INK, width=5)
draw.text((1032, 824), "All community-area land", font=font(14), fill=SOFT)
draw.line((1210, 835, 1250, 835), fill=COBALT, width=6)
draw.text((1262, 824), "Lowest-hardship quartile", font=font(14), fill=SOFT)
draw.line((980, 870, 1020, 870), fill=CORAL, width=6)
draw.text((1032, 859), "Highest-hardship quartile", font=font(14), fill=SOFT)
gap = low_coverage["1.0"] - high_coverage["1.0"]
equal_gap = low_equal_community_1km - high_equal_community_1km
draw.text((80, 824), f"Pooled land area (n=20 areas per group): low minus high hardship = {gap:.1f} percentage points at 1 km.", font=font(17, True), fill=INK)
draw.text((80, 865), f"Equal-community weighting narrows the contrast to {equal_gap:.1f} points; hardship is historical (2008–2012).", font=font(13), fill=SOFT)
image.save(OUT / "chicago_coverage_curve.png", optimize=True)

# Figure 3 — community-level distributional diagnostic.
W, H = 1600, 980
image = Image.new("RGB", (W, H), PAPER)
draw = ImageDraw.Draw(image)
draw.text((78, 58), "FIGURE 03 · DISTRIBUTIONAL DIAGNOSTIC", font=font(16, True), fill=COBALT)
draw.text((78, 98), "No monotonic citywide hardship gradient", font=font(44, serif=True), fill=INK)
draw.text((80, 166), "Each mark is one community area; distance is the mean across sampled land cells.", font=font(18), fill=SOFT)
draw.line((80, 218, 1520, 218), fill=INK, width=2)

scatter_left, scatter_top, scatter_right, scatter_bottom = 100, 300, 980, 810
x_min, x_max = min(x_hardship), max(x_hardship)
y_min, y_max = 0, max(10, math.ceil(max(y_access)))
for tick in (0, 20, 40, 60, 80, 100):
    x = scatter_left + (tick - x_min) / (x_max - x_min) * (scatter_right - scatter_left)
    if scatter_left <= x <= scatter_right:
        draw.line((x, scatter_top, x, scatter_bottom), fill=HAIR, width=1)
        draw.text((x - 10, scatter_bottom + 16), str(tick), font=font(12), fill=SOFT)
for km_tick in (0, 2, 4, 6, 8, 10):
    y = scatter_bottom - (km_tick - y_min) / (y_max - y_min) * (scatter_bottom - scatter_top)
    if scatter_top <= y <= scatter_bottom:
        draw.line((scatter_left, y, scatter_right, y), fill=HAIR, width=1)
        draw.text((55, y - 8), f"{km_tick:g}", font=font(12), fill=SOFT)

line_y1 = intercept + slope * x_min
line_y2 = intercept + slope * x_max
draw.line(
    (
        scatter_left,
        scatter_bottom - (line_y1 - y_min) / (y_max - y_min) * (scatter_bottom - scatter_top),
        scatter_right,
        scatter_bottom - (line_y2 - y_min) / (y_max - y_min) * (scatter_bottom - scatter_top),
    ),
    fill=CORAL,
    width=4,
)
label_set = {row["ca"] for row in access_gaps[:5]}
for row in community_metrics:
    x = scatter_left + (row["hardship"] - x_min) / (x_max - x_min) * (scatter_right - scatter_left)
    y = scatter_bottom - (row["mean_km"] - y_min) / (y_max - y_min) * (scatter_bottom - scatter_top)
    colour = CORAL if row["ca"] in high_ca else (COBALT if row["ca"] in low_ca else GREEN)
    draw.ellipse((x - 6, y - 6, x + 6, y + 6), fill=colour)
    if row["ca"] in label_set:
        draw.text((x + 8, y - 13), row["name"], font=font(11), fill=INK)
draw.text((scatter_left, 852), "Hardship index (2008–2012)", font=font(14), fill=SOFT)
draw.text((20, scatter_top - 28), "Mean km", font=font(13), fill=SOFT)
draw.ellipse((390, 852, 402, 864), fill=COBALT); draw.text((412, 848), "LOW Q1", font=font(11), fill=SOFT)
draw.ellipse((500, 852, 512, 864), fill=GREEN); draw.text((522, 848), "MIDDLE", font=font(11), fill=SOFT)
draw.ellipse((625, 852, 637, 864), fill=CORAL); draw.text((647, 848), "HIGH Q4", font=font(11), fill=SOFT)

table_x = 1060
draw.text((table_x, 286), "LARGEST AREA-WEIGHTED ACCESS GAPS", font=font(13, True), fill=COBALT)
draw.text((table_x, 322), "Community area", font=font(13), fill=SOFT)
draw.text((1455, 322), "Mean km", font=font(13), fill=SOFT)
draw.line((table_x, 350, 1520, 350), fill=INK, width=2)
for index, row in enumerate(access_gaps[:7]):
    y = 378 + index * 55
    draw.text((table_x, y), f"{index + 1:02d}", font=font(12), fill=SOFT)
    draw.text((table_x + 38, y - 2), row["name"], font=font(15), fill=INK)
    draw.text((1455, y - 2), f"{row['mean_km']:.2f}", font=font(15, True), fill=INK)
    draw.line((table_x, y + 30, 1520, y + 30), fill=HAIR, width=1)

draw.text((table_x, 790), f"Spearman ρ = {rho:.2f}", font=font(22, serif=True), fill=INK)
draw.text((table_x, 826), "Descriptive community-level rank statistic", font=font(14), fill=SOFT)
draw.text((table_x, 858), f"Linear slope = {slope * 10:.2f} km per +10 hardship points", font=font(14), fill=SOFT)
draw.text((80, 928), "No population weights, walking-network impedance or service frequency are included; community averages may conceal within-area inequality.", font=font(13), fill=SOFT)
image.save(OUT / "chicago_equity_diagnostics.png", optimize=True)

summary = {
    "status": "exploratory spatial analysis; not causal",
    "sources": [
        "City of Chicago / CTA 'L' Stops (8pix-ypme)",
        "City of Chicago community-area boundaries (igwz-8jzy)",
        "Selected socioeconomic indicators 2008–2012 (kn9c-c2s2)",
    ],
    "raw_stop_records": len(raw_stops),
    "unique_stations": len(stations),
    "station_centroids_outside_community_area_union": outside_station_count,
    "community_areas": len(areas),
    "sampled_city_cells": len(distances),
    "city_coverage_pct": {key: round(value, 2) for key, value in city_coverage.items()},
    "median_nearest_km": round(statistics.median(distances), 3),
    "p90_nearest_km": round(percentile(distances, 0.9), 3),
    "hardship_quartile_cutpoints": [round(hardship_q1, 2), round(hardship_q3, 2)],
    "low_hardship_coverage_1km_pct": round(low_coverage["1.0"], 2),
    "high_hardship_coverage_1km_pct": round(high_coverage["1.0"], 2),
    "coverage_gap_low_minus_high_pp": round(gap, 2),
    "low_hardship_equal_community_mean_coverage_1km_pct": round(low_equal_community_1km, 2),
    "high_hardship_equal_community_mean_coverage_1km_pct": round(high_equal_community_1km, 2),
    "equal_community_coverage_gap_low_minus_high_pp": round(equal_gap, 2),
    "spearman_hardship_mean_distance": round(rho, 4),
    "linear_slope_km_per_10_hardship_points": round(slope * 10, 4),
    "boundary_sensitivity_coverage_1km_pct": {
        "official_community_area_union": round(city_1km, 2),
        "station_convex_hull": round(hull_1km, 2),
    },
    "grid_sensitivity_coverage_1km_pct": {key: round(value, 2) for key, value in grid_sensitivity.items()},
    "largest_access_gaps": [
        {
            "community_area": row["name"],
            "mean_km": round(row["mean_km"], 3),
            "within_1km_pct": round(row["within_1km_pct"], 2),
            "hardship": row["hardship"],
        }
        for row in access_gaps
    ],
    "limitations": [
        "Straight-line distance is not walking-network distance or travel time.",
        "Area-weighted exposure is not population-weighted accessibility.",
        "The socioeconomic indicators are historical (2008–2012) while the station extract has no reproducible as-of field; this is a historical ecological overlay, not a current resident-equity estimate.",
        "Station proximity omits frequency, reliability, transfers, fares, disability access and destinations.",
        "Community-area averages conceal within-area variation.",
    ],
}
(ROOT / "chicago_summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")
print(json.dumps(summary, indent=2))