From 4e20d216a55949f8786f4deafe44411aa519e919 Mon Sep 17 00:00:00 2001 From: rlemke Date: Mon, 13 Jul 2026 19:15:32 -0700 Subject: [PATCH] feat(osm-equity): uniform square-tile analysis units + areal equity interpolation Adds an alternative to irregular census tracts: - build_tiles(bbox, sqmi=2.0): a regular grid of ~2-square-mile square tiles as tract-shaped units (uniform size => cleaner cross-area comparison, less modifiable-areal-unit-problem noise than variable tracts). - areal_interpolate_equity(tiles, tract_records): dasymetric-style areal interpolation that intersection-area-weights ACS tract equity onto each tile (tiles over no populated tract are omitted, not zeroed). The OSM-quality metrics and the MapLibre heat map work unchanged on tiles. On real SF data the 2-sq-mi tiling gives 45 tiles and raises Global Moran's I of attribute completeness from 0.33 (tracts) to 0.48 (tiles) - a cleaner spatial signal. 3 offline unit tests added; 16 offline tests pass; ruff clean. Co-Authored-By: Claude Opus 4.8 --- examples/osm-equity/README.md | 2 + .../handlers/shared/sources_real.py | 81 +++++++++++++++++++ examples/osm-equity/tests/test_osm_equity.py | 49 +++++++++++ 3 files changed, 132 insertions(+) diff --git a/examples/osm-equity/README.md b/examples/osm-equity/README.md index 8b1457ad..7430ab7b 100644 --- a/examples/osm-equity/README.md +++ b/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). diff --git a/examples/osm-equity/handlers/shared/sources_real.py b/examples/osm-equity/handlers/shared/sources_real.py index 45b7a901..0791a7a2 100644 --- a/examples/osm-equity/handlers/shared/sources_real.py +++ b/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 # --------------------------------------------------------------------------- diff --git a/examples/osm-equity/tests/test_osm_equity.py b/examples/osm-equity/tests/test_osm_equity.py index 9c1df43b..ae3de718 100644 --- a/examples/osm-equity/tests/test_osm_equity.py +++ b/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 # ---------------------------------------------------------------------------