BigBrainParking/tools/citymap/georef.py

#!/usr/bin/env python3
"""Georeference the parking-map page coordinates against OSM.

The base map is a north-up Web Mercator screenshot, so page -> mercator is a
uniform scale plus a translation (3 free params, not 6). Control points are
street centrelines identified by name: an avenue pins X, a street pins Y.

Streets in Sandpoint jog between the north and south halves of the grid, so each
control street's mercator coordinate is measured only over the span the page
segment actually covers, not over the whole way.
"""
import json
import math

R = 6378137.0


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


def unmerc(X, Y):
    return (math.degrees(2 * math.atan(math.exp(Y / R)) - math.pi / 2), math.degrees(X / R))


osm = json.load(open("osm.json"))
ways = {}
for w in osm["elements"]:
    n = w.get("tags", {}).get("name")
    if not n or "geometry" not in w:
        continue
    ways.setdefault(n, []).append([merc(p["lat"], p["lon"]) for p in w["geometry"]])


def centreline(name, axis, lo, hi):
    """Mean coordinate on `axis` of `name`, over the other axis' [lo,hi] window."""
    other = 1 - axis
    vals = []
    for g in ways.get(name, []):
        for (x0, y0), (x1, y1) in zip(g, g[1:]):
            p0, p1 = (x0, y0), (x1, y1)
            if not (lo <= p0[other] <= hi or lo <= p1[other] <= hi):
                continue
            vals.append((p0[axis] + p1[axis]) / 2)
    return sum(vals) / len(vals) if vals else None


# Page grid lines read off the rendered map, with the mercator window each spans.
# X window for E-W streets / Y window for N-S avenues, in mercator metres.
XW = (-12974700, -12974000)  # 5th Ave .. 1st Ave
YW = (6152200, 6153300)  # Lake St .. Poplar St

AVENUES = [  # page x, OSM name
    (83.1, "North 5th Avenue"),
    (120.0, "North 4th Avenue"),
    (163.8, "North 3rd Avenue"),
    (207.7, "North 2nd Avenue"),
    (235.7, "North 1st Avenue"),
    (164.3, "South 3rd Avenue"),
    (198.4, "South 2nd Avenue"),
]
STREETS = [  # page y, OSM name
    (184.4, "Poplar Street"),
    (228.3, "Alder Street"),
    (272.3, "Cedar Street"),
    (314.9, "Oak Street"),
    (359.1, "Church Street"),
    (398.0, "Pine Street"),
    (438.2, "Lake Street"),
    (492.1, "Superior Street"),
]


def fit(pairs, flip):
    """Least-squares v = s*p + t. Returns (s, t, residuals)."""
    n = len(pairs)
    sp = sum(p for p, v in pairs)
    sv = sum(v for p, v in pairs)
    spp = sum(p * p for p, v in pairs)
    spv = sum(p * v for p, v in pairs)
    s = (n * spv - sp * sv) / (n * spp - sp * sp)
    t = (sv - s * sp) / n
    return s, t, [(p, v, s * p + t - v) for p, v in pairs]


ax = [(px, centreline(n, 0, *YW)) for px, n in AVENUES]
ay = [(py, centreline(n, 1, *XW)) for py, n in STREETS]
print("control points (mercator metres):")
for (px, n), (_, v) in zip(AVENUES, ax):
    print(f"  x {px:7.1f}  {n:20s} {v if v is None else round(v,1)}")
for (py, n), (_, v) in zip(STREETS, ay):
    print(f"  y {py:7.1f}  {n:20s} {v if v is None else round(v,1)}")

ax = [(p, v) for p, v in ax if v is not None]
ay = [(p, v) for p, v in ay if v is not None]

COS = math.cos(math.radians(48.278))  # mercator metres -> ground metres here


def report(label, pairs, names):
    s, t, res = fit(pairs, False)
    print(f"\n{label}: scale={s:.4f} merc-m/pt  ({abs(s)*COS:.4f} ground-m/pt), offset={t:.1f}")
    for (p, v, r), nm in zip(res, names):
        print(f"   {nm:20s} page={p:7.1f} residual={r*COS:7.1f} ground-m")
    rms = math.sqrt(sum(r * r for _, _, r in res) / len(res)) * COS
    print(f"   RMS = {rms:.1f} ground-m")
    return s, t, rms


sx, tx, rx = report("X (avenues)", ax, [n for _, n in AVENUES])
sy, ty, ry = report("Y (streets)", ay, [n for _, n in STREETS])
print(f"\nscale ratio |sy/sx| = {abs(sy/sx):.4f}  (1.0 == truly uniform / north-up)")

json.dump({"sx": sx, "tx": tx, "sy": sy, "ty": ty}, open("fit_raw.json", "w"), indent=1)