"""Compute per-hex factor scores for San Francisco and emit data.json."""
import csv, json, math, os
from primitives import *

HERE = os.path.dirname(os.path.abspath(__file__))
P = lambda f: os.path.join(HERE, f)

LAT0 = 37.7600
MX = 111320.0 * math.cos(math.radians(LAT0))   # metres per degree longitude
MY = 110574.0                                   # metres per degree latitude

def xy(lat, lon):
    return (lon * MX, lat * MY)

def dist_m(lat1, lon1, lat2, lon2):
    return math.hypot((lon1 - lon2) * MX, (lat1 - lat2) * MY)

def dist_to_seg(lat, lon, alat, alon, blat, blon):
    px, py = xy(lat, lon); ax, ay = xy(alat, alon); bx, by = xy(blat, blon)
    dx, dy = bx - ax, by - ay
    L2 = dx * dx + dy * dy
    if L2 == 0:
        return math.hypot(px - ax, py - ay), 0.0
    t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / L2))
    return math.hypot(px - (ax + t * dx), py - (ay + t * dy)), t

def dist_to_line(lat, lon, pts):
    """Min distance in metres to a polyline, plus interpolation position."""
    best, bt, bi = 1e12, 0.0, 0
    for i in range(len(pts) - 1):
        a, b = pts[i], pts[i + 1]
        d, t = dist_to_seg(lat, lon, a[0], a[1], b[0], b[1])
        if d < best:
            best, bt, bi = d, t, i
    return best, bt, bi

def in_poly(lon, lat, poly):
    n = len(poly); c = False; j = n - 1
    for i in range(n):
        xi, yi = poly[i]; xj, yj = poly[j]
        if ((yi > lat) != (yj > lat)) and (lon < (xj - xi) * (lat - yi) / (yj - yi) + xi):
            c = not c
        j = i
    return c

# ---------------------------------------------------------------- load inputs
grid = json.load(open(P("grid.json")))
elev = {int(r["id"]): float(r["elev_m"]) for r in csv.DictReader(open(P("elevations.csv")))}

cells = []
for c in grid["cells"]:
    e = elev.get(c["id"])
    if e is None or e < -3.0:      # drop open-water cells the boundary over-covered
        continue
    c = dict(c); c["elev"] = max(0.0, e); c["elev_raw"] = e
    cells.append(c)
print("cells kept:", len(cells), "of", len(grid["cells"]))

def load(fn):
    return list(csv.DictReader(open(P(fn))))

parks = [(float(r["lat"]), float(r["lon"]), float(r["acres"]), r["type"]) for r in load("parks.csv")]
schools = [(float(r["lat"]), float(r["lon"]), r["low_grade"], r["high_grade"]) for r in load("schools.csv")]
libraries = [(float(r["lat"]), float(r["lon"])) for r in load("libraries.csv")]
biz = {"grocery": [], "convenience": [], "restaurant": [], "bar": []}
for r in load("businesses.csv"):
    if r["category"] in biz:
        biz[r["category"]].append((float(r["lat"]), float(r["lon"])))
print({k: len(v) for k, v in biz.items()})

playgrounds = [p for p in parks if "Playground" in p[3] or "Neighborhood Park" in p[3]]

# Snap curated supermarkets to the nearest MEASURED grocery point within 350 m
SUPER = []
for name, la, lo in SUPERMARKETS:
    best, bp = 1e12, None
    for g in biz["grocery"]:
        d = dist_m(la, lo, g[0], g[1])
        if d < best:
            best, bp = d, g
    SUPER.append((name, bp[0], bp[1]) if best <= 350 else (name, la, lo))
snapped = sum(1 for i, s in enumerate(SUPER) if (s[1], s[2]) != (SUPERMARKETS[i][1], SUPERMARKETS[i][2]))
print("supermarkets snapped to measured grocery points:", snapped, "/", len(SUPER))

