v0.6.0: city parking-map overlay + local time tracking, no IPS API
All checks were successful
build-apk / build (push) Successful in 10m38s
All checks were successful
build-apk / build (push) Successful in 10m38s
Adds the City of Sandpoint's printed "Downtown & Waterfront Public Parking" map as a georeferenced overlay, and lets you track your time on any of its areas without ever touching the ParkSmarter/IPS API. Georeferencing (tools/citymap/) - The PDF carries no geo metadata, so the page->WebMercator affine is recovered by fitting the drawing to OSM street centrelines. - pdftocairo writes stroked street segments with per-path matrix() transforms in local coords while filled lots are absolute; both are handled. The five legend swatches share the real geometry's colours and are identified by stroke-width and position, then dropped. - 49 areas, fitted to RMS 4.1 m (X) / 3.5 m (Y). On-street segments land a mean 4.0 m from the nearest OSM road. sp-039/040 sit further out because they are angled bays along the old rail corridor, on no named road at all. - Sandpoint's grid jogs 38 m between N 2nd Ave and S 2nd Ave; the page shows the same jog at the fitted scale, which independently confirms the fit. App - Map tab: "City map" layer in the legend's colours, tappable. - "Park here" pins the car from GPS and auto-detects the containing area (40 m snap). With no fix it asks you to tap the spot instead, so the pin never depends on GPS working. - The pin lives in its own storage key, not inside the session: pinning the car without starting a timer must survive backing out of the screen. - Durations cap at the posted limit — a 2-hour space is not offered a 4-hour timer. Lots and no-limit spots get the long options. - Reuses the existing foreground-service countdown. The second notification button reads "+1 hr" for a city area rather than "Extend": there is nothing to buy, so it edits the local timer and says so. - Account -> Align city map: nudge/scale/rotate the whole overlay against a live GPS fix. Save-on-phone needs no admin token, since the person who can see the misalignment is the one standing on the street. Server - parking_areas + map_overlay tables, public read, admin replace-all. The areas come from one source document, so replacement is wholesale rather than an upsert. Dropped geometryCenter from the geo module: on the real data it returns a point in the water for the crescent City Beach lot and mid-block for L-shaped runs. Nothing used it. Tests: 8 geometry tests in app/, 5 area/overlay tests in server/. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
parent
ad55559f55
commit
148e1635d3
23 changed files with 3027 additions and 21 deletions
72
tools/citymap/README.md
Normal file
72
tools/citymap/README.md
Normal file
|
|
@ -0,0 +1,72 @@
|
|||
# Georeferencing the city parking map
|
||||
|
||||
Turns the City of Sandpoint's printed **Downtown & Waterfront Public Parking** PDF into
|
||||
`app/src/features/citymap/parkingAreas.json` — the colour-coded overlay the app draws and
|
||||
hit-tests against.
|
||||
|
||||
Run this again when the city publishes a new edition of the map.
|
||||
|
||||
## Why it needs doing at all
|
||||
|
||||
The PDF carries **no** georeferencing metadata (no `/Measure`, `/GPTS`, `/LPTS`, `/GEO`,
|
||||
`/Viewport`, `/GCS`). It is a north-up Web Mercator screenshot of a slippy map with vector
|
||||
parking stripes drawn on top, so page coordinates relate to the world by a plain affine
|
||||
transform — which has to be recovered by fitting the drawing to something we already know
|
||||
the coordinates of. That something is OpenStreetMap's street centrelines.
|
||||
|
||||
Two structural details cost the most time, so they are worth knowing up front:
|
||||
|
||||
- `pdftocairo` writes each **stroked** street segment with its own `matrix()` transform and
|
||||
*local* coordinates, while **filled** lots are in absolute page coordinates. Both have to
|
||||
be handled or the streets land in a heap near the origin.
|
||||
- The legend swatches are drawn in the same five colours as the real geometry. They are
|
||||
identified by stroke-width 7 inside the legend card's x-band and dropped.
|
||||
|
||||
## Pipeline
|
||||
|
||||
```bash
|
||||
# 0. deps: poppler-utils (pdftocairo, pdftotext, pdfimages), python3, curl
|
||||
pdftocairo -svg downtown_and_waterfront_public_parking_map.pdf map.svg
|
||||
|
||||
# 1. vector geometry -> page coordinates, bucketed by the legend's five colours
|
||||
python3 extract_map.py map.svg map_page_coords.json
|
||||
|
||||
# 2. OSM street centrelines for downtown Sandpoint
|
||||
curl -s --data-binary @roads.overpass https://overpass-api.de/api/interpreter -o osm.json
|
||||
|
||||
# 3. fit page -> Web Mercator against named streets; prints per-street residuals
|
||||
python3 georef.py # writes fit_raw.json
|
||||
|
||||
# 4. apply the fit, name each area from OSM, verify, emit the GeoJSON
|
||||
python3 build_geojson.py # writes parking_areas.geojson
|
||||
|
||||
cp parking_areas.geojson ../../app/src/features/citymap/parkingAreas.json
|
||||
```
|
||||
|
||||
## What "good" looks like
|
||||
|
||||
`georef.py` prints a residual per control street and `build_geojson.py` prints how far each
|
||||
on-street segment sits from the nearest OSM road. The current edition fits to:
|
||||
|
||||
| Check | Result |
|
||||
| --- | --- |
|
||||
| X control residual (avenues) | **RMS 4.1 m** |
|
||||
| Y control residual (streets) | **RMS 3.5 m** |
|
||||
| On-street segments vs nearest OSM road | **mean 4.0 m**, 40 of 42 under 10 m |
|
||||
|
||||
The two segments over 10 m (`sp-039`, `sp-040`) are correct, not errors: they are angled
|
||||
bays along the old rail corridor that sit on no named road at all — the nearest way is a
|
||||
service alley 19 m off. Anything much worse than the table above means a control street was
|
||||
mis-identified; `georef.py`'s per-street residuals will say which.
|
||||
|
||||
Residual error is also correctable after the fact without re-running any of this — the app's
|
||||
**Account → Align city map** screen nudges the whole overlay against a live GPS fix and
|
||||
persists the correction.
|
||||
|
||||
## Control points
|
||||
|
||||
`georef.py` maps page grid lines to OSM street names by hand (`AVENUES` / `STREETS`). Sandpoint's
|
||||
grid **jogs** between its north and south halves — North 2nd Ave and South 2nd Ave are 38 m
|
||||
apart — so each control street is measured only over the span its page segment actually
|
||||
covers, and both halves are used as independent control points. That jog is a useful sanity
|
||||
check: the page shows the same 9.3 pt offset, which at the fitted scale is 38 m.
|
||||
190
tools/citymap/build_geojson.py
Normal file
190
tools/citymap/build_geojson.py
Normal file
|
|
@ -0,0 +1,190 @@
|
|||
#!/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")
|
||||
195
tools/citymap/extract_map.py
Normal file
195
tools/citymap/extract_map.py
Normal file
|
|
@ -0,0 +1,195 @@
|
|||
#!/usr/bin/env python3
|
||||
"""Extract the color-coded parking geometry from the city's parking-map PDF.
|
||||
|
||||
pdftocairo emits every street segment as a stroked <path> carrying its own
|
||||
matrix() transform with local coordinates, and every lot as an absolute filled
|
||||
<path>. So: parse the path, flatten curves, push through the path's own matrix,
|
||||
and bucket by the exact colour pdftocairo wrote.
|
||||
|
||||
Output is GeoJSON-shaped but still in SVG *page* coordinates (y down); the
|
||||
georeferencing step turns that into lat/lon.
|
||||
"""
|
||||
import json
|
||||
import re
|
||||
import sys
|
||||
from xml.etree import ElementTree as ET
|
||||
|
||||
SVG = sys.argv[1] if len(sys.argv) > 1 else "/tmp/map.svg"
|
||||
OUT = sys.argv[2] if len(sys.argv) > 2 else "map_page_coords.json"
|
||||
|
||||
# Colours exactly as pdftocairo writes them, mapped to the map legend.
|
||||
COLORS = {
|
||||
"rgb(76.499939%, 76.899719%, 76.098633%)": "no_limit",
|
||||
"rgb(96.899414%, 54.499817%, 3.09906%)": "limit_4h",
|
||||
"rgb(79.998779%, 77.2995%, 25.898743%)": "limit_3h",
|
||||
"rgb(82.699585%, 40.39917%, 79.998779%)": "free_2h",
|
||||
"rgb(45.899963%, 69.799805%, 34.899902%)": "green_lot",
|
||||
}
|
||||
|
||||
NUM = re.compile(r"[-+]?(?:\d*\.\d+|\d+\.?)(?:[eE][-+]?\d+)?")
|
||||
|
||||
|
||||
def parse_matrix(s):
|
||||
"""matrix(a,b,c,d,e,f) -> tuple. Identity when absent."""
|
||||
if not s:
|
||||
return (1.0, 0.0, 0.0, 1.0, 0.0, 0.0)
|
||||
m = re.search(r"matrix\s*\(([^)]*)\)", s)
|
||||
if not m:
|
||||
return (1.0, 0.0, 0.0, 1.0, 0.0, 0.0)
|
||||
v = [float(x) for x in NUM.findall(m.group(1))]
|
||||
return tuple(v[:6])
|
||||
|
||||
|
||||
def apply(mtx, x, y):
|
||||
a, b, c, d, e, f = mtx
|
||||
return (a * x + c * y + e, b * x + d * y + f)
|
||||
|
||||
|
||||
def bezier(p0, p1, p2, p3, steps=8):
|
||||
"""Flatten a cubic to points. The map's curves are gentle; 8 is plenty."""
|
||||
out = []
|
||||
for i in range(1, steps + 1):
|
||||
t = i / steps
|
||||
u = 1 - t
|
||||
out.append(
|
||||
(
|
||||
u * u * u * p0[0] + 3 * u * u * t * p1[0] + 3 * u * t * t * p2[0] + t * t * t * p3[0],
|
||||
u * u * u * p0[1] + 3 * u * u * t * p1[1] + 3 * u * t * t * p2[1] + t * t * t * p3[1],
|
||||
)
|
||||
)
|
||||
return out
|
||||
|
||||
|
||||
def parse_path(d):
|
||||
"""Return a list of subpaths [(points, closed)] in the path's local space."""
|
||||
tokens = re.findall(r"([MmLlHhVvCcSsZz])|([-+]?(?:\d*\.\d+|\d+\.?)(?:[eE][-+]?\d+)?)", d)
|
||||
subpaths, pts = [], []
|
||||
cur = (0.0, 0.0)
|
||||
start = (0.0, 0.0)
|
||||
cmd = None
|
||||
i = 0
|
||||
flat = [(c, n) for c, n in tokens]
|
||||
|
||||
def nums(k):
|
||||
nonlocal i
|
||||
vals = []
|
||||
while len(vals) < k and i < len(flat) and flat[i][1]:
|
||||
vals.append(float(flat[i][1]))
|
||||
i += 1
|
||||
return vals
|
||||
|
||||
while i < len(flat):
|
||||
c, n = flat[i]
|
||||
if c:
|
||||
cmd = c
|
||||
i += 1
|
||||
elif cmd is None:
|
||||
i += 1
|
||||
continue
|
||||
if cmd in "Zz":
|
||||
if pts:
|
||||
subpaths.append((pts, True))
|
||||
pts = []
|
||||
cur = start
|
||||
cmd = None
|
||||
continue
|
||||
if cmd in "Mm":
|
||||
v = nums(2)
|
||||
if len(v) < 2:
|
||||
break
|
||||
if pts:
|
||||
subpaths.append((pts, False))
|
||||
pts = []
|
||||
cur = (v[0], v[1]) if cmd == "M" else (cur[0] + v[0], cur[1] + v[1])
|
||||
start = cur
|
||||
pts = [cur]
|
||||
cmd = "L" if cmd == "M" else "l"
|
||||
elif cmd in "Ll":
|
||||
v = nums(2)
|
||||
if len(v) < 2:
|
||||
break
|
||||
cur = (v[0], v[1]) if cmd == "L" else (cur[0] + v[0], cur[1] + v[1])
|
||||
pts.append(cur)
|
||||
elif cmd in "Hh":
|
||||
v = nums(1)
|
||||
if not v:
|
||||
break
|
||||
cur = (v[0], cur[1]) if cmd == "H" else (cur[0] + v[0], cur[1])
|
||||
pts.append(cur)
|
||||
elif cmd in "Vv":
|
||||
v = nums(1)
|
||||
if not v:
|
||||
break
|
||||
cur = (cur[0], v[0]) if cmd == "V" else (cur[0], cur[1] + v[0])
|
||||
pts.append(cur)
|
||||
elif cmd in "Cc":
|
||||
v = nums(6)
|
||||
if len(v) < 6:
|
||||
break
|
||||
if cmd == "C":
|
||||
p1, p2, p3 = (v[0], v[1]), (v[2], v[3]), (v[4], v[5])
|
||||
else:
|
||||
p1 = (cur[0] + v[0], cur[1] + v[1])
|
||||
p2 = (cur[0] + v[2], cur[1] + v[3])
|
||||
p3 = (cur[0] + v[4], cur[1] + v[5])
|
||||
pts.extend(bezier(cur, p1, p2, p3))
|
||||
cur = p3
|
||||
else:
|
||||
i += 1
|
||||
if pts:
|
||||
subpaths.append((pts, False))
|
||||
return subpaths
|
||||
|
||||
|
||||
def main():
|
||||
tree = ET.parse(SVG)
|
||||
root = tree.getroot()
|
||||
ns = "{http://www.w3.org/2000/svg}"
|
||||
|
||||
features = []
|
||||
for el in root.iter(ns + "path"):
|
||||
stroke = (el.get("stroke") or "").strip()
|
||||
fill = (el.get("fill") or "").strip()
|
||||
kind = COLORS.get(stroke) or COLORS.get(fill)
|
||||
if not kind:
|
||||
continue
|
||||
is_stroke = stroke in COLORS
|
||||
mtx = parse_matrix(el.get("transform"))
|
||||
width = float(el.get("stroke-width") or 0)
|
||||
for local, closed in parse_path(el.get("d") or ""):
|
||||
world = [apply(mtx, x, y) for x, y in local]
|
||||
if len(world) < 2:
|
||||
continue
|
||||
xs = [p[0] for p in world]
|
||||
ys = [p[1] for p in world]
|
||||
features.append(
|
||||
{
|
||||
"kind": kind,
|
||||
"geom": "line" if is_stroke else "polygon",
|
||||
"closed": closed,
|
||||
"strokeWidth": width,
|
||||
"bbox": [min(xs), min(ys), max(xs), max(ys)],
|
||||
"points": [[round(x, 3), round(y, 3)] for x, y in world],
|
||||
}
|
||||
)
|
||||
|
||||
with open(OUT, "w") as fh:
|
||||
json.dump({"viewBox": [0, 0, 491.87, 529.043], "features": features}, fh, indent=1)
|
||||
|
||||
from collections import Counter
|
||||
print(f"{len(features)} features -> {OUT}")
|
||||
for (k, g), n in sorted(Counter((f["kind"], f["geom"]) for f in features).items()):
|
||||
print(f" {k:10s} {g:8s} {n}")
|
||||
# Where are they? Legend swatches cluster in one corner; real geometry spreads out.
|
||||
print("\nbbox spread by kind:")
|
||||
for k in COLORS.values():
|
||||
fs = [f for f in features if f["kind"] == k]
|
||||
if not fs:
|
||||
continue
|
||||
print(
|
||||
f" {k:10s} x[{min(f['bbox'][0] for f in fs):7.1f},{max(f['bbox'][2] for f in fs):7.1f}] "
|
||||
f"y[{min(f['bbox'][1] for f in fs):7.1f},{max(f['bbox'][3] for f in fs):7.1f}]"
|
||||
)
|
||||
|
||||
|
||||
main()
|
||||
114
tools/citymap/georef.py
Normal file
114
tools/citymap/georef.py
Normal file
|
|
@ -0,0 +1,114 @@
|
|||
#!/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)
|
||||
3
tools/citymap/roads.overpass
Normal file
3
tools/citymap/roads.overpass
Normal file
|
|
@ -0,0 +1,3 @@
|
|||
[out:json][timeout:60];
|
||||
way["highway"](48.2600,-116.5750,48.2900,-116.5300);
|
||||
out geom;
|
||||
Loading…
Add table
Add a link
Reference in a new issue