From 662a41445609005b583e57a3e7752a23ff506508 Mon Sep 17 00:00:00 2001 From: rlemke Date: Mon, 13 Jul 2026 18:38:51 -0700 Subject: [PATCH] feat(osm-equity): wire Microsoft Open Buildings as the extrinsic benchmark MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Completes the real-source set. FetchRegionFootprints in real mode now pulls Microsoft Global ML Building Footprints for the region bbox: - keyless quadkey-tiled GeoJSONL (dataset-links.csv index, ~7MB, cached); zoom-9 tiles covering the bbox are fetched + cached, buildings filtered to the bbox as centroid points. Capped at 16 tiles (city/metro scale); regions out of coverage keep has_reference=false (still gated, not fabricated). - ComputeExtrinsicQuality now counts reference buildings per tract and gets nearest-neighbour positional accuracy through a cached STRtree over the footprint layer (fast at ~300k footprints), same index path as the OSM clip. Honest domain caveat, documented: footprint_completeness = OSM/reference caps at 1.0 where OSM out-maps the ML layer (e.g. San Francisco, ~168k OSM buildings vs ~23k MS footprints in-bbox) — the metric discriminates best where OSM lags the reference. Verified live on SF: real completeness + 24-80m positional accuracy. Live footprints test added (opt-in). 13 offline tests pass; ruff clean. Co-Authored-By: Claude Opus 4.8 --- .../handlers/shared/equity_utils.py | 56 +++++---- .../handlers/shared/sources_real.py | 111 ++++++++++++++++++ examples/osm-equity/tests/test_osm_equity.py | 10 ++ 3 files changed, 154 insertions(+), 23 deletions(-) diff --git a/examples/osm-equity/handlers/shared/equity_utils.py b/examples/osm-equity/handlers/shared/equity_utils.py index 78fab3b1..65d67db9 100644 --- a/examples/osm-equity/handlers/shared/equity_utils.py +++ b/examples/osm-equity/handlers/shared/equity_utils.py @@ -5,10 +5,10 @@ selected by :func:`equity_mode` (env ``FW_EQUITY_SOURCE``): * ``real`` -> live sources in :mod:`sources_real`: Census TIGERweb tract - polygons, the Census ACS 5-year API (needs ``CENSUS_API_KEY``), and OSM - from Overpass. Realistic at city/metro scale. The extrinsic footprint - benchmark is not yet wired, so real mode honestly reports has_reference= - false (gated) rather than synthesising it. + polygons, the Census ACS 5-year API (needs ``CENSUS_API_KEY``), OSM from + Overpass, and Microsoft Global ML Building Footprints as the extrinsic + benchmark. Realistic at city/metro scale; where the region is out of + footprint coverage, has_reference stays false (gated). * ``offline`` (default) -> a deterministic generator that tiles the region into a grid of synthetic tracts whose income gradient is spatially smooth AND drives OSM density/recency/diversity (the "digital divide" signal), so @@ -297,13 +297,21 @@ def fetch_region_footprints( covers the region - the workflow then records has_reference=false per tract rather than fabricating an extrinsic score. - TODO(real): pull Microsoft/Google Open Buildings (or municipal open data) - for the region bbox. Not yet wired, so in real mode the extrinsic benchmark - is honestly reported as unavailable (has_reference=false) rather than - synthesised - the gating is the whole point of improvement #4. + Real mode -> Microsoft Global ML Building Footprints for the region bbox + (city/metro scale). If no tiles cover the region, has_reference stays false. + Offline -> a synthetic reference layer. """ if equity_mode() == "real": - return "", False + path = out_dir() / f"footprints_{_slug(region)}.geojson" + if path.exists() and not update: + return str(path), True + from . import sources_real as R + + feats = R.real_fetch_footprints(region_bbox(region), str(out_dir())) + if not feats: + return "", False + _write_geojson(path, feats) + return str(path), True # Offline: a deployment might have no benchmark for some regions/sources. available = source.strip().lower() in { "microsoft_open_buildings", @@ -493,23 +501,25 @@ def compute_extrinsic_quality( } poly = shapely_wkt.loads(geometry_wkt) osm_buildings = [ - f for f in _read_geojson(osm_path) if f["properties"].get("osm_kind") == "building" - ] - ref_buildings = [ - f - for f in _read_geojson(footprints_path) - if f["properties"].get("geoid") == geoid or poly.covers(shape(f["geometry"])) + shape(f["geometry"]) + for f in _read_geojson(osm_path) + if f["properties"].get("osm_kind") == "building" ] - ref_n = len(ref_buildings) + # reference buildings inside this tract, via the cached STRtree over the + # whole footprint layer (fast even at ~300k Microsoft footprints) + ref_feats, ref_reps, ref_tree = _region_index(footprints_path) + ref_idx = [i for i in ref_tree.query(poly) if poly.covers(ref_reps[i])] + ref_pts = [ref_reps[i] for i in ref_idx] + ref_n = len(ref_pts) completeness = (len(osm_buildings) / ref_n) if ref_n else 0.0 - # positional accuracy: mean nearest-neighbour distance OSM->reference (m) - if osm_buildings and ref_buildings: - ref_pts = [shape(f["geometry"]) for f in ref_buildings] + # positional accuracy: mean nearest-neighbour distance OSM->reference (m), + # nearest looked up through the footprint STRtree + if osm_buildings and ref_n: dists = [] - for b in osm_buildings[:80]: # cap for speed - p = shape(b["geometry"]) - nn = min(_haversine_km((p.x, p.y), (r.x, r.y)) for r in ref_pts[:120]) - dists.append(nn * 1000.0) + for p in osm_buildings[:120]: # cap for speed + nearest = ref_tree.nearest(p) + r = ref_reps[nearest] + dists.append(_haversine_km((p.x, p.y), (r.x, r.y)) * 1000.0) positional = float(np.mean(dists)) else: positional = -1.0 diff --git a/examples/osm-equity/handlers/shared/sources_real.py b/examples/osm-equity/handlers/shared/sources_real.py index ab07a1a6..0c90064d 100644 --- a/examples/osm-equity/handlers/shared/sources_real.py +++ b/examples/osm-equity/handlers/shared/sources_real.py @@ -5,6 +5,10 @@ * census tract geometry -> Census TIGERweb ArcGIS REST (keyless) * equity variables -> Census ACS 5-year API (needs CENSUS_API_KEY) * OSM features -> Overpass API (buildings/highways/amenities+meta) + * extrinsic footprints -> Microsoft Global ML Building Footprints (keyless, + quadkey-tiled GeoJSONL) as the authoritative + reference. NB: where OSM out-maps the ML layer + (e.g. well-mapped cities), OSM/reference caps at 1. Scope: real mode targets a CITY / METRO bounding box. Overpass and TIGERweb will not return a whole continent in one request, and ACS is fetched per county @@ -298,3 +302,110 @@ def bbox_area_km2(bbox: tuple[float, float, float, float]) -> float: xmin, ymin, xmax, ymax = bbox kpl, kpo = _km_per_deg((ymin + ymax) / 2) return abs((xmax - xmin) * kpo) * abs((ymax - ymin) * kpl) + + +# --------------------------------------------------------------------------- +# Microsoft Global ML Building Footprints — authoritative extrinsic benchmark +# --------------------------------------------------------------------------- + +_MS_INDEX_URL = "https://minedbuildings.z5.web.core.windows.net/global-buildings/dataset-links.csv" +_MS_ZOOM = 9 +_MS_MAX_TILES = 16 # cap: real footprints are a city/metro benchmark +_ms_index: dict[str, list[str]] | None = None + + +def _tile_xy(lon: float, lat: float, z: int) -> tuple[int, int]: + sinlat = math.sin(math.radians(lat)) + n = 1 << z + x = int((lon + 180.0) / 360.0 * n) + y = int((0.5 - math.log((1 + sinlat) / (1 - sinlat)) / (4 * math.pi)) * n) + return max(0, min(n - 1, x)), max(0, min(n - 1, y)) + + +def _xy_to_quadkey(x: int, y: int, z: int) -> str: + qk = "" + for i in range(z, 0, -1): + d, mask = 0, 1 << (i - 1) + if x & mask: + d += 1 + if y & mask: + d += 2 + qk += str(d) + return qk + + +def _bbox_quadkeys(bbox: tuple[float, float, float, float], z: int = _MS_ZOOM) -> set[str]: + xmin, ymin, xmax, ymax = bbox + x0, y0 = _tile_xy(xmin, ymax, z) + x1, y1 = _tile_xy(xmax, ymin, z) + return { + _xy_to_quadkey(x, y, z) + for x in range(min(x0, x1), max(x0, x1) + 1) + for y in range(min(y0, y1), max(y0, y1) + 1) + } + + +def _load_ms_index() -> dict[str, list[str]]: + global _ms_index + if _ms_index is not None: + return _ms_index + import csv + import io + + txt = requests.get(_MS_INDEX_URL, headers=_UA, timeout=90).text + idx: dict[str, list[str]] = {} + for row in csv.DictReader(io.StringIO(txt)): + idx.setdefault(row["QuadKey"], []).append(row["Url"]) + _ms_index = idx + return idx + + +def _poly_centroid(coords: list) -> tuple[float, float]: + ring = coords[0] + return sum(p[0] for p in ring) / len(ring), sum(p[1] for p in ring) / len(ring) + + +def real_fetch_footprints(bbox: tuple[float, float, float, float], cache_dir: str) -> list[dict]: + """Microsoft Global ML Building Footprints within the bbox, as centroid + Point features. Returns [] when no tiles cover the region or the region is + too large (the caller then reports has_reference=false).""" + import gzip + import json as _json + from pathlib import Path + + qks = _bbox_quadkeys(bbox) + if len(qks) > _MS_MAX_TILES: + return [] # region too large for the footprint benchmark + index = _load_ms_index() + urls: list[str] = [] + for qk in qks: + urls.extend(index.get(qk, [])) + if not urls: + return [] + xmin, ymin, xmax, ymax = bbox + tdir = Path(cache_dir) / "ms_tiles" + tdir.mkdir(parents=True, exist_ok=True) + feats: list[dict] = [] + for url in urls: + name = url.rsplit("/", 1)[-1] + local = tdir / name + if not local.exists(): + r = requests.get(url, headers=_UA, timeout=240) + r.raise_for_status() + local.write_bytes(r.content) + for line in gzip.decompress(local.read_bytes()).decode().splitlines(): + if not line.strip(): + continue + g = _json.loads(line).get("geometry") + if not g or g.get("type") != "Polygon": + continue + cx, cy = _poly_centroid(g["coordinates"]) + if xmin <= cx <= xmax and ymin <= cy <= ymax: + feats.append( + { + "type": "Feature", + "properties": {"ref": "building"}, + "geometry": {"type": "Point", "coordinates": [cx, cy]}, + } + ) + return feats diff --git a/examples/osm-equity/tests/test_osm_equity.py b/examples/osm-equity/tests/test_osm_equity.py index 1368279d..9c1df43b 100644 --- a/examples/osm-equity/tests/test_osm_equity.py +++ b/examples/osm-equity/tests/test_osm_equity.py @@ -242,3 +242,13 @@ def test_tiger_and_acs_and_overpass(self, tmp_path, monkeypatch): feats = U._read_geojson(path) kinds = {f["properties"]["osm_kind"] for f in feats} assert {"building", "road", "poi"} & kinds and len(feats) > 1000 + + def test_microsoft_open_buildings(self, tmp_path, monkeypatch): + monkeypatch.setenv("FW_CACHE_ROOT", str(tmp_path)) + monkeypatch.setenv("FW_EQUITY_SOURCE", "real") + from handlers.shared import equity_utils as U + + path, available = U.fetch_region_footprints("San Francisco, CA", "microsoft", False, []) + assert available and path + refs = U._read_geojson(path) + assert len(refs) > 1000 and all(r["properties"]["ref"] == "building" for r in refs[:5])