# ------------------------------------------------- neighbourhood assignment
PARK_NBHDS = {"Golden Gate Park", "Presidio", "McLaren Park", "Lincoln Park"}
for c in cells:
    best, bn = 1e12, None
    for nb in NEIGHBORHOODS:
        d = dist_m(c["lat"], c["lon"], nb[1], nb[2])
        if d < best:
            best, bn = d, nb
    c["nb"] = bn[0]; c["nbdata"] = bn
    c["residential"] = 0 if bn[0] in PARK_NBHDS else 1

# Self-consistent area normalisation: neighbourhood area = its share of the grid
nb_cells = {}
for c in cells:
    nb_cells[c["nb"]] = nb_cells.get(c["nb"], 0) + 1
HEX_KM2 = (grid["spacing_m"] / 1000.0) ** 2 * math.sqrt(3) / 2
nb_density = {}
for nb in NEIGHBORHOODS:
    n = nb_cells.get(nb[0], 0)
    area = max(0.25, n * HEX_KM2)
    nb_density[nb[0]] = (nb[3] / area, nb[4] / area, nb[5] / area)

# ------------------------------------------------------------- slope (MEASURED)
by_id = {c["id"]: c for c in cells}
pos = {(round(c["lat"], 5), round(c["lon"], 5)): c for c in cells}
for c in cells:
    grads = []
    for o in cells:
        if o["id"] == c["id"]:
            continue
        d = dist_m(c["lat"], c["lon"], o["lat"], o["lon"])
        if d < 620:
            grads.append(abs(o["elev"] - c["elev"]) / d)
    c["slope_pct"] = (sorted(grads)[-2] if len(grads) >= 2 else (grads[0] if grads else 0.0)) * 100

# --------------------------------------------------------- distance to shore
shore_segs = []
B = SF_BOUNDARY
for i in range(len(B) - 1):
    if COUNTY_LINE_RANGE[0] <= i <= COUNTY_LINE_RANGE[1]:
        continue                      # inland county line, not a shoreline
    shore_segs.append((B[i], B[i + 1]))
for i in range(len(TREASURE_ISLAND) - 1):
    shore_segs.append((TREASURE_ISLAND[i], TREASURE_ISLAND[i + 1]))

for c in cells:
    best = 1e12
    for (a, b) in shore_segs:
        d, _ = dist_to_seg(c["lat"], c["lon"], a[1], a[0], b[1], b[0])
        best = min(best, d)
    c["shore_m"] = best

# ------------------------------------------------- near-road exposure surface
def road_exposure(c, lines, decay_near=150.0, decay_far=450.0):
    """Near-road decay: steep inside ~150 m, tail out to ~450 m (HEI 2010/2022).

    Returns (summed exposure, single worst road). The sum captures cumulative
    burden in places ringed by arterials; the max captures the distinct penalty
    of fronting one very big road, which the sum alone understates — a house on
    19th Ave should score badly even though only one arterial is near it.
    """
    tot = 0.0; mx = 0.0
    for entry in lines:
        w, pts = entry[1], entry[-1]
        d, _, _ = dist_to_line(c["lat"], c["lon"], [(p[0], p[1]) for p in pts])
        e = w * (0.72 * math.exp(-d / decay_near) + 0.28 * math.exp(-d / decay_far))
        tot += e
        if e > mx: mx = e
    return tot, mx

for c in cells:
    c["fw_exp"], c["fw_max"] = road_exposure(c, FREEWAYS, 170.0, 520.0)
    c["art_exp"], c["art_max"] = road_exposure(c, ARTERIALS, 110.0, 330.0)
    hin = [(a[0], min(1.0, a[2] / 8.0), a[3]) for a in ARTERIALS if a[2] > 0]
    c["hin_exp"], c["hin_max"] = road_exposure(c, hin, 90.0, 260.0)
    c["comm_exp"], _ = road_exposure(c, COMMERCIAL, 150.0, 400.0)

