feat(osm-equity): uniform 2-sq-mi tile units + areal equity interpolation

examples/osm-equity/README.md

@@ -48,6 +48,8 @@ is selected by `FW_EQUITY_SOURCE`:
   grid of synthetic tracts whose income gradient drives OSM richness, so the
   pipeline runs with **no network** and the analysis recovers a known signal
   (used by the test suite).
+> **Analysis units.** The default unit is the census tract. `sources_real.build_tiles(bbox, sqmi=2.0)` produces an alternative **uniform square-tile grid** (e.g. 2 sq mi), and `areal_interpolate_equity` dasymetrically area-weights ACS tract values onto the tiles — uniform units give a cleaner heat map and reduce modifiable-areal-unit noise (SF Moran's I 0.33 tract → 0.48 tile).
+
 - **`real`** — live sources (`handlers/shared/sources_real.py`): Census
   **TIGERweb** tract polygons, the Census **ACS 5-year API** (`CENSUS_API_KEY`),
   and **OSM via Overpass** (buildings/highways/amenities + edit metadata).

examples/osm-equity/handlers/shared/sources_real.py

@@ -26,6 +26,7 @@
 from typing import Any
 
 import requests
+from shapely import wkt as shapely_wkt
 from shapely.geometry import shape
 
 _UA = {"User-Agent": "facetwork-osm-equity/1.0 (research; github.com/rlemke/facetwork)"}
@@ -315,6 +316,86 @@ def bbox_area_km2(bbox: tuple[float, float, float, float]) -> float:
     return abs((xmax - xmin) * kpo) * abs((ymax - ymin) * kpl)
 
 
+# ---------------------------------------------------------------------------
+# Regular square-tile analysis units (alternative to irregular census tracts)
+# ---------------------------------------------------------------------------
+
+_SQMI_TO_KM2 = 2.589988
+
+
+def build_tiles(bbox: tuple[float, float, float, float], sqmi: float = 2.0) -> list[dict[str, Any]]:
+    """Break a bbox into a regular grid of ~`sqmi`-square-mile tiles, as
+    tract-shaped analysis units (uniform size makes cross-area OSM-quality
+    comparison cleaner than variable census tracts)."""
+    xmin, ymin, xmax, ymax = bbox
+    side_km = math.sqrt(sqmi * _SQMI_TO_KM2)
+    kpl, kpo = _km_per_deg((ymin + ymax) / 2)
+    dlat, dlon = side_km / kpl, side_km / kpo
+    nx = max(1, math.ceil((xmax - xmin) / dlon))
+    ny = max(1, math.ceil((ymax - ymin) / dlat))
+    tiles: list[dict[str, Any]] = []
+    for i in range(nx):
+        for j in range(ny):
+            x0, y0 = xmin + i * dlon, ymin + j * dlat
+            x1, y1 = min(x0 + dlon, xmax), min(y0 + dlat, ymax)
+            kpl2, kpo2 = _km_per_deg((y0 + y1) / 2)
+            area = abs((x1 - x0) * kpo2) * abs((y1 - y0) * kpl2)
+            tiles.append(
+                {
+                    "geoid": f"TILE_{i:03d}_{j:03d}",
+                    "name": f"Tile {i}.{j}",
+                    "geometry_wkt": f"POLYGON (({x0} {y0}, {x1} {y0}, {x1} {y1}, {x0} {y1}, {x0} {y0}))",
+                    "area_km2": round(area, 4),
+                }
+            )
+    return tiles
+
+
+_EQ_FIELDS = [
+    "median_income",
+    "pct_rent_burdened",
+    "pct_bipoc",
+    "pct_no_hs_diploma",
+    "pct_limited_english",
+    "pct_zero_vehicle",
+    "pct_internet_subscription",
+]
+
+
+def areal_interpolate_equity(tiles: list[dict], tract_records: list[dict]) -> dict[str, dict]:
+    """Assign each tile an equity profile by intersection-area-weighting the
+    ACS values of the census tracts it overlaps (dasymetric-style areal
+    interpolation). Tracts with no income data are excluded. Returns
+    {tile_geoid: equity_dict}; tiles covering no populated tract are omitted."""
+    from shapely.strtree import STRtree
+
+    geoms = [shapely_wkt.loads(t["geometry_wkt"]) for t in tract_records]
+    tree = STRtree(geoms)
+    out: dict[str, dict] = {}
+    for tile in tiles:
+        poly = shapely_wkt.loads(tile["geometry_wkt"])
+        acc = dict.fromkeys(_EQ_FIELDS, 0.0)
+        wsum = 0.0
+        for idx in tree.query(poly):
+            eq = tract_records[idx]["equity"]
+            if eq["median_income"] <= 0:
+                continue
+            inter = poly.intersection(geoms[idx])
+            w = inter.area
+            if w <= 0:
+                continue
+            for f in _EQ_FIELDS:
+                acc[f] += eq[f] * w
+            wsum += w
+        if wsum <= 0:
+            continue
+        ev = {"geoid": tile["geoid"]}
+        for f in _EQ_FIELDS:
+            ev[f] = round(acc[f] / wsum, 1 if f == "median_income" else 4)
+        out[tile["geoid"]] = ev
+    return out
+
+
 # ---------------------------------------------------------------------------
 # Microsoft Global ML Building Footprints — authoritative extrinsic benchmark
 # ---------------------------------------------------------------------------

examples/osm-equity/tests/test_osm_equity.py

@@ -105,6 +105,55 @@ def test_footprint_source_unknown_is_unavailable(self, tmp_path, monkeypatch):
         assert avail is False and path == ""
 
 
+# ---------------------------------------------------------------------------
+# 2b. Tile analysis units (offline-pure: no network)
+# ---------------------------------------------------------------------------
+class TestTiling:
+    def test_build_tiles_are_two_sqmi(self):
+        from handlers.shared import sources_real as R
+
+        tiles = R.build_tiles((-122.52, 37.70, -122.36, 37.83), sqmi=2.0)
+        assert len(tiles) > 20
+        full = [t["area_km2"] for t in tiles if t["area_km2"] > 5.0]
+        assert full and all(abs(a - 2 * 2.589988) < 0.05 for a in full)  # ~5.18 km2
+        assert len({t["geoid"] for t in tiles}) == len(tiles)  # unique ids
+
+    def test_areal_interpolation_area_weights(self):
+        from handlers.shared import sources_real as R
+
+        # two side-by-side tracts, a tile straddling both 50/50 -> mean income
+        tracts = [
+            {"geometry_wkt": "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))", "equity": _equity(100000)},
+            {"geometry_wkt": "POLYGON ((1 0, 2 0, 2 1, 1 1, 1 0))", "equity": _equity(50000)},
+        ]
+        tile = {
+            "geoid": "TILE_000_000",
+            "geometry_wkt": "POLYGON ((0.5 0, 1.5 0, 1.5 1, 0.5 1, 0.5 0))",
+        }
+        out = R.areal_interpolate_equity([tile], tracts)
+        assert abs(out["TILE_000_000"]["median_income"] - 75000) < 1  # 50/50 area weight
+
+    def test_areal_interpolation_skips_no_income(self):
+        from handlers.shared import sources_real as R
+
+        tracts = [{"geometry_wkt": "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))", "equity": _equity(0)}]
+        tile = {"geoid": "T", "geometry_wkt": "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))"}
+        assert R.areal_interpolate_equity([tile], tracts) == {}  # no populated tract -> omitted
+
+
+def _equity(income):
+    return {
+        "geoid": "x",
+        "median_income": float(income),
+        "pct_rent_burdened": 0.3,
+        "pct_bipoc": 0.5,
+        "pct_no_hs_diploma": 0.1,
+        "pct_limited_english": 0.1,
+        "pct_zero_vehicle": 0.2,
+        "pct_internet_subscription": 0.8,
+    }
+
+
 # ---------------------------------------------------------------------------
 # 3. End-to-end: drive every handler in the workflow's data order
 # ---------------------------------------------------------------------------