"""Build a local OSM feature extract for the configured tag groups. Why this exists. The pipeline's OSM need is not a query workload: the twelve tag groups in ``config/overpass_tag_groups.yaml`` are fixed, and ``out center tags`` means the only data wanted is a point and its tags. Asking a shared query service that question per address, per bbox, cost 94% of a cold run's wall clock, failed a pilot user outright, and -- because a 429 sends a tag group to a different mirror -- silently varied the data vintage between groups of the same run. A committed extract answers the same question locally, at one pinned vintage. Governance follows the Google Trends importer, deliberately. This writes ``manifest.yaml``, which nothing reads, carrying the digest of the SQLite artifact it just built. An owner reviews it and copies it to ``manifest.v1.yaml``; that is what the adapter opens, and it refuses an artifact whose digest does not match the approved manifest. The bytes are too large for Git, so the manifest is the committed record and the artifact is shipped beside it -- which is why the digest binding is the whole of the integrity story. Each extract is read three times: once for the relations and their membership, once for the node refs of the member ways, and once for the nodes and ways. That order exists because relations come last in a PBF while their coordinates come from objects that came earlier, and because the node-location table the third pass builds is the memory ceiling, so it is built once and both way geometries and relation members are resolved from it. Built per state rather than from a single national PBF because the node-location index is the memory ceiling: California alone peaks near 2.7 GB, and a national pass would want 15-25 GB. Per state it is bounded by the largest state, and each state's replication timestamp is recorded separately, so the manifest can state a vintage per source rather than one figure covering all of them. Usage: python scripts/build_osm_extract.py --pbf-dir DIR --output DIR """ from __future__ import annotations import argparse import hashlib import json import sqlite3 import sys import time from datetime import UTC, datetime from pathlib import Path import yaml ROOT = Path(__file__).resolve().parents[1] sys.path.insert(0, str(ROOT / "src")) from market_research.adapters.osm_selectors import ( # noqa: E402 candidate_keys, load_tag_groups, ) ARTIFACT_NAME = "osm-extract.v1.sqlite" # One relation-only pass per level of nesting, and a bound because relation # membership can be cyclic: recording that a limit was hit is better than either # looping forever or resolving an unknown fraction of a relation's extent. MAX_RELATION_DEPTH = 4 SCHEMA = """ CREATE TABLE feature ( id INTEGER PRIMARY KEY, osm_type TEXT NOT NULL, osm_id INTEGER NOT NULL, lat REAL NOT NULL, lon REAL NOT NULL, tags TEXT NOT NULL, -- How many coordinates the bbox centre was computed from. Nodes carry 0. For -- ways this resolves cross-state duplicates: Geofabrik clips with complete -- ways, so a way spanning a state line appears in both files, and the copy -- assembled from more nodes has the truer bbox centre. Relations carry the -- number of member coordinates resolved, which is not a completeness claim -- -- relations are NOT clipped complete, so their extents are unioned across -- extracts instead of arbitrated between them (see merge_relation). node_count INTEGER NOT NULL, UNIQUE (osm_type, osm_id) ); CREATE TABLE feature_group ( feature_id INTEGER NOT NULL REFERENCES feature(id) ON DELETE CASCADE, group_id TEXT NOT NULL, PRIMARY KEY (feature_id, group_id) ); CREATE INDEX feature_group_by_group ON feature_group(group_id); -- Degenerate rectangles: every feature is the point Overpass would have reported -- as its centre, and filter_to_polygon tests exactly that point downstream. CREATE VIRTUAL TABLE feature_rtree USING rtree(id, minlon, maxlon, minlat, maxlat); """ def digest_file(path: Path) -> str: hasher = hashlib.sha256() with path.open("rb") as stream: for chunk in iter(lambda: stream.read(1 << 20), b""): hasher.update(chunk) return hasher.hexdigest() def replication_timestamp(pbf: Path) -> str | None: """The OSM vintage the extract was cut from, per its file header.""" import osmium header = osmium.io.Reader(str(pbf)).header() return header.get("osmosis_replication_timestamp") or None def scan_relations(pbf: Path, groups, keys) -> tuple[dict, dict, dict]: """Matched relations, and the membership needed to place them. Two things make relations a separate pass rather than another entity bit on the feature pass. Their coordinates are not in the object: a relation's ``out center`` is the centre of the bounding box over its members' geometry, recursively, so the member ids have to be known before the node and way passes can collect what those members need. And relations come last in a PBF, so learning the ids from the same pass that resolves them is not possible. Returns the matched relations, the member lists of every relation reachable from them, and counters for what could not be reached. """ import osmium matched: dict[int, tuple[dict, list[str]]] = {} members: dict[int, list[tuple[str, int]]] = {} for obj in (osmium.FileProcessor(str(pbf), osmium.osm.RELATION) .with_filter(osmium.filter.KeyFilter(*keys))): tags = dict(obj.tags) hits = sorted(gid for gid, group in groups.items() if group.matches(tags)) if not hits: continue matched[obj.id] = (tags, hits) members[obj.id] = [(member.type, member.ref) for member in obj.members] # A matched relation may hold sub-relations -- ``site=climbing`` areas made # of crag relations are the shape this config actually produces -- and their # geometry is part of the parent's extent. One extra relation-only pass per # nesting level; ids already asked for are never asked again, so a member # clipped out of this extract does not loop. asked = set(members) depth = 0 while True: wanted = sorted({ref for member in members.values() for kind, ref in member if kind == "r" and ref not in asked}) if not wanted: truncated = 0 break if depth >= MAX_RELATION_DEPTH: truncated = len(wanted) break depth += 1 asked.update(wanted) for obj in (osmium.FileProcessor(str(pbf), osmium.osm.RELATION) .with_filter(osmium.filter.IdFilter(wanted))): members[obj.id] = [(member.type, member.ref) for member in obj.members] stats = { "relations_matched": len(matched), "nesting_depth": depth, "nesting_truncated": truncated, # A sub-relation this extract does not contain: the parent crosses the # extract's boundary, and its extent is therefore partial here. Union # across extracts covers a state line; nothing covers a national border. "member_relations_absent": len(asked - set(members)), } return matched, members, stats def member_geometry(rel_id: int, members: dict) -> tuple[set[int], set[int]]: """The node and way ids reachable from one relation. Iterative with a visited set because relation membership can be cyclic, and a cycle here would otherwise be an unbounded recursion in a build that takes forty minutes to reach the state it fails on. """ nodes: set[int] = set() ways: set[int] = set() seen: set[int] = set() stack = [rel_id] while stack: current = stack.pop() if current in seen: continue seen.add(current) for kind, ref in members.get(current, ()): if kind == "n": nodes.add(ref) elif kind == "w": ways.add(ref) else: stack.append(ref) return nodes, ways def member_way_nodes(pbf: Path, way_ids: set[int]) -> dict[int, list[int]]: """The node refs of the member ways, read by id in a way-only pass. Cheap because osmium filters on id in C++ and no other entity is decoded; the refs are resolved to coordinates from the feature pass's location table, so this does not need its own. """ import osmium refs: dict[int, list[int]] = {} if not way_ids: return refs for obj in (osmium.FileProcessor(str(pbf), osmium.osm.WAY) .with_filter(osmium.filter.IdFilter(sorted(way_ids)))): refs[obj.id] = [node.ref for node in obj.nodes] return refs def extract_state(pbf: Path, groups, keys) -> tuple[list[tuple], dict, dict]: """Every feature in one PBF matching any group, with its bbox centre. Three reads of the file: the relations and their membership, the node refs of the member ways, then the nodes and ways themselves. The last one builds the location table, which is the memory ceiling, so it runs once and both the way geometries and the relation members are resolved from it. """ import osmium matched_relations, members, stats = scan_relations(pbf, groups, keys) wanted_nodes: set[int] = set() wanted_ways: set[int] = set() reach: dict[int, tuple[set[int], set[int]]] = {} for rel_id in matched_relations: nodes, ways = member_geometry(rel_id, members) reach[rel_id] = (nodes, ways) wanted_nodes |= nodes wanted_ways |= ways way_refs = member_way_nodes(pbf, wanted_ways) stats["member_ways_absent"] = len(wanted_ways) - len(way_refs) for refs in way_refs.values(): wanted_nodes.update(refs) processor = (osmium.FileProcessor(str(pbf), osmium.osm.NODE | osmium.osm.WAY) .with_filter(osmium.filter.KeyFilter(*keys)) .with_locations()) rows: list[tuple] = [] unlocatable = 0 for obj in processor: tags = dict(obj.tags) matched = sorted(gid for gid, group in groups.items() if group.matches(tags)) if not matched: continue if obj.is_node(): kind, lat, lon, node_count = "node", obj.location.lat, obj.location.lon, 0 else: kind = "way" lats, lons = [], [] for node in obj.nodes: if node.location.valid(): lats.append(node.location.lat) lons.append(node.location.lon) if not lats: unlocatable += 1 continue node_count = len(lats) lat = (min(lats) + max(lats)) / 2 lon = (min(lons) + max(lons)) / 2 rows.append((kind, obj.id, round(lat, 7), round(lon, 7), json.dumps(tags, sort_keys=True, ensure_ascii=False), node_count, ",".join(matched))) stats["unlocatable"] = unlocatable # The location table is filtered-object-independent: pyosmium caches every # node's location before the key filter runs, so member nodes carrying no # tags at all are in it. locations = processor.node_location_storage extents: dict[int, tuple] = {} for rel_id, (nodes, ways) in reach.items(): coordinates = set(nodes) for way_id in ways: coordinates.update(way_refs.get(way_id, ())) lats, lons = [], [] for node_id in coordinates: try: location = locations.get(node_id) except (KeyError, osmium.InvalidLocationError): continue if location.valid(): lats.append(location.lat) lons.append(location.lon) if not lats: continue tags, hits = matched_relations[rel_id] extents[rel_id] = (min(lats), max(lats), min(lons), max(lons), len(lats), json.dumps(tags, sort_keys=True, ensure_ascii=False), ",".join(hits)) stats["relations_placed"] = len(extents) stats["relations_unplaced"] = len(matched_relations) - len(extents) return rows, extents, stats def merge_relation(current: tuple | None, incoming: tuple) -> tuple: """Combine one relation's extent across state extracts. Ways are arbitrated -- Geofabrik clips with complete ways, so each extract holds the whole geometry and the copies agree. Relations are not clipped complete: a route crossing a state line has some members in one file and the rest in another, and picking either copy would report the centre of a truncated extent as though it were the object's. So the extents are unioned, and the tags come from whichever copy resolved more of the geometry. """ if current is None: return incoming tags, hits = (incoming[5:7] if incoming[4] > current[4] else current[5:7]) return (min(current[0], incoming[0]), max(current[1], incoming[1]), min(current[2], incoming[2]), max(current[3], incoming[3]), max(current[4], incoming[4]), tags, hits) def relation_rows(extents: dict[int, tuple]) -> list[tuple]: """The merged relation extents as feature rows, centre of the bbox.""" return [("relation", rel_id, round((min_lat + max_lat) / 2, 7), round((min_lon + max_lon) / 2, 7), tags, count, hits) for rel_id, (min_lat, max_lat, min_lon, max_lon, count, tags, hits) in sorted(extents.items())] def insert(connection: sqlite3.Connection, rows: list[tuple]) -> tuple[int, int, int]: """Insert, resolving cross-state duplicates. Returns (new, replaced, conflicts).""" new = replaced = disagreed = 0 for kind, osm_id, lat, lon, tags, node_count, matched in rows: existing = connection.execute( "SELECT id, lat, lon, node_count FROM feature WHERE osm_type=? AND osm_id=?", (kind, osm_id)).fetchone() if existing is None: cursor = connection.execute( "INSERT INTO feature (osm_type, osm_id, lat, lon, tags, node_count) " "VALUES (?,?,?,?,?,?)", (kind, osm_id, lat, lon, tags, node_count)) feature_id = cursor.lastrowid connection.execute( "INSERT INTO feature_rtree (id, minlon, maxlon, minlat, maxlat) " "VALUES (?,?,?,?,?)", (feature_id, lon, lon, lat, lat)) new += 1 else: feature_id, old_lat, old_lon, old_count = existing # A materially different centre for the same object means one copy was # assembled from a partial node set. Counted, not hidden: silent # disagreement here is a feature placed in the wrong neighbourhood. if abs(old_lat - lat) > 1e-4 or abs(old_lon - lon) > 1e-4: disagreed += 1 if node_count > old_count: connection.execute( "UPDATE feature SET lat=?, lon=?, tags=?, node_count=? WHERE id=?", (lat, lon, tags, node_count, feature_id)) connection.execute( "UPDATE feature_rtree SET minlon=?, maxlon=?, minlat=?, maxlat=? " "WHERE id=?", (lon, lon, lat, lat, feature_id)) replaced += 1 connection.executemany( "INSERT OR IGNORE INTO feature_group (feature_id, group_id) VALUES (?,?)", [(feature_id, gid) for gid in matched.split(",")]) return new, replaced, disagreed def main() -> None: parser = argparse.ArgumentParser() parser.add_argument("--pbf-dir", type=Path, required=True) parser.add_argument("--output", type=Path, required=True) parser.add_argument("--config", type=Path, default=ROOT / "config" / "overpass_tag_groups.yaml") args = parser.parse_args() config = yaml.safe_load(args.config.read_text(encoding="utf-8")) groups = load_tag_groups(config) keys = sorted(candidate_keys(groups)) pbfs = sorted(args.pbf_dir.glob("*-latest.osm.pbf")) if not pbfs: raise SystemExit(f"no *-latest.osm.pbf under {args.pbf_dir}") args.output.mkdir(parents=True, exist_ok=True) artifact = args.output / ARTIFACT_NAME artifact.unlink(missing_ok=True) connection = sqlite3.connect(artifact) connection.executescript(SCHEMA) print(f"{len(groups)} groups, {len(keys)} prefilter keys, {len(pbfs)} extracts") sources: list[dict] = [] totals = {"new": 0, "replaced": 0, "disagreed": 0, "unlocatable": 0, "relations_matched": 0, "relations_unplaced": 0, "member_ways_absent": 0, "member_relations_absent": 0, "nesting_truncated": 0, "relations_multi_state": 0, "relations_extent_unioned": 0} # Relations are held until every extract has been read, because a relation # split across a state line is only whole once its extents are unioned. relations: dict[int, tuple] = {} started = time.monotonic() for index, pbf in enumerate(pbfs, start=1): slug = pbf.name.removesuffix("-latest.osm.pbf") at = time.monotonic() vintage = replication_timestamp(pbf) rows, extents, stats = extract_state(pbf, groups, keys) new, replaced, disagreed = insert(connection, rows) connection.commit() totals["new"] += new totals["replaced"] += replaced totals["disagreed"] += disagreed for key in ("unlocatable", "relations_matched", "relations_unplaced", "member_ways_absent", "member_relations_absent", "nesting_truncated"): totals[key] += stats[key] for rel_id, incoming in extents.items(): current = relations.get(rel_id) merged = merge_relation(current, incoming) if current is not None: totals["relations_multi_state"] += 1 if merged[:4] != current[:4]: totals["relations_extent_unioned"] += 1 relations[rel_id] = merged sources.append({ "slug": slug, "pbf_sha256": digest_file(pbf), "replication_timestamp": vintage, "matched_features": len(rows), "matched_relations": stats["relations_matched"], }) print(f"[{index:2}/{len(pbfs)}] {slug:22} {len(rows):7} matched " f"{new:7} new {replaced:5} replaced {disagreed:4} centre-conflicts " f"{stats['relations_placed']:5} relations " f"{time.monotonic() - at:6.1f}s vintage={vintage}") new, _, _ = insert(connection, relation_rows(relations)) print(f"\n{new} relations written, {totals['relations_multi_state']} seen in more " f"than one extract ({totals['relations_extent_unioned']} with an extent the " "union extended)") connection.execute("ANALYZE") connection.commit() features = connection.execute("SELECT count(*) FROM feature").fetchone()[0] per_group = dict(connection.execute( "SELECT group_id, count(*) FROM feature_group GROUP BY group_id " "ORDER BY group_id").fetchall()) per_type = dict(connection.execute( "SELECT osm_type, count(*) FROM feature GROUP BY osm_type " "ORDER BY osm_type").fetchall()) connection.close() vintages = sorted({s["replication_timestamp"] for s in sources if s["replication_timestamp"]}) manifest = { "version": "v1", "status": "extract_unapproved", "approval_state": "provisional", "source_id": "osm_local_extract", "artifact": { "path": ARTIFACT_NAME, "sha256": digest_file(artifact), "bytes": artifact.stat().st_size, }, "tag_config": { "path": str(args.config.relative_to(ROOT)), "sha256": digest_file(args.config), }, "built_at": datetime.now(UTC).isoformat().replace("+00:00", "Z"), "attribution": "OpenStreetMap contributors / ODbL", # One vintage per source, and the range stated outright. The Overpass path # could not do this: a 429 moved a tag group to another mirror and its # vintage moved with it, while oldest_base_timestamp reported the minimum # as though it described the whole fetch. "vintage_earliest": vintages[0] if vintages else None, "vintage_latest": vintages[-1] if vintages else None, "feature_count": features, "features_per_group": per_group, "features_per_type": per_type, # Stated because a reader cannot tell a relation that was never matched # from one that could not be placed. The absent members are the extract's # boundary showing up in the data: a relation reaching outside the 51 # geographies is placed from the part that is inside. "relation_pass": { "matched": totals["relations_matched"], "placed": per_type.get("relation", 0), "unplaced": totals["relations_unplaced"], "multi_state": totals["relations_multi_state"], "extent_unioned": totals["relations_extent_unioned"], "member_ways_absent": totals["member_ways_absent"], "member_relations_absent": totals["member_relations_absent"], "nesting_truncated": totals["nesting_truncated"], }, "sources": sources, } (args.output / "manifest.yaml").write_text( yaml.safe_dump(manifest, sort_keys=False, allow_unicode=True), encoding="utf-8") print(f"\n{features} features in {artifact.stat().st_size / 1e6:.1f} MB " f"({time.monotonic() - started:.0f}s)") print(f"vintages {vintages[0] if vintages else '?'} .. " f"{vintages[-1] if vintages else '?'}") print(f"cross-state: {totals['replaced']} replaced by a more complete copy, " f"{totals['disagreed']} centre disagreements, " f"{totals['unlocatable']} ways with no usable location") print(f"relations: {totals['relations_matched']} matched, " f"{per_type.get('relation', 0)} placed, " f"{totals['relations_unplaced']} with no locatable member, " f"{totals['member_ways_absent']} member ways and " f"{totals['member_relations_absent']} member relations outside the extracts" + (f", NESTING TRUNCATED at depth {MAX_RELATION_DEPTH} " f"({totals['nesting_truncated']} unresolved)" if totals["nesting_truncated"] else "")) print(f"\nwrote {args.output / 'manifest.yaml'} (unapproved). Review it and copy " "to manifest.v1.yaml to make it loadable.") if __name__ == "__main__": main()