#!/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()