diff --git a/examples/osm-equity/README.md b/examples/osm-equity/README.md index 4a888797..8b1457ad 100644 --- a/examples/osm-equity/README.md +++ b/examples/osm-equity/README.md @@ -42,12 +42,22 @@ tests/test_osm_equity.py FFL-compile + unit + end-to-end signal recovery Handlers do **real** geometry (shapely) and statistics (scipy Spearman, a pure-numpy permutation Moran's I, a Gaussian-kernel local GWR). The data source -defaults to a deterministic offline generator that tiles the region into a grid -of synthetic tracts whose income gradient is spatially smooth and drives OSM -richness — so the whole pipeline runs with **no network** and the analysis -recovers a real signal. Points marked `TODO(real)` in `equity_utils.py` are -where a deployment swaps in Census TIGER, the ACS API (`CENSUS_API_KEY`), an -Overpass/ohsome pull, and Microsoft/Google Open Buildings. +is selected by `FW_EQUITY_SOURCE`: + +- **`offline`** (default) — a deterministic generator tiles the region into a + 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). +- **`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). + Realistic at **city/metro scale**. The extrinsic footprint benchmark isn't + wired, so real mode honestly reports `has_reference=false` (gated). + +```bash +# real run (needs network + CENSUS_API_KEY): +FW_EQUITY_SOURCE=real python -c "..." # or a runner with the env set +``` ## Run the tests diff --git a/examples/osm-equity/handlers/shared/equity_utils.py b/examples/osm-equity/handlers/shared/equity_utils.py index d2054363..78fab3b1 100644 --- a/examples/osm-equity/handlers/shared/equity_utils.py +++ b/examples/osm-equity/handlers/shared/equity_utils.py @@ -1,21 +1,20 @@ """Shared computation for the osm.equity mapping-equity study. Everything here is real geometry (shapely) + real statistics (scipy + -pure-numpy Moran's I and locally-weighted regression). The *data source* -has two paths: - - * a deterministic offline generator (default) that tiles the study - region into a grid of synthetic census tracts whose income gradient - is spatially smooth AND drives OSM feature density/recency/diversity - (the "digital divide" signal) - so the downstream correlation, - spatial-autocorrelation and GWR steps find a real relationship and - the whole workflow runs with no network; and - * hooks where a real deployment would swap in TIGER tracts, an Overpass - /ohsome OSM pull, an ACS API call, and an authoritative footprint - layer (Microsoft/Google Open Buildings). Those are marked TODO(real). - -Determinism: all randomness is seeded from the geoid, so a re-run with the -same inputs is byte-identical (no os.urandom / unseeded RNG). +pure-numpy Moran's I and locally-weighted regression). The *data source* is +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. + * ``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 + the whole workflow runs with no network and the analysis recovers a known + relationship. Used by the test-suite; seeded from the geoid, so re-runs are + byte-identical (no os.urandom / unseeded RNG). """ from __future__ import annotations @@ -105,6 +104,13 @@ def grid_size() -> int: return 6 +def equity_mode() -> str: + """'real' -> live TIGER/ACS/Overpass sources; 'offline' (default) -> the + deterministic generator. Opt in with FW_EQUITY_SOURCE=real (needs network + + CENSUS_API_KEY; realistic at city/metro scale).""" + return os.environ.get("FW_EQUITY_SOURCE", "offline").strip().lower() + + # --------------------------------------------------------------------------- # Phase 1 / 2 - study area, tracts, region OSM, footprints # --------------------------------------------------------------------------- @@ -145,10 +151,12 @@ def tract_ij(geoid: str) -> tuple[int, int]: def make_tracts(region: str) -> list[dict[str, Any]]: - """Tile the region bbox into a grid of tract polygons (Phase 2 units). + """Phase 2 analysis units. Real mode -> Census TIGERweb tract polygons for + the region bbox; offline -> a synthetic grid tiling the bbox.""" + if equity_mode() == "real": + from . import sources_real as R - TODO(real): replace with a Census TIGER/Line tract fetch for the region. - """ + return R.real_resolve_tracts(region_bbox(region)) n = grid_size() lon0, lat0, lon1, lat1 = region_bbox(region) dlon = (lon1 - lon0) / n @@ -266,10 +274,16 @@ def fetch_region_osm(region: str, update: bool, tracts: list[dict]) -> tuple[str snap_path = out_dir() / f"osm_{_slug(region)}.snapshot" if path.exists() and snap_path.exists() and not update: return str(path), snap_path.read_text().strip(), True - features = generate_region_osm(region, tracts) - _write_geojson(path, features) - # snapshot stamp is derived from content, not wall-clock (determinism) - snapshot = "osm-" + hashlib.sha256(path.read_bytes()).hexdigest()[:12] + if equity_mode() == "real": + from . import sources_real as R + + features, snapshot = R.real_fetch_region_osm(region_bbox(region)) + _write_geojson(path, features) + else: + features = generate_region_osm(region, tracts) + _write_geojson(path, features) + # snapshot stamp is derived from content, not wall-clock (determinism) + snapshot = "osm-" + hashlib.sha256(path.read_bytes()).hexdigest()[:12] snap_path.write_text(snapshot) return str(path), snapshot, False @@ -284,9 +298,13 @@ def fetch_region_footprints( tract rather than fabricating an extrinsic score. TODO(real): pull Microsoft/Google Open Buildings (or municipal open data) - for the region bbox. + 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. """ - # A deployment might have no benchmark for some regions/sources. + if equity_mode() == "real": + return "", False + # Offline: a deployment might have no benchmark for some regions/sources. available = source.strip().lower() in { "microsoft_open_buildings", "google_open_buildings", @@ -328,11 +346,16 @@ def fetch_region_footprints( def fetch_census_equity(geoid: str, acs_year: int, region: str | None = None) -> dict[str, Any]: """Per-tract ACS equity variables (Phase 2 key variables). - TODO(real): call the Census ACS API (needs CENSUS_API_KEY) for the given - year and tract, pulling B19013 (income), B25070 (rent burden), B03002 - (race), B15003 (education), B16004 (English), B08201 (vehicles), - B28002 (internet). + Real mode -> Census ACS 5-year API (B19013 income, B25070 rent burden, + B03002 race, B06009 education, C16002 English, B08201 vehicles, B28002 + internet), fetched per county and cached. Offline -> synthetic profile. """ + if equity_mode() == "real": + from . import sources_real as R + + return R.real_fetch_census_equity( + geoid, int(acs_year), os.environ.get("CENSUS_API_KEY", "") + ) prof = tract_profile(geoid, region) wealth = prof["wealth"] rng = _rng(geoid, "acs", acs_year) @@ -353,15 +376,38 @@ def _pct(base, span, lo=0.0, hi=1.0): } +# region-extract STRtree cache: {region_osm_path: (features, rep_points, STRtree)} +# built once, reused across all per-tract clips (so a 254k-feature real extract +# is indexed a single time rather than rescanned per tract). +_clip_index: dict[str, tuple[list[dict], list, Any]] = {} + + +def _region_index(region_osm_path: str): + cached = _clip_index.get(region_osm_path) + if cached is not None: + return cached + from shapely.strtree import STRtree + + feats = _read_geojson(region_osm_path) + reps = [] + for f in feats: + g = shape(f["geometry"]) + reps.append(g if g.geom_type == "Point" else g.representative_point()) + tree = STRtree(reps) + entry = (feats, reps, tree) + _clip_index[region_osm_path] = entry + return entry + + def clip_tract_osm(region_osm_path: str, geometry_wkt: str) -> str: - """Clip the shared region extract to a tract polygon - LOCAL, no network.""" + """Clip the shared region extract to a tract polygon - LOCAL, no network. + Uses a cached STRtree over the region so each tract is a fast index query.""" poly = shapely_wkt.loads(geometry_wkt) + feats, reps, tree = _region_index(region_osm_path) kept = [] - for f in _read_geojson(region_osm_path): - geom = shape(f["geometry"]) - rep = geom if geom.geom_type == "Point" else geom.representative_point() - if poly.covers(rep): - kept.append(f) + for i in tree.query(poly): # bbox candidates + if poly.covers(reps[i]): + kept.append(feats[i]) tag = hashlib.sha256((region_osm_path + geometry_wkt).encode()).hexdigest()[:10] path = out_dir() / f"tract_{tag}.geojson" _write_geojson(path, kept) @@ -659,13 +705,15 @@ def render_maplibre_map( inline GeoJSON, CARTO Voyager raster basemap, click popups, a legend, and a findings banner. Data-desert tracts get a red outline.""" fc = json.dumps({"type": "FeatureCollection", "features": feats}) - xs, ys = [], [] + # total bounds via shapely (handles Polygon + MultiPolygon uniformly) + xmin = ymin = float("inf") + xmax = ymax = float("-inf") for f in feats: - for ring in f["geometry"]["coordinates"]: - for x, y in ring: - xs.append(x) - ys.append(y) - bounds = [[min(xs), min(ys)], [max(xs), max(ys)]] if xs else [[-124.4, 32.5], [-114.1, 42.0]] + a, b, c, d = shape(f["geometry"]).bounds + xmin, ymin, xmax, ymax = min(xmin, a), min(ymin, b), max(xmax, c), max(ymax, d) + bounds = ( + [[xmin, ymin], [xmax, ymax]] if xmin != float("inf") else [[-124.4, 32.5], [-114.1, 42.0]] + ) rho = stats_res.get("spearman_rho", 0.0) moran = stats_res.get("morans_i", 0.0) trend = stats_res.get("temporal_trend", "n/a") @@ -780,7 +828,9 @@ def build_report(records: list[dict], stats_res: dict, title: str) -> tuple[str, # the interactive heat map (this is the returned map_path) map_path = out_dir() / "equity_map.html" - map_path.write_text(render_maplibre_map(feats, title, stats_res, deserts)) + map_path.write_text( + render_maplibre_map(feats, title, stats_res, deserts, synthetic=equity_mode() != "real") + ) rho = stats_res.get("spearman_rho", 0.0) moran = stats_res.get("morans_i", 0.0) diff --git a/examples/osm-equity/handlers/shared/sources_real.py b/examples/osm-equity/handlers/shared/sources_real.py new file mode 100644 index 00000000..ab07a1a6 --- /dev/null +++ b/examples/osm-equity/handlers/shared/sources_real.py @@ -0,0 +1,300 @@ +"""Live data sources for the osm.equity study (opt-in via FW_EQUITY_SOURCE=real). + +Replaces the deterministic offline generator with real APIs: + + * 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) + +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 +(one HTTP round-trip per county, cached), so continent-scale runs should use the +offline generator. All fetches retry with fallback endpoints; failures raise +explicitly rather than silently degrading (the caller decides mode). +""" + +from __future__ import annotations + +import math +import time +from datetime import UTC, datetime +from typing import Any + +import requests +from shapely.geometry import shape + +_UA = {"User-Agent": "facetwork-osm-equity/1.0 (research; github.com/rlemke/facetwork)"} + +_TIGER_TRACTS = ( + "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Tracts_Blocks/MapServer/4/query" +) +_OVERPASS_EPS = [ + "https://maps.mail.ru/osm/tools/overpass/api/interpreter", + "https://overpass-api.de/api/interpreter", + "https://overpass.private.coffee/api/interpreter", +] + +# ACS 5-year detailed-table variables (one request per county returns all). +_ACS_VARS = [ + "B19013_001E", # median household income + "B25070_001E", + "B25070_007E", + "B25070_008E", + "B25070_009E", + "B25070_010E", # rent burden >=30% + "B03002_001E", + "B03002_003E", # total pop, white-non-hispanic (BIPOC = 1 - w/total) + "B06009_001E", + "B06009_002E", # education total, less-than-HS + "C16002_001E", + "C16002_004E", + "C16002_007E", + "C16002_010E", + "C16002_013E", # limited-English HH + "B08201_001E", + "B08201_002E", # households, zero-vehicle + "B28002_001E", + "B28002_013E", # internet total, no-internet +] + +_NULL = -666666600 # ACS uses -666666666 etc. as missing-value sentinels + +_census_cache: dict[tuple[str, str, int], dict[str, dict]] = {} + + +def _km_per_deg(lat: float) -> tuple[float, float]: + return 110.574, 111.320 * math.cos(math.radians(lat)) + + +# --------------------------------------------------------------------------- +# TIGERweb — real census-tract polygons intersecting the region bbox +# --------------------------------------------------------------------------- + + +def real_resolve_tracts(bbox: tuple[float, float, float, float]) -> list[dict[str, Any]]: + xmin, ymin, xmax, ymax = bbox + tracts: list[dict[str, Any]] = [] + offset = 0 + while True: + params = { + "geometry": f"{xmin},{ymin},{xmax},{ymax}", + "geometryType": "esriGeometryEnvelope", + "inSR": "4326", + "spatialRel": "esriSpatialRelIntersects", + "outFields": "GEOID,NAME", + "returnGeometry": "true", + "outSR": "4326", + "f": "geojson", + "resultOffset": offset, + "resultRecordCount": 800, + } + r = requests.get(_TIGER_TRACTS, params=params, headers=_UA, timeout=60) + r.raise_for_status() + feats = r.json().get("features", []) + if not feats: + break + for f in feats: + if not f.get("geometry"): + continue + geom = shape(f["geometry"]) + c = geom.centroid + kpl, kpo = _km_per_deg(c.y) + p = f.get("properties", {}) + tracts.append( + { + "geoid": p.get("GEOID", ""), + "name": p.get("NAME", ""), + "geometry_wkt": geom.wkt, + "area_km2": round(geom.area * kpl * kpo, 4), + } + ) + if len(feats) < 800: + break + offset += 800 + if not tracts: + raise RuntimeError("TIGERweb returned no tracts for the region bbox") + return tracts + + +# --------------------------------------------------------------------------- +# Census ACS — real equity variables (per county, cached) +# --------------------------------------------------------------------------- + + +def _num(v: Any) -> float | None: + try: + f = float(v) + except (TypeError, ValueError): + return None + return f if f > _NULL else None + + +def _parse_acs_row(row: list, idx: dict[str, int]) -> dict[str, Any]: + """Turn one ACS API row into an EquityVars dict (derived ratios, nulls dropped).""" + + def g(k: str) -> float | None: + return _num(row[idx[k]]) + + def ratio(numer_keys: list[str], denom_key: str) -> float: + d = g(denom_key) or 0 + if not d: + return 0.0 + n = sum(x for x in (g(k) for k in numer_keys) if x) + return round(n / d, 4) + + geoid = row[idx["state"]] + row[idx["county"]] + row[idx["tract"]] + inc = g("B19013_001E") + bt, wnh = g("B03002_001E") or 0, g("B03002_003E") or 0 + it, noi = g("B28002_001E") or 0, g("B28002_013E") or 0 + return { + "geoid": geoid, + "median_income": float(inc) if inc is not None else 0.0, + "pct_rent_burdened": ratio( + ["B25070_007E", "B25070_008E", "B25070_009E", "B25070_010E"], "B25070_001E" + ), + "pct_bipoc": round(1 - wnh / bt, 4) if bt else 0.0, + "pct_no_hs_diploma": ratio(["B06009_002E"], "B06009_001E"), + "pct_limited_english": ratio( + ["C16002_004E", "C16002_007E", "C16002_010E", "C16002_013E"], "C16002_001E" + ), + "pct_zero_vehicle": ratio(["B08201_002E"], "B08201_001E"), + "pct_internet_subscription": round(1 - noi / it, 4) if it else 0.0, + } + + +def real_region_census(state: str, county: str, year: int, key: str) -> dict[str, dict]: + ck = (state, county, year) + if ck in _census_cache: + return _census_cache[ck] + url = f"https://api.census.gov/data/{year}/acs/acs5" + params = { + "get": "NAME," + ",".join(_ACS_VARS), + "for": "tract:*", + "in": f"state:{state} county:{county}", + "key": key, + } + r = requests.get(url, params=params, headers=_UA, timeout=60) + r.raise_for_status() + rows = r.json() + idx = {h: i for i, h in enumerate(rows[0])} + out: dict[str, dict] = {} + for row in rows[1:]: + rec = _parse_acs_row(row, idx) + out[rec["geoid"]] = rec + _census_cache[ck] = out + return out + + +def real_fetch_census_equity(geoid: str, year: int, key: str) -> dict[str, Any]: + if not key: + raise RuntimeError("real mode needs CENSUS_API_KEY for ACS equity data") + state, county = geoid[:2], geoid[2:5] + reg = real_region_census(state, county, year, key) + v = reg.get(geoid) + if v is None: + # tract with no ACS row (e.g. water/institutional) — explicit, not fabricated + raise RuntimeError(f"no ACS data for tract {geoid}") + return v + + +# --------------------------------------------------------------------------- +# Overpass — real OSM features for the region bbox +# --------------------------------------------------------------------------- + + +def _overpass(query: str, retries: int = 3) -> list[dict]: + last = "" + for _ in range(retries): + for ep in _OVERPASS_EPS: + try: + r = requests.post(ep, data={"data": query}, headers=_UA, timeout=180) + if r.status_code == 200 and r.headers.get("content-type", "").startswith( + "application/json" + ): + return r.json().get("elements", []) + last = f"{ep} -> HTTP {r.status_code}" + except Exception as e: # noqa: BLE001 + last = f"{ep} -> {str(e)[:60]}" + time.sleep(3) + raise RuntimeError(f"Overpass unavailable after retries ({last})") + + +def _recency_days(ts: str | None, now: datetime) -> int: + if not ts: + return 0 + try: + dt = datetime.strptime(ts, "%Y-%m-%dT%H:%M:%SZ").replace(tzinfo=UTC) + return max(0, (now - dt).days) + except ValueError: + return 0 + + +def real_fetch_region_osm(bbox: tuple[float, float, float, float]) -> tuple[list[dict], str]: + xmin, ymin, xmax, ymax = bbox + b = f"{ymin},{xmin},{ymax},{xmax}" # Overpass order: S,W,N,E + now = datetime.now(UTC) + + # buildings + amenities: centroids are enough (count/density/tags) + q_bp = ( + f"[out:json][timeout:180];(" + f'way["building"]({b});node["amenity"]({b});way["amenity"]({b});' + f");out center meta tags;" + ) + # highways need geometry for length + q_hw = f'[out:json][timeout:180];(way["highway"]({b}););out geom meta tags;' + + feats: list[dict] = [] + for el in _overpass(q_bp): + tags = el.get("tags", {}) or {} + rec = _recency_days(el.get("timestamp"), now) + if el["type"] == "node": + lon, lat = el.get("lon"), el.get("lat") + else: + c = el.get("center") or {} + lon, lat = c.get("lon"), c.get("lat") + if lon is None or lat is None: + continue + if "building" in tags: + props = {"osm_kind": "building", "building": tags["building"], "recency_days": rec} + else: + props = { + "osm_kind": "poi", + "amenity": tags.get("amenity", "other"), + "recency_days": rec, + } + for k in ("name", "addr:street", "opening_hours"): + if tags.get(k): + props[k] = tags[k] + feats.append( + { + "type": "Feature", + "properties": props, + "geometry": {"type": "Point", "coordinates": [lon, lat]}, + } + ) + + for el in _overpass(q_hw): + geom = el.get("geometry") or [] + if len(geom) < 2: + continue + coords = [[p["lon"], p["lat"]] for p in geom] + feats.append( + { + "type": "Feature", + "properties": { + "osm_kind": "road", + "highway": (el.get("tags") or {}).get("highway", "road"), + "recency_days": _recency_days(el.get("timestamp"), now), + }, + "geometry": {"type": "LineString", "coordinates": coords}, + } + ) + if not feats: + raise RuntimeError("Overpass returned no OSM features for the region bbox") + return feats, "osm-live-" + now.strftime("%Y%m%d") + + +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) diff --git a/examples/osm-equity/tests/test_osm_equity.py b/examples/osm-equity/tests/test_osm_equity.py index b77e0e9a..1368279d 100644 --- a/examples/osm-equity/tests/test_osm_equity.py +++ b/examples/osm-equity/tests/test_osm_equity.py @@ -14,6 +14,8 @@ import subprocess import sys +import pytest + _EXAMPLE_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) _REPO_ROOT = os.path.abspath(os.path.join(_EXAMPLE_ROOT, "..", "..")) _FFL = os.path.join(_EXAMPLE_ROOT, "ffl", "osm_equity.ffl") @@ -217,3 +219,26 @@ def test_extrinsic_na_without_reference(self, tmp_path, monkeypatch): records, _, _ = self._run(tmp_path, monkeypatch, benchmark="no_such_source") assert all(r["extrinsic"]["has_reference"] is False for r in records) assert all(r["extrinsic"]["footprint_completeness"] == -1.0 for r in records) + + +# --------------------------------------------------------------------------- +# 4. Live real sources (opt-in; skipped unless FW_EQUITY_LIVE_TEST=1 + key) +# --------------------------------------------------------------------------- +@pytest.mark.skipif( + os.environ.get("FW_EQUITY_LIVE_TEST") != "1" or not os.environ.get("CENSUS_API_KEY"), + reason="live test needs FW_EQUITY_LIVE_TEST=1 and CENSUS_API_KEY (network)", +) +class TestRealSources: + def test_tiger_and_acs_and_overpass(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 + + tracts = U.make_tracts("San Francisco, CA") # TIGERweb + assert len(tracts) > 50 and all(len(t["geoid"]) == 11 for t in tracts) + eq = U.fetch_census_equity(tracts[0]["geoid"], 2022) # ACS + assert eq["median_income"] >= 0 and 0 <= eq["pct_bipoc"] <= 1 + path, snap, cached = U.fetch_region_osm("San Francisco, CA", False, tracts) # Overpass + feats = U._read_geojson(path) + kinds = {f["properties"]["osm_kind"] for f in feats} + assert {"building", "road", "poi"} & kinds and len(feats) > 1000