"""Compute per-hex factor scores for San Francisco at 200 m and emit data2.json."""
import base64, csv, json, math, os
import numpy as np
from primitives import *

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

MX = 111320.0 * math.cos(math.radians(37.76))
MY = 110574.0

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

# ---------------------------------------------------------------- cells
grid = json.load(open(P("grid200.json")))
elev = {int(r["id"]): float(r["elev_m"]) for r in load("elevations200.csv")}
cells = []
for c in grid["cells"]:
    e = elev.get(c["id"])
    if e is None or e < -3.0:
        continue
    c = dict(c); c["elev"] = max(0.0, e)
    cells.append(c)
NC = len(cells)
print("cells kept: %d of %d" % (NC, len(grid["cells"])))

clat = np.array([c["lat"] for c in cells]); clon = np.array([c["lon"] for c in cells])
cx = clon * MX; cy = clat * MY
celev = np.array([c["elev"] for c in cells])

def pt_dists(lats, lons):
    """Distance in metres from every cell to every given point -> (NC, k)."""
    px = np.asarray(lons) * MX; py = np.asarray(lats) * MY
    return np.hypot(cx[:, None] - px[None, :], cy[:, None] - py[None, :])

def seg_dist(alat, alon, blat, blon):
    """Distance from every cell to one segment, plus interpolation position."""
    ax, ay = alon * MX, alat * MY; bx, by = blon * MX, blat * MY
    dx, dy = bx - ax, by - ay; L2 = dx * dx + dy * dy
    if L2 == 0:
        return np.hypot(cx - ax, cy - ay), np.zeros(NC)
    t = np.clip(((cx - ax) * dx + (cy - ay) * dy) / L2, 0.0, 1.0)
    return np.hypot(cx - (ax + t * dx), cy - (ay + t * dy)), t

def line_dist(pts):
    """Min distance to a polyline; returns (dist, seg_index, t)."""
    best = np.full(NC, 1e12); bi = np.zeros(NC, int); bt = np.zeros(NC)
    for i in range(len(pts) - 1):
        d, t = seg_dist(pts[i][0], pts[i][1], pts[i + 1][0], pts[i + 1][1])
        m = d < best
        best = np.where(m, d, best); bt = np.where(m, t, bt); bi = np.where(m, i, bi)
    return best, bi, bt

def in_rings(rings):
    """Boolean per cell: inside any of the given lon/lat rings (ray casting)."""
    inside = np.zeros(NC, bool)
    for ring in rings:
        r = np.asarray(ring)
        if len(r) < 3: continue
        xi, yi = r[:, 0], r[:, 1]
        xj, yj = np.roll(xi, 1), np.roll(yi, 1)
        cnt = np.zeros(NC, int)
        for k in range(len(r)):
            a, b = yi[k], yj[k]
            if a == b: continue
            straddle = ((a > clat) != (b > clat))
            xint = (xj[k] - xi[k]) * (clat - a) / (b - a) + xi[k]
            cnt += (straddle & (clon < xint)).astype(int)
        inside |= (cnt % 2 == 1)
    return inside

# ---------------------------------------------------------------- inputs
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"])) 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"])))
stops = [(r["mode"], float(r["lat"]), float(r["lon"])) for r in load("transit_stops.csv")]
playgrounds = [p for p in parks if "Playground" in p[3] or "Neighborhood Park" in p[3]]
print({k: len(v) for k, v in biz.items()}, "stops:", len(stops))

haz_liq = json.load(open(P("haz_liquefaction.json")))
haz_tsu = json.load(open(P("haz_tsunami.json")))
shore = json.load(open(P("base_shore.json")))
bstreets = json.load(open(P("base_streets.json")))
# base_parks2: 219 parks / 309 polygons / 5 holes, ~7 m tolerance, winding
# normalised (outers CCW, holes CW), MultiPolygon structure preserved.
bparks = json.load(open(P("base_parks2.json")))

# Federal / state / other non-city open space, from the California Protected
# Areas Database. Rec & Parks covers ONLY city-managed property, so without this
# the Presidio (1,314 ac), Ocean Beach, Fort Funston, Crissy Field, Lands End,
# Fort Mason, Mount Sutro and Baker Beach were all invisible to the park score —
# homes bordering a national park were scoring in the bottom third.
bfed = json.load(open(P("base_parks_fed.json")))
FED_SKIP_ACCESS = {"Restricted Access", "No Public Access"}