# ------------------------------------------------------------------- transit
def transit_downtown(c):
    best = 1e9; via = ""
    for name, mode, stops in RAIL_LINES:
        d, t, i = dist_to_line(c["lat"], c["lon"], [(s[0], s[1]) for s in stops])
        ride = stops[i][2] + t * (stops[i + 1][2] - stops[i][2])
        walk = d / 80.0                       # 80 m/min walking
        pen = 0.0 if mode == "bart" else 1.5  # Muni Metro reliability penalty
        tt = walk + ride + pen
        if walk > 14: tt += (walk - 14) * 0.7  # long walks get less attractive
        if tt < best: best, via = tt, name
    for name, q, stops in BUS_CORRIDORS:
        d, t, i = dist_to_line(c["lat"], c["lon"], [(s[0], s[1]) for s in stops])
        ride = stops[i][2] + t * (stops[i + 1][2] - stops[i][2])
        walk = d / 80.0
        tt = walk + ride + (1.0 - q) * 9.0 + 3.0
        if tt < best: best, via = tt, name
    return best, via

def transit_sfo(c):
    best = 1e9
    for (la, lo, m) in BART_TO_SFO:
        walk = dist_m(c["lat"], c["lon"], la, lo) / 80.0
        if walk <= 20:
            best = min(best, walk + m)
    dt, _ = transit_downtown(c)
    best = min(best, dt + 30 + 6)     # ride to the Market corridor, then BART
    return best

def drive_sfo(c):
    best = 1e9
    for (_n, la, lo, m) in SFO_RAMPS:
        d = dist_m(c["lat"], c["lon"], la, lo)
        best = min(best, d / 420.0 + 2.0 + m)   # ~25 km/h surface streets
    return best

for c in cells:
    c["t_dt"], c["t_via"] = transit_downtown(c)
    c["t_sfo"] = transit_sfo(c)
    c["d_sfo"] = drive_sfo(c)
    stop_q = 0.0
    for name, mode, stops in RAIL_LINES:
        d, _, _ = dist_to_line(c["lat"], c["lon"], [(s[0], s[1]) for s in stops])
        stop_q = max(stop_q, (1.0 if mode == "bart" else 0.85) * math.exp(-d / 500.0))
    for name, q, stops in BUS_CORRIDORS:
        d, _, _ = dist_to_line(c["lat"], c["lon"], [(s[0], s[1]) for s in stops])
        stop_q = max(stop_q, q * 0.62 * math.exp(-d / 350.0))
    c["stop_q"] = stop_q

# ------------------------------------------------------- amenity accessibility
def decayed_count(c, pts, scale, cap=None):
    s = 0.0
    for p in pts:
        d = dist_m(c["lat"], c["lon"], p[0], p[1])
        if d < scale * 5:
            s += math.exp(-d / scale)
    return min(s, cap) if cap else s

def nearest(c, pts):
    best = 1e12
    for p in pts:
        d = dist_m(c["lat"], c["lon"], p[0], p[1])
        if d < best: best = d
    return best

for c in cells:
    c["d_super"] = nearest(c, [(s[1], s[2]) for s in SUPER])
    c["d_park"] = 1e12; c["park_acc"] = 0.0
    for (la, lo, ac, ty) in parks:
        d = dist_m(c["lat"], c["lon"], la, lo)
        w = min(1.0, math.sqrt(max(ac, 0.1)) / 7.0)     # acreage saturates
        c["park_acc"] += w * math.exp(-d / 300.0)       # fast decay, per the brief
        if d < c["d_park"]: c["d_park"] = d
    c["d_play"] = nearest(c, [(p[0], p[1]) for p in playgrounds])
    c["d_school"] = nearest(c, [(s[0], s[1]) for s in schools])
    c["d_lib"] = nearest(c, [(l[0], l[1]) for l in libraries])
    c["d_hosp"] = nearest(c, [(h[1], h[2]) for h in HOSPITALS])
    c["corner"] = decayed_count(c, biz["convenience"] + biz["grocery"], 180.0, 8.0)
    c["food"] = decayed_count(c, biz["restaurant"], 350.0, 60.0)
    c["night"] = decayed_count(c, biz["bar"], 400.0, 14.0)
    c["walk_poi"] = (decayed_count(c, biz["restaurant"], 420.0, 70.0) * 1.0 +
                     decayed_count(c, biz["grocery"], 420.0, 12.0) * 3.0 +
                     decayed_count(c, biz["convenience"], 420.0, 8.0) * 2.0 +
                     decayed_count(c, biz["bar"], 420.0, 16.0) * 1.0)

