tokyo-wangan-racing/reference/ramps/docs/osm_daikoku.py

#!/usr/bin/env python3
"""Pull real Daikoku JCT / PA geometry from OpenStreetMap via Overpass.

Outputs measured, real-world numbers:
  - PA service-area polygon (extent, area)
  - Each ramp/link way: length, turn direction, fitted circle radius
  - Bridge levels present

All distances in metres. Local ENU projection about the PA centre so we can
fit circles in a flat plane without map-projection distortion.
"""
import json, math, urllib.request, urllib.parse, sys

LAT0, LON0 = 35.46246657, 139.68149001   # user's Google Earth point
R_EARTH = 6378137.0

# Generous bbox around the whole interchange
S, W, N, E = 35.4545, 139.6700, 35.4705, 139.6935

Q = f"""
[out:json][timeout:120];
(
  way["highway"]({S},{W},{N},{E});
  way["amenity"="parking"]({S},{W},{N},{E});
  way["highway"="services"]({S},{W},{N},{E});
  way["highway"="rest_area"]({S},{W},{N},{E});
);
out body geom;
"""

ENDPOINTS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
]


def fetch():
    last = None
    for ep in ENDPOINTS:
        try:
            data = urllib.parse.urlencode({"data": Q}).encode()
            req = urllib.request.Request(ep, data=data,
                                         headers={"User-Agent": "daikoku-ref/1.0"})
            with urllib.request.urlopen(req, timeout=150) as r:
                return json.loads(r.read().decode())
        except Exception as e:
            last = f"{ep}: {e}"
            print(f"  endpoint failed {last}", file=sys.stderr)
    raise SystemExit(f"all overpass endpoints failed: {last}")


def enu(lat, lon):
    """Local east/north metres relative to (LAT0, LON0)."""
    x = math.radians(lon - LON0) * R_EARTH * math.cos(math.radians(LAT0))
    y = math.radians(lat - LAT0) * R_EARTH
    return x, y


def polyline_len(pts):
    return sum(math.dist(pts[i], pts[i + 1]) for i in range(len(pts) - 1))


def fit_circle(pts):
    """Algebraic (Kasa) circle fit. Returns (cx, cy, R, rms_residual)."""
    n = len(pts)
    if n < 3:
        return None
    sx = sum(p[0] for p in pts) / n
    sy = sum(p[1] for p in pts) / n
    u = [p[0] - sx for p in pts]
    v = [p[1] - sy for p in pts]
    Suu = sum(a * a for a in u); Svv = sum(b * b for b in v)
    Suv = sum(u[i] * v[i] for i in range(n))
    Suuu = sum(a ** 3 for a in u); Svvv = sum(b ** 3 for b in v)
    Suvv = sum(u[i] * v[i] * v[i] for i in range(n))
    Svuu = sum(v[i] * u[i] * u[i] for i in range(n))
    d = 2 * (Suu * Svv - Suv * Suv)
    if abs(d) < 1e-9:
        return None
    uc = (Svv * (Suuu + Suvv) - Suv * (Svvv + Svuu)) / d
    vc = (Suu * (Svvv + Svuu) - Suv * (Suuu + Suvv)) / d
    cx, cy = uc + sx, vc + sy
    rs = [math.dist(p, (cx, cy)) for p in pts]
    R = sum(rs) / len(rs)
    rms = math.sqrt(sum((r - R) ** 2 for r in rs) / len(rs))
    return cx, cy, R, rms


def total_turn(pts):
    """Signed cumulative heading change in degrees. +ccw / -cw."""
    tot = 0.0
    for i in range(1, len(pts) - 1):
        a1 = math.atan2(pts[i][1] - pts[i - 1][1], pts[i][0] - pts[i - 1][0])
        a2 = math.atan2(pts[i + 1][1] - pts[i][1], pts[i + 1][0] - pts[i][0])
        d = math.degrees(a2 - a1)
        while d > 180:  d -= 360
        while d < -180: d += 360
        tot += d
    return tot