ALL_PARKS = []
for pk in bparks["parks"]:
    ALL_PARKS.append({"n": pk["n"], "acres": float(pk.get("acres") or 0.0),
                      "type": pk.get("type") or "", "polys": pk["polys"], "src": "city"})
for pk in bfed["parks"]:
    # A fenced reservoir is not a neighbourhood park; drop restricted access.
    if (pk.get("access") or "") in FED_SKIP_ACCESS:
        continue
    ALL_PARKS.append({"n": pk["n"], "acres": float(pk.get("acres") or 0.0),
                      "type": "", "polys": pk["polys"], "src": "fed"})
print("parks: %d city + %d non-city = %d (%.0f acres non-city)"
      % (len(bparks["parks"]), len(ALL_PARKS) - len(bparks["parks"]), len(ALL_PARKS),
         sum(p["acres"] for p in ALL_PARKS if p["src"] == "fed")))

def poly_dist(outer):
    """Distance from every cell to a polygon: 0 inside, else to nearest edge."""
    r = np.asarray(outer)
    if len(r) < 3:
        return np.full(NC, 1e12)
    edge = np.full(NC, 1e12)
    for k in range(len(r)):
        a, b = r[k], r[(k + 1) % len(r)]
        d, _ = seg_dist(a[1], a[0], b[1], b[0])
        edge = np.minimum(edge, d)
    xi, yi = r[:, 0], r[:, 1]
    xj, yj = np.roll(xi, 1), np.roll(yi, 1)
    cnt = np.zeros(NC, int)
    for k in range(len(r)):
        a, b = yi[k], yj[k]
        if a == b: continue
        straddle = ((a > clat) != (b > clat))
        xint = (xj[k] - xi[k]) * (clat - a) / (b - a) + xi[k]
        cnt += (straddle & (clon < xint)).astype(int)
    return np.where(cnt % 2 == 1, 0.0, edge)

# Snap curated supermarkets to the nearest measured grocery point within 350 m
gl = np.array([g[0] for g in biz["grocery"]]); gn = np.array([g[1] for g in biz["grocery"]])
SUPER = []
for name, la, lo in SUPERMARKETS:
    d = np.hypot((gn - lo) * MX, (gl - la) * MY)
    j = int(np.argmin(d))
    SUPER.append((name, gl[j], gn[j]) if d[j] <= 350 else (name, la, lo))

# ---------------------------------------------------------------- neighbourhoods
nbd = pt_dists([n[1] for n in NEIGHBORHOODS], [n[2] for n in NEIGHBORHOODS])
nb_i = np.argmin(nbd, axis=1)
PARK_NBHDS = {"Golden Gate Park", "Presidio", "McLaren Park", "Lincoln Park"}
nb_name = [NEIGHBORHOODS[i][0] for i in nb_i]
residential = np.array([0 if n in PARK_NBHDS else 1 for n in nb_name])

HEX_KM2 = (grid["spacing_m"] / 1000.0) ** 2 * math.sqrt(3) / 2
counts = np.bincount(nb_i, minlength=len(NEIGHBORHOODS))
dens = np.zeros((len(NEIGHBORHOODS), 3))
for i, nb in enumerate(NEIGHBORHOODS):
    area = max(0.25, counts[i] * HEX_KM2)
    dens[i] = (nb[3] / area, nb[4] / area, nb[5] / area)

# Crime and 311 are native to 41 neighbourhoods, not 200 m. Hard Voronoi
# assignment would stamp visible polygon edges onto a gradient map, so blend
# the 5 nearest centroids by inverse distance instead. The value is still
# neighbourhood-resolution — the blending only avoids implying a false edge.
K = 5
near = np.argsort(nbd, axis=1)[:, :K]
w = 1.0 / np.maximum(np.take_along_axis(nbd, near, axis=1), 250.0) ** 2
w /= w.sum(axis=1, keepdims=True)
viol_d = (dens[near, 0] * w).sum(axis=1)
prop_d = (dens[near, 1] * w).sum(axis=1)
c311_d = (dens[near, 2] * w).sum(axis=1)

