#!/usr/bin/env python3 """Build street_limits.bin: the offline posted-speed-limit network for the Speed Cams app. Run by a developer before a release, never on the phone and never as an Xcode build phase. Standard library only, Python 3.9+ (run it with `python3 -I`). The app bundles the output and matches GPS fixes against it on the device; no location ever leaves the phone. python3 -I build_speed_limits.py build # all sources, OSM fills the gaps python3 -I build_speed_limits.py build --no-osm # official sources only (no ODbL data) python3 -I build_speed_limits.py fetch # download / refresh the cache only python3 -I build_speed_limits.py verify out/street_limits.bin python3 -I build_speed_limits.py lookup out/street_limits.bin 34.1649 -118.3900 --course 90 python3 -I build_speed_limits.py selftest SOURCES AND PRECEDENCE (speed limits are treated as public facts; the LADOT layers are the primary source, credited as "Speed limit data: City of Los Angeles Department of Transportation"; we store our own normalized representation, one limit per road piece, and never ship the raw files). A lower layer only contributes the parts of its lines that no higher layer already covers ("carving": 5 m samples, a sample is covered when an accepted line lies within the carve distance at a bearing within 30 deg). Every LADOT layer is clipped to the City of Los Angeles boundary (LA County eGIS city boundaries; a line is kept up to 15 m outside it, so boundary roads stay), and the base-map rows LADOT itself marks SECT_ID 'outside' are also used as a mask for the LADOT add steps. 0. Official documents (street_limits_overrides.json, each entry with a citation and a review date): City of Glendale GMC 10.24.020 values (source 8, County road centerline geometry) and LADOT Speed Safety System program documents (source 9, applied to the LADOT records of the named street between two cross streets). 1. LADOT SpeedLimits_LAMC_35mph (Jan 2025; LA City streets posted <= 35 mph). Overrides the base map per 5 m sample: every part of a base segment that runs within 15 m of a LAMC line at a bearing within 30 deg takes the LAMC value (the segment is split at the ends of the LAMC coverage). LAMC lines that match no base segment are added. 2. LADOT City_of_Los_Angeles_Speed_Limit_Map (LAMC 80.81 listed streets, 2018 base), patched by SECT_ID with the AB 43 change layer (Progress_Status = 'Complete'); AB 43 rows whose SECT_ID is not in the base also lower the base parts they run along. AB 43 never raises a limit. Segments whose SECT_ID is 'outside' are dropped. Rows without a SECT_ID are carved by the rows that have one. A base part at 35 mph or less that no LAMC line covers (also after a 40 m retry with the same street name) is a probable stale 2018 value: it is removed where an OSM way with a higher value runs along it (the OSM value fills it), and kept elsewhere (see item 6). 3. LADOT SpeedLimits_CVC: LA City streets not in the LAMC, 25 mph prima facie (LADOT's own labelling). Rows whose SECT_ID is in the base map are dropped, and so are walk streets, vacated streets, Rec & Parks roads, Housing Authority roads, private and state-highway rows (ST_TYPE WK, VA, RP, HA, PA, ST). A CVC part is also dropped where a faster official record of the same street runs alongside within 25 m (frontage roads and duplicate centerlines), and where, on the city edge (within 30 m of the boundary), a County line with a posted value other than 25 runs alongside within 15 m. 4. LA County DPW Speed Network (County roads, mostly unincorporated areas; --no-county leaves it out). Rows with conflicting values ("25 and 35") carry no limit. 5. Caltrans D7 SHS posted speed limits: freeway mainlines only (route 2 and route 105 only inside their freeway sections). Its surface-highway values are not used by default (--caltrans-surface uses them anyway). Where a Caltrans line is not used, an OSM motorway or trunk running along it may not carry a value above the Caltrans one. 6. OpenStreetMap maxspeed ("N mph" only) fills the remaining gaps, (c) OpenStreetMap contributors, ODbL 1.0. Untagged OSM roads in the gaps are kept with no limit (source 0) so the matcher knows the driver is on a road with no data instead of snapping to a nearby road that has one. Every OSM-derived record has source 0 or 1. Safety exception (off with --strict-precedence): LADOT's CVC 25 is a prima facie default, and a lower posted limit governs; where an OSM way with a limit below 25 runs along a CVC piece, that part of the CVC piece is removed and the OSM value is used; likewise a probable stale base part (item 2) is removed where OSM has a higher value. Because these rules let OSM decide which official parts are kept, the official records of a default build are not fully independent of OSM: to remove OSM, rebuild with --no-osm (which also turns the rules off); do not just filter codes 0 and 1 out. Official records that still overlap within 8 m with different limits are cut in favour of the higher-ranked source (or the longer record) before OSM is added. A suppress list (street_limits_suppress.json) can only remove OSM values. 7. Nothing else: no invented defaults (the only 25 mph defaults are LADOT's CVC layer). School-zone and time-of-day limits are not in any source and are not represented (maxspeed:conditional is only flagged). OUTPUTS (in --out): street_limits.bin, street_limits.manifest.json, street_limits.NOTICE.txt, street_limits.stats.json / .stats.txt (the report, with every gate result) and street_limits.review.json (manual-review lists that change nothing in the file: Caltrans vs OSM freeway disagreement sections, LADOT CVC 25 segments where OSM says 35 or more, and OSM records with a limit in one direction only). BUILD GATES (the build fails, and keeps the previous output, when one fails; the stats report lists every finding; --allow-gate-failures writes the file anyway for diagnosis): - overlapping official records: two official records running within 8 m of each other at a bearing within 30 deg with different limits, over 20 m or more; - LADOT records outside the city: a LADOT record running more than 20 m outside the City of Los Angeles boundary for 10 m or more; - camera segments: every camera in the app's cameras.json (read only) with a limit is compared with the dataset at points along its segment, snapped to the named street; any disagreement fails unless street_limits_camera_waivers.json documents it; - --release: also fails while the NOTICE has no publish URL (--publish-url) or a waiver or override is past its review date. OUTPUT FILE FORMAT (street_limits.bin, formatVersion 2, all integers little-endian) Byte layout is identical to the OSM prototype's version 1 (96-byte header, 12-byte records, UInt16 fixed-point points, uniform grid index, name table). Version 2 changes the meaning of the source nibble and of flag bit 1, so a version-1 reader must refuse it. Coordinates are quantized to Q units of 1e-5 degree (quantumE6 = 10 micro-degrees), measured from the bounding-box origin (south-west corner): lat = (originLatE6 + latQ * quantumE6) / 1e6 lon = (originLonE6 + lonQ * quantumE6) / 1e6 The bounding box 33.70..34.35 N, -118.70..-118.10 E is 65000 x 60000 Q units, so a point is two UInt16. Resolution ~1.1 m N-S and ~0.92 m E-W. Header, 96 bytes: off type field 0 char[4] magic "LASL" 4 u16 formatVersion (2) 6 u16 headerSize (96) 8 i32 originLatE6 (33700000) 12 i32 originLonE6 (-118700000) 16 u32 quantumE6 (10) 20 u16 cellSizeQ (512 Q units = 0.00512 deg, ~568 m N-S x ~472 m E-W) 22 u16 gridCols (118, cells along longitude) 24 u16 gridRows (127, cells along latitude) 26 u16 indexMarginM (30: every record is listed in each cell it passes within this many meters of, so a lookup with radius <= 30 m only reads the single cell that contains the fix) 28 u32 segmentCount (number of records; "segment" = one polyline) 32 u32 pointCount 36 u32 cellRefCount 40 u32 nameCount 44 u32 dataDate (YYYYMMDD, UTC date the sources were retrieved: "as of") 48 u32 buildDate (YYYYMMDD, UTC) 52 u32 segmentsOffset 56 u32 pointsOffset 60 u32 cellStartOffset 64 u32 cellRefsOffset 68 u32 nameOffsetsOffset 72 u32 nameBytesOffset 76 u32 nameBytesLength 80 u32 fileLength 84 u32 crc32 (zlib CRC-32 of bytes [headerSize, fileLength)) 88 u32 sourcesMask (bit n set when any record has source n) 92 u32 osmDataDate (YYYYMMDD of the OSM snapshot; 0 when built --no-osm) Sections (every offset is a multiple of 4), in this order: segments segmentCount x 12 bytes: +0 u32 firstPoint index into points +4 u16 pointCount >= 2 +6 u16 nameIndex 0xFFFF = no name +8 u8 speedFwdMph limit when travelling in point order (0 = none) +9 u8 speedBwdMph limit when travelling against point order +10 u8 classSource low nibble = road class, high nibble = source +11 u8 flags road class: 0 motorway, 1 motorway_link, 2 trunk, 3 trunk_link, 4 primary, 5 primary_link, 6 secondary, 7 secondary_link, 8 tertiary, 9 tertiary_link, 10 unclassified, 11 residential, 12 living_street. For OSM records it is the highway tag. For official records it is derived: LADOT ST_DESIG for the base map, 11 for CVC, from the limit for LAMC <= 35 and County lines, 0 for Caltrans freeways. Only "is a freeway" (0 or 1) matters to the matcher. source (high nibble): 0 none: no limit known. OSM road geometry kept as "known road" evidence 1 osm: OpenStreetMap maxspeed (ODbL) 2 ladot_cvc_25: LADOT SpeedLimits_CVC, 25 mph prima facie 3 la_county_dpw: LA County Public Works Speed Network 4 ladot_lamc_35: LADOT SpeedLimits_LAMC_35mph (Jan 2025) 5 ladot_speed_limit_map: LADOT City_of_Los_Angeles_Speed_Limit_Map 6 ladot_ab43: base map segment corrected by the AB 43 change layer 7 caltrans_d7: Caltrans District 7 SHS speed limit (freeway mainline) 8 glendale: City of Glendale GMC 10.24.020 (Ord. 6040, 2025-08-19), on LA County road centerline geometry (street_limits_overrides.json) 9 ladot_document: LADOT program documents (Speed Safety System camera maps, 2026), on LADOT geometry (street_limits_overrides.json) OSM-derived records are exactly those with source 0 or 1. Readers must accept sources 0..9 (a version-2 file never uses 10..15). flags: bit0 oneway (legal only in point order; OSM oneway=-1 ways are stored reversed; official records are two-way centerlines or carriageways), bit1 prima facie (set exactly when source == 2), bit2 bridge, bit3 tunnel, bit4 has maxspeed:conditional (OSM; the stored speed is the unconditional one), bits5-7 OSM layer + 3, clamped to 0..7 (3 = ground level for all official records) points pointCount x 4 bytes: u16 latQ, u16 lonQ. A record's points are points[firstPoint ..< firstPoint + pointCount]. cellStart (gridCols * gridRows + 1) x u32. Cell index = row * gridCols + col, row = latQ / cellSizeQ, col = lonQ / cellSizeQ (integer division). The cell's records are cellRefs[cellStart[i] ..< cellStart[i + 1]]. cellRefs cellRefCount x u32 record indices, ascending within a cell. nameOffsets (nameCount + 1) x u32 byte offsets into nameBytes. nameBytes UTF-8. Name i = nameBytes[nameOffsets[i] ..< nameOffsets[i + 1]]. "name" or "nameref" (ref like "I 405"); either part may be empty. Zero padding up to the next multiple of 4. Records are sorted in Morton (Z) order of their midpoint, so one cell's records sit close together and a memory-mapped file touches few pages per lookup. """ from __future__ import annotations import argparse import array import collections import datetime as dt import glob import hashlib import itertools import json import math import os import pickle import random import re import struct import sys import time import urllib.error import urllib.parse import urllib.request import zlib # --------------------------------------------------------------------------- constants SOUTH, WEST, NORTH, EAST = 33.70, -118.70, 34.35, -118.10 MAGIC = b"LASL" FORMAT_VERSION = 2 HEADER_SIZE = 96 ORIGIN_LAT_E6 = 33_700_000 ORIGIN_LON_E6 = -118_700_000 QUANTUM_E6 = 10 Q_PER_DEG = 1_000_000 // QUANTUM_E6 LAT_SPAN_Q = round((NORTH - SOUTH) * Q_PER_DEG) # 65000 LON_SPAN_Q = round((EAST - WEST) * Q_PER_DEG) # 60000 CELL_Q = 512 SEGMENT_STRUCT = struct.Struct(" None: print(f"[{time.strftime('%H:%M:%S')}] {message}", file=sys.stderr, flush=True) def km(meters: float) -> float: return round(meters / 1000.0, 1) # --------------------------------------------------------------------------- HTTP def http_get(url: str, timeout: int = 120, retries: int = 6, first_delay: float = 4.0) -> bytes: """GET with exponential backoff on 429/5xx and network errors (honours Retry-After).""" delay = first_delay for attempt in range(1, retries + 1): request = urllib.request.Request(url, headers={"User-Agent": USER_AGENT}) try: with urllib.request.urlopen(request, timeout=timeout) as response: return response.read() except urllib.error.HTTPError as error: if error.code not in (429, 500, 502, 503, 504) or attempt == retries: raise retry_after = error.headers.get("Retry-After") if error.headers else None if retry_after and retry_after.isdigit(): delay = max(delay, float(retry_after)) reason = f"HTTP {error.code}" except (urllib.error.URLError, TimeoutError, ConnectionError, OSError) as error: if attempt == retries: raise reason = str(error) log(f" {reason}; retry {attempt}/{retries - 1} in {delay:.0f}s") time.sleep(delay) delay = min(delay * 2, 120) raise RuntimeError("unreachable") def fetch_arcgis(source_id: str, cache_root: str, refresh: bool, pause: float, page_size: int = 2000) -> dict: """Page through an ArcGIS FeatureServer query (f=geojson, outSR 4326). Every page is cached raw under cache/arcgis//, so an interrupted run resumes; meta.json records the server count, retrieval time and the layer's last data edit.""" spec = ARCGIS_SOURCES[source_id] folder = os.path.join(cache_root, source_id) meta_path = os.path.join(folder, "meta.json") if os.path.exists(meta_path) and not refresh: with open(meta_path, encoding="utf-8") as f: meta = json.load(f) if meta.get("complete"): log(f"cached {source_id}: {meta['featureCount']} features in {len(meta['pages'])} pages") return meta if refresh and os.path.isdir(folder): for old in glob.glob(os.path.join(folder, "page_*.geojson")): os.remove(old) os.makedirs(folder, exist_ok=True) base = {"where": spec["where"], "outFields": spec["outFields"], "outSR": 4326, "returnGeometry": "true"} if spec["bbox"]: w, s_, e, n = spec["bbox"] if isinstance(spec["bbox"], tuple) else (WEST, SOUTH, EAST, NORTH) base.update(geometry=f"{w},{s_},{e},{n}", geometryType="esriGeometryEnvelope", inSR=4326, spatialRel="esriSpatialRelIntersects") layer = json.loads(http_get(spec["url"] + "?f=json")) with open(os.path.join(folder, "layer.json"), "w", encoding="utf-8") as f: json.dump(layer, f) time.sleep(pause) count_q = {k: v for k, v in base.items() if k not in ("outFields", "outSR", "returnGeometry")} count_q.update(returnCountOnly="true", f="json") expected = json.loads(http_get(spec["url"] + "/query?" + urllib.parse.urlencode(count_q)))["count"] time.sleep(pause) pages, total, offset, page_no = [], 0, 0, 0 while True: name = f"page_{page_no:04d}.geojson" path = os.path.join(folder, name) if os.path.exists(path): with open(path, "rb") as f: payload = f.read() else: query = dict(base, orderByFields=spec["orderBy"], resultOffset=offset, resultRecordCount=page_size, geometryPrecision=7, f="geojson") log(f"{source_id}: page {page_no} (offset {offset} of {expected})") payload = http_get(spec["url"] + "/query?" + urllib.parse.urlencode(query)) doc = json.loads(payload) if "error" in doc: raise RuntimeError(f"{source_id}: ArcGIS error {doc['error']}") with open(path + ".part", "wb") as f: f.write(payload) os.replace(path + ".part", path) time.sleep(pause) doc = json.loads(payload) n = len(doc.get("features") or []) pages.append(name) total += n exceeded = (doc.get("properties") or {}).get("exceededTransferLimit") or doc.get("exceededTransferLimit") if n == 0 or total >= expected or (n < page_size and not exceeded): break offset += n page_no += 1 if total != expected: raise RuntimeError(f"{source_id}: got {total} features, the server counts {expected}") info = layer.get("editingInfo") or {} edit = info.get("dataLastEditDate") or info.get("lastEditDate") meta = { "source": source_id, "layerUrl": spec["url"], "itemUrl": spec.get("item"), "layerName": layer.get("name"), "where": spec["where"], "outFields": spec["outFields"], "bboxFilter": list(spec["bbox"]) if isinstance(spec["bbox"], tuple) else spec["bbox"], "featureCount": total, "pages": pages, "retrievedUtc": dt.datetime.now(dt.timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"), "dataLastEditUtc": (dt.datetime.fromtimestamp(edit / 1000, dt.timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") if edit else None), "complete": True, } with open(meta_path, "w", encoding="utf-8") as f: json.dump(meta, f, indent=1) log(f"{source_id}: {total} features") return meta def load_features(cache_root: str, source_id: str) -> list: folder = os.path.join(cache_root, source_id) with open(os.path.join(folder, "meta.json"), encoding="utf-8") as f: meta = json.load(f) features = [] for name in meta["pages"]: with open(os.path.join(folder, name), encoding="utf-8") as f: features.extend(json.load(f).get("features") or []) return features def fetch_md5(url: str): candidates = [url + ".md5"] if "-latest." in url: candidates.append(url.replace("-latest.", ".") + ".md5") for candidate in candidates: try: text = http_get(candidate, timeout=30, retries=2).decode("ascii", "replace") except Exception: # noqa: BLE001 - md5 files are optional continue match = re.match(r"\s*([0-9a-f]{32})\b", text) if match: return match.group(1) return None def download_file(url: str, dest: str, refresh: bool) -> str: if os.path.exists(dest) and not refresh: log(f"cached {os.path.basename(dest)} ({os.path.getsize(dest) / 1e6:.1f} MB)") return dest os.makedirs(os.path.dirname(dest), exist_ok=True) for attempt in (1, 2): log(f"downloading {url}") tmp = dest + ".part" digest = hashlib.md5() received, next_report = 0, 50_000_000 request = urllib.request.Request(url, headers={"User-Agent": USER_AGENT}) with urllib.request.urlopen(request, timeout=120) as response, open(tmp, "wb") as out: total = int(response.headers.get("Content-Length") or 0) while True: chunk = response.read(1 << 20) if not chunk: break out.write(chunk) digest.update(chunk) received += len(chunk) if received >= next_report: log(f" {received / 1e6:.0f} / {total / 1e6:.0f} MB") next_report += 50_000_000 expected = fetch_md5(url) if expected is None or expected == digest.hexdigest(): os.replace(tmp, dest) log(f" done, {received / 1e6:.1f} MB, md5 {'ok' if expected else 'not published'}") return dest log(f" md5 mismatch (got {digest.hexdigest()}, expected {expected})") os.remove(tmp) raise RuntimeError(f"md5 mismatch for {url}") # --------------------------------------------------------------------------- OSM PBF reader # Minimal decoder for the OSM PBF format (fileformat.proto / osmformat.proto): enough to # read ways with their tags and node references, then the coordinates of those nodes. _VARINT_RE = re.compile(rb"[\x80-\xff]*[\x00-\x7f]") def _uvarint(buf: bytes, pos: int): byte = buf[pos] if byte < 0x80: return byte, pos + 1 result, shift = byte & 0x7F, 7 pos += 1 while True: byte = buf[pos] pos += 1 result |= (byte & 0x7F) << shift if byte < 0x80: return result, pos shift += 7 def _fields(buf: bytes, pos: int, end: int): """Yield (field, wire_type, value): an int for varints, (start, end) for bytes.""" while pos < end: key, pos = _uvarint(buf, pos) wire = key & 7 if wire == 0: value, pos = _uvarint(buf, pos) yield key >> 3, 0, value elif wire == 2: length, pos = _uvarint(buf, pos) yield key >> 3, 2, (pos, pos + length) pos += length elif wire == 1: yield key >> 3, 1, buf[pos:pos + 8] pos += 8 elif wire == 5: yield key >> 3, 5, buf[pos:pos + 4] pos += 4 else: raise ValueError(f"unsupported protobuf wire type {wire}") def _packed_uvarints(buf: bytes, start: int, end: int) -> list: out = [] append = out.append for m in _VARINT_RE.findall(buf, start, end): n = len(m) if n == 1: append(m[0]) elif n == 2: append((m[0] & 127) | (m[1] << 7)) elif n == 3: append((m[0] & 127) | ((m[1] & 127) << 7) | (m[2] << 14)) else: value = shift = 0 for byte in m: value |= (byte & 127) << shift shift += 7 append(value) return out def _packed_sint_deltas(buf: bytes, start: int, end: int) -> list: """Packed zigzag sint64 deltas -> absolute values.""" out = [] append = out.append for m in _VARINT_RE.findall(buf, start, end): n = len(m) if n == 1: v = m[0] elif n == 2: v = (m[0] & 127) | (m[1] << 7) elif n == 3: v = (m[0] & 127) | ((m[1] & 127) << 7) | (m[2] << 14) elif n == 4: v = (m[0] & 127) | ((m[1] & 127) << 7) | ((m[2] & 127) << 14) | (m[3] << 21) else: v = shift = 0 for byte in m: v |= (byte & 127) << shift shift += 7 append((v >> 1) ^ -(v & 1)) return list(itertools.accumulate(out)) def _signed64(value: int) -> int: return value - (1 << 64) if value >= 1 << 63 else value def _iter_blobs(path: str): with open(path, "rb") as f: while True: head = f.read(4) if len(head) < 4: return (header_len,) = struct.unpack(">I", head) header = f.read(header_len) blob_type, size = "", 0 for field, _, value in _fields(header, 0, len(header)): if field == 1: blob_type = header[value[0]:value[1]].decode() elif field == 3: size = value offset = f.tell() f.seek(size, os.SEEK_CUR) yield blob_type, offset, size def _blob_payload(blob: bytes) -> bytes: for field, _, value in _fields(blob, 0, len(blob)): if field == 1: return blob[value[0]:value[1]] if field == 3: return zlib.decompress(blob[value[0]:value[1]]) if field == 4: import lzma return lzma.decompress(blob[value[0]:value[1]]) if field in (5, 6, 7): raise ValueError("PBF blob compression (bzip2/lz4/zstd) not supported; use a zlib PBF") raise ValueError("empty PBF blob") def _block_parts(data: bytes): strings: list = [] groups: list = [] granularity, lat_offset, lon_offset = 100, 0, 0 for field, _, value in _fields(data, 0, len(data)): if field == 1: for f2, _, v2 in _fields(data, value[0], value[1]): if f2 == 1: strings.append(data[v2[0]:v2[1]]) elif field == 2: groups.append(value) elif field == 17: granularity = value elif field == 19: lat_offset = _signed64(value) elif field == 20: lon_offset = _signed64(value) return strings, groups, granularity, lat_offset, lon_offset def _group_kind(data: bytes, group) -> int: start, end = group if start >= end: return 0 key, _ = _uvarint(data, start) return key >> 3 # 1 nodes, 2 dense nodes, 3 ways, 4 relations WANTED_TAGS = ( "highway", "maxspeed", "maxspeed:forward", "maxspeed:backward", "maxspeed:conditional", "name", "ref", "oneway", "junction", "bridge", "tunnel", "layer", "access", "motor_vehicle", "motorcar", "area", ) _WANTED_BYTES = {t.encode(): t for t in WANTED_TAGS} _DRIVABLE_BYTES = {c.encode() for c in ROAD_CLASSES} def read_pbf_ways(paths: list) -> tuple: """Return (ways, osm_timestamp) where ways = [(way_id, tags, [(lat, lon), ...]), ...].""" started = time.time() timestamp = None way_tags: dict = {} way_refs: dict = {} node_blobs: dict = {path: [] for path in paths} for path in paths: log(f"PBF pass 1 (ways): {os.path.basename(path)}") with open(path, "rb") as f: for blob_type, offset, size in _iter_blobs(path): f.seek(offset) data = _blob_payload(f.read(size)) if blob_type == "OSMHeader": for field, _, value in _fields(data, 0, len(data)): if field == 32: stamp = dt.datetime.fromtimestamp(_signed64(value), dt.timezone.utc) timestamp = stamp if timestamp is None else min(timestamp, stamp) continue if blob_type != "OSMData": continue strings, groups, *_ = _block_parts(data) kinds = {_group_kind(data, g) for g in groups} if kinds & {1, 2}: node_blobs[path].append((offset, size)) if 3 not in kinds: continue try: hw_key = strings.index(b"highway") except ValueError: continue wanted = {i: _WANTED_BYTES[s] for i, s in enumerate(strings) if s in _WANTED_BYTES} drivable = {i for i, s in enumerate(strings) if s in _DRIVABLE_BYTES} if not drivable: continue for group in groups: if _group_kind(data, group) != 3: continue for field, _, value in _fields(data, group[0], group[1]): if field != 3: continue way_id = 0 keys_range = vals_range = refs_range = None for f2, _, v2 in _fields(data, value[0], value[1]): if f2 == 1: way_id = v2 elif f2 == 2: keys_range = v2 elif f2 == 3: vals_range = v2 elif f2 == 8: refs_range = v2 if keys_range is None or refs_range is None: continue keys = _packed_uvarints(data, *keys_range) if hw_key not in keys: continue vals = _packed_uvarints(data, *vals_range) if vals[keys.index(hw_key)] not in drivable or way_id in way_tags: continue way_tags[way_id] = {wanted[k]: strings[v].decode("utf-8", "replace") for k, v in zip(keys, vals) if k in wanted} way_refs[way_id] = array.array("q", _packed_sint_deltas(data, *refs_range)) log(f" {len(way_tags)} drivable ways so far ({time.time() - started:.0f}s)") needed = set() for refs in way_refs.values(): needed.update(refs) log(f"PBF pass 2 (nodes): {len(needed)} node coordinates needed") coords: dict = {} for path in paths: with open(path, "rb") as f: for offset, size in node_blobs[path]: f.seek(offset) data = _blob_payload(f.read(size)) _, groups, granularity, lat_offset, lon_offset = _block_parts(data) for group in groups: kind = _group_kind(data, group) if kind == 2: for field, _, value in _fields(data, group[0], group[1]): if field != 2: continue id_r = lat_r = lon_r = None for f2, _, v2 in _fields(data, value[0], value[1]): if f2 == 1: id_r = v2 elif f2 == 8: lat_r = v2 elif f2 == 9: lon_r = v2 if id_r is None: continue ids = _packed_sint_deltas(data, *id_r) hits = needed.intersection(ids) if not hits: continue lats = _packed_sint_deltas(data, *lat_r) lons = _packed_sint_deltas(data, *lon_r) for k in itertools.compress(range(len(ids)), map(hits.__contains__, ids)): lat_e7 = (lat_offset + granularity * lats[k] + 50) // 100 lon_e7 = (lon_offset + granularity * lons[k] + 50) // 100 coords[ids[k]] = ((lat_e7 + 900_000_000) << 32) | (lon_e7 + 1_800_000_000) elif kind == 1: for field, _, value in _fields(data, group[0], group[1]): if field != 1: continue node_id = lat = lon = 0 for f2, _, v2 in _fields(data, value[0], value[1]): if f2 == 1: node_id = (v2 >> 1) ^ -(v2 & 1) elif f2 == 8: lat = (v2 >> 1) ^ -(v2 & 1) elif f2 == 9: lon = (v2 >> 1) ^ -(v2 & 1) if node_id in needed: lat_e7 = (lat_offset + granularity * lat + 50) // 100 lon_e7 = (lon_offset + granularity * lon + 50) // 100 coords[node_id] = ((lat_e7 + 900_000_000) << 32) | (lon_e7 + 1_800_000_000) log(f" {os.path.basename(path)}: {len(coords)} coordinates found ({time.time() - started:.0f}s)") del needed ways = [] for way_id, refs in way_refs.items(): run: list = [] for ref in refs: packed = coords.get(ref) if packed is None: if len(run) >= 2: ways.append((way_id, way_tags[way_id], run)) run = [] continue run.append((((packed >> 32) - 900_000_000) / 1e7, ((packed & 0xFFFFFFFF) - 1_800_000_000) / 1e7)) if len(run) >= 2: ways.append((way_id, way_tags[way_id], run)) log(f"PBF read: {len(ways)} way runs, {time.time() - started:.0f}s") return ways, timestamp # --------------------------------------------------------------------------- normalization _MPH_RE = re.compile(r"^\s*(\d{1,3})\s*mph\s*$", re.IGNORECASE) _KMH_RE = re.compile(r"^\s*(\d{1,3})\s*(km/h|kmh|kph)\s*$", re.IGNORECASE) _TRUE = {"yes", "true", "1"} def valid_mph(mph: int) -> bool: """California posts limits in 5 mph steps; 5..75 is the plausible range.""" return mph % 5 == 0 and 5 <= mph <= 75 def parse_mph(value, issues: collections.Counter) -> int: """OSM maxspeed: '35 mph' -> 35. 0 when absent or unusable (counted in `issues`).""" if value is None: return 0 match = _MPH_RE.match(value) if match: mph = int(match.group(1)) elif value.strip().isdigit(): issues["osm_unitless_dropped"] += 1 # km/h by OSM rules, mph by intent: unsafe return 0 else: match = _KMH_RE.match(value) if not match: issues["osm_unparsed:" + value.strip()[:40]] += 1 return 0 mph = int(5 * round(int(match.group(1)) / 1.609344 / 5)) issues["osm_kmh_converted"] += 1 if not valid_mph(mph): issues["osm_implausible_dropped"] += 1 return 0 return mph def parse_speed_field(value) -> int: """Official layers: 35, '35', '35 mph' -> 35; '25 and 35', 0, None -> 0.""" if value is None: return 0 if isinstance(value, (int, float)): mph = int(value) else: match = re.fullmatch(r"\s*(\d{1,2})\s*(mph)?\s*", str(value), re.IGNORECASE) if not match: return 0 mph = int(match.group(1)) return mph if valid_mph(mph) else 0 def normalize_osm_way(tags: dict, issues: collections.Counter): """Return (key, reverse geometry?) or None for non-drivable ways. key = (class_code, speed_fwd, speed_bwd, source, flags, label). No defaults are invented: an untagged way gets source none and speed 0.""" highway = tags.get("highway") if highway not in CLASS_CODE or tags.get("area") == "yes": return None motor = tags.get("motor_vehicle") or tags.get("motorcar") if motor == "no" or (tags.get("access") == "no" and motor not in ("yes", "designated", "permissive", "destination")): issues["osm_dropped_no_motor_access"] += 1 return None general = parse_mph(tags.get("maxspeed"), issues) fwd = parse_mph(tags.get("maxspeed:forward"), issues) or general bwd = parse_mph(tags.get("maxspeed:backward"), issues) or general oneway_tag = (tags.get("oneway") or "").lower() reverse = oneway_tag in ("-1", "reverse") if oneway_tag in _TRUE or reverse: oneway = True elif oneway_tag in ("no", "false", "0"): oneway = False else: oneway = highway in ("motorway", "motorway_link") or tags.get("junction") in ("roundabout", "circular") if reverse: fwd, bwd = bwd, fwd flags = FLAG_ONEWAY if oneway else 0 if tags.get("bridge") not in (None, "no"): flags |= FLAG_BRIDGE if tags.get("tunnel") not in (None, "no"): flags |= FLAG_TUNNEL if tags.get("maxspeed:conditional"): flags |= FLAG_CONDITIONAL try: layer = int(float(tags.get("layer", "0"))) except ValueError: layer = 0 flags |= (max(-3, min(4, layer)) + 3) << LAYER_SHIFT if fwd or bwd: source = SRC_OSM else: source = SRC_NONE fwd = bwd = 0 name = " ".join((tags.get("name") or "").split()) ref = " ".join((tags.get("ref") or "").split()) label = f"{name}\t{ref}" if ref else name return (CLASS_CODE[highway], fwd, bwd, source, flags, label), reverse def reversed_key(key: tuple) -> tuple: cls, fwd, bwd, source, flags, label = key return (cls, bwd, fwd, source, flags, label) _SUFFIXES = { "AV": "Avenue", "AVE": "Avenue", "AVENUE": "Avenue", "BL": "Boulevard", "BLVD": "Boulevard", "ST": "Street", "DR": "Drive", "RD": "Road", "ROAD": "Road", "PL": "Place", "WY": "Way", "WAY": "Way", "CT": "Court", "LN": "Lane", "LANE": "Lane", "TR": "Terrace", "TER": "Terrace", "TERR": "Terrace", "CIR": "Circle", "CI": "Circle", "PKWY": "Parkway", "PY": "Parkway", "HWY": "Highway", "HY": "Highway", "SQ": "Square", "TRL": "Trail", "TL": "Trail", "LP": "Loop", "CRES": "Crescent", "PZ": "Plaza", "PLZ": "Plaza", "CYN": "Canyon", "FWY": "Freeway", "FY": "Freeway", "EXPY": "Expressway", "RDWY": "Roadway", "MALL": "Mall", "VIS": "Vista", "PT": "Point", "WK": "Walk", "AL": "Alley", "ALY": "Alley", "GRN": "Green", "HTS": "Heights", } _DIRECTIONS = {"N": "North", "S": "South", "E": "East", "W": "West", "NORTH": "North", "SOUTH": "South", "EAST": "East", "WEST": "West"} _ORDINAL_RE = re.compile(r"^(\d+)(ST|ND|RD|TH)$", re.IGNORECASE) def _title_token(token: str) -> str: match = _ORDINAL_RE.match(token) if match: return match.group(1) + match.group(2).lower() if token.isdigit() or "/" in token: return token if token.upper().startswith("MC") and len(token) > 2: return "Mc" + token[2:].capitalize() return "-".join(part.capitalize() for part in token.split("-")) def street_name(prefix_dir, name, suffix=None, suffix_dir=None) -> str: """BOE / County style 'W', 'SUNSET', 'BLVD' -> 'West Sunset Boulevard'. When `suffix` is None the street type is taken from the last word of `name`.""" tokens = " ".join(str(name or "").replace("\t", " ").split()).split(" ") tokens = [t for t in tokens if t] if not tokens: return "" if suffix is None and len(tokens) > 1 and tokens[-1].upper() in _SUFFIXES: suffix = tokens.pop() words = [_title_token(t) for t in tokens] suffix = (suffix or "").strip().upper() if suffix: words.append(_SUFFIXES.get(suffix, _title_token(suffix))) prefix = _DIRECTIONS.get((prefix_dir or "").strip().upper()) if prefix and words and words[0] not in ("North", "South", "East", "West"): words.insert(0, prefix) tail = _DIRECTIONS.get((suffix_dir or "").strip().upper()) if tail: words.append(tail) return " ".join(words) def caltrans_ref(route: str) -> str: route = str(route).strip() if route in INTERSTATES: return f"I {route}" if route in US_ROUTES: return f"US {route}" return f"CA {route}" # --------------------------------------------------------------------------- geometry def to_xy(lat: float, lon: float): return (lon - WEST) * M_PER_DEG_LON, (lat - SOUTH) * M_PER_DEG_LAT def distance_m(a, b) -> float: return math.hypot((b[1] - a[1]) * M_PER_DEG_LON, (b[0] - a[0]) * M_PER_DEG_LAT) def polyline_length_m(coords) -> float: return sum(distance_m(coords[i], coords[i + 1]) for i in range(len(coords) - 1)) def _clip_segment(x0, y0, x1, y1, xmin, ymin, xmax, ymax): """Liang-Barsky. Returns (t0, t1) of the part inside the rectangle, or None.""" dx, dy = x1 - x0, y1 - y0 t0, t1 = 0.0, 1.0 for p, q in ((-dx, x0 - xmin), (dx, xmax - x0), (-dy, y0 - ymin), (dy, ymax - y0)): if p == 0: if q < 0: return None continue t = q / p if p < 0: if t > t1: return None if t > t0: t0 = t else: if t < t0: return None if t < t1: t1 = t return t0, t1 def clip_to_bbox(coords) -> list: """Split a polyline into the runs that lie inside the service bounding box.""" if all(SOUTH <= lat <= NORTH and WEST <= lon <= EAST for lat, lon in coords): return [coords] runs, run = [], [] for (lat0, lon0), (lat1, lon1) in zip(coords, coords[1:]): hit = _clip_segment(lon0, lat0, lon1, lat1, WEST, SOUTH, EAST, NORTH) if hit is None: if len(run) >= 2: runs.append(run) run = [] continue t0, t1 = hit a = (lat0 + (lat1 - lat0) * t0, lon0 + (lon1 - lon0) * t0) b = (lat0 + (lat1 - lat0) * t1, lon0 + (lon1 - lon0) * t1) if t0 > 0 or not run: if len(run) >= 2: runs.append(run) run = [a] run.append(b) if t1 < 1: if len(run) >= 2: runs.append(run) run = [] if len(run) >= 2: runs.append(run) return runs def simplify(coords: list, tolerance_m: float) -> list: """Douglas-Peucker on a local metric projection (distance to segment, so rings work).""" n = len(coords) if n <= 2: return coords xs = [(lon - WEST) * M_PER_DEG_LON for _, lon in coords] ys = [(lat - SOUTH) * M_PER_DEG_LAT for lat, _ in coords] keep = bytearray(n) keep[0] = keep[n - 1] = 1 tol2 = tolerance_m * tolerance_m stack = [(0, n - 1)] while stack: i, j = stack.pop() if j <= i + 1: continue ax, ay, bx, by = xs[i], ys[i], xs[j], ys[j] dx, dy = bx - ax, by - ay length2 = dx * dx + dy * dy best, best_k = -1.0, -1 for k in range(i + 1, j): px, py = xs[k] - ax, ys[k] - ay if length2 == 0.0: d2 = px * px + py * py else: t = (px * dx + py * dy) / length2 if t <= 0.0: d2 = px * px + py * py elif t >= 1.0: qx, qy = xs[k] - bx, ys[k] - by d2 = qx * qx + qy * qy else: cross = px * dy - py * dx d2 = cross * cross / length2 if d2 > best: best, best_k = d2, k if best > tol2: keep[best_k] = 1 stack.append((i, best_k)) stack.append((best_k, j)) return [coords[k] for k in range(n) if keep[k]] def split_by_length(coords: list, max_m: float) -> list: chunks, current, length = [], [coords[0]], 0.0 for k in range(1, len(coords)): step = distance_m(coords[k - 1], coords[k]) if length + step > max_m and len(current) >= 2: chunks.append(current) current, length = [coords[k - 1]], 0.0 current.append(coords[k]) length += step chunks.append(current) return chunks def merge_chains(pieces: list) -> list: """Join pieces with identical attributes end to end where exactly one continuation exists (two-way pieces may be joined reversed).""" def node(p): return (round(p[0] * 1e7), round(p[1] * 1e7)) ends = collections.defaultdict(list) for i, (key, coords) in enumerate(pieces): ends[node(coords[0])].append((i, True)) ends[node(coords[-1])].append((i, False)) used = bytearray(len(pieces)) merged = [] def continuation(at, key, oneway, seed, want_start): found = None for j, is_start in ends[at]: if used[j] or j == seed: continue other_key, _ = pieces[j] if is_start == want_start and other_key == key: candidate = (j, False) elif (not oneway and is_start != want_start and not (other_key[4] & FLAG_ONEWAY) and reversed_key(other_key) == key): candidate = (j, True) else: continue if found is not None: return None found = candidate return found for seed in range(len(pieces)): if used[seed]: continue used[seed] = 1 key, coords = pieces[seed] line = list(coords) oneway = bool(key[4] & FLAG_ONEWAY) while True: nxt = continuation(node(line[-1]), key, oneway, seed, True) if nxt is None: break j, flip = nxt used[j] = 1 extra = pieces[j][1][::-1] if flip else pieces[j][1] line.extend(extra[1:]) while True: prv = continuation(node(line[0]), key, oneway, seed, False) if prv is None: break j, flip = prv used[j] = 1 extra = pieces[j][1][::-1] if flip else pieces[j][1] line[:0] = extra[:-1] merged.append((key, line)) return merged def dedupe_pieces(pieces: list) -> tuple: """Drop pieces that repeat another piece exactly, or exactly reversed with the reversed attributes (two-way only: merge_chains would join such a pair into a line that goes out and back). Returns (pieces, number dropped).""" seen = set() out = [] dropped = 0 for key, coords in pieces: sig = tuple((round(lat * 1e6), round(lon * 1e6)) for lat, lon in coords) if (sig, key) in seen or (not key[4] & FLAG_ONEWAY and (sig[::-1], reversed_key(key)) in seen): dropped += 1 continue seen.add((sig, key)) out.append((key, coords)) return out, dropped def reversal_vertices(xy: list, min_angle: float = 165.0) -> list: """Indices of interior vertices where a polyline turns back on itself.""" limit = math.cos(math.radians(min_angle)) out = [] for k in range(1, len(xy) - 1): ax, ay = xy[k][0] - xy[k - 1][0], xy[k][1] - xy[k - 1][1] bx, by = xy[k + 1][0] - xy[k][0], xy[k + 1][1] - xy[k][1] la, lb = math.hypot(ax, ay), math.hypot(bx, by) if la > 0 and lb > 0 and (ax * bx + ay * by) / (la * lb) < limit: out.append(k) return out def iter_samples(coords: list, step: float): """(x, y, bearing mod 180) every `step` m along a polyline (at the step centres).""" xy = [to_xy(lat, lon) for lat, lon in coords] cum = [0.0] for a, b in zip(xy, xy[1:]): cum.append(cum[-1] + math.hypot(b[0] - a[0], b[1] - a[1])) total = cum[-1] if total <= 0.0: return n = max(1, int(math.ceil(total / step))) spacing = total / n seg_i = 0 for k in range(n): s = (k + 0.5) * spacing while seg_i < len(xy) - 2 and cum[seg_i + 1] < s: seg_i += 1 (ax, ay), (bx, by) = xy[seg_i], xy[seg_i + 1] length = cum[seg_i + 1] - cum[seg_i] t = (s - cum[seg_i]) / length if length > 0 else 0.0 yield ax + t * (bx - ax), ay + t * (by - ay), math.degrees(math.atan2(bx - ax, by - ay)) % 180.0, spacing def feature_lines(feature: dict) -> list: """GeoJSON LineString / MultiLineString -> [[(lat, lon), ...], ...], with the parts of a MultiLineString chained where one part starts at the previous part's end.""" geometry = feature.get("geometry") or {} if geometry.get("type") == "LineString": parts = [geometry.get("coordinates") or []] elif geometry.get("type") == "MultiLineString": parts = geometry.get("coordinates") or [] else: return [] lines: list = [] for part in parts: pts = [(pt[1], pt[0]) for pt in part if pt and len(pt) >= 2] dedup = [p for i, p in enumerate(pts) if i == 0 or p != pts[i - 1]] if len(dedup) < 2: continue if lines and distance_m(lines[-1][-1], dedup[0]) < 0.05: lines[-1].extend(dedup[1:]) else: lines.append(dedup) return lines def midpoint_and_bearing(coords: list): """Point halfway along the polyline (by length) and the bearing there (mod 180).""" xy = [to_xy(lat, lon) for lat, lon in coords] lengths = [math.hypot(b[0] - a[0], b[1] - a[1]) for a, b in zip(xy, xy[1:])] half = sum(lengths) / 2.0 for (a, b), step in zip(zip(xy, xy[1:]), lengths): if step > 0 and half <= step: t = half / step return (a[0] + t * (b[0] - a[0]), a[1] + t * (b[1] - a[1]), math.degrees(math.atan2(b[0] - a[0], b[1] - a[1])) % 180.0) half -= step a, b = xy[0], xy[-1] return a[0], a[1], math.degrees(math.atan2(b[0] - a[0], b[1] - a[1])) % 180.0 class LineIndex: """Uniform grid of line segments in the local metric projection. Each segment is inserted into every cell within `pad` meters, so a query within `pad` reads one cell.""" def __init__(self, cell: float = 50.0, pad: float = 30.0): self.cell = cell self.pad = pad self.grid: dict = {} self.segs: list = [] def add(self, coords: list, payload, family: int) -> None: cell, pad, grid, segs = self.cell, self.pad, self.grid, self.segs xy = [to_xy(lat, lon) for lat, lon in coords] for (x0, y0), (x1, y1) in zip(xy, xy[1:]): dx, dy = x1 - x0, y1 - y0 l2 = dx * dx + dy * dy if l2 == 0.0: continue index = len(segs) segs.append((x0, y0, dx, dy, l2, math.degrees(math.atan2(dx, dy)) % 180.0, payload, family, math.sqrt(l2))) reach = pad + cell * 0.7072 for gx in range(int((min(x0, x1) - pad) // cell), int((max(x0, x1) + pad) // cell) + 1): cx = (gx + 0.5) * cell for gy in range(int((min(y0, y1) - pad) // cell), int((max(y0, y1) + pad) // cell) + 1): cy = (gy + 0.5) * cell t = max(0.0, min(1.0, ((cx - x0) * dx + (cy - y0) * dy) / l2)) if math.hypot(x0 + t * dx - cx, y0 + t * dy - cy) > reach: continue bucket = grid.get((gx, gy)) if bucket is None: grid[(gx, gy)] = [index] else: bucket.append(index) def nearest(self, x: float, y: float, bearing: float, max_d: float, max_angle: float, families: int = FAM_SURFACE | FAM_FREEWAY, end_tol: float = 3.0, accept=None): """(segment tuple, distance) of the nearest aligned segment within max_d whose perpendicular foot lies on the segment or at most `end_tol` m past its ends (so a line that merely ends near the point, like the previous block of the same street, does not count as running alongside it), or None. `accept(payload)` can veto candidates (for example: same street name only).""" assert max_d <= self.pad, (max_d, self.pad) bucket = self.grid.get((int(x // self.cell), int(y // self.cell))) if not bucket: return None best, hit = max_d * max_d, None segs = self.segs for i in bucket: seg = segs[i] if not seg[7] & families: continue diff = abs(seg[5] - bearing) % 180.0 if diff > 90.0: diff = 180.0 - diff if diff > max_angle: continue x0, y0, dx, dy, l2 = seg[0], seg[1], seg[2], seg[3], seg[4] t = ((x - x0) * dx + (y - y0) * dy) / l2 if t < 0.0: if -t * seg[8] > end_tol: continue t = 0.0 elif t > 1.0: if (t - 1.0) * seg[8] > end_tol: continue t = 1.0 ex, ey = x0 + t * dx - x, y0 + t * dy - y d2 = ex * ex + ey * ey if d2 <= best and (accept is None or accept(seg[6])): best, hit = d2, seg return None if hit is None else (hit, math.sqrt(best)) def carve(coords: list, index: LineIndex, max_d: float, families: int, tally=None, step: float = CARVE_STEP_M, max_angle: float = CARVE_ANGLE, accept=None, point_filter=None) -> list: """Parts of `coords` that no indexed line covers (see the module docstring). `tally(payload, meters, distance, x, y)` is called for every covered sample. `accept` filters candidate payloads; `point_filter(x, y)` limits carving to some samples.""" xy = [to_xy(lat, lon) for lat, lon in coords] cum = [0.0] for a, b in zip(xy, xy[1:]): cum.append(cum[-1] + math.hypot(b[0] - a[0], b[1] - a[1])) total = cum[-1] if total <= 0.0: return [] n = max(1, int(math.ceil(total / step))) spacing = total / n covered = bytearray(n) seg_i = 0 any_cover = False for k in range(n): s = (k + 0.5) * spacing while seg_i < len(xy) - 2 and cum[seg_i + 1] < s: seg_i += 1 (ax, ay), (bx, by) = xy[seg_i], xy[seg_i + 1] length = cum[seg_i + 1] - cum[seg_i] t = (s - cum[seg_i]) / length if length > 0 else 0.0 x, y = ax + t * (bx - ax), ay + t * (by - ay) bearing = math.degrees(math.atan2(bx - ax, by - ay)) % 180.0 if point_filter is not None and not point_filter(x, y): continue hit = index.nearest(x, y, bearing, max_d, max_angle, families, accept=accept) if hit is not None: covered[k] = 1 any_cover = True if tally is not None: tally(hit[0][6], spacing, hit[1], x, y) if not any_cover: return [coords] # Ignore short covered blips inside uncovered stretches, then collect uncovered runs. runs = [] k = 0 while k < n: j = k while j < n and covered[j] == covered[k]: j += 1 runs.append([covered[k], k, j]) k = j if len(runs) > 1: # short covered blips (also at either end: junction overlaps) are ignored for run in runs: if run[0] and (run[2] - run[1]) * spacing < MAX_GAP_FILL_M: run[0] = 0 merged_runs = [] for run in runs: if merged_runs and merged_runs[-1][0] == run[0]: merged_runs[-1][2] = run[2] else: merged_runs.append(run) out = [] for flag, start, end in merged_runs: if flag: continue a, b = start * spacing, end * spacing if b - a < MIN_KEEP_M: continue out.append(sub_polyline(coords, cum, a, b)) return [part for part in out if len(part) >= 2] def sub_polyline(coords: list, cum: list, a: float, b: float) -> list: """The part of the polyline between arc lengths a and b (meters).""" def point_at(s): for i in range(len(cum) - 1): if cum[i + 1] >= s: length = cum[i + 1] - cum[i] t = (s - cum[i]) / length if length > 0 else 0.0 (la0, lo0), (la1, lo1) = coords[i], coords[i + 1] return (la0 + t * (la1 - la0), lo0 + t * (lo1 - lo0)), i return coords[-1], len(coords) - 2 start, i0 = point_at(a) end, i1 = point_at(b) pts = [start] + [coords[i] for i in range(i0 + 1, i1 + 1)] + [end] return [p for i, p in enumerate(pts) if i == 0 or distance_m(p, pts[i - 1]) > 0.01] def label_runs(coords: list, label_at, step: float = CARVE_STEP_M, min_run_m: float = MIN_RUN_M) -> list: """Split a polyline where a per-sample label changes. `label_at(x, y, bearing)` is called every `step` m; runs shorter than `min_run_m` join their longer neighbour (so a label flickering at a junction does not cut the line). Returns [(label, part)] in order; the parts are contiguous and together cover the whole polyline.""" xy = [to_xy(lat, lon) for lat, lon in coords] cum = [0.0] for a, b in zip(xy, xy[1:]): cum.append(cum[-1] + math.hypot(b[0] - a[0], b[1] - a[1])) total = cum[-1] if total <= 0.0: return [] n = max(1, int(math.ceil(total / step))) spacing = total / n labels = [] seg_i = 0 for k in range(n): s = (k + 0.5) * spacing while seg_i < len(xy) - 2 and cum[seg_i + 1] < s: seg_i += 1 (ax, ay), (bx, by) = xy[seg_i], xy[seg_i + 1] length = cum[seg_i + 1] - cum[seg_i] t = (s - cum[seg_i]) / length if length > 0 else 0.0 labels.append(label_at(ax + t * (bx - ax), ay + t * (by - ay), math.degrees(math.atan2(bx - ax, by - ay)) % 180.0)) runs = [] for k, label in enumerate(labels): if runs and runs[-1][0] == label: runs[-1][2] = k + 1 else: runs.append([label, k, k + 1]) while len(runs) > 1: short = min(range(len(runs)), key=lambda i: runs[i][2] - runs[i][1]) if (runs[short][2] - runs[short][1]) * spacing >= min_run_m: break left = runs[short - 1] if short > 0 else None right = runs[short + 1] if short + 1 < len(runs) else None if right is None or (left is not None and left[2] - left[1] >= right[2] - right[1]): left[2] = runs[short][2] else: right[1] = runs[short][1] del runs[short] merged = [] for run in runs: if merged and merged[-1][0] == run[0]: merged[-1][2] = run[2] else: merged.append(run) runs = merged if len(runs) == 1: return [(runs[0][0], coords)] out = [] for label, start, end in runs: part = sub_polyline(coords, cum, start * spacing, end * spacing) if len(part) >= 2: out.append((label, part)) return out class CityBoundary: """Polygon(s) of one city in the local metric projection: inside test (even-odd per feature, any feature) and distance to the boundary, with a coarse raster cache so most samples need no polygon arithmetic.""" def __init__(self, features: list, cell: float = 100.0, pad: float = 40.0): self.cell, self.pad = cell, pad self.polys = [] # per feature: list of rings [(x, y), ...] self.edges = [] # (x0, y0, x1, y1) for feature in features: geometry = feature.get("geometry") or {} if geometry.get("type") == "Polygon": polygons = [geometry.get("coordinates") or []] elif geometry.get("type") == "MultiPolygon": polygons = geometry.get("coordinates") or [] else: continue rings = [] for polygon in polygons: for ring in polygon: pts = [to_xy(pt[1], pt[0]) for pt in ring] if len(pts) >= 3: rings.append(pts) self.edges.extend((a[0], a[1], b[0], b[1]) for a, b in zip(pts, pts[1:] + pts[:1]) if a != b) if rings: self.polys.append(rings) self.band = 50.0 self.bands = collections.defaultdict(list) # y band -> [(poly id, edge)] poly_edges = [] for pid, rings in enumerate(self.polys): for pts in rings: for a, b in zip(pts, pts[1:] + pts[:1]): if a[1] != b[1]: poly_edges.append((pid, a[0], a[1], b[0], b[1])) for edge in poly_edges: for band in range(int(min(edge[2], edge[4]) // self.band), int(max(edge[2], edge[4]) // self.band) + 1): self.bands[band].append(edge) self.grid = collections.defaultdict(list) for i, (x0, y0, x1, y1) in enumerate(self.edges): for gx in range(int((min(x0, x1) - pad) // cell), int((max(x0, x1) + pad) // cell) + 1): for gy in range(int((min(y0, y1) - pad) // cell), int((max(y0, y1) + pad) // cell) + 1): self.grid[(gx, gy)].append(i) self.uniform: dict = {} def _inside(self, x: float, y: float) -> bool: crossings = collections.Counter() for pid, x0, y0, x1, y1 in self.bands.get(int(y // self.band), ()): if (y0 > y) != (y1 > y) and x < x0 + (y - y0) * (x1 - x0) / (y1 - y0): crossings[pid] += 1 return any(c % 2 for c in crossings.values()) def _edge_distance(self, x: float, y: float) -> float: best = float("inf") for i in self.grid.get((int(x // self.cell), int(y // self.cell)), ()): x0, y0, x1, y1 = self.edges[i] dx, dy = x1 - x0, y1 - y0 l2 = dx * dx + dy * dy t = 0.0 if l2 == 0 else max(0.0, min(1.0, ((x - x0) * dx + (y - y0) * dy) / l2)) best = min(best, math.hypot(x0 + t * dx - x, y0 + t * dy - y)) return best # inf when farther than `pad` from every edge def signed_distance(self, x: float, y: float) -> float: """Negative inside, positive outside: meters to the boundary, capped at +-pad.""" key = (int(x // self.cell), int(y // self.cell)) if key not in self.grid: # no edge within `pad` of this cell state = self.uniform.get(key) if state is None: state = self.uniform[key] = self._inside((key[0] + 0.5) * self.cell, (key[1] + 0.5) * self.cell) return -self.pad if state else self.pad d = min(self._edge_distance(x, y), self.pad) return -d if self._inside(x, y) else d def outside_m(self, x: float, y: float) -> float: return max(0.0, self.signed_distance(x, y)) def clip(self, coords: list, keep_m: float = CITY_KEEP_M) -> list: """The parts of a polyline that lie inside the city or at most keep_m outside it.""" runs = label_runs(coords, lambda x, y, b: self.signed_distance(x, y) <= keep_m) return [part for keep, part in runs if keep and polyline_length_m(part) >= 1.0] # --------------------------------------------------------------------------- official layers def ladot_designation_class(desig, mph: int) -> int: """BOE street designation (ST_DESIG) -> road class hint.""" if desig in (1, 2, 12, 13, 17, 23): return C_PRIMARY if desig in (3, 11, 15, 16, 21, 22, 52, 53, 54, 55, 66, 67): return C_SECONDARY if desig in (4, 5, 6, 19): return C_TERTIARY return C_PRIMARY if mph >= 40 else C_SECONDARY def class_from_mph(mph: int) -> int: if mph >= 40: return C_PRIMARY if mph >= 30: return C_TERTIARY return C_RESIDENTIAL def name_key(label: str) -> str: """Comparable street name: 'West Imperial Highway', 'Imperial Hwy' -> 'imperial highway'.""" tokens = (label or "").split("\t")[0].replace(".", " ").split() if len(tokens) > 1 and tokens[0].upper() in _DIRECTIONS: tokens = tokens[1:] if len(tokens) > 1 and tokens[-1].upper() in _DIRECTIONS: tokens = tokens[:-1] if len(tokens) > 1 and tokens[-1].upper() in _SUFFIXES: tokens[-1] = _SUFFIXES[tokens[-1].upper()] return " ".join(t.lower() for t in tokens) def street_name_full(full: str) -> str: """'W Glenoaks Blvd' (direction inside the name) -> 'West Glenoaks Boulevard'.""" tokens = (full or "").split() if len(tokens) > 1 and tokens[0].upper() in ("N", "S", "E", "W"): return street_name(tokens[0], " ".join(tokens[1:])) return street_name(None, full) class Builder: """Collects accepted records layer by layer. A layer is carved only against the layers committed before it (never against itself, so consecutive blocks of one street do not trim each other), except where a layer says otherwise.""" PAD = 40.0 def __init__(self): self.records: list = [] # (key, coords) self.families: list = [] # carving family of each official record self.index = LineIndex(cell=50.0, pad=self.PAD) self.payloads: list = [] # payload id -> (source, value, class, label, name key) self.pending: list = [] self.stats: dict = collections.OrderedDict() self.stale_candidates: list = [] # base parts <= 35 that no LAMC line covers (coords objects) def accept(self, key: tuple, coords: list, family: int) -> None: self.records.append((key, coords)) self.families.append(family) self.pending.append((key, coords, family)) def reindex(self) -> None: """Rebuild the carving index from the (official) records accepted so far.""" self.index = LineIndex(cell=50.0, pad=self.PAD) self.payloads = [] self.pending = [(key, coords, fam) for (key, coords), fam in zip(self.records, self.families)] self.commit() def commit(self) -> None: for key, coords, family in self.pending: payload = len(self.payloads) self.payloads.append((key[3], max(key[1], key[2]), key[0], key[5], name_key(key[5]))) self.index.add(coords, payload, family) self.pending = [] def carve(self, coords, max_d, families, tally=None, accept=None): return carve(coords, self.index, max_d, families, tally, accept=accept) def _length(parts) -> float: return sum(polyline_length_m(p) for p in parts) def build_ladot(builder: Builder, cache: str, city: CityBoundary, county_index) -> list: """All four LADOT layers (see the module docstring). Returns the accepted CVC pieces [(coords, sect_id, st_type, name)] for the review report.""" base = load_features(cache, "ladot_map") lamc = load_features(cache, "ladot_lamc35") ab43 = load_features(cache, "ladot_ab43") st = collections.Counter() transitions = collections.Counter() samples = collections.Counter() def in_city(coords, layer): out = [] for run in clip_to_bbox(coords): parts = city.clip(run) st[f"{layer}_km_clipped_outside_city"] += polyline_length_m(run) - _length(parts) out.extend(parts) return out # LAMC <= 35 lines: the per-sample override index (15 m) and the name retry (40 m). lamc_lines = [] lamc_index = LineIndex(cell=50.0, pad=LAMC_NAME_RETRY_M) for feature in lamc: props = feature.get("properties") or {} mph = parse_speed_field(props.get("CurrentSp")) name = street_name(None, props.get("STNAME_A")) for coords in feature_lines(feature): lamc_lines.append((coords, mph, name, name_key(name))) if mph: lamc_index.add(coords, len(lamc_lines) - 1, FAM_SURFACE) st["lamc_features"] += 1 st["lamc_features_without_limit" if not mph else "lamc_features_with_limit"] += 1 # Base pre-pass: the SECT_IDs it lists, and the 'outside' rows as a mask. base_sects = set() outside_mask = LineIndex(cell=50.0, pad=30.0) for feature in base: sect = ((feature.get("properties") or {}).get("SECT_ID") or "").strip() if sect.lower() == "outside": for coords in feature_lines(feature): outside_mask.add(coords, 0, FAM_SURFACE) elif sect: base_sects.add(sect) ab43_by_sect = collections.defaultdict(list) for feature in ab43: props = feature.get("properties") or {} proposed = parse_speed_field(props.get("Proposed_Speed_Limit_mph")) if proposed: ab43_by_sect[(props.get("SECT_ID") or "").strip()].append((proposed, feature)) # AB 43 rows that are not in the base by SECT_ID: they also lower the base parts they # run along (geometric join), and are added where no base line covers them. ab43_extra = LineIndex(cell=50.0, pad=OVERRIDE_MATCH_M) ab43_extra_rows = [] for sect, rows in ab43_by_sect.items(): if sect and sect in base_sects: continue for proposed, feature in rows: props = feature.get("properties") or {} name = street_name(props.get("STR_DIR"), props.get("STR_NAME")) for coords in feature_lines(feature): ab43_extra_rows.append((coords, proposed, name)) ab43_extra.add(coords, proposed, FAM_SURFACE) st["ab43_rows_not_in_base"] += 1 def base_pieces(feature, carve_first: bool) -> None: props = feature.get("properties") or {} sect = (props.get("SECT_ID") or "").strip() mph = parse_speed_field(props.get("SpeedLimit")) source = SRC_LADOT_MAP rows = ab43_by_sect.get(sect) if sect else None if rows: proposed = min(p for p, _ in rows) new = proposed if not mph else min(mph, proposed) # never raises a limit if new != mph: st["base_patched_by_ab43"] += 1 transitions[f"ab43 {mph}->{new} (records)"] += 1 mph, source = new, SRC_LADOT_AB43 else: st["base_ab43_no_change"] += 1 name = street_name(props.get("TDIR"), props.get("STNAME"), props.get("STSFX") or "", props.get("SFXDIR")) key_name = name_key(name) def same_street(i): return lamc_lines[i][3] == key_name def lower_ab43(value): return 0 < value < mph def label_at(x, y, bearing): hit = lamc_index.nearest(x, y, bearing, OVERRIDE_MATCH_M, CARVE_ANGLE) if hit is not None: return (SRC_LADOT_LAMC35, lamc_lines[hit[0][6]][1]) if mph: hit = ab43_extra.nearest(x, y, bearing, OVERRIDE_MATCH_M, CARVE_ANGLE, accept=lower_ab43) if hit is not None: return (SRC_LADOT_AB43, hit[0][6]) if source == SRC_LADOT_MAP and 0 < mph <= 35: hit = lamc_index.nearest(x, y, bearing, LAMC_NAME_RETRY_M, CARVE_ANGLE, accept=same_street) if hit is not None: samples["base_km_lamc35_matched_by_name_retry_40m"] += 1 return (SRC_LADOT_LAMC35, lamc_lines[hit[0][6]][1]) return ("stale", 0) return (source, mph) for coords in feature_lines(feature): for run in in_city(coords, "base"): parts = [run] if carve_first: parts = builder.carve(run, OFFICIAL_CARVE_M, FAM_SURFACE) st["base_no_sect_km_carved_by_sect_rows"] += polyline_length_m(run) - _length(parts) for part in parts: for (src, value), piece in label_runs(part, label_at): length = polyline_length_m(piece) if src == "stale": # kept for now; see stale_base_yields_to_higher_osm st["base_le35_no_lamc_match_km"] += length st["base_le35_no_lamc_match_pieces"] += 1 cls = ladot_designation_class(props.get("ST_DESIG"), mph) builder.accept((cls, mph, mph, SRC_LADOT_MAP, GROUND_LAYER_FLAGS, name), piece, FAM_SURFACE) builder.stale_candidates.append(piece) continue if not value: st["base_km_dropped_no_limit"] += length continue if src == SRC_LADOT_LAMC35: st["base_km_lamc35_" + ("same" if value == mph else "changed")] += length if value != mph: transitions[f"lamc35 {mph}->{value} (km)"] += length elif src == SRC_LADOT_AB43 and source != SRC_LADOT_AB43: st["base_km_lowered_by_ab43_geometric"] += length transitions[f"ab43 geometric {mph}->{value} (km)"] += length else: st["base_km_kept_" + SOURCE_NAMES[src]] += length cls = ladot_designation_class(props.get("ST_DESIG"), value) builder.accept((cls, value, value, src, GROUND_LAYER_FLAGS, name), piece, FAM_SURFACE) st["base_pieces_accepted"] += 1 for feature in base: sect = ((feature.get("properties") or {}).get("SECT_ID") or "").strip() st["base_features"] += 1 if sect.lower() == "outside": st["base_dropped_outside_city"] += 1 elif sect: base_pieces(feature, carve_first=False) builder.commit() for feature in base: # rows without a SECT_ID only fill what the SECT_ID rows leave if not ((feature.get("properties") or {}).get("SECT_ID") or "").strip(): st["base_features_without_sect_id"] += 1 base_pieces(feature, carve_first=True) builder.commit() def add_lines(lines, layer, make_key, name_retry=False): for coords, mph, name, *rest in lines: if not mph: continue key_name = name_key(name) for run in in_city(coords, layer): parts = builder.carve(run, OFFICIAL_CARVE_M, FAM_SURFACE) if name_retry: # a LAMC line offset from its base street parts = [q for p in parts for q in builder.carve(p, LAMC_NAME_RETRY_M, FAM_SURFACE, accept=lambda pl: builder.payloads[pl][4] == key_name)] kept = [q for p in parts for q in carve(p, outside_mask, OFFICIAL_CARVE_M, FAM_SURFACE)] st[f"{layer}_km_masked_by_outside_rows"] += _length(parts) - _length(kept) for part in kept: builder.accept(make_key(mph, name), part, FAM_SURFACE) st[f"{layer}_added_pieces"] += 1 st[f"{layer}_added_km"] += polyline_length_m(part) builder.commit() add_lines(lamc_lines, "lamc35", lambda mph, name: (class_from_mph(mph), mph, mph, SRC_LADOT_LAMC35, GROUND_LAYER_FLAGS, name), name_retry=True) add_lines(ab43_extra_rows, "ab43", lambda mph, name: (class_from_mph(mph), mph, mph, SRC_LADOT_AB43, GROUND_LAYER_FLAGS, name)) # CVC: LA City streets not in the LAMC, 25 mph prima facie. cvc = load_features(cache, "ladot_cvc") excluded_types = {"WK": "walk street", "VA": "vacated", "RP": "Rec & Parks road", "HA": "Housing Authority road", "PA": "private", "ST": "state highway"} cvc_pieces = [] def on_city_edge(x, y): return city.signed_distance(x, y) >= -COUNTY_EDGE_M for feature in cvc: props = feature.get("properties") or {} st["cvc_features"] += 1 sect = (props.get("SECT_ID") or "").strip() st_type = (props.get("ST_TYPE") or "").strip().upper() if sect and sect in base_sects: st["cvc_dropped_sect_in_lamc_base"] += 1 continue if st_type in excluded_types: st[f"cvc_dropped_type_{st_type}_{excluded_types[st_type].replace(' ', '_')}"] += 1 continue name = street_name(props.get("ST_DIR"), props.get("ST_NAME")) key_name = name_key(name) def faster_same_street(pl, key_name=key_name): return builder.payloads[pl][4] == key_name and builder.payloads[pl][1] > 25 for coords in feature_lines(feature): for run in in_city(coords, "cvc"): parts = builder.carve(run, OFFICIAL_CARVE_M, FAM_SURFACE) st["cvc_km_carved_by_lamc_layers"] += polyline_length_m(run) - _length(parts) kept = [q for p in parts for q in carve(p, outside_mask, OFFICIAL_CARVE_M, FAM_SURFACE)] st["cvc_km_masked_by_outside_rows"] += _length(parts) - _length(kept) parts = kept kept = [q for p in parts for q in builder.carve(p, FRONTAGE_CARVE_M, FAM_SURFACE, accept=faster_same_street)] st["cvc_km_dropped_faster_same_street_within_25m"] += _length(parts) - _length(kept) parts = kept if county_index is not None: kept = [q for p in parts for q in carve(p, county_index, OFFICIAL_CARVE_M, FAM_SURFACE, accept=lambda v: v != 25, point_filter=on_city_edge)] st["cvc_km_dropped_county_posted_on_city_edge"] += _length(parts) - _length(kept) parts = kept for part in parts: builder.accept((C_RESIDENTIAL, 25, 25, SRC_LADOT_CVC, GROUND_LAYER_FLAGS | FLAG_PRIMA_FACIE, name), part, FAM_SURFACE) cvc_pieces.append((part, sect, st_type, name)) st["cvc_pieces_accepted"] += 1 builder.commit() for key, count in samples.items(): st[key] += count * CARVE_STEP_M st["transitions"] = {k: (km(v) if "(km)" in k else v) for k, v in transitions.most_common()} builder.stats["ladot"] = {k: (km(v) if "_km" in k and not isinstance(v, dict) else v) for k, v in st.items()} return cvc_pieces def county_lines(cache: str) -> list: out = [] for feature in load_features(cache, "la_county"): props = feature.get("properties") or {} mph = parse_speed_field(props.get("Speed")) out.append((feature, props, mph)) return out def county_posted_index(lines: list) -> LineIndex: index = LineIndex(cell=50.0, pad=30.0) for feature, _, mph in lines: if mph: for coords in feature_lines(feature): index.add(coords, mph, FAM_SURFACE) return index def build_county(builder: Builder, lines: list) -> None: st = collections.Counter() for feature, props, mph in lines: st["features"] += 1 raw = props.get("Speed") if not mph: st["dropped_conflicting_or_invalid" if raw else "dropped_no_speed"] += 1 continue name = street_name(None, props.get("Rd_Dissolve_SplitLineAtPoint_STREET")) for coords in feature_lines(feature): for run in clip_to_bbox(coords): parts = builder.carve(run, OFFICIAL_CARVE_M, FAM_SURFACE) st["km_carved_by_ladot"] += polyline_length_m(run) - _length(parts) for part in parts: builder.accept((class_from_mph(mph), mph, mph, SRC_COUNTY, GROUND_LAYER_FLAGS, name), part, FAM_SURFACE) st["pieces_accepted"] += 1 builder.commit() builder.stats["la_county"] = {k: (km(v) if k.startswith("km_") else v) for k, v in st.items()} # --------------------------------------------------------------------------- document overrides def load_overrides(path) -> list: if not path: return [] with open(path, encoding="utf-8") as f: doc = json.load(f) entries = doc["overrides"] if isinstance(doc, dict) else doc for entry in entries: assert entry["source"] in ("glendale", "ladot_document"), entry assert valid_mph(int(entry["mph"])), entry for field in ("id", "street", "from", "to", "citation", "reviewBy"): assert field in entry, (field, entry) return entries def _span_test(entry: dict): """(x, y) -> True when the point projects between the entry's from and to points (cross-street intersections) and lies within `maxOffsetM` of that chord.""" ax, ay = to_xy(*entry["from"]["at"]) bx, by = to_xy(*entry["to"]["at"]) dx, dy = bx - ax, by - ay l2 = dx * dx + dy * dy length = math.sqrt(l2) max_off = float(entry.get("maxOffsetM", 40.0)) def inside(x, y): t = ((x - ax) * dx + (y - ay) * dy) / l2 if t < 0.0 or t > 1.0: return False return abs((x - ax) * dy - (y - ay) * dx) / length <= max_off return inside def apply_overrides(builder: Builder, cache: str, entries: list) -> dict: """Official-document values (see the module docstring, item 0). 'ladot_document' entries replace the value of the LADOT records of the named street inside the span; 'glendale' entries add County centerline geometry of the named street inside the span, carved by the official records already accepted (the OSM layer is carved by them).""" report = [] pattern_cache = {} cams = None for entry in entries: source = SOURCE_CODES[entry["source"]] mph = int(entry["mph"]) inside = _span_test(entry) street = pattern_cache.setdefault(entry["street"], re.compile(entry["street"], re.IGNORECASE)) meters = collections.Counter() if source == SRC_LADOT_DOC: records, families = [], [] for (key, coords), fam in zip(builder.records, builder.families): cls, fwd, bwd, src, flags, label = key if src not in LADOT_SOURCES or not street.search(label.split("\t")[0]): records.append((key, coords)) families.append(fam) continue for hit, part in label_runs(coords, lambda x, y, b: inside(x, y)): if hit: meters[f"{SOURCE_NAMES[src]} {max(fwd, bwd)} -> {mph}"] += polyline_length_m(part) part_key = (cls, mph, mph, source, flags & ~FLAG_PRIMA_FACIE, label) else: part_key = key records.append((part_key, part)) families.append(fam) builder.records, builder.families = records, families builder.reindex() else: if cams is None: cams = load_features(cache, OVERRIDE_GEOMETRY_SOURCES["glendale"]) for feature in cams: props = feature.get("properties") or {} if (props.get("Type") not in ("Primary", "Secondary") or props.get("Elevation") not in (None, "Surface") or not street.search(props.get("FullName") or "")): continue name = street_name_full(props.get("FullName")) cls = C_PRIMARY if props.get("Type") == "Primary" else C_SECONDARY for coords in feature_lines(feature): for run in clip_to_bbox(coords): for hit, part in label_runs(run, lambda x, y, b: inside(x, y)): if not hit: continue for piece in builder.carve(part, OFFICIAL_CARVE_M, FAM_SURFACE): builder.accept((cls, mph, mph, source, GROUND_LAYER_FLAGS, name), piece, FAM_SURFACE) meters[f"added {mph}"] += polyline_length_m(piece) builder.commit() total = sum(meters.values()) report.append({"id": entry["id"], "source": entry["source"], "mph": mph, "km": km(total), "detailKm": {k: km(v) for k, v in meters.items()}, "citation": entry["citation"], "reviewBy": entry["reviewBy"]}) if total < 50.0: raise RuntimeError(f"override {entry['id']}: only {total:.0f} m matched; check its street and span") return {"entries": report} def caltrans_freeway_parts(route: str, mph: int, coords: list) -> list: """[(is_freeway, part)]: routes in CALTRANS_FREEWAY_ROUTES are freeway everywhere, a route with sections only inside its section rectangles, and never below the minimum.""" if mph < CALTRANS_FREEWAY_MIN_MPH: return [(False, coords)] if route in CALTRANS_FREEWAY_ROUTES: return [(True, coords)] rects = CALTRANS_FREEWAY_SECTIONS.get(route) if not rects: return [(False, coords)] xy_rects = [(to_xy(s_, w), to_xy(n, e)) for w, s_, e, n in rects] def in_section(x, y, bearing): return any(lo[0] <= x <= hi[0] and lo[1] <= y <= hi[1] for lo, hi in xy_rects) return label_runs(coords, in_section) def build_caltrans(builder: Builder, cache: str, use_surface: bool) -> tuple: """Freeway mainlines become records. Returns (surface lines for QA, cap index): the cap index holds every Caltrans line that did not become a record, so an OSM motorway or trunk along it cannot carry a higher value (see build_osm).""" features = load_features(cache, "caltrans_d7") st = collections.Counter() surface = [] cap = LineIndex(cell=50.0, pad=30.0) for feature in features: props = feature.get("properties") or {} route = str(props.get("route") or "").strip() mph = parse_speed_field(props.get("speed_limit")) st["features"] += 1 if not mph: st["dropped_no_speed"] += 1 continue for coords in feature_lines(feature): for run in clip_to_bbox(coords): for freeway, piece in caltrans_freeway_parts(route, mph, run): length = polyline_length_m(piece) if not freeway: surface.append((route, mph, piece)) st["surface_km"] += length if (route in CALTRANS_FREEWAY_ROUTES or route in CALTRANS_FREEWAY_SECTIONS) \ and mph < CALTRANS_FREEWAY_MIN_MPH: st[f"freeway_route_value_rejected_below_{CALTRANS_FREEWAY_MIN_MPH}_km"] += length if not use_surface: cap.add(piece, (route, mph), FAM_SURFACE | FAM_FREEWAY) continue ref = caltrans_ref(route) cls = C_MOTORWAY if freeway else C_PRIMARY family = FAM_FREEWAY if freeway else FAM_SURFACE parts = [piece] if freeway else builder.carve(piece, OSM_CARVE_M, FAM_SURFACE) for part in parts: builder.accept((cls, mph, mph, SRC_CALTRANS, GROUND_LAYER_FLAGS, f"\t{ref}"), part, family) st["freeway_km" if freeway else "surface_km_used"] += polyline_length_m(part) builder.stats["caltrans"] = {k: (km(v) if k.endswith("_km") or "_km_" in k else v) for k, v in st.items()} return surface, cap def caltrans_surface_adjacency(builder: Builder) -> dict: """QA only: Caltrans freeway records running within 8 m of an aligned official surface record (a matcher on the surface street can snap to the freeway there; the class tells the two apart). Freeway and surface lines never carve each other.""" tally = collections.Counter() where = {} for key, coords in builder.records: if key[3] != SRC_CALTRANS or key[0] != C_MOTORWAY: continue def record(pl, meters, distance, x, y, key=key): other = builder.payloads[pl] name = f"{key[5].strip()} vs {other[3]} ({other[1]})" tally[name] += meters where.setdefault(name, [round(SOUTH + y / M_PER_DEG_LAT, 5), round(WEST + x / M_PER_DEG_LON, 5)]) carve(coords, builder.index, OVERLAP_GATE_M, FAM_SURFACE, record, accept=lambda pl: builder.payloads[pl][0] != SRC_CALTRANS) return {"km": km(sum(tally.values())), "top": [{"pair": k, "km": km(v), "at": where[k]} for k, v in tally.most_common(10)]} def caltrans_surface_qa(builder: Builder, surface: list) -> dict: """Compare Caltrans surface-highway values with the accepted LADOT / County lines.""" tally = collections.Counter() diffs = collections.Counter() for route, mph, coords in surface: def record(payload, meters, distance, x, y, route=route, mph=mph): source, value = builder.payloads[payload][:2] group = "ladot" if source in LADOT_SOURCES else SOURCE_NAMES[source] tally[f"{group}_matched_m"] += meters if value == mph: tally[f"{group}_agree_m"] += meters else: tally[f"{group}_caltrans_{'higher' if mph > value else 'lower'}_m"] += meters diffs[f"SR {route}: caltrans {mph} vs {group} {value}"] += meters carve(coords, builder.index, OSM_CARVE_M, FAM_SURFACE, record) out = {k[:-2] + "_km": km(v) for k, v in sorted(tally.items())} out["top_differences_km"] = {k: km(v) for k, v in diffs.most_common(12)} return out OSM_CHAINS_VERSION = 3 # 3: exact and reversed duplicate pieces removed before chaining def osm_chains(cache: str, pbf_paths: list) -> tuple: """Normalized, clipped and chained OSM ways, cached as a pickle keyed on the PBFs.""" stamp = [(os.path.basename(p), os.path.getsize(p), int(os.path.getmtime(p))) for p in pbf_paths] derived = os.path.join(cache, "derived") path = os.path.join(derived, "osm_chains.pickle") if os.path.exists(path): with open(path, "rb") as f: saved = pickle.load(f) if saved.get("stamp") == stamp and saved.get("version") == OSM_CHAINS_VERSION: log(f"cached OSM chains ({len(saved['chains'])})") return saved["chains"], saved["timestamp"], collections.Counter(saved["issues"]) raw_ways, osm_time = read_pbf_ways(pbf_paths) issues: collections.Counter = collections.Counter() pieces = [] for _, tags, coords in raw_ways: normalized = normalize_osm_way(tags, issues) if normalized is None: continue key, reverse = normalized if reverse: coords = coords[::-1] for run in clip_to_bbox(coords): pieces.append((key, run)) del raw_ways pieces, dropped = dedupe_pieces(pieces) issues["osm_duplicate_way_pieces_dropped"] += dropped chains = merge_chains(pieces) log(f"OSM: {len(pieces)} pieces in the box -> {len(chains)} chains") os.makedirs(derived, exist_ok=True) with open(path + ".part", "wb") as f: pickle.dump({"version": OSM_CHAINS_VERSION, "stamp": stamp, "chains": chains, "timestamp": osm_time, "issues": dict(issues)}, f, protocol=pickle.HIGHEST_PROTOCOL) os.replace(path + ".part", path) return chains, osm_time, issues def stale_base_yields_to_higher_osm(builder: Builder, chains: list) -> dict: """A base-map part at 35 mph or less that no LAMC_35 line covers is probably a stale 2018 value (the Jan 2025 LAMC_35 layer lists the City streets posted 35 or less), but not always (a 2026 camera document confirms one such 35). Where an OSM way with a higher value runs along it, the part is removed so the OSM value fills it; elsewhere the base value stays. Off with --strict-precedence and in --no-osm builds.""" stale = {id(c) for c in builder.stale_candidates} if not stale: return {} high = LineIndex(cell=50.0, pad=30.0) for key, coords in chains: cls, fwd, bwd, source, flags, label = key if source == SRC_OSM and cls not in (C_MOTORWAY, C_MOTORWAY_LINK): high.add(coords, max(fwd, bwd), FAM_SURFACE) st = collections.Counter() records, families = [], [] for (key, coords), fam in zip(builder.records, builder.families): if id(coords) not in stale: records.append((key, coords)) families.append(fam) continue value = max(key[1], key[2]) def label_at(x, y, bearing, value=value): hit = high.nearest(x, y, bearing, OSM_CARVE_M, CARVE_ANGLE) return hit[0][6] if hit is not None and hit[0][6] > value else None for osm_value, part in label_runs(coords, label_at): length = polyline_length_m(part) if osm_value: st[f"base {value} -> osm {osm_value}"] += length continue st["kept_base_value"] += length records.append((key, part)) families.append(fam) builder.records, builder.families = records, families builder.reindex() return {"km": {k: km(v) for k, v in st.most_common()}} CONFLICT_RANK = {SRC_LADOT_DOC: 0, SRC_GLENDALE: 0, SRC_LADOT_LAMC35: 1, SRC_LADOT_AB43: 2, SRC_LADOT_MAP: 3, SRC_LADOT_CVC: 4, SRC_COUNTY: 5, SRC_CALTRANS: 6} def resolve_official_conflicts(records: list) -> tuple: """Official records that still run on top of each other (within OVERLAP_GATE_M, aligned, same family) with different limits, typically two rows of one layer (layers are not carved against themselves) or a shallow junction: the record of the higher-ranked source, or of the same source the longer one, keeps the overlap; the other is cut there. Runs on the final, simplified records, so the overlap gate checks the same geometry.""" official = [i for i, (key, _) in enumerate(records) if key[3] in OFFICIAL_SOURCES and max(key[1], key[2])] order = sorted(official, key=lambda i: (CONFLICT_RANK[records[i][0][3]], -polyline_length_m(records[i][1]))) kept_index = LineIndex(cell=50.0, pad=OVERLAP_GATE_M) out = {} st = collections.Counter() for i in order: key, coords = records[i] fam = _family_of(key[0]) value = max(key[1], key[2]) parts = carve(coords, kept_index, OVERLAP_GATE_M, fam, accept=lambda v: v != value) removed = polyline_length_m(coords) - _length(parts) if removed > 0.5: st[f"{SOURCE_NAMES[key[3]]} km cut"] += removed st[f"{SOURCE_NAMES[key[3]]} records cut"] += 1 out[i] = parts for part in parts: kept_index.add(part, value, fam) result = [] for i, item in enumerate(records): if i in out: result.extend((item[0], part) for part in out[i]) else: result.append(item) return result, {k: (km(v) if "km" in k else v) for k, v in st.most_common()} def cvc_yields_to_lower_osm(builder: Builder, chains: list) -> dict: """Safety rule (on unless --strict-precedence): the CVC 25 mph value is a prima facie default, which a lower posted limit overrides. Where an OSM way with a limit below 25 runs along a CVC piece, that part of the CVC piece is removed, so the OSM value fills the gap (as source osm) instead of a 25 that may overstate the posted limit.""" low = LineIndex(cell=50.0, pad=30.0) count = 0 for key, coords in chains: cls, fwd, bwd, source, flags, label = key if source == SRC_OSM and cls not in (C_MOTORWAY, C_MOTORWAY_LINK) and 0 < max(fwd, bwd) < 25: low.add(coords, max(fwd, bwd), FAM_SURFACE) count += 1 removed = collections.Counter() records, families = [], [] for (key, coords), fam in zip(builder.records, builder.families): if key[3] != SRC_LADOT_CVC: records.append((key, coords)) families.append(fam) continue tally = collections.Counter() parts = carve(coords, low, OSM_CARVE_M, FAM_SURFACE, lambda p, m, d, x, y: tally.update({p: m})) if not tally: records.append((key, coords)) families.append(fam) continue removed.update(tally) for part in parts: records.append((key, part)) families.append(fam) builder.records, builder.families = records, families builder.reindex() return {"osm_ways_below_25": count, "cvc_km_yielded_by_osm_value": {f"{mph} mph": km(m) for mph, m in sorted(removed.items())}} def cap_osm_freeways(chains: list, cap: LineIndex) -> tuple: """An OSM motorway or trunk that runs along a Caltrans line that did not become a record (a value below the freeway minimum, or a surface section) may not carry a value above the Caltrans one there: those parts keep their geometry with no limit.""" out = [] capped = collections.Counter() for key, coords in chains: cls, fwd, bwd, source, flags, label = key value = max(fwd, bwd) if source != SRC_OSM or cls not in (C_MOTORWAY, C_TRUNK): out.append((key, coords)) continue def label_at(x, y, bearing, value=value): hit = cap.nearest(x, y, bearing, OSM_FREEWAY_CARVE_M, CARVE_ANGLE, accept=lambda pl: pl[1] < value) return None if hit is None else hit[0][6] for hit, part in label_runs(coords, label_at): if hit is None: out.append((key, part)) else: out.append(((cls, 0, 0, SRC_NONE, flags & ~FLAG_CONDITIONAL, label), part)) capped[f"route {hit[0]}: osm {value} above caltrans {hit[1]} ({label.replace(chr(9), ' | ')})"] += \ polyline_length_m(part) return out, {k: km(v) for k, v in capped.most_common()} def _sections(samples: list, gap_m: float = 60.0) -> list: """Group (chain, x, y, meters) samples into contiguous sections.""" sections = [] for chain, x, y, meters in samples: last = sections[-1] if sections else None if last and last["chain"] == chain and math.hypot(x - last["_x"], y - last["_y"]) <= gap_m: last["_m"] += meters last["_x"], last["_y"] = x, y else: sections.append({"chain": chain, "_x0": x, "_y0": y, "_x": x, "_y": y, "_m": meters}) out = [] for sec in sections: def ll(x, y): return [round(SOUTH + y / M_PER_DEG_LAT, 5), round(WEST + x / M_PER_DEG_LON, 5)] out.append({"from": ll(sec["_x0"], sec["_y0"]), "to": ll(sec["_x"], sec["_y"]), "km": round(sec["_m"] / 1000, 2)}) return out def build_osm(builder: Builder, chains: list, cap=None) -> dict: """Adds the OSM fill. Returns review lists (Caltrans vs OSM disagreement sections).""" st = collections.Counter() agreement = collections.Counter() diffs = collections.Counter() offsets = collections.Counter() by_route = collections.defaultdict(collections.Counter) official_higher = collections.Counter() disagreement = collections.defaultdict(list) if cap is not None: chains, capped = cap_osm_freeways(chains, cap) st_capped = capped else: st_capped = {} for chain_no, (key, coords) in enumerate(chains): cls, fwd, bwd, source, flags, label = key osm_value = max(fwd, bwd) def record(payload, meters, distance, x, y, source=source, fwd=fwd, bwd=bwd, cls=cls, label=label, chain_no=chain_no): off_source, value, _, off_label = builder.payloads[payload][:4] group = "ladot" if off_source in LADOT_SOURCES else SOURCE_NAMES[off_source] fam = "freeway" if cls in (C_MOTORWAY, C_TRUNK) and off_source == SRC_CALTRANS else "surface" agreement[f"{group}_{fam}_matched_m"] += meters route = off_label.split("\t")[-1] if off_source == SRC_CALTRANS else None if route: offsets[min(int(distance // 2) * 2, 18)] += meters if source == SRC_OSM: if value in (fwd, bwd): outcome = "same" else: outcome = "osm_higher" if max(fwd, bwd) > value else "osm_lower" diffs[f"{group}: osm {max(fwd, bwd)} vs {value}"] += meters if outcome == "osm_lower": street = off_label.split("\t")[0] or off_label.strip() official_higher[f"{street} ({group} {value}, osm {max(fwd, bwd)})"] += meters else: outcome = "osm_no_limit" agreement[f"{group}_{fam}_{outcome}_m"] += meters if route: by_route[route][f"caltrans {value} / {outcome}"] += meters if source == SRC_OSM and value not in (fwd, bwd) and fam == "freeway": disagreement[(route, value, max(fwd, bwd))].append((chain_no, x, y, meters)) length = polyline_length_m(coords) if cls in (C_MOTORWAY_LINK, C_TRUNK_LINK): parts = [coords] # ramps: no official data covers them elif cls == C_MOTORWAY: parts = builder.carve(coords, OSM_FREEWAY_CARVE_M, FAM_FREEWAY, record) elif cls == C_TRUNK: parts = builder.carve(coords, OSM_CARVE_M, FAM_SURFACE | FAM_FREEWAY, record) else: parts = builder.carve(coords, OSM_CARVE_M, FAM_SURFACE, record) kept = sum(polyline_length_m(p) for p in parts) bucket = "osm" if source == SRC_OSM else "none" st[f"input_km_{bucket}"] += length st[f"kept_km_{bucket}"] += kept for part in parts: builder.records.append((key, part)) # not indexed: nothing below OSM builder.stats["osm"] = {k: km(v) for k, v in sorted(st.items())} builder.stats["osm"]["agreement_with_official_km"] = {k[:-2] + "_km": km(v) for k, v in sorted(agreement.items())} builder.stats["osm"]["official_higher_than_osm_top_streets_km"] = { k: km(v) for k, v in official_higher.most_common(20)} builder.stats["osm"]["caltrans_vs_osm_by_route_km"] = { route: {k: km(v) for k, v in sorted(c.items()) if v >= 50} for route, c in sorted(by_route.items(), key=lambda kv: kv[0])} builder.stats["osm"]["top_differences_km"] = {k: km(v) for k, v in diffs.most_common(15)} builder.stats["osm"]["caltrans_vs_osm_carriageway_offset_km_by_2m_bin"] = { f"{b}-{b + 2}m": km(v) for b, v in sorted(offsets.items())} builder.stats["osm"]["osm_freeway_values_capped_by_unused_caltrans_lines_km"] = st_capped review = [] for (route, cal, osm), samples in sorted(disagreement.items(), key=lambda kv: -sum(m for *_, m in kv[1])): sections = [sec for sec in _sections(samples) if sec["km"] >= 0.2] review.append({"route": route, "caltrans": cal, "osm": osm, "km": round(sum(m for *_, m in samples) / 1000, 1), "direction": "caltrans lower (false red if OSM is right)" if cal < osm else "caltrans higher (under-warning if OSM is right)", "sections": [{k: v for k, v in sec.items() if k != "chain"} for sec in sections]}) return {"caltransVsOsmFreeway": review} def apply_suppressions(records: list, path: str) -> tuple: """Suppress list entries only ever REMOVE OSM limits: matching records keep their geometry with source none and speed 0. Entry: {"name": regex, "bbox": [w, s, e, n], "aboveMph": optional int, "reason", "addedOn", "reviewBy"}.""" with open(path, encoding="utf-8") as f: entries = json.load(f) report = [] out = list(records) for entry in entries: pattern = re.compile(entry["name"], re.IGNORECASE) w, s, e, n = entry["bbox"] above = int(entry.get("aboveMph", 0)) hits = 0 meters = 0.0 for i, (key, coords) in enumerate(out): cls, fwd, bwd, source, flags, label = key if source != SRC_OSM or max(fwd, bwd) <= above: continue if not pattern.search(label.split("\t")[0]): continue lat, lon = coords[len(coords) // 2] if not (w <= lon <= e and s <= lat <= n): continue out[i] = ((cls, 0, 0, SRC_NONE, flags & ~FLAG_CONDITIONAL, label), coords) hits += 1 meters += polyline_length_m(coords) report.append({"name": entry["name"], "reason": entry.get("reason"), "records": hits, "km": km(meters), "reviewBy": entry.get("reviewBy")}) return out, report # --------------------------------------------------------------------------- writer def quantize(lat: float, lon: float): return (min(LAT_SPAN_Q, max(0, round((lat - SOUTH) * Q_PER_DEG))), min(LON_SPAN_Q, max(0, round((lon - WEST) * Q_PER_DEG)))) def _morton(a: int, b: int) -> int: code = 0 for bit in range(16): code |= ((a >> bit) & 1) << (2 * bit) | ((b >> bit) & 1) << (2 * bit + 1) return code def write_resource(records: list, out_path: str, data_date: int, osm_date: int, margin_m: int) -> dict: """records: [(key, coords)] -> binary file. Returns layout numbers.""" cols = math.ceil(LON_SPAN_Q / CELL_Q) + (1 if LON_SPAN_Q % CELL_Q == 0 else 0) rows = math.ceil(LAT_SPAN_Q / CELL_Q) + (1 if LAT_SPAN_Q % CELL_Q == 0 else 0) quantized = [] for key, coords in records: points = [] for lat, lon in coords: q = quantize(lat, lon) if not points or points[-1] != q: points.append(q) if len(points) >= 2: quantized.append((key, points)) quantized.sort(key=lambda item: (_morton(item[1][len(item[1]) // 2][0] >> 4, item[1][len(item[1]) // 2][1] >> 4), item[1])) labels = sorted({key[5] for key, _ in quantized if key[5]}) if len(labels) >= NO_NAME: raise RuntimeError(f"{len(labels)} distinct names exceed the u16 name index") label_index = {label: i for i, label in enumerate(labels)} seg_bytes = bytearray(SEGMENT_STRUCT.size * len(quantized)) points_arr = array.array("H") cell_sets = [None] * (rows * cols) margin_lat = math.ceil(margin_m / (M_PER_DEG_LAT / Q_PER_DEG)) margin_lon = math.ceil(margin_m / (M_PER_DEG_LON / Q_PER_DEG)) sources_mask = 0 for rec, (key, points) in enumerate(quantized): cls, fwd, bwd, source, flags, label = key sources_mask |= 1 << source SEGMENT_STRUCT.pack_into(seg_bytes, rec * SEGMENT_STRUCT.size, len(points_arr) // 2, len(points), label_index.get(label, NO_NAME), fwd, bwd, cls | (source << 4), flags) for lat_q, lon_q in points: points_arr.append(lat_q) points_arr.append(lon_q) for (la0, lo0), (la1, lo1) in zip(points, points[1:]): r0 = max(0, (min(la0, la1) - margin_lat) // CELL_Q) r1 = min(rows - 1, (max(la0, la1) + margin_lat) // CELL_Q) c0 = max(0, (min(lo0, lo1) - margin_lon) // CELL_Q) c1 = min(cols - 1, (max(lo0, lo1) + margin_lon) // CELL_Q) single = r0 == r1 and c0 == c1 for r in range(r0, r1 + 1): for c in range(c0, c1 + 1): if not single and _clip_segment( lo0, la0, lo1, la1, c * CELL_Q - margin_lon, r * CELL_Q - margin_lat, (c + 1) * CELL_Q + margin_lon, (r + 1) * CELL_Q + margin_lat) is None: continue cell = r * cols + c bucket = cell_sets[cell] if bucket is None: cell_sets[cell] = [rec] elif bucket[-1] != rec: bucket.append(rec) cell_start = array.array("I", [0]) cell_refs = array.array("I") for bucket in cell_sets: if bucket: cell_refs.extend(bucket) cell_start.append(len(cell_refs)) name_offsets = array.array("I", [0]) name_blob = bytearray() for label in labels: name_blob += label.encode("utf-8") name_offsets.append(len(name_blob)) name_length = len(name_blob) name_blob += b"\0" * (-len(name_blob) % 4) if sys.byteorder != "little": for arr in (points_arr, cell_start, cell_refs, name_offsets): arr.byteswap() segments_offset = HEADER_SIZE points_offset = segments_offset + len(seg_bytes) cell_start_offset = points_offset + len(points_arr) * 2 cell_refs_offset = cell_start_offset + len(cell_start) * 4 name_offsets_offset = cell_refs_offset + len(cell_refs) * 4 name_bytes_offset = name_offsets_offset + len(name_offsets) * 4 file_length = name_bytes_offset + len(name_blob) body = b"".join([bytes(seg_bytes), points_arr.tobytes(), cell_start.tobytes(), cell_refs.tobytes(), name_offsets.tobytes(), bytes(name_blob)]) assert HEADER_SIZE + len(body) == file_length build_date = int(dt.datetime.now(dt.timezone.utc).strftime("%Y%m%d")) header = HEADER_STRUCT.pack( MAGIC, FORMAT_VERSION, HEADER_SIZE, ORIGIN_LAT_E6, ORIGIN_LON_E6, QUANTUM_E6, CELL_Q, cols, rows, margin_m, len(quantized), len(points_arr) // 2, len(cell_refs), len(labels), data_date, build_date, segments_offset, points_offset, cell_start_offset, cell_refs_offset, name_offsets_offset, name_bytes_offset, name_length, file_length, zlib.crc32(body) & 0xFFFFFFFF, sources_mask, osm_date) assert len(header) == HEADER_SIZE os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True) with open(out_path + ".part", "wb") as f: f.write(header) f.write(body) os.replace(out_path + ".part", out_path) return { "records": len(quantized), "points": len(points_arr) // 2, "cellRefs": len(cell_refs), "names": len(labels), "gridCols": cols, "gridRows": rows, "sourcesMask": sources_mask, "bytes": {"header": HEADER_SIZE, "segments": len(seg_bytes), "points": len(points_arr) * 2, "cellStart": len(cell_start) * 4, "cellRefs": len(cell_refs) * 4, "nameOffsets": len(name_offsets) * 4, "nameBytes": len(name_blob), "total": file_length}, } # --------------------------------------------------------------------------- reader class StreetLimitFile: """Reference reader. The Swift reader in the app must behave the same way.""" def __init__(self, path: str, check_crc: bool = True): with open(path, "rb") as f: self.buf = f.read() if len(self.buf) < HEADER_SIZE: raise ValueError("file shorter than the header") (magic, version, header_size, self.origin_lat_e6, self.origin_lon_e6, self.quantum_e6, self.cell_q, self.cols, self.rows, self.margin_m, self.count, self.point_count, self.ref_count, self.name_count, self.data_date, self.build_date, self.seg_off, self.pts_off, self.cs_off, self.cr_off, self.no_off, self.nb_off, self.nb_len, self.file_length, self.crc, self.sources_mask, self.osm_date) = HEADER_STRUCT.unpack_from(self.buf, 0) if magic != MAGIC or version != FORMAT_VERSION or header_size != HEADER_SIZE: raise ValueError(f"not a street_limits.bin v{FORMAT_VERSION} file (magic {magic!r}, version {version})") if self.file_length != len(self.buf): raise ValueError("file length mismatch") if check_crc and zlib.crc32(self.buf[HEADER_SIZE:]) & 0xFFFFFFFF != self.crc: raise ValueError("CRC mismatch") sections = [(self.seg_off, self.count * 12), (self.pts_off, self.point_count * 4), (self.cs_off, (self.cols * self.rows + 1) * 4), (self.cr_off, self.ref_count * 4), (self.no_off, (self.name_count + 1) * 4), (self.nb_off, self.nb_len)] previous_end = HEADER_SIZE for offset, length in sections: if offset % 4 or offset < previous_end or offset + length > len(self.buf): raise ValueError(f"section at {offset} out of order, misaligned or out of bounds") previous_end = offset + length self.points = array.array("H", self.buf[self.pts_off:self.pts_off + self.point_count * 4]) self.cell_start = array.array("I", self.buf[self.cs_off:self.cs_off + (self.cols * self.rows + 1) * 4]) self.cell_refs = array.array("I", self.buf[self.cr_off:self.cr_off + self.ref_count * 4]) self.name_offsets = array.array("I", self.buf[self.no_off:self.no_off + (self.name_count + 1) * 4]) if sys.byteorder != "little": for arr in (self.points, self.cell_start, self.cell_refs, self.name_offsets): arr.byteswap() self.q_per_deg = 1e6 / self.quantum_e6 def raw_record(self, i: int) -> tuple: return SEGMENT_STRUCT.unpack_from(self.buf, self.seg_off + i * SEGMENT_STRUCT.size) def record(self, i: int) -> dict: first, count, name, fwd, bwd, class_source, flags = self.raw_record(i) source = class_source >> 4 return {"index": i, "first": first, "count": count, "name": self.name(name), "fwd": fwd, "bwd": bwd, "class": ROAD_CLASSES[class_source & 15] if class_source & 15 < 13 else "?", "source": SOURCE_NAMES.get(source, "?"), "sourceCode": source, "flags": flags, "oneway": bool(flags & FLAG_ONEWAY), "primaFacie": bool(flags & FLAG_PRIMA_FACIE), "osmDerived": source in OSM_DERIVED, "layer": (flags >> LAYER_SHIFT) - 3} def name(self, index: int): if index == NO_NAME: return "" if index >= self.name_count: return None start = self.nb_off + self.name_offsets[index] end = self.nb_off + self.name_offsets[index + 1] return self.buf[start:end].decode("utf-8") def record_points(self, rec: dict) -> list: base = rec["first"] * 2 return [(self.points[base + 2 * k], self.points[base + 2 * k + 1]) for k in range(rec["count"])] def to_q(self, lat: float, lon: float): return ((lat * 1e6 - self.origin_lat_e6) / self.quantum_e6, (lon * 1e6 - self.origin_lon_e6) / self.quantum_e6) def _distance(self, rec: dict, qlat: float, qlon: float): """(meters, segment bearing in degrees) from a Q-space point to a record.""" my = M_PER_DEG_LAT / self.q_per_deg lat = (self.origin_lat_e6 + qlat * self.quantum_e6) / 1e6 mx = 111_320.0 * math.cos(math.radians(lat)) / self.q_per_deg pts = self.record_points(rec) best = (float("inf"), 0.0) for (la0, lo0), (la1, lo1) in zip(pts, pts[1:]): ax, ay = (lo0 - qlon) * mx, (la0 - qlat) * my bx, by = (lo1 - qlon) * mx, (la1 - qlat) * my dx, dy = bx - ax, by - ay length2 = dx * dx + dy * dy t = 0.0 if length2 == 0 else max(0.0, min(1.0, -(ax * dx + ay * dy) / length2)) d = math.hypot(ax + t * dx, ay + t * dy) if d < best[0]: best = (d, math.degrees(math.atan2(dx, dy)) % 360.0) return best def candidates(self, lat: float, lon: float, radius_m: float, brute_force: bool = False) -> list: radius_m = min(radius_m, float(self.margin_m)) qlat, qlon = self.to_q(lat, lon) if brute_force: indices = range(self.count) else: if not (0 <= qlat < self.rows * self.cell_q and 0 <= qlon < self.cols * self.cell_q): return [] cell = int(qlat // self.cell_q) * self.cols + int(qlon // self.cell_q) indices = self.cell_refs[self.cell_start[cell]:self.cell_start[cell + 1]] found = [] for i in indices: rec = self.record(i) d, bearing = self._distance(rec, qlat, qlon) if d <= radius_m: found.append((d, bearing, rec)) found.sort(key=lambda item: (item[0], item[2]["index"])) return found def lookup(self, lat: float, lon: float, course=None, radius_m: float = 20.0, max_angle: float = 35.0): """Nearest record within radius whose bearing agrees with the course (either direction unless oneway). Mirrors the prototype's reference rule.""" for d, bearing, rec in self.candidates(lat, lon, radius_m): forward = True if course is not None: diff = abs((course - bearing + 180.0) % 360.0 - 180.0) if diff <= max_angle: forward = True elif diff >= 180.0 - max_angle and not rec["oneway"]: forward = False else: continue speed = rec["fwd"] if forward else rec["bwd"] return dict(rec, distance_m=round(d, 2), bearing=round(bearing, 2), speed_mph=speed, forward=forward) return None def verify(path: str, samples: int = 1500, seed: int = 7) -> dict: """Structure, per-record and index checks; the index is compared with an independent coarse-grid brute-force search at random points near roads.""" store = StreetLimitFile(path) problems = [] if (store.origin_lat_e6, store.origin_lon_e6, store.quantum_e6, store.cell_q) != \ (ORIGIN_LAT_E6, ORIGIN_LON_E6, QUANTUM_E6, CELL_Q): problems.append("origin / quantum / cell size differ from the app constants") if store.cols * store.cell_q < LON_SPAN_Q or store.rows * store.cell_q < LAT_SPAN_Q: problems.append("grid does not cover the bounding box") for i in range(len(store.cell_start) - 1): if store.cell_start[i] > store.cell_start[i + 1]: problems.append(f"cellStart not monotonic at {i}") break refs = store.cell_refs[store.cell_start[i]:store.cell_start[i + 1]] if any(refs[k] >= refs[k + 1] for k in range(len(refs) - 1)) or any(r >= store.count for r in refs): problems.append(f"cell {i} refs not ascending or out of range") break if store.cell_start[-1] != store.ref_count: problems.append("last cellStart != cellRefCount") for k in range(store.name_count): if store.name_offsets[k] > store.name_offsets[k + 1] or store.name_offsets[k + 1] > store.nb_len: problems.append(f"name offsets broken at {k}") break try: if store.name(k).count("\t") > 1: problems.append(f"name {k} has more than one TAB") except UnicodeDecodeError: problems.append(f"name {k} is not UTF-8") mask = 0 coarse = collections.defaultdict(list) coarse_q, pad = 2000, 40 for i in range(store.count): first, count, name, fwd, bwd, class_source, flags = store.raw_record(i) cls, source = class_source & 15, class_source >> 4 mask |= 1 << source if count < 2 or first + count > store.point_count: problems.append(f"record {i}: point range") if name != NO_NAME and name >= store.name_count: problems.append(f"record {i}: name index") if cls > 12 or source > MAX_SOURCE: problems.append(f"record {i}: class {cls} / source {source}") if bool(flags & FLAG_PRIMA_FACIE) != (source == SRC_LADOT_CVC): problems.append(f"record {i}: prima facie flag vs source {source}") for mph in (fwd, bwd): if mph and not valid_mph(mph): problems.append(f"record {i}: speed {mph}") if (source == SRC_NONE) != (fwd == 0 and bwd == 0): problems.append(f"record {i}: source {source} with speeds {fwd}/{bwd}") if len(problems) > 20: break rec = {"first": first, "count": count} pts = store.record_points(rec) las = [p[0] for p in pts] los = [p[1] for p in pts] for gy in range((min(las) - pad) // coarse_q, (max(las) + pad) // coarse_q + 1): for gx in range((min(los) - pad) // coarse_q, (max(los) + pad) // coarse_q + 1): coarse[(gy, gx)].append(i) if mask != store.sources_mask: problems.append(f"sourcesMask {store.sources_mask:#x} != records {mask:#x}") rng = random.Random(seed) checked = mismatches = with_hits = 0 radius = float(store.margin_m) - 1.0 for _ in range(samples): rec = store.record(rng.randrange(store.count)) pts = store.record_points(rec) la, lo = pts[rng.randrange(len(pts))] lat = (store.origin_lat_e6 + la * store.quantum_e6) / 1e6 + rng.uniform(-0.0003, 0.0003) lon = (store.origin_lon_e6 + lo * store.quantum_e6) / 1e6 + rng.uniform(-0.0003, 0.0003) qlat, qlon = store.to_q(lat, lon) if not (0 <= qlat <= LAT_SPAN_Q and 0 <= qlon <= LON_SPAN_Q): continue fast = {c[2]["index"] for c in store.candidates(lat, lon, radius)} slow = set() for i in coarse.get((int(qlat // coarse_q), int(qlon // coarse_q)), ()): if store._distance(store.record(i), qlat, qlon)[0] <= radius: slow.add(i) checked += 1 with_hits += bool(slow) mismatches += fast != slow if mismatches: problems.append(f"{mismatches} index mismatches") return {"ok": not problems, "problems": problems, "crc": "ok", "records": store.count, "formatVersion": FORMAT_VERSION, "indexChecks": checked, "checksWithCandidates": with_hits, "indexMismatches": mismatches, "radiusM": radius} # --------------------------------------------------------------------------- stats def coverage_report(records: list) -> dict: by_class = collections.defaultdict(collections.Counter) by_area = {name: collections.defaultdict(collections.Counter) for name in SUB_AREAS} for key, coords in records: cls, fwd, bwd, source, flags, label = key length = polyline_length_m(coords) name = ROAD_CLASSES[cls] bucket = SOURCE_NAMES[source] by_class[name][bucket] += length by_class[name]["total"] += length mid = coords[len(coords) // 2] group = "major" if name in MAJOR_CLASSES else ("links" if name in LINK_CLASSES else "local") for area, (w, s, e, n) in SUB_AREAS.items(): if w <= mid[1] <= e and s <= mid[0] <= n: by_area[area][group][bucket] += length by_area[area][group]["total"] += length def summarize(counter): total = counter.get("total", 0.0) with_limit = total - counter.get("none", 0.0) out = {"km_total": km(total), "km_with_limit": km(with_limit), "pct_with_limit": round(100 * with_limit / total, 1) if total else None} for src in SOURCE_NAMES.values(): if counter.get(src): out[f"km_{src}"] = km(counter[src]) return out def combine(names): total = collections.Counter() for name in names: total.update(by_class.get(name, {})) return total return { "by_class": {name: summarize(by_class[name]) for name in ROAD_CLASSES if name in by_class}, "major_roads_motorway_to_tertiary": summarize(combine(MAJOR_CLASSES)), "arterials_primary_secondary_tertiary": summarize(combine(("primary", "secondary", "tertiary"))), "links": summarize(combine(LINK_CLASSES)), "all_roads": summarize(combine(ROAD_CLASSES)), "sub_areas_major_roads": {area: summarize(groups["major"]) for area, groups in by_area.items()}, "sub_areas_local_roads": {area: summarize(groups["local"]) for area, groups in by_area.items()}, } def counts_by_source_and_limit(records: list) -> tuple: by_source = collections.OrderedDict() by_limit = collections.Counter() by_limit_km = collections.Counter() by_source_limit = collections.defaultdict(collections.Counter) for key, coords in records: cls, fwd, bwd, source, flags, label = key name = SOURCE_NAMES[source] entry = by_source.setdefault(name, {"code": source, "records": 0, "km": 0.0}) length = polyline_length_m(coords) entry["records"] += 1 entry["km"] += length for mph in sorted({fwd, bwd}): by_limit[str(mph)] += 1 by_limit_km[str(mph)] += length by_source_limit[name][str(mph)] += length for entry in by_source.values(): entry["km"] = km(entry["km"]) order = sorted(by_limit, key=int) limits = {mph: {"records": by_limit[mph], "km": km(by_limit_km[mph])} for mph in order} detail = {src: {mph: km(v) for mph, v in sorted(c.items(), key=lambda kv: int(kv[0]))} for src, c in by_source_limit.items()} return dict(sorted(by_source.items(), key=lambda kv: kv[1]["code"])), limits, detail # --------------------------------------------------------------------------- finalize and gates def finalize_records(records: list) -> tuple: """Quantize every record exactly as the writer will, split it where the quantized line turns back on itself, and drop what the file must not carry: fewer than 2 points, shorter than MIN_RECORD_M, or an exact (two-way: also reversed) duplicate. The records returned carry the dequantized coordinates, so every count and length in the manifest and the stats describes the file as written.""" st = collections.Counter() seen = set() out = [] my, mx = M_PER_DEG_LAT / Q_PER_DEG, M_PER_DEG_LON / Q_PER_DEG for key, coords in records: q = [] for lat, lon in coords: point = quantize(lat, lon) if not q or q[-1] != point: q.append(point) pieces = [q] if len(q) >= 3: cuts = reversal_vertices([(b * mx, a * my) for a, b in q]) if cuts: st["split_at_reversal"] += len(cuts) bounds = [0] + cuts + [len(q) - 1] pieces = [q[i:j + 1] for i, j in zip(bounds, bounds[1:])] for piece in pieces: if len(piece) < 2: st["dropped_under_2_points"] += 1 continue length = sum(math.hypot((b[0] - a[0]) * my, (b[1] - a[1]) * mx) for a, b in zip(piece, piece[1:])) if length < MIN_RECORD_M: st[f"dropped_shorter_than_{MIN_RECORD_M:g}m"] += 1 continue sig = tuple(piece) if (sig, key) in seen or (not key[4] & FLAG_ONEWAY and (sig[::-1], reversed_key(key)) in seen): st["dropped_duplicate"] += 1 continue seen.add((sig, key)) out.append((key, [(SOUTH + a / Q_PER_DEG, WEST + b / Q_PER_DEG) for a, b in piece])) return out, dict(st) def _family_of(cls: int) -> int: return FAM_FREEWAY if cls in (C_MOTORWAY, C_MOTORWAY_LINK) else FAM_SURFACE def gate_official_overlap(records: list) -> dict: """Official records that run on top of each other with different limits.""" index = LineIndex(cell=50.0, pad=OVERLAP_GATE_M) official = [i for i, (key, _) in enumerate(records) if key[3] in OFFICIAL_SOURCES and max(key[1], key[2])] value = {i: max(records[i][0][1], records[i][0][2]) for i in official} for i in official: index.add(records[i][1], i, _family_of(records[i][0][0])) pairs = collections.Counter() where = {} for i in official: key, coords = records[i] def other(pl, i=i): return pl != i and value[pl] != value[i] def tally(pl, meters, distance, x, y, i=i): pair = (min(i, pl), max(i, pl)) pairs[pair] += meters where.setdefault(pair, (x, y)) carve(coords, index, OVERLAP_GATE_M, _family_of(key[0]), tally, accept=other) findings = [] for (a, b), meters in pairs.most_common(): meters /= 2.0 # every overlap is seen from both records if meters < OVERLAP_GATE_RUN_M: continue x, y = where[(a, b)] findings.append({"family": "freeway" if _family_of(records[a][0][0]) == FAM_FREEWAY else "surface", "m": round(meters), "at": [round(SOUTH + y / M_PER_DEG_LAT, 6), round(WEST + x / M_PER_DEG_LON, 6)], "records": [{"name": records[r][0][5].replace("\t", " | "), "mph": value[r], "source": SOURCE_NAMES[records[r][0][3]]} for r in (a, b)]}) return {"ok": not findings, "rule": f"official records within {OVERLAP_GATE_M:g} m, aligned within " f"{CARVE_ANGLE:g} deg, different limits, over {OVERLAP_GATE_RUN_M:g} m or more", "findings": len(findings), "examples": findings[:40], "belowThresholdM": round(sum(m / 2.0 for m in pairs.values() if m / 2.0 < OVERLAP_GATE_RUN_M))} def gate_ladot_outside(records: list, city: CityBoundary) -> dict: findings = [] checked = 0.0 for key, coords in records: if key[3] not in LADOT_SOURCES: continue run = worst = 0.0 first = None for x, y, _, spacing in iter_samples(coords, CARVE_STEP_M): checked += spacing out = city.outside_m(x, y) if out > CITY_GATE_M: run += spacing worst = max(worst, out) first = first or (x, y) elif run: break if run >= 2 * CARVE_STEP_M: x, y = first findings.append({"name": key[5], "mph": max(key[1], key[2]), "source": SOURCE_NAMES[key[3]], "outsideRunM": round(run), "maxOutsideM": round(worst), "at": [round(SOUTH + y / M_PER_DEG_LAT, 6), round(WEST + x / M_PER_DEG_LON, 6)]}) return {"ok": not findings, "rule": f"no LADOT record more than {CITY_GATE_M:g} m outside the City of Los " f"Angeles boundary for {2 * CARVE_STEP_M:g} m or more", "kmChecked": km(checked), "findings": len(findings), "examples": findings[:40]} def _nearest_on(records: list, candidates: list, lat: float, lon: float, radius: float, boxes: dict): """(distance m, record index, (lat, lon) of the nearest point) over candidate records.""" px, py = to_xy(lat, lon) dlat, dlon = radius / M_PER_DEG_LAT, radius / M_PER_DEG_LON best = None for r in candidates: box = boxes.get(r) if box is None: coords = records[r][1] box = boxes[r] = (min(c[0] for c in coords), min(c[1] for c in coords), max(c[0] for c in coords), max(c[1] for c in coords)) if lat < box[0] - dlat or lat > box[2] + dlat or lon < box[1] - dlon or lon > box[3] + dlon: continue coords = records[r][1] for (la0, lo0), (la1, lo1) in zip(coords, coords[1:]): ax, ay = to_xy(la0, lo0) bx, by = to_xy(la1, lo1) dx, dy = bx - ax, by - ay l2 = dx * dx + dy * dy t = 0.0 if l2 == 0 else max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / l2)) d = math.hypot(ax + t * dx - px, ay + t * dy - py) if d <= radius and (best is None or d < best[0]): best = (d, r, (la0 + t * (la1 - la0), lo0 + t * (lo1 - lo0))) return best def load_waivers(path) -> dict: if not path or not os.path.exists(path): return {} with open(path, encoding="utf-8") as f: doc = json.load(f) out = {} for entry in doc.get("waivers", []): for field in ("camera", "datasetMph", "reason", "evidence", "reviewBy"): assert field in entry, (field, entry) out[entry["camera"]] = entry return out def gate_cameras(records: list, cameras_path, waivers: dict, today: str) -> dict: """Every camera with a posted limit is compared with the dataset along its segment. The camera's end points are snapped to the named street first (some are a few hundred meters off), then points at fixed fractions of the snapped segment are snapped to the street again, so a curving street is followed. A camera disagrees when any sampled point carries a different limit. Freeway records are never used. Only cameras whose known end points lie inside the dataset's bounding box are checked (cameras.json is statewide; SF, San Jose, Oakland and Malibu sit outside the box and are only counted).""" if not cameras_path or not os.path.exists(cameras_path): return {"ok": True, "skipped": f"no cameras file at {cameras_path}"} with open(cameras_path, encoding="utf-8") as f: cameras = json.load(f).get("cameras", []) def in_box(point) -> bool: return bool(point) and SOUTH <= point["lat"] <= NORTH and WEST <= point["lon"] <= EAST outside = [cam for cam in cameras if cam.get("speedLimitMph") and not (in_box(cam.get("start")) and (not cam.get("end") or in_box(cam.get("end"))))] outside_ids = {cam["id"] for cam in outside} cameras = [cam for cam in cameras if cam["id"] not in outside_ids] by_name = collections.defaultdict(list) for i, (key, _) in enumerate(records): if key[0] not in (C_MOTORWAY, C_MOTORWAY_LINK) and key[5]: by_name[name_key(key[5].replace(",", " "))].append(i) # "Martin Luther King, Jr Boulevard" results = [] counts = collections.Counter() boxes: dict = {} for cam in cameras: mph = cam.get("speedLimitMph") if not mph: continue cands = by_name.get(name_key(cam["street"].replace(",", " ")), []) # A Glendale camera is never on a LADOT (LA City) record, an LA camera never on a # Glendale one: keeps a camera from snapping to a same-named street across the line. foreign = LADOT_SOURCES if cam.get("city") != "Los Angeles" else (SRC_GLENDALE,) cands = [r for r in cands if records[r][0][3] not in foreign] entry = {"camera": cam["id"], "street": cam["street"], "cameraMph": mph, "quality": cam.get("coordinateQuality")} start = _nearest_on(records, cands, cam["start"]["lat"], cam["start"]["lon"], 400.0, boxes) end = None if cam.get("end") and cam.get("coordinateQuality") != "one_end_only": end = _nearest_on(records, cands, cam["end"]["lat"], cam["end"]["lon"], 400.0, boxes) if start is None: entry["status"] = "street_not_found" counts["street_not_found"] += 1 results.append(entry) continue entry["snapM"] = [round(start[0])] + ([round(end[0])] if end else []) one_end = end is None if one_end: # Only the "from" end is known, so the direction of the segment is not: every # named-street value within 120 m of it is collected, and the camera agrees # when its limit is among them. points = [] for r in cands: hit = _nearest_on(records, [r], start[2][0], start[2][1], 120.0, boxes) if hit is not None: points.append(hit) else: fractions = (0.25, 0.375, 0.5, 0.625, 0.75) if cam.get("coordinateQuality") == "approximate" \ else (0.15, 0.3, 0.45, 0.6, 0.75, 0.85) (la0, lo0), (la1, lo1) = start[2], end[2] points = [_nearest_on(records, cands, la0 + t * (la1 - la0), lo0 + t * (lo1 - lo0), 150.0, boxes) for t in fractions] points = [p for p in points if p is not None] values = [] for d, r, at in points: key = records[r][0] values.append({"mph": max(key[1], key[2]), "source": SOURCE_NAMES[key[3]], "at": [round(at[0], 6), round(at[1], 6)]}) known = [v for v in values if v["mph"]] wrong = [v for v in known if v["mph"] != mph] if one_end and any(v["mph"] == mph for v in known): wrong = [] entry["values"] = sorted({v["mph"] for v in values}) entry["sources"] = sorted({v["source"] for v in values}) if wrong: entry["status"] = "mismatch" entry["mismatches"] = wrong elif known: entry["status"] = "agree" if len(known) == len(values) else "agree_partial_no_data" else: entry["status"] = "no_data" waiver = waivers.get(cam["id"]) if entry["status"] == "mismatch" and waiver: if sorted(waiver["datasetMph"]) == sorted({v["mph"] for v in wrong}): entry["status"] = "mismatch_waived" entry["waiver"] = {"reason": waiver["reason"], "reviewBy": waiver["reviewBy"], "expired": waiver["reviewBy"] < today} else: entry["waiverStale"] = f"waiver expects {waiver['datasetMph']}, dataset has " \ f"{sorted({v['mph'] for v in wrong})}" counts[entry["status"]] += 1 results.append(entry) unused = sorted(set(waivers) - {r["camera"] for r in results if r.get("status") == "mismatch_waived"}) matched = counts["agree"] + counts["agree_partial_no_data"] + counts["mismatch"] + counts["mismatch_waived"] return {"ok": counts["mismatch"] == 0 and not unused, "counts": dict(counts), "agreementPct": round(100.0 * (counts["agree"] + counts["agree_partial_no_data"]) / matched, 1) if matched else None, "outsideAreaNotChecked": len(outside), "unusedWaivers": unused, "expiredWaivers": sorted(r["camera"] for r in results if r.get("waiver", {}).get("expired")), "results": results} def review_cvc_vs_osm(cvc_pieces: list, chains: list) -> dict: """LADOT CVC 25 pieces along which OSM has a higher value (for LADOT / manual review; nothing in the file changes).""" index = LineIndex(cell=50.0, pad=OSM_CARVE_M) for i, (coords, _, _, _) in enumerate(cvc_pieces): index.add(coords, i, FAM_SURFACE) tally = collections.defaultdict(collections.Counter) for key, coords in chains: cls, fwd, bwd, source, flags, label = key value = max(fwd, bwd) if source != SRC_OSM or value <= 25: continue carve(coords, index, OSM_CARVE_M, FAM_SURFACE, lambda pl, m, d, x, y, v=value: tally[pl].update({v: m})) by_type = collections.defaultdict(collections.Counter) listed = [] for pl, values in tally.items(): coords, sect, st_type, name = cvc_pieces[pl] for v, meters in values.items(): by_type[st_type or "?"][v] += meters top, meters = values.most_common(1)[0] if top >= 35 and meters >= 50: mid = coords[len(coords) // 2] listed.append({"sectId": sect, "stType": st_type, "name": name, "osmMph": top, "m": round(meters), "at": [round(mid[0], 6), round(mid[1], 6)]}) listed.sort(key=lambda e: (e["stType"] != "SE", -e["m"])) return {"kmByStTypeAndOsmValue": {t: {str(v): km(m) for v, m in sorted(c.items())} for t, c in sorted(by_type.items())}, "selectSegmentsOsm35OrMore": [e for e in listed if e["stType"] == "SE"], "otherSegmentsOsm35OrMore": [e for e in listed if e["stType"] != "SE"][:200]} def review_one_direction(records: list) -> list: out = [] for key, coords in records: cls, fwd, bwd, source, flags, label = key if not flags & FLAG_ONEWAY and (fwd == 0) != (bwd == 0): mid = coords[len(coords) // 2] out.append({"name": label, "fwd": fwd, "bwd": bwd, "source": SOURCE_NAMES[source], "at": [round(mid[0], 6), round(mid[1], 6)], "why": "OSM tags maxspeed for one direction only and no maxspeed; the other " "direction is left unknown (no invented value)"}) return out # --------------------------------------------------------------------------- build def resolve_pbfs(args, cache: str) -> list: if args.osm_pbf: return [os.path.abspath(p) for p in args.osm_pbf] urls = args.pbf_url or DEFAULT_PBF_URLS return [download_file(url, os.path.join(cache, "osm", url.rsplit("/", 1)[-1]), args.refresh) for url in urls] LADOT_LAYER_TITLES = { "ladot_lamc35": "SpeedLimits_LAMC_35mph", "ladot_map": "City_of_Los_Angeles_Speed_Limit_Map", "ladot_ab43": "SpeedLimitChange011 (AB 43 speed limit changes)", "ladot_cvc": "SpeedLimits_CVC", } SOURCE_NOTES = { "ladot_cvc": {"itemUrl": "the layer is published on LADOT's own ArcGIS Server only; no ArcGIS Online item " "was found for it", "dataLastEditUtc": "the server does not expose editingInfo for this layer"}, "ladot_map": {"dataLastEditUtc": "the server does not expose editingInfo for this layer; the base is the " "2018 LAMC 80.81 list, patched by AB 43 and LAMC_35"}, "la_city_boundary": {"itemUrl": "queried directly from the County's public MapServer (layer 19, City " "Boundaries); not recorded as a hub item", "dataLastEditUtc": "a MapServer layer: it exposes no editingInfo"}, "glendale_cams": {"itemUrl": "queried directly from the County eGIS FeatureServer (CAMS_ADDRESS_LINES); " "not recorded as a hub item"}, } def build(args) -> None: started = time.time() today = dt.datetime.now(dt.timezone.utc).strftime("%Y-%m-%d") out_dir = os.path.abspath(args.out) cache = os.path.abspath(args.cache) arcgis_cache = os.path.join(cache, "arcgis") os.makedirs(out_dir, exist_ok=True) wanted = [sid for sid in ARCGIS_SOURCES if not (args.no_county and sid == "la_county")] metas = {sid: fetch_arcgis(sid, arcgis_cache, args.refresh, args.pause) for sid in wanted} retrieved = min(m["retrievedUtc"] for m in metas.values()) data_date = int(retrieved[:10].replace("-", "")) city = CityBoundary(load_features(arcgis_cache, "la_city_boundary")) county = None if args.no_county else county_lines(arcgis_cache) county_index = county_posted_index(county) if county else None builder = Builder() cvc_pieces = build_ladot(builder, arcgis_cache, city, county_index) log(f"LADOT: {len(builder.records)} records accepted") overrides = load_overrides(args.overrides) if overrides: builder.stats["overrides"] = apply_overrides(builder, arcgis_cache, overrides) log(f"overrides: {len(overrides)} entries applied") if county: build_county(builder, county) log(f"County: {len(builder.records)} records accepted so far") surface, cap = build_caltrans(builder, arcgis_cache, args.caltrans_surface) builder.stats["caltrans"]["surface_qa_vs_official"] = caltrans_surface_qa(builder, surface) builder.commit() builder.stats["caltrans"]["freeway_next_to_official_surface"] = caltrans_surface_adjacency(builder) log(f"Caltrans: {len(builder.records)} records accepted so far") official_count = len(builder.records) osm_time = None osm_issues: collections.Counter = collections.Counter() pbf_paths = [] review = {} if not args.no_osm: pbf_paths = resolve_pbfs(args, cache) chains, osm_time, osm_issues = osm_chains(cache, pbf_paths) if not args.strict_precedence: builder.stats["cvc_yields_to_lower_osm"] = cvc_yields_to_lower_osm(builder, chains) log(f"CVC yield rule: {builder.stats['cvc_yields_to_lower_osm']}") builder.stats["stale_base_yields_to_higher_osm"] = stale_base_yields_to_higher_osm(builder, chains) log(f"stale base rule: {builder.stats['stale_base_yields_to_higher_osm']}") review.update(build_osm(builder, chains, cap)) review["cvcVsOsm"] = review_cvc_vs_osm(cvc_pieces, chains) log(f"OSM: {len(builder.records) - official_count} records added") osm_date = int(osm_time.strftime("%Y%m%d")) if osm_time else 0 # Chain identical neighbours, simplify, split long lines. pieces, duplicates = dedupe_pieces(builder.records) chains = merge_chains(pieces) records = [] points_before = sum(len(c) for _, c in builder.records) for key, coords in chains: for chunk in split_by_length(simplify(coords, args.simplify_m), args.max_chunk_m): if len(chunk) >= 2: records.append((key, chunk)) log(f"{len(builder.records)} pieces ({duplicates} duplicates dropped) -> {len(chains)} chains -> " f"{len(records)} records ({points_before} points before simplification)") suppressed = [] if args.suppress: records, suppressed = apply_suppressions(records, args.suppress) records, builder.stats["official_conflicts_resolved"] = resolve_official_conflicts(records) log(f"official overlaps with different limits cut: {builder.stats['official_conflicts_resolved']}") records, hygiene = finalize_records(records) hygiene["duplicate_pieces_dropped_before_chaining"] = duplicates hygiene["retracing_records_after_finalize"] = sum( 1 for _, c in records if reversal_vertices([to_xy(*p) for p in c])) review["oneDirectionOnlyRecords"] = review_one_direction(records) log(f"finalized: {len(records)} records {hygiene}") # Gates. gates = { "officialOverlap": gate_official_overlap(records), "ladotOutsideCity": gate_ladot_outside(records, city), "cameras": gate_cameras(records, args.cameras, load_waivers(args.camera_waivers), today), } release = [] if any(not g["ok"] for g in gates.values()): release.append("a build gate failed") if not args.publish_url: release.append("no --publish-url: the NOTICE still has the placeholder (ODbL 4.6/4.7)") expired = [e["id"] for e in overrides if e["reviewBy"] < today] expired += gates["cameras"].get("expiredWaivers", []) if expired: release.append(f"past their review date: {expired}") gates["release"] = {"ok": not release, "blockers": release, "required": bool(args.release)} for name, gate in gates.items(): log(f"gate {name}: {'ok' if gate['ok'] else 'FAILED'}" + (f" ({gate.get('findings')} findings)" if "findings" in gate else "") + (f" {gate.get('counts')}" if "counts" in gate else "") + (f" {gate.get('blockers')}" if gate.get("blockers") else "")) failed = [name for name, gate in gates.items() if not gate["ok"] and (name != "release" or args.release)] if failed and not args.allow_gate_failures: path = os.path.join(out_dir, "street_limits.gates.FAILED.json") with open(path, "w", encoding="utf-8") as f: json.dump(gates, f, indent=1) raise RuntimeError(f"build gates failed: {failed}; the previous output is unchanged; details in {path}") out_path = os.path.join(out_dir, "street_limits.bin") layout = write_resource(records, out_path, data_date, osm_date, args.index_margin_m) assert layout["records"] == len(records), (layout["records"], len(records)) with open(out_path, "rb") as f: blob = f.read() zlib_size = len(zlib.compress(blob, 9)) sha256 = hashlib.sha256(blob).hexdigest() check = verify(out_path) if not check["ok"]: raise RuntimeError(f"verification failed: {check['problems']}") coverage = coverage_report(records) by_source, by_limit, by_source_limit = counts_by_source_and_limit(records) assert sum(v["records"] for v in by_source.values()) == layout["records"] generated = dt.datetime.now(dt.timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") sources = [] for sid, meta in metas.items(): entry = { "id": sid, "layerName": meta.get("layerName"), "url": meta["layerUrl"], "itemUrl": ARCGIS_SOURCES[sid].get("item"), "query": {"where": meta["where"], "outFields": meta["outFields"], "bboxFilter": meta["bboxFilter"]}, "featureCount": meta["featureCount"], "retrievedUtc": meta["retrievedUtc"], "dataLastEditUtc": meta.get("dataLastEditUtc"), } notes = {k: v for k, v in SOURCE_NOTES.get(sid, {}).items() if entry.get(k) is None} if notes: entry["notes"] = notes sources.append(entry) if overrides: sources.append({"id": "overrides", "file": os.path.basename(args.overrides), "entries": [{k: e.get(k) for k in ("id", "source", "street", "mph", "citation", "url", "addedOn", "reviewBy")} for e in overrides]}) if not args.no_osm: sources.append({ "id": "osm", "layerName": "OpenStreetMap highway ways (maxspeed)", "url": ", ".join(args.pbf_url or DEFAULT_PBF_URLS) if not args.osm_pbf else "local PBF extract(s)", "files": [os.path.basename(p) for p in pbf_paths], "snapshotUtc": osm_time.strftime("%Y-%m-%dT%H:%M:%SZ") if osm_time else None, }) source_codes = [] for code, info in SOURCE_INFO.items(): if not layout["sourcesMask"] & (1 << code): continue source_codes.append(dict({"code": code, "name": SOURCE_NAMES[code], "osmDerived": code in OSM_DERIVED}, **info)) attribution = [] mask = layout["sourcesMask"] if any(mask & (1 << c) for c in LADOT_SOURCES): attribution.append(SOURCE_INFO[SRC_LADOT_LAMC35]["attribution"]) if mask & (1 << SRC_GLENDALE): attribution.append(SOURCE_INFO[SRC_GLENDALE]["attribution"]) county_parts = (["Department of Public Works, Speed Network"] if mask & (1 << SRC_COUNTY) else []) + \ (["eGIS road centerlines"] if mask & (1 << SRC_GLENDALE) else []) if county_parts: attribution.append("Contains data from County of Los Angeles, " + "; ".join(county_parts) + ". Not endorsed by the County.") if mask & (1 << SRC_CALTRANS): attribution.append(SOURCE_INFO[SRC_CALTRANS]["attribution"]) if any(mask & (1 << c) for c in OSM_DERIVED): attribution.append(SOURCE_INFO[SRC_OSM]["attribution"]) yielded = builder.stats.get("cvc_yields_to_lower_osm", {}).get("cvc_km_yielded_by_osm_value", {}) stale_yield = {k: v for k, v in builder.stats.get("stale_base_yields_to_higher_osm", {}).get("km", {}).items() if k != "kept_base_value"} manifest = { "dataset": "street_limits", "version": f"{FORMAT_VERSION}.{data_date}", "formatVersion": FORMAT_VERSION, "generatedUtc": generated, "dataDate": data_date, "osmDataDate": osm_date or None, "file": {"name": "street_limits.bin", "bytes": layout["bytes"]["total"], "zlibBytes": zlib_size, "sha256": sha256, "crc32": f"{StreetLimitFile(out_path).crc:08x}"}, "bbox": {"south": SOUTH, "west": WEST, "north": NORTH, "east": EAST}, "grid": {"originLatE6": ORIGIN_LAT_E6, "originLonE6": ORIGIN_LON_E6, "quantumE6": QUANTUM_E6, "cellSizeQ": CELL_Q, "cols": layout["gridCols"], "rows": layout["gridRows"], "indexMarginM": args.index_margin_m}, "attribution": attribution, "disclaimer": "Speed limits may be missing, outdated or wrong. Posted signs always govern.", "sourceCodes": source_codes, "sources": sources, "counts": {"records": layout["records"], "points": layout["points"], "names": layout["names"], "osmDerivedRecords": sum(v["records"] for k, v in by_source.items() if v["code"] in OSM_DERIVED), "bySource": by_source, "byLimitMph": by_limit}, "osmInfluenceOnOfficialRecords": { "cvcKmRemovedWhereOsmHasALowerValue": yielded, "staleBaseKmRemovedWhereOsmHasAHigherValue": stale_yield, "note": ("Official records never take geometry or values from OpenStreetMap. In this build OSM " "decided where two official values are left out: LADOT CVC 25 parts where OSM has a value " "below 25 mph, and 2018 base-map parts at 35 mph or less with no LAMC_35 match where OSM has " "a higher value. Dropping source codes 0 and 1 from this file therefore leaves those streets " "with no limit. A file free of OSM influence is built with --no-osm.") if (yielded or stale_yield) else "None: this build does not use OSM to shape official records.", }, "gates": {name: {k: v for k, v in gate.items() if k in ("ok", "rule", "findings", "counts", "agreementPct", "blockers", "required", "skipped")} for name, gate in gates.items()}, "releaseReady": gates["release"]["ok"], "notIncluded": ["school-zone and time-of-day limits", "advisory speeds", "invented defaults (the only 25 mph defaults are LADOT's CVC layer)"] + ([] if args.caltrans_surface else ["Caltrans surface-highway values (failed QA)"]) + (["LA County Speed Network (--no-county)"] if args.no_county else []), } with open(os.path.join(out_dir, "street_limits.manifest.json"), "w", encoding="utf-8") as f: json.dump(manifest, f, indent=1) f.write("\n") stats = { "generatedUtc": generated, "output": out_path, "options": {"noOsm": args.no_osm, "noCounty": args.no_county, "caltransSurface": args.caltrans_surface, "suppress": args.suppress, "overrides": args.overrides, "cameras": args.cameras, "cameraWaivers": args.camera_waivers, "strictPrecedence": args.strict_precedence, "publishUrl": args.publish_url, "release": args.release, "simplifyM": args.simplify_m, "maxChunkM": args.max_chunk_m, "indexMarginM": args.index_margin_m}, "conflation": {"overrideMatchM": OVERRIDE_MATCH_M, "lamcNameRetryM": LAMC_NAME_RETRY_M, "officialCarveM": OFFICIAL_CARVE_M, "frontageCarveM": FRONTAGE_CARVE_M, "countyEdgeM": COUNTY_EDGE_M, "cityKeepM": CITY_KEEP_M, "osmCarveM": OSM_CARVE_M, "osmFreewayCarveM": OSM_FREEWAY_CARVE_M, "carveAngleDeg": CARVE_ANGLE, "carveStepM": CARVE_STEP_M, "minKeepM": MIN_KEEP_M, "minRunM": MIN_RUN_M, "minRecordM": MIN_RECORD_M}, "layout": layout, "zlibCompressedBytes": zlib_size, "verify": check, "bySource": by_source, "byLimitMph": by_limit, "kmBySourceAndLimit": by_source_limit, "coverage": coverage, "layers": builder.stats, "suppressed": suppressed, "geometryHygiene": hygiene, "gates": gates, "osmTagIssues": dict(osm_issues.most_common(25)), "seconds": round(time.time() - started, 1), } with open(os.path.join(out_dir, "street_limits.stats.json"), "w", encoding="utf-8") as f: json.dump(stats, f, indent=1) f.write("\n") with open(os.path.join(out_dir, "street_limits.review.json"), "w", encoding="utf-8") as f: json.dump(dict({"generatedUtc": generated, "purpose": "manual review lists; nothing here changes the file"}, **review), f, indent=1) f.write("\n") with open(os.path.join(out_dir, "street_limits.stats.txt"), "w", encoding="utf-8") as f: f.write(render_stats_text(manifest, stats, review)) write_notice(os.path.join(out_dir, "street_limits.NOTICE.txt"), manifest, metas, osm_time, args.publish_url, overrides) failed_path = os.path.join(out_dir, "street_limits.gates.FAILED.json") if os.path.exists(failed_path) and not failed: os.remove(failed_path) major = coverage["major_roads_motorway_to_tertiary"] log(f"wrote {out_path}: {layout['bytes']['total'] / 1e6:.2f} MB ({zlib_size / 1e6:.2f} MB zlib), " f"{layout['records']} records; major roads {major['pct_with_limit']}% of {major['km_total']} km " f"have a limit; {time.time() - started:.0f}s") def render_stats_text(manifest: dict, stats: dict, review: dict) -> str: lines = [f"street_limits.bin {manifest['version']} generated {manifest['generatedUtc']}", f"{manifest['file']['bytes']:,} bytes ({manifest['file']['zlibBytes']:,} zlib), " f"{manifest['counts']['records']:,} records, {manifest['counts']['points']:,} points, " f"{manifest['counts']['names']:,} names; verify ok={stats['verify']['ok']}, " f"{stats['verify']['indexMismatches']} index mismatches in {stats['verify']['indexChecks']} checks", f"data as of {manifest['dataDate']}; OSM snapshot {manifest['osmDataDate']}; " f"release ready: {manifest['releaseReady']}", "", "Gates:"] for name, gate in stats["gates"].items(): detail = {k: gate[k] for k in ("findings", "counts", "agreementPct", "blockers", "unusedWaivers") if gate.get(k)} lines.append(f" {name:<18} {'ok' if gate['ok'] else 'FAILED'} {json.dumps(detail)}") cams = stats["gates"]["cameras"].get("results", []) for cam in cams: if cam.get("status") not in ("agree",): lines.append(f" camera {cam['camera']}: {cam.get('status')} camera {cam['cameraMph']} dataset " f"{cam.get('values')} {cam.get('sources')}" + (f" (waived: {cam['waiver']['reason'][:90]})" if cam.get("waiver") else "")) lines += ["", "Records and km by source:"] for name, entry in stats["bySource"].items(): lines.append(f" {entry['code']} {name:<24} {entry['records']:>7,} records {entry['km']:>9,.1f} km") lines += ["", "km by limit (each direction value counted once per record):"] for mph, entry in stats["byLimitMph"].items(): lines.append(f" {mph:>3} mph {entry['records']:>7,} records {entry['km']:>9,.1f} km") cov = stats["coverage"] lines += ["", "Share of road length with a limit:"] for label, key in (("all roads", "all_roads"), ("major roads (motorway..tertiary)", "major_roads_motorway_to_tertiary"), ("arterials (primary..tertiary)", "arterials_primary_secondary_tertiary"), ("links / ramps", "links")): c = cov[key] lines.append(f" {label:<34} {c['pct_with_limit']}% of {c['km_total']:,} km") for name, c in cov["by_class"].items(): lines.append(f" {name:<16} {c['pct_with_limit']}% of {c['km_total']:,} km") lines.append(" sub-areas, major roads:") for name, c in cov["sub_areas_major_roads"].items(): lines.append(f" {name:<34} {c['pct_with_limit']}% of {c['km_total']:,} km") lines.append(" sub-areas, local roads:") for name, c in cov["sub_areas_local_roads"].items(): lines.append(f" {name:<34} {c['pct_with_limit']}% of {c['km_total']:,} km") lines += ["", "Geometry hygiene:", json.dumps(stats["geometryHygiene"])] lines += ["", "Layer details:", json.dumps(stats["layers"], indent=1)] if stats.get("suppressed"): lines += ["", "Suppressed OSM values:", json.dumps(stats["suppressed"], indent=1)] if review.get("caltransVsOsmFreeway"): lines += ["", "Caltrans vs OSM freeway disagreements (km; sections in street_limits.review.json):"] for entry in review["caltransVsOsmFreeway"]: if entry["km"] >= 1.0: lines.append(f" {entry['route']:<8} caltrans {entry['caltrans']} / osm {entry['osm']}: " f"{entry['km']} km, {entry['direction']}") if review.get("cvcVsOsm"): lines += ["", f"CVC 25 select (SE) segments where OSM says 35+: " f"{len(review['cvcVsOsm']['selectSegmentsOsm35OrMore'])} (list in street_limits.review.json)"] return "\n".join(lines) + "\n" def write_notice(path: str, manifest: dict, metas: dict, osm_time, publish_url, overrides) -> None: mask = sum(1 << s["code"] for s in manifest["sourceCodes"]) lines = ["street_limits.bin: data sources and credits", ""] if any(mask & (1 << c) for c in LADOT_SOURCES): lines += ["Speed limit data: City of Los Angeles Department of Transportation.", "Layers retrieved " + metas["ladot_map"]["retrievedUtc"][:10] + ":"] for sid, title in LADOT_LAYER_TITLES.items(): lines += [f" {title}: {ARCGIS_SOURCES[sid]['url']}"] if mask & (1 << SRC_LADOT_DOC): lines += [" LADOT Speed Safety System program documents (2026), per street_limits_overrides.json:"] lines += [f" {e['citation']}: {e.get('url', '')}" for e in overrides if e["source"] == "ladot_document"] lines += ["Stored as a normalized derived representation; the raw files are not redistributed.", ""] if mask & (1 << SRC_GLENDALE): lines += ["Speed limits: City of Glendale Municipal Code 10.24.020 (as amended by Ord. 6040, 2025-08-19),", "transcribed for the street segments listed in street_limits_overrides.json.", ""] county_parts = (["Department of Public Works, Speed Network"] if mask & (1 << SRC_COUNTY) else []) + \ (["eGIS road centerlines (CAMS)"] if mask & (1 << SRC_GLENDALE) else []) + ["eGIS city boundaries"] lines += ["Contains data from: County of Los Angeles, " + "; ".join(county_parts) + ".", "County of Los Angeles Enterprise GIS. Accessed " + retrieved_day(metas) + "."] lines += [ARCGIS_SOURCES[sid]["url"] for sid, code in (("la_county", SRC_COUNTY), ("glendale_cams", SRC_GLENDALE)) if mask & (1 << code)] lines += [ARCGIS_SOURCES["la_city_boundary"]["url"] + " (used to clip the LADOT layers; not redistributed)", "Used under the County of Los Angeles eGIS Terms of Use. Not endorsed by the County.", ""] if mask & (1 << SRC_CALTRANS): lines += ["Source: California Department of Transportation (Caltrans), District 7,", "SHS existing conditions: Speed limit (mph). " + ARCGIS_SOURCES["caltrans_d7"]["url"], ""] if any(mask & (1 << c) for c in OSM_DERIVED): influence = dict(manifest["osmInfluenceOnOfficialRecords"]["cvcKmRemovedWhereOsmHasALowerValue"]) stale_km = sum(manifest["osmInfluenceOnOfficialRecords"]["staleBaseKmRemovedWhereOsmHasAHigherValue"].values()) lines += ["Speed limit data (c) OpenStreetMap contributors, available under the Open Database", "License 1.0 (https://opendatacommons.org/licenses/odbl/1-0/).", "See https://www.openstreetmap.org/copyright.", f"OSM snapshot: {osm_time.isoformat() if osm_time else 'unknown'}.", "The records with source code 0 or 1 are derived from OpenStreetMap and form a", "Derivative Database under the ODbL 1.0. A machine-readable copy of them and the", f"script that builds this file: {publish_url or ''}."] if influence or stale_km: lines += ["Official records take no geometry or values from OpenStreetMap. OSM values only", "decided where two official values are left out: LADOT's 25 mph prima facie value", "where OSM has a lower limit (" + ", ".join(f"{v} km at {k}" for k, v in influence.items()) + "),", f"and {stale_km:.1f} km of 2018 LADOT base-map values of 35 mph or less that the 2025", "LAMC_35 layer no longer lists, where OSM has a higher limit. A build with --no-osm", "has no OSM content and no OSM influence."] lines += [""] lines += ["LA City local streets show the 25 mph prima facie limit (California Vehicle Code", "22352) as labelled by LADOT. School-zone and time-of-day limits are not included.", "", "Speed limits may be missing, outdated or wrong. Posted signs always govern."] with open(path, "w", encoding="utf-8") as f: f.write("\n".join(lines) + "\n") def retrieved_day(metas: dict) -> str: days = sorted({m["retrievedUtc"][:10] for sid, m in metas.items() if sid in ("la_county", "glendale_cams", "la_city_boundary")}) return ", ".join(days) if days else "unknown" # --------------------------------------------------------------------------- self test def selftest() -> None: issues = collections.Counter() assert parse_mph("35 mph", issues) == 35 assert parse_mph("35", issues) == 0 assert parse_mph("50 km/h", issues) == 30 assert parse_mph("51 mph", issues) == 0 for bad in ("none", "signals", "walk", "CA:urban"): assert parse_mph(bad, issues) == 0 assert parse_speed_field("25 and 35") == 0 and parse_speed_field("35") == 35 and parse_speed_field(0) == 0 assert street_name("W", "SUNSET", "BLVD", " ") == "West Sunset Boulevard" assert street_name("N", "102ND ST") == "North 102nd Street" assert street_name(None, "MCCADDEN PL") == "McCadden Place" assert caltrans_ref("405") == "I 405" and caltrans_ref("101") == "US 101" and caltrans_ref("134") == "CA 134" # Carving: a 400 m line, the middle 200 m covered by a parallel line 8 m away. lat0, lon0 = 34.0, -118.4 dlon = 1.0 / M_PER_DEG_LON dlat = 1.0 / M_PER_DEG_LAT line = [(lat0, lon0 + k * 10 * dlon) for k in range(41)] index = LineIndex() index.add([(lat0 + 8 * dlat, lon0 + 100 * dlon), (lat0 + 8 * dlat, lon0 + 300 * dlon)], 0, FAM_SURFACE) parts = carve(line, index, 15.0, FAM_SURFACE) lengths = sorted(round(polyline_length_m(p)) for p in parts) assert len(parts) == 2 and all(93 <= x <= 107 for x in lengths), lengths assert carve(line, index, 5.0, FAM_SURFACE) == [line] # too far: nothing carved assert carve(line, index, 15.0, FAM_FREEWAY) == [line] # wrong family perpendicular = LineIndex() perpendicular.add([(lat0 - 100 * dlat, lon0 + 200 * dlon), (lat0 + 100 * dlat, lon0 + 200 * dlon)], 0, FAM_SURFACE) assert carve(line, perpendicular, 15.0, FAM_SURFACE) == [line] # crossing street: kept previous_block = LineIndex() previous_block.add([(lat0, lon0 - 200 * dlon), (lat0, lon0)], 0, FAM_SURFACE) kept = carve(line, previous_block, 15.0, FAM_SURFACE) # continues end to end assert len(kept) == 1 and abs(polyline_length_m(kept[0]) - 400) < 0.5, kept # Writer / reader round trip. import tempfile records = [((C_PRIMARY, 35, 35, SRC_LADOT_LAMC35, GROUND_LAYER_FLAGS, "Magnolia Boulevard"), line), ((C_RESIDENTIAL, 25, 25, SRC_LADOT_CVC, GROUND_LAYER_FLAGS | FLAG_PRIMA_FACIE, "Elm Street"), [(lat0 + 50 * dlat, lon0), (lat0 + 50 * dlat, lon0 + 300 * dlon)]), ((C_MOTORWAY, 65, 65, SRC_CALTRANS, GROUND_LAYER_FLAGS, "\tI 5"), [(lat0 - 60 * dlat, lon0), (lat0 - 60 * dlat, lon0 + 300 * dlon)]), ((C_RESIDENTIAL, 0, 0, SRC_NONE, GROUND_LAYER_FLAGS, ""), [(lat0 + 25 * dlat, lon0), (lat0 + 25 * dlat, lon0 + 300 * dlon)])] with tempfile.TemporaryDirectory() as tmp: path = os.path.join(tmp, "t.bin") write_resource(records, path, 20261027, 20261027, 30) result = verify(path, samples=50) assert result["ok"], result store = StreetLimitFile(path) hit = store.lookup(lat0 + 3 * dlat, lon0 + 150 * dlon, course=90.0) assert hit and hit["speed_mph"] == 35 and hit["source"] == "ladot_lamc_35", hit hit = store.lookup(lat0 - 58 * dlat, lon0 + 150 * dlon, course=270.0) assert hit and hit["speed_mph"] == 65 and hit["name"] == "\tI 5", hit assert store.lookup(33.69, -118.71) is None # Names. assert name_key("West Imperial Highway") == name_key("Imperial Hwy") == "imperial highway" assert name_key("Vista Del Mar") == "vista del mar" and name_key("North Broadway") == "broadway" assert street_name_full("W Glenoaks Blvd") == "West Glenoaks Boulevard" # Per-sample labels: a 400 m line, 30 mph between 100 and 300 m, a 4 m flicker ignored. def label_at(x, y, bearing): d = x - to_xy(lat0, lon0)[0] return 30 if 100 <= d <= 300 or 350 <= d <= 354 else 40 runs = label_runs(line, label_at) assert [lab for lab, _ in runs] == [40, 30, 40], [lab for lab, _ in runs] assert all(abs(polyline_length_m(p) - e) < 6 for (_, p), e in zip(runs, (100, 200, 100))), runs assert abs(sum(polyline_length_m(p) for _, p in runs) - 400) < 0.5 # City boundary: a 1 km square; inside, outside and the clip with its 15 m tolerance. square = {"geometry": {"type": "Polygon", "coordinates": [[ [lon0, lat0], [lon0 + 1000 * dlon, lat0], [lon0 + 1000 * dlon, lat0 + 1000 * dlat], [lon0, lat0 + 1000 * dlat], [lon0, lat0]]]}} city = CityBoundary([square]) cx, cy = to_xy(lat0 + 500 * dlat, lon0 + 500 * dlon) assert city.signed_distance(cx, cy) < 0 and city.outside_m(*to_xy(lat0 - 30 * dlat, lon0 + 500 * dlon)) > 25 assert abs(city.outside_m(*to_xy(lat0 - 10 * dlat, lon0 + 500 * dlon)) - 10) < 0.5 crossing = [(lat0 + 500 * dlat, lon0 + k * 20 * dlon) for k in range(60)] # 0 .. 1180 m east kept = city.clip(crossing) assert len(kept) == 1 and 1005 <= polyline_length_m(kept[0]) <= 1025, polyline_length_m(kept[0]) # Finalize: an out-and-back two-way record collapses to one, a 1 m record is dropped, # a oneway U is split into its two directions. out_back = line[:5] + line[3::-1] tiny = [(lat0, lon0), (lat0, lon0 + 1.0 * dlon)] u_turn = [(lat0, lon0), (lat0, lon0 + 100 * dlon), (lat0, lon0 + 10 * dlon)] two_way = (C_RESIDENTIAL, 25, 25, SRC_LADOT_CVC, GROUND_LAYER_FLAGS | FLAG_PRIMA_FACIE, "A") one_way = (C_RESIDENTIAL, 25, 0, SRC_OSM, GROUND_LAYER_FLAGS | FLAG_ONEWAY, "B") done, st = finalize_records([(two_way, out_back), (two_way, tiny), (one_way, u_turn)]) assert len(done) == 3 and st.get("dropped_duplicate") == 1 and st.get("split_at_reversal") == 2, (len(done), st) pieces, dropped = dedupe_pieces([(two_way, line), (two_way, line[::-1]), (one_way, line), (one_way, line[::-1])]) assert dropped == 1 and len(pieces) == 3 print("selftest ok") # --------------------------------------------------------------------------- CLI def main(argv=None) -> int: here = os.path.dirname(os.path.abspath(__file__)) parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) sub = parser.add_subparsers(dest="command", required=True) def common(p): p.add_argument("--cache", default=os.path.join(here, "cache"), help="download cache (raw pages, PBFs)") p.add_argument("--pause", type=float, default=1.5, help="seconds between ArcGIS requests") p.add_argument("--refresh", action="store_true", help="ignore cached downloads") b = sub.add_parser("build", help="download sources and write street_limits.bin") common(b) b.add_argument("--out", default=os.path.join(here, "out")) b.add_argument("--no-osm", action="store_true", help="official sources only (no ODbL content)") b.add_argument("--osm-pbf", action="append", help="use local OSM PBF extract(s) instead of downloading") b.add_argument("--pbf-url", action="append", help=f"OSM PBF extract URL(s) covering the box (default: OSM France LA + Ventura; " f"alternative: {GEOFABRIK_SOCAL})") b.add_argument("--strict-precedence", action="store_true", help="do not let a lower OSM limit replace LADOT's CVC 25 mph prima facie default") b.add_argument("--caltrans-surface", action="store_true", help="also use Caltrans surface-highway values (off: they fail QA)") default_suppress = os.path.join(here, "street_limits_suppress.json") b.add_argument("--suppress", default=default_suppress if os.path.exists(default_suppress) else None, help="JSON list of OSM values to remove (see apply_suppressions); default: " "street_limits_suppress.json next to this script when present") b.add_argument("--no-suppress", dest="suppress", action="store_const", const=None) b.add_argument("--no-county", action="store_true", help="leave the LA County Speed Network out (spec section 7's licensing fallback)") default_overrides = os.path.join(here, "street_limits_overrides.json") b.add_argument("--overrides", default=default_overrides if os.path.exists(default_overrides) else None, help="official-document values (sources 8 and 9); default: street_limits_overrides.json") b.add_argument("--no-overrides", dest="overrides", action="store_const", const=None) b.add_argument("--cameras", default=DEFAULT_CAMERAS, help="the app's cameras.json, read only, for the camera-segment gate") default_waivers = os.path.join(here, "street_limits_camera_waivers.json") b.add_argument("--camera-waivers", default=default_waivers, help="documented camera disagreements that do not fail the camera gate") b.add_argument("--publish-url", help="public URL of the ODbL copy, written into the NOTICE") b.add_argument("--release", action="store_true", help="fail unless the build is release ready (publish URL, nothing past its review date)") b.add_argument("--allow-gate-failures", action="store_true", help="write the outputs even when a gate fails (diagnosis only; the manifest records it)") b.add_argument("--simplify-m", type=float, default=5.0) b.add_argument("--max-chunk-m", type=float, default=1000.0) b.add_argument("--index-margin-m", type=int, default=30) f = sub.add_parser("fetch", help="download / refresh the ArcGIS layers into the cache only") common(f) v = sub.add_parser("verify", help="check a street_limits.bin file") v.add_argument("path") q = sub.add_parser("lookup", help="match a coordinate against a street_limits.bin file") q.add_argument("path") q.add_argument("lat", type=float) q.add_argument("lon", type=float) q.add_argument("--course", type=float) q.add_argument("--radius", type=float, default=20.0) sub.add_parser("selftest", help="unit checks of parsing, carving and the writer/reader") args = parser.parse_args(argv) if args.command == "build": build(args) elif args.command == "fetch": for sid in ARCGIS_SOURCES: fetch_arcgis(sid, os.path.join(os.path.abspath(args.cache), "arcgis"), args.refresh, args.pause) elif args.command == "verify": result = verify(args.path) print(json.dumps(result, indent=2)) return 0 if result["ok"] else 1 elif args.command == "lookup": store = StreetLimitFile(args.path) print(json.dumps(store.lookup(args.lat, args.lon, args.course, args.radius), indent=2)) else: selftest() return 0 if __name__ == "__main__": sys.exit(main())