# ------------------------------------------------------------------ hazards
for c in cells:
    liq = 0.0
    for (_n, sev, poly) in LIQUEFACTION:
        if in_poly(c["lon"], c["lat"], poly):
            liq = max(liq, sev)
        else:                                   # soften the polygon edge
            dmin = 1e12
            for i in range(len(poly)):
                a, b = poly[i], poly[(i + 1) % len(poly)]
                d, _ = dist_to_seg(c["lat"], c["lon"], a[1], a[0], b[1], b[0])
                dmin = min(dmin, d)
            if dmin < 350:
                liq = max(liq, sev * (1.0 - dmin / 350.0) * 0.6)
    # Low elevation right at the shore is a liquefaction tell even off-polygon
    if c["elev"] < 4 and c["shore_m"] < 500:
        liq = max(liq, 0.55)
    c["liq"] = liq

    e, sh = c["elev"], c["shore_m"]
    # Cal OES tsunami inundation for SF tracks low elevation near open water
    c["tsu"] = max(0.0, min(1.0, (12.0 - e) / 12.0)) * math.exp(-sh / 700.0)
    # Sea level rise: NOAA-style intermediate-high plus storm surge
    c["slr"] = max(0.0, min(1.0, (4.5 - e) / 4.5)) * math.exp(-sh / 900.0)
    # Slope over a 400-600 m baseline smooths out local steepness, so the observed
    # distribution tops out near 29% (p90 = 12%). Threshold at 10% and scale by 16
    # to spread the real range instead of flagging only the top 5% of cells.
    c["landslide"] = max(0.0, min(1.0, (c["slope_pct"] - 10.0) / 16.0))

# Pluvial flood: DERIVED topographic sink + ENCODED known hotspots.
# A true sink has NO lower neighbour, so use min(neighbours) - own elevation.
# Using the MEAN instead would falsely flag every hillside, where uphill
# neighbours drag the average up even though water drains away downhill.
for c in cells:
    around = [o["elev"] for o in cells
              if o["id"] != c["id"] and dist_m(c["lat"], c["lon"], o["lat"], o["lon"]) < 850]
    bowl = 0.0
    if around:
        sink = min(around) - c["elev"]          # >0 only in a genuine depression
        bowl = max(0.0, min(1.0, sink / 7.0))
    hot = 0.0
    for (_n, sev, la, lo, rad) in FLOOD_HOTSPOTS:
        d = dist_m(c["lat"], c["lon"], la, lo)
        hot = max(hot, sev * math.exp(-(d / rad) ** 2))
    # Flat low-lying ground near the shore also ponds in heavy rain
    lowflat = 1.0 if (c["elev"] < 5 and c["slope_pct"] < 3) else 0.0
    c["flood"] = min(1.0, 0.60 * bowl + 0.80 * hot + 0.25 * lowflat)