# ---------------------------------------------------------------- slope
IDX = {}
for i in range(NC):
    IDX.setdefault((int(cx[i] // 400), int(cy[i] // 400)), []).append(i)
def neigh(i, radius):
    r = int(radius // 400) + 1
    gxi, gyi = int(cx[i] // 400), int(cy[i] // 400)
    out = []
    for a in range(gxi - r, gxi + r + 1):
        for b in range(gyi - r, gyi + r + 1):
            out.extend(IDX.get((a, b), []))
    o = np.array(out)
    d = np.hypot(cx[o] - cx[i], cy[o] - cy[i])
    return o[(d < radius) & (d > 1)], d[(d < radius) & (d > 1)]

slope = np.zeros(NC)
for i in range(NC):
    o, d = neigh(i, 320)
    if len(o):
        g = np.abs(celev[o] - celev[i]) / d
        slope[i] = (np.sort(g)[-2] if len(g) >= 2 else g[0]) * 100
print("slope p50 %.1f%%  p90 %.1f%%  max %.1f%%" % tuple(np.percentile(slope, [50, 90, 100])))

# ---------------------------------------------------------------- shore distance
shore_m = np.full(NC, 1e12)
for ring in shore["rings"]:
    for k in range(len(ring)):
        a, b = ring[k], ring[(k + 1) % len(ring)]
        mlat = (a[1] + b[1]) / 2; mlon = (a[0] + b[0]) / 2
        # skip the inland San Mateo county line: it is a land border, not shore
        if mlat < 37.7125 and -122.5000 < mlon < -122.3980:
            continue
        d, _ = seg_dist(a[1], a[0], b[1], b[0])
        shore_m = np.minimum(shore_m, d)

# ---------------------------------------------------------------- road exposure
def exposure(entries, near_m, far_m):
    """entries = [(strength, polyline)]. Returns (summed, max single road)."""
    tot = np.zeros(NC); mx = np.zeros(NC)
    for strength, pts in entries:
        d, _, _ = line_dist(pts)
        e = strength * (0.72 * np.exp(-d / near_m) + 0.28 * np.exp(-d / far_m))
        tot += e; mx = np.maximum(mx, e)
    return tot, mx

def aadt_of(name, fallback):
    v = ROAD_AADT.get(name)
    return (v[0] if v else fallback) / 100000.0

fw = [(aadt_of(n, 150000), [(p[0], p[1]) for p in pts]) for (n, _w, pts) in FREEWAYS]
ar = [(aadt_of(a[0], 15000), [(p[0], p[1]) for p in a[3]]) for a in ARTERIALS]
hn = [(min(1.0, a[2] / 8.0), [(p[0], p[1]) for p in a[3]]) for a in ARTERIALS if a[2] > 0]
cm = [(i, [(p[0], p[1]) for p in pts]) for (_n, i, pts) in COMMERCIAL]

fw_exp, fw_max = exposure(fw, 170.0, 520.0)
ar_exp, ar_max = exposure(ar, 110.0, 330.0)
hin_exp, hin_max = exposure(hn, 90.0, 260.0)
cm_exp, _ = exposure(cm, 150.0, 400.0)

# ---------------------------------------------------------------- transit
t_dt = np.full(NC, 1e9); t_via = [""] * NC
for name, mode, st in RAIL_LINES:
    d, bi, bt = line_dist([(s[0], s[1]) for s in st])
    ride = np.array([st[bi[i]][2] + bt[i] * (st[bi[i] + 1][2] - st[bi[i]][2]) for i in range(NC)])
    walk = d / 80.0
    tt = walk + ride + (0.0 if mode == "bart" else 1.5)
    tt = tt + np.where(walk > 14, (walk - 14) * 0.7, 0.0)
    m = tt < t_dt
    for i in np.where(m)[0]: t_via[i] = name
    t_dt = np.where(m, tt, t_dt)
for name, q, st in BUS_CORRIDORS:
    d, bi, bt = line_dist([(s[0], s[1]) for s in st])
    ride = np.array([st[bi[i]][2] + bt[i] * (st[bi[i] + 1][2] - st[bi[i]][2]) for i in range(NC)])
    tt = d / 80.0 + ride + (1.0 - q) * 9.0 + 3.0
    m = tt < t_dt
    for i in np.where(m)[0]: t_via[i] = name
    t_dt = np.where(m, tt, t_dt)

bd = pt_dists([b[0] for b in BART_TO_SFO], [b[1] for b in BART_TO_SFO])
bmin = np.array([b[2] for b in BART_TO_SFO])
walkable = bd / 80.0
t_sfo = np.min(np.where(walkable <= 20, walkable + bmin, 1e9), axis=1)
t_sfo = np.minimum(t_sfo, t_dt + 36)

rd = pt_dists([r[1] for r in SFO_RAMPS], [r[2] for r in SFO_RAMPS])
d_sfo = np.min(rd / 420.0 + 2.0 + np.array([r[3] for r in SFO_RAMPS]), axis=1)

# Transit stop access, now from 3,260 real SFMTA stops + 8 BART stations.
sd = pt_dists([s[1] for s in stops], [s[2] for s in stops])
is_bart = np.array([1.0 if s[0] == "bart" else 0.42 for s in stops])
stop_q = np.max(is_bart[None, :] * np.exp(-sd / 420.0), axis=1) + \
         0.05 * (np.exp(-sd / 300.0)).sum(axis=1).clip(0, 12)

# ---------------------------------------------------------------- amenities
def decayed(pts, scale, cap):
    if not pts: return np.zeros(NC)
    d = pt_dists([p[0] for p in pts], [p[1] for p in pts])
    return np.minimum(np.exp(-d / scale).sum(axis=1), cap)

def nearest(pts):
    return np.min(pt_dists([p[0] for p in pts], [p[1] for p in pts]), axis=1)

d_super = nearest([(s[1], s[2]) for s in SUPER])
# Two grocery questions, deliberately separate factors:
#   walk  — usable on foot, so decay hard past a ~10 min walk
#   drive — inside an estimated 5 min drive at city surface speeds (~2 km)
# Walk: logistic centred on ~900 m, which is the distance people actually stop
# being willing to carry groceries (roughly an 11 min walk). An exponential
# decay tuned tighter than this scored a comfortable 8 min walk at 26/100.
#   200 m -> 97 | 500 m -> 86 | 800 m -> 62 | 1.2 km -> 20 | 1.6 km -> 4
super_walk = 1.0 / (1.0 + np.exp((d_super - 900.0) / 220.0))
# Drive: city surface streets average ~25 km/h with signals, plus parking time.
drive_min = d_super / 420.0 + 1.5
super_drive = 1.0 / (1.0 + np.exp((drive_min - 5.0) * 1.5))

# Park access measured to the polygon BOUNDARY, not the centroid. Golden Gate
# Park's centroid sits ~2.5 km from its own edge, so centroid distance scored a
# house across the street from it as if the park were miles away. Distance is 0
# for a cell inside a park.
PLAY_TYPES = ("Playground", "Neighborhood Park", "Mini Park")
d_park = np.full(NC, 1e12)
d_play = np.full(NC, 1e12)
park_acc = np.zeros(NC)
park_dmin = {}
for pk in ALL_PARKS:
    wgt = min(1.0, math.sqrt(max(pk["acres"], 0.1)) / 7.0)
    dmin = np.full(NC, 1e12)
    for poly in pk["polys"]:
        dmin = np.minimum(dmin, poly_dist(poly["outer"]))
    park_dmin[pk["n"]] = dmin
    d_park = np.minimum(d_park, dmin)
    park_acc += wgt * np.exp(-dmin / 300.0)
    # Playgrounds stay city-only: a national seashore is not a play structure.
    if pk["src"] == "city" and any(t in pk["type"] for t in PLAY_TYPES):
        d_play = np.minimum(d_play, dmin)
print("park boundary distance: median %.0f m, p90 %.0f m, inside a park: %d cells"
      % (np.median(d_park), np.percentile(d_park, 90), (d_park == 0).sum()))

# Parkland flag, now driven by real geometry rather than a neighbourhood guess.
# A cell inside a LARGE park is open space, not housing stock, so it is excluded
# from the ranking. The 20-acre floor matters: a 200 m hex centred in a 2-acre
# plaza still has houses around its edges, one deep inside Golden Gate Park does
# not.
#
# The old rule flagged any cell whose nearest neighbourhood centroid was Golden
# Gate Park / McLaren / Lincoln Park, which is a crude Voronoi that wrongly
# excluded real housing along those parks' edges. Real polygons supersede it —
# and now that CPAD supplies the Presidio's true boundary, the last hand-rolled
# neighbourhood fallback is gone too. Every exclusion is geometric.
big_park = np.zeros(NC, bool)
for pk in ALL_PARKS:
    if pk["acres"] < 20.0:
        continue
    big_park |= (park_dmin[pk["n"]] == 0.0)
residential = np.where(big_park, 0, 1)
print("parkland excluded: %d cells, all from park geometry" % (residential == 0).sum())
d_school = nearest(schools)
d_lib = nearest(libraries)
d_hosp = nearest([(h[1], h[2]) for h in HOSPITALS])
corner = decayed(biz["convenience"] + biz["grocery"], 180.0, 8.0)
food = decayed(biz["restaurant"], 350.0, 60.0)
night = decayed(biz["bar"], 400.0, 14.0)
walk_poi = (decayed(biz["restaurant"], 420.0, 70.0)
            + 3.0 * decayed(biz["grocery"], 420.0, 12.0)
            + 2.0 * decayed(biz["convenience"], 420.0, 8.0)
            + decayed(biz["bar"], 420.0, 16.0))

# ---------------------------------------------------------------- hazards
liq = in_rings([p["ring"] for p in haz_liq["polygons"]]).astype(float)
tsu = in_rings([p["ring"] for p in haz_tsu["polygons"]]).astype(float)
print("in CGS liquefaction zone: %d cells (%.0f%%)" % (liq.sum(), 100 * liq.mean()))
print("in CGS tsunami hazard area: %d cells (%.0f%%)" % (tsu.sum(), 100 * tsu.mean()))

slr = np.clip((4.5 - celev) / 4.5, 0, 1) * np.exp(-shore_m / 900.0)
landslide = np.clip((slope - 10.0) / 16.0, 0, 1)

bowl = np.zeros(NC)
for i in range(NC):
    o, _d = neigh(i, 850)
    if len(o): bowl[i] = np.clip((celev[o].min() - celev[i]) / 7.0, 0, 1)
hot = np.zeros(NC)
for (_n, sev, la, lo, rad) in FLOOD_HOTSPOTS:
    d = np.hypot((clon - lo) * MX, (clat - la) * MY)
    hot = np.maximum(hot, sev * np.exp(-(d / rad) ** 2))
lowflat = ((celev < 5) & (slope < 3)).astype(float)
flood = np.minimum(1.0, 0.60 * bowl + 0.80 * hot + 0.25 * lowflat)

# ---------------------------------------------------------------- microclimate
F = FOG
sun = 1.0 - np.exp(-(clon - F["west_edge_lon"]) / F["decay_deg"])
lee = (clon > F["ridge_lon"]) & (clat >= F["lee_lat_range"][0]) & (clat <= F["lee_lat_range"][1])
sun += np.where(lee, F["lee_bonus"] * np.clip((clon - F["ridge_lon"]) / 0.020, 0, 1), 0)
sun -= np.where(clat > F["gate_lat"], F["gate_penalty"] * np.clip((clat - F["gate_lat"]) / 0.012, 0, 1), 0)
sun -= np.where(clon < F["coast_lon"], F["coast_penalty"] * np.clip((F["coast_lon"] - clon) / 0.014, 0, 1), 0)
sun = np.clip(sun, 0, 1.15)

# ---------------------------------------------------------------- factors
RAW = {
    "air":            -(fw_exp * 2.2 + ar_exp * 0.75 + ar_max * 2.6),
    "noise":          -(fw_exp * 1.8 + ar_exp * 0.9 + ar_max * 2.8 + cm_exp * 0.5),
    "microclimate":   sun,
    "slope":          -slope,
    "liquefaction":   -liq,
    "tsunami":        -tsu,
    "slr":            -slr,
    "flood":          -flood,
    "landslide":      -landslide,
    "transit_dt":     -t_dt,
    "transit_sfo":    -t_sfo,
    "drive_sfo":      -d_sfo,
    "walkability":    walk_poi,
    "bikeability":    -slope * 1.9 + walk_poi * 0.16 - hin_max * 22.0,
    "transit_score":  stop_q,
    "street_calm":    -(ar_exp * 0.9 + ar_max * 3.4 + fw_exp * 1.2),
    "ped_safety":     -(hin_exp * 0.7 + hin_max * 3.0),
    "supermarket_walk":  super_walk,
    "supermarket_drive": super_drive,
    "corner_store":   corner,
    # Distance dominates, size supplements. A pure acreage-weighted sum ranked a
    # block 110 m from a playground below median, because its neighbours are
    # small — but "is there a park I can walk to" is the primary question, and
    # the brief asked for fast distance decay. 65% proximity to the nearest park
    # of any size, 35% the weighted abundance of parkland around it.
    "parks":          0.65 * np.exp(-d_park / 250.0) + 0.35 * np.minimum(park_acc, 3.0) / 3.0,
    "food_nightlife": food + night * 2.2,
    "crime_violent":  -viol_d,
    "crime_property": -prop_d,
    "street_cond":    -c311_d,
    "schools":        -d_school,
    "playgrounds":    -d_play,
    "library":        -d_lib,
    "hospital":       -d_hosp,
}
# Hazard layers keep a linear 0-1 scale so 100 means "no modelled risk" rather
# than "best in the city"; rank-scaling would tie the majority of the city at
# an arbitrary mid value. Liquefaction and tsunami are statutory in/out zones,
# so they are genuinely binary — that is the data, not an artefact.
# Two different unit conventions, and mixing them up pins a factor at 100:
#   UNIT_NEG — raw is -severity in [-1, 0], so score = 100 * (1 + raw)
#   UNIT_POS — raw is already goodness in [0, 1], so score = 100 * raw
UNIT_NEG = {"liquefaction", "tsunami", "slr", "flood", "landslide"}
UNIT_POS = {"supermarket_walk", "supermarket_drive"}

def pct(v):
    order = np.argsort(v, kind="mergesort")
    ranks = np.empty(NC)
    r = np.arange(NC, dtype=float)
    ranks[order] = r
    sv = v[order]
    # average ranks within ties so equal values get equal scores
    i = 0
    while i < NC:
        j = i
        while j + 1 < NC and sv[j + 1] == sv[i]: j += 1
        if j > i: ranks[order[i:j + 1]] = (i + j) / 2.0
        i = j + 1
    return 100.0 * ranks / max(1, NC - 1)

scores = {}
for k, v in RAW.items():
    v = np.asarray(v, dtype=float)
    if k in UNIT_NEG:   s = 100.0 * (1.0 + v)
    elif k in UNIT_POS: s = 100.0 * v
    else:               s = pct(v)
    scores[k] = np.clip(np.rint(s), 0, 100).astype(np.uint8)

# ---------------------------------------------------------------- output
def b64(arr):
    return base64.b64encode(np.asarray(arr, dtype=np.uint8).tobytes()).decode()

def q(arr, mul=1):
    """Quantise a raw series to uint8 for compact transport."""
    a = np.clip(np.rint(np.asarray(arr, float) * mul), 0, 255)
    return b64(a)

via_names = sorted(set(t_via))
out = {
    "meta": {"spacing_m": grid["spacing_m"], "n": NC,
             "powell": POWELL_ST, "sfo": SFO,
             "sources": {"liq": haz_liq["source"], "liq_type": haz_liq["layer_type"],
                         "tsu": haz_tsu["source"], "tsu_scenario": haz_tsu["scenario"]}},
    "factors": [{"key": k, "label": l, "group": g, "prov": p, "w": wt, "res": rs, "note": nt}
                for (k, l, g, p, wt, rs, nt) in FACTORS],
    "cells": {
        "lat": [round(float(x), 5) for x in clat],
        "lon": [round(float(x), 5) for x in clon],
        "nbi": [int(i) for i in nb_i],
        "res": b64(residential),
    },
    "nbnames": [n[0] for n in NEIGHBORHOODS],
    "scores": {k: b64(v) for k, v in scores.items()},
    "raw": {
        "elev":   q(np.clip(celev, 0, 280) * (255 / 280)),
        "slope":  q(np.clip(slope, 0, 40) * (255 / 40)),
        "t_dt":   q(np.clip(t_dt, 0, 80) * (255 / 80)),
        "t_sfo":  q(np.clip(t_sfo, 0, 120) * (255 / 120)),
        "d_sfo":  q(np.clip(d_sfo, 0, 60) * (255 / 60)),
        "d_super": q(np.clip(d_super, 0, 2550) / 10),
        "d_park": q(np.clip(d_park, 0, 1275) / 5),
        "d_play": q(np.clip(d_play, 0, 2550) / 10),
        "d_school": q(np.clip(d_school, 0, 2550) / 10),
        "d_lib":  q(np.clip(d_lib, 0, 5100) / 20),
        "d_hosp": q(np.clip(d_hosp, 0, 5100) / 20),
        "food":   q(np.clip(food, 0, 60) * (255 / 60)),
        "night":  q(np.clip(night, 0, 14) * (255 / 14)),
        "corner": q(np.clip(corner, 0, 8) * (255 / 8)),
        "sun":    q(np.clip(sun, 0, 1.15) * (255 / 1.15)),
        "shore":  q(np.clip(shore_m, 0, 5100) / 20),
        "drivemin": q(np.clip(drive_min, 0, 25.5) * 10),
        "liq":    b64(liq * 255), "tsu": b64(tsu * 255),
        "slr":    q(slr * 255), "flood": q(flood * 255),
        "viol":   q(np.clip(viol_d, 0, 2550) / 10),
        "prop":   q(np.clip(prop_d, 0, 12750) / 50),
        "c311":   q(np.clip(c311_d, 0, 51000) / 200),
    },
    "rawscale": {"elev": 280 / 255, "slope": 40 / 255, "t_dt": 80 / 255, "t_sfo": 120 / 255,
                 "d_sfo": 60 / 255, "d_super": 10, "d_park": 5, "d_play": 10, "d_school": 10,
                 "d_lib": 20, "d_hosp": 20, "food": 60 / 255, "night": 14 / 255,
                 "corner": 8 / 255, "sun": 1.15 / 255, "shore": 20, "drivemin": 0.1,
                 "liq": 1 / 255, "tsu": 1 / 255, "slr": 1 / 255, "flood": 1 / 255,
                 "viol": 10, "prop": 50, "c311": 200},
    "via": {"names": via_names, "idx": b64([via_names.index(v) for v in t_via])},
    "base": {
        "shore": [[[round(p[0], 5), round(p[1], 5)] for p in r] for r in shore["rings"]],
        "streets": bstreets,
        # One entry per park, each keeping its own polygons and holes, so the
        # renderer can draw parks as separate paths. Merging every ring into a
        # single evenodd path made overlapping parks cancel into holes.
        # City and non-city parks both render, so the Presidio, Ocean Beach,
        # Lands End and Fort Funston stop showing as blank land on the basemap.
        "parks": [{"n": p["n"], "f": 1 if p["src"] == "fed" else 0,
                   "polys": [{"o": poly["outer"], "h": poly.get("holes", [])}
                             for poly in p["polys"]]}
                  for p in ALL_PARKS],
        "haz_liq": [p["ring"] for p in haz_liq["polygons"]],
        "haz_tsu": [p["ring"] for p in haz_tsu["polygons"]],
        "freeways": [{"n": n, "p": [[p[0], p[1]] for p in pts]} for (n, _w, pts) in FREEWAYS],
        "rail": [{"n": n, "m": m, "p": [[s[0], s[1]] for s in st]} for (n, m, st) in RAIL_LINES],
        "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},
        ],
    },
}
json.dump(out, open(P("data2.json"), "w"), separators=(",", ":"))
print("wrote data2.json: %d KB" % (os.path.getsize(P("data2.json")) // 1024))

# ---------------------------------------------------------------- sanity
def at(lat, lon):
    return int(np.argmin(np.hypot((clon - lon) * MX, (clat - lat) * MY)))
print("\n%-22s %4s %4s %4s %4s %4s %4s %4s %4s %4s" % (
    "probe", "air", "noiz", "sun", "liq", "tsu", "slr", "swlk", "sdrv", "tDT"))
for n, y, x in [("19th Ave & Noriega", 37.7537, -122.4757), ("2 blk E of 19th", 37.7537, -122.4700),
                ("Outer Sunset 45th", 37.7560, -122.5000), ("Valencia & 20th", 37.7585, -122.4210),
                ("Mission Bay", 37.7700, -122.3930), ("Marina", 37.8035, -122.4370),
                ("Twin Peaks", 37.7530, -122.4470), ("Bernal top", 37.7430, -122.4160),
                ("Dogpatch", 37.7570, -122.3885), ("Presidio Hts", 37.7880, -122.4530)]:
    i = at(y, x)
    print("%-22s %4d %4d %4d %4d %4d %4d %4d %4d %4.0f" % (
        n, scores["air"][i], scores["noise"][i], scores["microclimate"][i],
        scores["liquefaction"][i], scores["tsunami"][i], scores["slr"][i],
        scores["supermarket_walk"][i], scores["supermarket_drive"][i], t_dt[i]))