def main():
    data = fetch()
    els = data.get("elements", [])
    print(f"overpass returned {len(els)} ways\n")

    ramps, areas = [], []
    for el in els:
        if el.get("type") != "way" or "geometry" not in el:
            continue
        t = el.get("tags", {})
        pts = [enu(g["lat"], g["lon"]) for g in el["geometry"]]
        if len(pts) < 2:
            continue
        rec = {"id": el["id"], "tags": t, "pts": pts,
               "len": polyline_len(pts)}
        hw = t.get("highway", "")
        if hw in ("services", "rest_area") or t.get("amenity") == "parking":
            areas.append(rec)
        elif hw.endswith("_link") or t.get("oneway") == "yes":
            ramps.append(rec)

    # ---- PA area ----
    print("=" * 70)
    print("PARKING / SERVICE AREA POLYGONS")
    print("=" * 70)
    for a in sorted(areas, key=lambda r: -r["len"])[:6]:
        xs = [p[0] for p in a["pts"]]; ys = [p[1] for p in a["pts"]]
        w = max(xs) - min(xs); h = max(ys) - min(ys)
        nm = a["tags"].get("name", a["tags"].get("name:en", "-"))
        print(f"  way {a['id']:>12}  {a['tags'].get('highway') or a['tags'].get('amenity'):<10}"
              f"  extent {w:6.0f} x {h:6.0f} m   perim {a['len']:7.0f} m   {nm}")
        print(f"       E {min(xs):7.0f}..{max(xs):7.0f}   N {min(ys):7.0f}..{max(ys):7.0f}")

    # ---- curved ramps: the loops ----
    print()
    print("=" * 70)
    print("CURVED RAMPS  (|turn| >= 120 deg)  -- the helix loops")
    print("=" * 70)
    print(f"  {'way id':>12} {'len(m)':>7} {'turn':>8} {'radius(m)':>10} {'rms':>6}  layer/bridge")
    loops = []
    for r in ramps:
        turn = total_turn(r["pts"])
        if abs(turn) < 120 or r["len"] < 60:
            continue
        fc = fit_circle(r["pts"])
        if not fc:
            continue
        cx, cy, R, rms = fc
        if R > 400:      # not a tight loop
            continue
        loops.append((r, turn, cx, cy, R, rms))
    for r, turn, cx, cy, R, rms in sorted(loops, key=lambda z: -abs(z[1])):
        t = r["tags"]
        lb = f"layer={t.get('layer','-')} bridge={t.get('bridge','-')}"
        print(f"  {r['id']:>12} {r['len']:7.0f} {turn:7.0f}\u00b0 {R:10.1f} {rms:6.1f}  {lb}")
        print(f"               centre E {cx:7.0f}  N {cy:7.0f}   "
              f"{'CCW/left' if turn>0 else 'CW/right'}")

    if loops:
        rr = [z[4] for z in loops]
        print(f"\n  loop radius: min {min(rr):.0f} m  max {max(rr):.0f} m  "
              f"mean {sum(rr)/len(rr):.0f} m")
        print(f"  in studs (1 stud = 0.28 m): mean {sum(rr)/len(rr)/0.28:.0f} studs")

    # ---- layers present ----
    print()
    print("=" * 70)
    print("DECK LAYERS (bridge ways)")
    print("=" * 70)
    lay = {}
    for r in ramps:
        if r["tags"].get("bridge"):
            L = r["tags"].get("layer", "0")
            lay.setdefault(L, []).append(r["len"])
    for L in sorted(lay, key=lambda s: int(s) if s.lstrip("-").isdigit() else 0):
        v = lay[L]
        print(f"  layer {L:>3}:  {len(v):3d} bridge ways, total {sum(v):8.0f} m")


if __name__ == "__main__":
    main()