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

[rlemke]

Completes the real-source set (the last TODO(real)). FetchRegionFootprints in real mode now pulls Microsoft Global ML Building Footprints for the region bbox — keyless quadkey-tiled GeoJSONL (7 MB index cached; zoom-9 tiles fetched+cached, buildings filtered to bbox). Capped at 16 tiles (city/metro scale); out-of-coverage regions keep has_reference=false (gated, not fabricated).

ComputeExtrinsicQuality now counts reference buildings per tract and computes nearest-neighbour positional accuracy through a cached STRtree over the footprint layer (fast at ~300k footprints).

Honest caveat (documented): footprint_completeness = OSM/reference caps at 1.0 where OSM out-maps the ML layer — e.g. SF has ~168k OSM buildings vs ~23k MS footprints in-bbox, so the metric discriminates best where OSM lags the reference. Verified live on SF (real completeness + 24–80 m positional accuracy). Opt-in live footprints test added. 13 offline tests pass; ruff clean.

🤖 Generated with Claude Code

[gemini-code-assist]

Code Review

This pull request integrates Microsoft Global ML Building Footprints as an authoritative extrinsic benchmark for the 'real' mode of the OSM equity tool, implementing tile-based downloading, caching, and an STRtree-based spatial index for faster quality calculations. The review feedback highlights several critical improvements, including optimizing memory usage and preventing cache corruption during tile downloads, handling non-point geometries defensively to avoid attribute errors, clamping latitudes to prevent math errors at the poles, caching the dataset index on disk to avoid redundant network requests, and correcting a centroid calculation bias caused by double-counting the closing vertex of polygons.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

[gemini-code-assist]

This block has two main issues:

1. gzip.decompress(local.read_bytes()).decode().splitlines() loads and decompresses the entire tile into memory at once, which can consume a significant amount of RAM for large footprint files. Streaming the file line-by-line using gzip.open is much more memory-efficient.
2. Writing directly to local can leave a corrupted file if the download or write is interrupted. Writing to a temporary file and atomically replacing it prevents cache corruption.

        if not local.exists():
            temp_local = local.with_suffix(".tmp")
            try:
                r = requests.get(url, headers=_UA, timeout=240)
                r.raise_for_status()
                temp_local.write_bytes(r.content)
                temp_local.replace(local)
            finally:
                if temp_local.exists():
                    try:
                        temp_local.unlink()
                    except Exception:
                        pass
        with gzip.open(local, "rt", encoding="utf-8") as f:
            for line in f:
                if not line.strip():

[gemini-code-assist]

If any OSM building geometry in osm_buildings is a Polygon or MultiPolygon (rather than a Point), accessing p.x or p.y will raise an AttributeError. Defensively extracting the representative point or centroid coordinates prevents crashes on non-point geometries.

Suggested change

	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))
	for p in osm_buildings[:120]: # cap for speed
	nearest = ref_tree.nearest(p)
	r = ref_reps[nearest]
	p_pt = p if p.geom_type == "Point" else p.representative_point()
	dists.append(_haversine_km((p_pt.x, p_pt.y), (r.x, r.y)) * 1000.0)
	positional = float(np.mean(dists))

[gemini-code-assist]

The latitude calculation in _tile_xy can raise a ZeroDivisionError or ValueError (math domain error) if the latitude is exactly at the poles (90 or -90 degrees). Clamping the latitude to the valid Web Mercator range [-85.05112878, 85.05112878] prevents potential crashes in high-latitude regions.

Suggested change

	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 _tile_xy(lon: float, lat: float, z: int) -> tuple[int, int]:
	lat = max(-85.05112878, min(85.05112878, lat))
	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))

[gemini-code-assist]

The Microsoft Global ML Building Footprints index is about 7 MB. Currently, _load_ms_index only caches it in memory, meaning it is downloaded over the network on every execution. Adding disk caching to cache_dir avoids redundant network requests, improving performance and reliability.

def _load_ms_index(cache_dir: str | None = None) -> dict[str, list[str]]:
    global _ms_index
    if _ms_index is not None:
        return _ms_index
    import csv
    import io
    from pathlib import Path

    txt = None
    if cache_dir:
        local_csv = Path(cache_dir) / "ms_dataset_links.csv"
        if local_csv.exists():
            try:
                txt = local_csv.read_text(encoding="utf-8")
            except Exception:
                pass

    if txt is None:
        txt = requests.get(_MS_INDEX_URL, headers=_UA, timeout=90).text
        if cache_dir:
            local_csv = Path(cache_dir) / "ms_dataset_links.csv"
            try:
                local_csv.parent.mkdir(parents=True, exist_ok=True)
                local_csv.write_text(txt, encoding="utf-8")
            except Exception:
                pass

    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

[gemini-code-assist]

Pass cache_dir to _load_ms_index to enable disk caching of the dataset links index.

Suggested change

	index = _load_ms_index()
	index = _load_ms_index(cache_dir)

[gemini-code-assist]

In closed GeoJSON Polygons, the first and last coordinates are identical. Averaging all coordinates directly in _poly_centroid double-counts the closing vertex, biasing the centroid. Excluding the last coordinate resolves this issue.

Suggested change

	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 _poly_centroid(coords: list) -> tuple[float, float]:
	ring = coords[0]
	pts = ring[:-1] if len(ring) > 1 and ring[0] == ring[-1] else ring
	if not pts:
	return 0.0, 0.0
	return sum(p[0] for p in pts) / len(pts), sum(p[1] for p in pts) / len(pts)