feat(osm-equity): wire Microsoft Open Buildings extrinsic benchmark

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

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

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])