# -------------------------------------------------------------- microclimate
F = FOG
for c in cells:
    sun = 1.0 - math.exp(-(c["lon"] - F["west_edge_lon"]) / F["decay_deg"])
    if c["lon"] > F["ridge_lon"] and F["lee_lat_range"][0] <= c["lat"] <= F["lee_lat_range"][1]:
        prox = min(1.0, (c["lon"] - F["ridge_lon"]) / 0.020)
        sun += F["lee_bonus"] * prox
    if c["lat"] > F["gate_lat"]:
        sun -= F["gate_penalty"] * min(1.0, (c["lat"] - F["gate_lat"]) / 0.012)
    if c["lon"] < F["coast_lon"]:
        sun -= F["coast_penalty"] * min(1.0, (F["coast_lon"] - c["lon"]) / 0.014)
    c["sun"] = max(0.0, min(1.15, sun))

# -------------------------------------------------- neighbourhood-linked data
for c in cells:
    v, p, t = nb_density[c["nb"]]
    c["viol_d"], c["prop_d"], c["c311_d"] = v, p, t

# ------------------------------------------------------------ raw -> factors
def inv(x):  return -x

RAW = {
    # Air: freeways dominate, but fronting one big arterial carries its own penalty
    "air":            lambda c: -(c["fw_exp"] * 2.2 + c["art_exp"] * 0.75 + c["art_max"] * 2.6),
    "noise":          lambda c: -(c["fw_exp"] * 1.8 + c["art_exp"] * 0.9 +
                                  c["art_max"] * 2.8 + c["comm_exp"] * 0.5),
    "microclimate":   lambda c: c["sun"],
    "slope":          lambda c: -c["slope_pct"],
    "liquefaction":   lambda c: -c["liq"],
    "tsunami":        lambda c: -c["tsu"],
    "slr":            lambda c: -c["slr"],
    "flood":          lambda c: -c["flood"],
    "landslide":      lambda c: -c["landslide"],
    "transit_dt":     lambda c: -c["t_dt"],
    "transit_sfo":    lambda c: -c["t_sfo"],
    "drive_sfo":      lambda c: -c["d_sfo"],
    "walkability":    lambda c: c["walk_poi"],
    "bikeability":    lambda c: -c["slope_pct"] * 1.9 + c["walk_poi"] * 0.16 - c["hin_max"] * 22.0,
    "transit_score":  lambda c: c["stop_q"],
    "street_calm":    lambda c: -(c["art_exp"] * 0.9 + c["art_max"] * 3.4 + c["fw_exp"] * 1.2),
    "ped_safety":     lambda c: -(c["hin_exp"] * 0.7 + c["hin_max"] * 3.0),
    "supermarket":    lambda c: -c["d_super"],
    "corner_store":   lambda c: c["corner"],
    "parks":          lambda c: c["park_acc"],
    "food_nightlife": lambda c: c["food"] + c["night"] * 2.2,
    "crime_violent":  lambda c: -c["viol_d"],
    "crime_property": lambda c: -c["prop_d"],
    "street_cond":    lambda c: -c["c311_d"],
    "schools":        lambda c: -c["d_school"],
    "playgrounds":    lambda c: -c["d_play"],
    "library":        lambda c: -c["d_lib"],
    "hospital":       lambda c: -c["d_hosp"],
}

# Hazard factors carry a large mass of genuinely zero-risk cells. Rank-scaling
# them would tie ~60% of the city at an arbitrary mid value and make bedrock
# read the same as sand. Scale these linearly on their true 0-1 severity so
# "100" means "no modelled risk" and the number is directly interpretable.
UNIT_SCALED = {"liquefaction", "tsunami", "slr", "flood", "landslide"}

def percentile_scale(vals):
    """Rank-based 0-100 so weights are comparable across differently-shaped factors."""
    order = sorted(range(len(vals)), key=lambda i: vals[i])
    out = [0.0] * len(vals)
    n = len(vals)
    i = 0
    while i < n:
        j = i
        while j + 1 < n and vals[order[j + 1]] == vals[order[i]]:
            j += 1
        mid = (i + j) / 2.0
        for k in range(i, j + 1):
            out[order[k]] = 100.0 * mid / max(1, n - 1)
        i = j + 1
    return out

