BigBrainParking/tools/citymap/build_geojson.py

#!/usr/bin/env python3
"""Apply the fitted transform and emit the final parking-areas GeoJSON.

Also verifies the result the only way that matters: every stroked segment should
land on an actual road, so measure each one's distance to the nearest OSM road
centreline. Lots are skipped in that check — they are off-street by definition.

Each feature gets a human name from OSM (the street it runs along, plus the two
cross streets it lies between) so the app can list areas without the map.
"""
import json
import math

R = 6378137.0
COS = math.cos(math.radians(48.278))

fit = json.load(open("fit_raw.json"))
SX, TX, SY, TY = fit["sx"], fit["tx"], fit["sy"], fit["ty"]


def to_merc(x, y):
    return (SX * x + TX, SY * y + TY)


def to_lonlat(x, y):
    X, Y = to_merc(x, y)
    return (round(math.degrees(X / R), 7), round(math.degrees(2 * math.atan(math.exp(Y / R)) - math.pi / 2), 7))


def merc(lat, lon):
    return (math.radians(lon) * R, math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) * R)


def dist_to_seg(p, a, b):
    px, py = p
    ax, ay = a
    bx, by = b
    dx, dy = bx - ax, by - ay
    L = dx * dx + dy * dy
    t = 0.0 if L == 0 else max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / L))
    return math.hypot(px - (ax + t * dx), py - (ay + t * dy))


# ---------------------------------------------------------------- OSM roads
osm = json.load(open("osm.json"))
SKIP = {"footway", "path", "cycleway", "steps", "track", "service"}
roads = []  # (a, b, name) in mercator
for w in osm["elements"]:
    t = w.get("tags", {})
    if t.get("highway") in SKIP or "geometry" not in w:
        continue
    g = [merc(p["lat"], p["lon"]) for p in w["geometry"]]
    nm = t.get("name")
    for a, b in zip(g, g[1:]):
        roads.append((a, b, nm))
named = [r for r in roads if r[2]]

SHORT = [
    ("North ", "N "), ("South ", "S "), ("East ", "E "), ("West ", "W "),
    (" Street", " St"), (" Avenue", " Ave"), (" Boulevard", " Blvd"),
    (" Road", " Rd"), (" Drive", " Dr"), (" Lane", " Ln"), (" Bridge", " Brg"),
]


def short(n):
    for a, b in SHORT:
        n = n.replace(a, b)
    return n


def nearest_name(p, exclude=None, limit=60.0):
    best, bestd = None, limit
    for a, b, nm in named:
        if nm == exclude:
            continue
        d = dist_to_seg(p, a, b)
        if d < bestd:
            best, bestd = nm, d
    return best


def describe(pts_merc, is_line):
    """'N 3rd Ave · Cedar St to Oak St' for a segment, or the nearest road for a lot."""
    mid = pts_merc[len(pts_merc) // 2]
    on = nearest_name(mid) if is_line else None
    if not is_line:
        near = nearest_name(mid, limit=200.0)
        return f"Lot off {short(near)}" if near else "City lot"
    ends = [pts_merc[0], pts_merc[-1]]
    cross = []
    for e in ends:
        c = nearest_name(e, exclude=on, limit=45.0)
        if c and short(c) not in cross:
            cross.append(short(c))
    if not on:
        # Bays along the old rail corridor sit on no named road — describe them
        # by what they are near rather than inventing a street.
        near = nearest_name(mid, limit=150.0)
        return f"Off-street bays near {short(near)}" if near else "Off-street bays"
    base = short(on)
    if len(cross) == 2:
        return f"{base} · {cross[0]} to {cross[1]}"
    if len(cross) == 1:
        return f"{base} · at {cross[0]}"
    return base


# --------------------------------------------------------------- build output
page = json.load(open("map_page_coords.json"))


def is_legend(f):
    """The five legend swatches: stroke-width 7 sitting in the legend card's x-band."""
    x0, y0, x1, y1 = f["bbox"]
    return f["strokeWidth"] > 5 and 25 < x0 < 28 and 60 < y0 < 130


# kind -> (short label, legend text, default tracked hours, colour)
KINDS = {
    "green_lot": ("City lot", "Paid hourly or permit", 2, "#75b259"),
    "free_2h": ("2-hour free", "Permits not valid", 2, "#d367cc"),
    "limit_3h": ("3-hour", "3-hour or permit", 3, "#ccc542"),
    "limit_4h": ("4-hour", "4-hour or permit", 4, "#f78b08"),
    "no_limit": ("No time limit", "No posted time limit", 0, "#c3c4c2"),
}

features = []
n_legend = 0
for i, f in enumerate(page["features"]):
    if is_legend(f):
        n_legend += 1
        continue
    is_line = f["geom"] == "line"
    pm = [to_merc(x, y) for x, y in f["points"]]
    coords = [to_lonlat(x, y) for x, y in f["points"]]
    if is_line:
        geom = {"type": "LineString", "coordinates": coords}
    else:
        if coords[0] != coords[-1]:
            coords.append(coords[0])
        geom = {"type": "Polygon", "coordinates": [coords]}
    label, legend, hours, color = KINDS[f["kind"]]
    features.append(
        {
            "type": "Feature",
            "id": f"sp-{i:03d}",
            "geometry": geom,
            "properties": {
                "id": f"sp-{i:03d}",
                "kind": f["kind"],
                "label": label,
                "legend": legend,
                "hours": hours,
                "color": color,
                "shape": "line" if is_line else "polygon",
                "name": describe(pm, is_line),
            },
        }
    )

fc = {
    "type": "FeatureCollection",
    "features": features,
    "metadata": {
        "source": "City of Sandpoint — Downtown & Waterfront Public Parking map",
        "generated": "from downtown_and_waterfront_public_parking_map.pdf",
        "georeference": "affine page->WebMercator fitted to OSM street centrelines",
    },
}
json.dump(fc, open("parking_areas.geojson", "w"), indent=1)
print(f"{len(features)} features written ({n_legend} legend swatches dropped)\n")
for f in features:
    p = f["properties"]
    print(f"  {p['id']} {p['kind']:10s} {p['shape']:7s} {p['name']}")

# ---- verification: distance from each on-street segment to the nearest road
worst = []
for f in features:
    if f["properties"]["shape"] != "line":
        continue
    ds = [min(dist_to_seg(merc(lat, lon), a, b) for a, b, _ in roads) * COS
          for lon, lat in f["geometry"]["coordinates"]]
    worst.append((max(ds), sum(ds) / len(ds), f["properties"]["id"], f["properties"]["kind"]))

worst.sort(reverse=True)
print(f"\non-street segments: {len(worst)}, mean offset from nearest road = "
      f"{sum(m for _, m, _, _ in worst)/len(worst):.1f} m")
print(f"segments with mean offset > 10 m: {sum(1 for _, m, _, _ in worst if m > 10)}")
for mx, mn, fid, kind in worst[:4]:
    print(f"  worst: {fid} {kind:10s} max={mx:6.1f} m mean={mn:6.1f} m")