scores = {}
for key, fn in RAW.items():
    vals = [fn(c) for c in cells]
    if key in UNIT_SCALED:
        scaled = [100.0 * (1.0 + v) for v in vals]      # v = -severity in [-1, 0]
    else:
        scaled = percentile_scale(vals)
    scores[key] = [round(max(0.0, min(100.0, v)), 1) for v in scaled]

# ------------------------------------------------------------------- output
out = {
    "meta": {
        "spacing_m": grid["spacing_m"],
        "dx": grid["dx"], "dy": grid["dy"],
        "n": len(cells),
        "powell": POWELL_ST, "sfo": SFO,
    },
    "factors": [
        {"key": k, "label": l, "group": g, "prov": p, "w": w}
        for (k, l, g, p, w) in FACTORS
    ],
    "cells": {
        "lat":  [round(c["lat"], 5) for c in cells],
        "lon":  [round(c["lon"], 5) for c in cells],
        "row":  [c["row"] for c in cells],
        "nb":   [c["nb"] for c in cells],
        "res":  [c["residential"] for c in cells],
    },
    "scores": scores,
    # Reference geography so the hex field reads as San Francisco rather than
    # an abstract lattice: coastline, the roads that drive the penalty surfaces,
    # the rail lines that drive transit access, and orientation labels.
    "geo": {
        "boundary":  [[round(p[1], 5), round(p[0], 5)] for p in SF_BOUNDARY],
        "ti":        [[round(p[1], 5), round(p[0], 5)] for p in TREASURE_ISLAND],
        "freeways":  [{"n": n, "w": w, "p": [[p[0], p[1]] for p in pts]}
                      for (n, w, pts) in FREEWAYS],
        "arterials": [{"n": n, "w": w, "p": [[p[0], p[1]] for p in pts]}
                      for (n, w, _h, pts) in
                      [(a[0], a[1], a[2], a[3]) for a in ARTERIALS] if w >= 0.60],
        "rail":      [{"n": n, "m": m, "p": [[s[0], s[1]] for s in stops]}
                      for (n, m, stops) in RAIL_LINES],
        "commercial":[{"n": n, "i": i, "p": [[p[0], p[1]] for p in pts]}
                      for (n, i, pts) in COMMERCIAL if i >= 0.70],
        "labels": [
            {"n": "Outer Sunset",   "lat": 37.7530, "lon": -122.4990},
            {"n": "Inner Sunset",   "lat": 37.7620, "lon": -122.4690},
            {"n": "Richmond",       "lat": 37.7800, "lon": -122.4790},
            {"n": "Presidio",       "lat": 37.7980, "lon": -122.4660},
            {"n": "Marina",         "lat": 37.8035, "lon": -122.4370},
            {"n": "Pacific Hts",    "lat": 37.7925, "lon": -122.4345},
            {"n": "North Beach",    "lat": 37.8025, "lon": -122.4105},
            {"n": "Downtown",       "lat": 37.7900, "lon": -122.4030},
            {"n": "SoMa",           "lat": 37.7780, "lon": -122.4085},
            {"n": "Hayes Valley",   "lat": 37.7760, "lon": -122.4260},
            {"n": "Haight",         "lat": 37.7700, "lon": -122.4470},
            {"n": "Castro",         "lat": 37.7615, "lon": -122.4350},
            {"n": "Mission",        "lat": 37.7595, "lon": -122.4160},
            {"n": "Mission Bay",    "lat": 37.7705, "lon": -122.3935},
            {"n": "Potrero Hill",   "lat": 37.7580, "lon": -122.4005},
            {"n": "Noe Valley",     "lat": 37.7500, "lon": -122.4325},
            {"n": "Twin Peaks",     "lat": 37.7520, "lon": -122.4475},
            {"n": "Bernal Hts",     "lat": 37.7395, "lon": -122.4155},
            {"n": "Dogpatch",       "lat": 37.7570, "lon": -122.3885},
            {"n": "Glen Park",      "lat": 37.7340, "lon": -122.4340},
            {"n": "W of Twin Peaks","lat": 37.7345, "lon": -122.4605},
            {"n": "Lakeshore",      "lat": 37.7250, "lon": -122.4845},
            {"n": "Excelsior",      "lat": 37.7240, "lon": -122.4270},
            {"n": "Bayview",        "lat": 37.7305, "lon": -122.3905},
            {"n": "Visitacion",     "lat": 37.7140, "lon": -122.4070},
            {"n": "Treasure Is.",   "lat": 37.8230, "lon": -122.3700},
        ],
    },
    "raw": {
        "elev":     [round(c["elev"], 1) for c in cells],
        "slope":    [round(c["slope_pct"], 1) for c in cells],
        "t_dt":     [round(c["t_dt"], 1) for c in cells],
        "t_via":    [c["t_via"] for c in cells],
        "t_sfo":    [round(c["t_sfo"], 1) for c in cells],
        "d_sfo":    [round(c["d_sfo"], 1) for c in cells],
        "d_super":  [int(c["d_super"]) for c in cells],
        "d_park":   [int(c["d_park"]) for c in cells],
        "d_play":   [int(c["d_play"]) for c in cells],
        "d_school": [int(c["d_school"]) for c in cells],
        "d_lib":    [int(c["d_lib"]) for c in cells],
        "d_hosp":   [int(c["d_hosp"]) for c in cells],
        "corner":   [round(c["corner"], 2) for c in cells],
        "food":     [round(c["food"], 1) for c in cells],
        "night":    [round(c["night"], 1) for c in cells],
        "liq":      [round(c["liq"], 2) for c in cells],
        "tsu":      [round(c["tsu"], 2) for c in cells],
        "slr":      [round(c["slr"], 2) for c in cells],
        "flood":    [round(c["flood"], 2) for c in cells],
        "sun":      [round(c["sun"], 2) for c in cells],
        "shore":    [int(c["shore_m"]) for c in cells],
        "fw":       [round(c["fw_exp"], 3) for c in cells],
        "art":      [round(c["art_exp"], 3) for c in cells],
        "viol_d":   [round(c["viol_d"], 0) for c in cells],
        "prop_d":   [round(c["prop_d"], 0) for c in cells],
        "c311_d":   [round(c["c311_d"], 0) for c in cells],
    },
}
json.dump(out, open(P("data.json"), "w"), separators=(",", ":"))
print("wrote data.json:", os.path.getsize(P("data.json")) // 1024, "KB")

# ------------------------------------------------------------ sanity output
def at(lat, lon):
    best, bi = 1e12, 0
    for i, c in enumerate(cells):
        d = dist_m(lat, lon, c["lat"], c["lon"])
        if d < best: best, bi = d, i
    return bi

probes = [
    ("19th Ave & Noriega",   37.7537, -122.4757),
    ("Outer Sunset 45th Av", 37.7560, -122.5000),
    ("Valencia & 20th",      37.7585, -122.4210),
    ("Mission Bay",          37.7700, -122.3930),
    ("Marina Blvd",          37.8040, -122.4370),
    ("Bernal Heights top",   37.7430, -122.4160),
    ("Twin Peaks slope",     37.7530, -122.4470),
    ("Noe Valley 24th",      37.7513, -122.4330),
    ("Bayview 3rd St",       37.7300, -122.3910),
    ("Inner Richmond",       37.7828, -122.4650),
]
print("\n%-22s %5s %5s %5s %5s %5s %5s %5s %5s" %
      ("probe", "air", "noise", "sun", "liq", "slr", "tDT", "park", "food"))
for (n, la, lo) in probes:
    i = at(la, lo)
    print("%-22s %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f" % (
        n, scores["air"][i], scores["noise"][i], scores["microclimate"][i],
        scores["liquefaction"][i], scores["slr"][i], out["raw"]["t_dt"][i],
        scores["parks"][i], scores["food_nightlife"][i]))