From b990f0abedffc83948e822027d3733da9875787d Mon Sep 17 00:00:00 2001 From: rlemke Date: Mon, 13 Jul 2026 17:07:59 -0700 Subject: [PATCH] =?UTF-8?q?feat(examples):=20osm-equity=20=E2=80=94=20Digi?= =?UTF-8?q?tal=20Divide=20&=20OSM=20Mapping=20Equity=20study?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A 5-phase FFL workflow (namespace osm.equity) encoding a tract-level study of how socioeconomic disparity shapes OSM data completeness/quality, with real handler implementations. Design: - Region OSM fetched ONCE and clipped per tract LOCALLY (no per-tract remote calls / rate limits); region fetches are UPDATE-GATED (reuse cache unless update=true). - Tract geometry carried into Phase 4 so Moran's I / GWR build real spatial weights; POI diversity (Shannon) measured and correlated with income. - Extrinsic footprint benchmark is pluggable and GATED — has_reference=false with an explicit N/A sentinel where no authoritative layer exists. - Provenance (acs_year + content-derived OSM snapshot) for reproducibility. Handlers do real geometry (shapely) + statistics (scipy Spearman, pure-numpy permutation Moran's I, Gaussian-kernel local GWR) over a deterministic offline data generator, so the pipeline runs with no network; TODO(real) markers show where TIGER / ACS API / Overpass / Open Buildings plug in. 13 tests pass: FFL compile+emit, util determinism, update-gating, extrinsic gating, and an end-to-end run that recovers the built-in digital-divide signal (Spearman rho>0.3 p<0.05, Moran's I>0) and flags low-income data deserts. ruff format + check clean. Co-Authored-By: Claude Opus 4.8 --- examples/osm-equity/README.md | 61 ++ examples/osm-equity/ffl/osm_equity.ffl | 197 +++++ examples/osm-equity/handlers/__init__.py | 33 + .../handlers/acquisition/__init__.py | 0 .../acquisition/acquisition_handlers.py | 84 ++ .../osm-equity/handlers/analysis/__init__.py | 0 .../handlers/analysis/analysis_handlers.py | 72 ++ .../osm-equity/handlers/design/__init__.py | 0 .../handlers/design/design_handlers.py | 48 ++ .../osm-equity/handlers/metrics/__init__.py | 0 .../handlers/metrics/metrics_handlers.py | 88 +++ .../osm-equity/handlers/reporting/__init__.py | 0 .../handlers/reporting/reporting_handlers.py | 55 ++ .../osm-equity/handlers/shared/__init__.py | 0 .../handlers/shared/equity_utils.py | 718 ++++++++++++++++++ examples/osm-equity/tests/__init__.py | 0 examples/osm-equity/tests/conftest.py | 28 + examples/osm-equity/tests/test_osm_equity.py | 218 ++++++ 18 files changed, 1602 insertions(+) create mode 100644 examples/osm-equity/README.md create mode 100644 examples/osm-equity/ffl/osm_equity.ffl create mode 100644 examples/osm-equity/handlers/__init__.py create mode 100644 examples/osm-equity/handlers/acquisition/__init__.py create mode 100644 examples/osm-equity/handlers/acquisition/acquisition_handlers.py create mode 100644 examples/osm-equity/handlers/analysis/__init__.py create mode 100644 examples/osm-equity/handlers/analysis/analysis_handlers.py create mode 100644 examples/osm-equity/handlers/design/__init__.py create mode 100644 examples/osm-equity/handlers/design/design_handlers.py create mode 100644 examples/osm-equity/handlers/metrics/__init__.py create mode 100644 examples/osm-equity/handlers/metrics/metrics_handlers.py create mode 100644 examples/osm-equity/handlers/reporting/__init__.py create mode 100644 examples/osm-equity/handlers/reporting/reporting_handlers.py create mode 100644 examples/osm-equity/handlers/shared/__init__.py create mode 100644 examples/osm-equity/handlers/shared/equity_utils.py create mode 100644 examples/osm-equity/tests/__init__.py create mode 100644 examples/osm-equity/tests/conftest.py create mode 100644 examples/osm-equity/tests/test_osm_equity.py diff --git a/examples/osm-equity/README.md b/examples/osm-equity/README.md new file mode 100644 index 00000000..4a734fbd --- /dev/null +++ b/examples/osm-equity/README.md @@ -0,0 +1,61 @@ +# osm-equity — Digital Divide & OSM Mapping Equity + +A 5-phase FFL study of how socioeconomic disparity shapes the completeness/ +quality of OpenStreetMap data, at census-tract scale. The FFL **is** the +methodology; the runtime schedules it by data flow. + +``` +ffl/osm_equity.ffl the workflow (namespace osm.equity) +handlers/ event-facet implementations (real geometry + stats) + shared/equity_utils.py computation + deterministic offline data source + design/ acquisition/ metrics/ analysis/ reporting/ +tests/test_osm_equity.py FFL-compile + unit + end-to-end signal recovery +``` + +## Phases → FFL + +| Phase | Plan section | FFL | +|-------|--------------|-----| +| 1 | Study design | `DefineStudyArea` + workflow params | +| 2 | Data acquisition | `ResolveStudyTracts`, `FetchRegionOSM`, `FetchRegionFootprints` (region-level, **update-gated**) | +| 3 | Metric calculation | per-tract `foreach`: `ClipTractOSM` (local), `FetchCensusEquity`, `ComputeAttributeQuality` (intrinsic), `ComputeExtrinsicQuality` (extrinsic, gated) | +| 4 | Statistical & spatial | `SpearmanCorrelation`, `MoransI`, `GeographicallyWeightedRegression`, `TemporalEvolution` — run concurrently over the full tract set | +| 5 | Actionable reporting | `BuildEquityReport` → HTML + choropleth GeoJSON + ranked data deserts | + +## Design choices (folded-in improvements) + +1. **Fetch OSM once, clip per tract locally.** Region OSM is pulled a single + time and each tract is obtained by a local geometry clip — no per-tract + remote calls (Overpass/ohsome rate limits). `ClipTractOSM` is pure-local. +2. **Update-gated fetch.** `FetchRegionOSM` / `FetchRegionFootprints` reuse the + cached extract and **only re-fetch when `update = true`**. +3. **Geometry carried into Phase 4** so Moran's I / GWR build a real spatial- + weights matrix (`TractQuality.geometry_wkt`). +4. **POI diversity measured** (Shannon entropy) and correlated with income. +5. **Extrinsic benchmark is pluggable + gated:** `ExtrinsicQuality.has_reference` + is `false` where no authoritative footprint layer exists — metrics are an + explicit N/A sentinel (`-1.0`), never silently zeroed. +6. **Provenance:** `acs_year` param + a content-derived OSM `snapshot` for + reproducible re-runs. + +## Data source + +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. + +## Run the tests + +```bash +pytest examples/osm-equity/tests -q +``` + +The end-to-end test drives every handler in the workflow's data order and +asserts the built-in digital-divide signal is recovered (Spearman ρ > 0.3, +p < 0.05; Moran's I > 0) and that data deserts are the low-income tracts. +Set `FW_EQUITY_GRID` to change the grid resolution (default 6 → 36 tracts). diff --git a/examples/osm-equity/ffl/osm_equity.ffl b/examples/osm-equity/ffl/osm_equity.ffl new file mode 100644 index 00000000..624bbade --- /dev/null +++ b/examples/osm-equity/ffl/osm_equity.ffl @@ -0,0 +1,197 @@ +namespace osm.equity { + + // =================================================================== + // Digital Divide & OSM Mapping Equity + // + // A 5-phase, census-tract-level study of how socioeconomic disparity + // shapes the completeness/quality of volunteered geographic + // information (OSM). The data flow *is* the methodology. + // + // Improvements folded in vs. the first sketch: + // 1. Region OSM is fetched ONCE (update-gated cache) and clipped + // per tract LOCALLY - no per-tract remote calls (rate limits). + // 2. Tract geometry (WKT) is carried into Phase 4 so Moran's I / GWR + // can build a real spatial-weights matrix. + // 3. POI diversity is measured (Phase 4.1 correlates it with income). + // 4. The extrinsic benchmark is a pluggable source and is GATED: + // ExtrinsicQuality.has_reference is false where no authoritative + // footprint layer exists (no silent zeroing). + // 5. Provenance (acs_year, osm snapshot) is carried for reproducibility. + // =================================================================== + + // ---- Analysis unit + data layers ---- + schema StudyArea { + region: String, + scale: String, + hypothesis: String, + acs_year: Long + } + + schema TractGeo { + geoid: String, + name: String, + // tract boundary as WKT (kept so spatial stats can build weights) + geometry_wkt: String, + area_km2: Double + } + + // Phase 2 "Key Equity Variables" - census/ACS + digital-access layers. + schema EquityVars { + geoid: String, + median_income: Double, + pct_rent_burdened: Double, + pct_bipoc: Double, + pct_no_hs_diploma: Double, + pct_limited_english: Double, + pct_zero_vehicle: Double, + pct_internet_subscription: Double + } + + // Phase 3.1 Intrinsic quality (self-referential OSM attributes). + schema AttributeQuality { + geoid: String, + building_density: Double, + road_km_per_km2: Double, + attr_completeness: Double, + // Shannon diversity of amenity/shop types (POI diversity) + poi_diversity: Double, + median_recency_days: Double + } + + // Phase 3.2 Extrinsic quality (OSM vs authoritative references). + schema ExtrinsicQuality { + geoid: String, + footprint_completeness: Double, + positional_accuracy_m: Double, + // false => no authoritative reference for this tract (metrics are N/A) + has_reference: Boolean + } + + // Joined per-tract record. Carries geometry for the spatial analysis. + schema TractQuality { + geoid: String, + geometry_wkt: String, + equity: EquityVars, + intrinsic: AttributeQuality, + extrinsic: ExtrinsicQuality + } + + schema StatResults { + spearman_rho: Double, + spearman_pvalue: Double, + morans_i: Double, + morans_pvalue: Double, + gwr_summary_path: String, + gwr_mean_r2: Double, + temporal_trend: String, + temporal_gap_change: Double + } + + // =================================================================== + // Phase 1 - Study design. + // =================================================================== + event facet DefineStudyArea(region: String, scale: String, hypothesis: String, acs_year: Long) => (area: StudyArea) + + // =================================================================== + // Phase 2 - Data acquisition. Region-level fetches happen ONCE and are + // update-gated: if the extract already exists it is reused unless + // `update` is set. Per-tract OSM is obtained by LOCAL clipping. + // =================================================================== + event facet ResolveStudyTracts(area: StudyArea) => (tracts: [TractGeo]) + event facet FetchRegionOSM(region: String, update: Boolean) => (osm_path: String, snapshot: String) + event facet FetchRegionFootprints(region: String, benchmark_source: String, update: Boolean) => (footprints_path: String, available: Boolean) + + // =================================================================== + // Phase 3 - Per-tract metric calculation (local clip + compute). + // =================================================================== + event facet ClipTractOSM(region_osm_path: String, geometry_wkt: String) => (osm_path: String) + event facet FetchCensusEquity(geoid: String, acs_year: Long) => (vars: EquityVars) + event facet ComputeAttributeQuality(osm_path: String, area_km2: Double, geoid: String) => (quality: AttributeQuality) + event facet ComputeExtrinsicQuality(osm_path: String, geometry_wkt: String, footprints_path: String, reference_available: Boolean) => (quality: ExtrinsicQuality) + event facet AssembleTractQuality(geoid: String, geometry_wkt: String, equity: EquityVars, intrinsic: AttributeQuality, extrinsic: ExtrinsicQuality) => (record: TractQuality) + + // Fan-out: one parallel branch per census tract. Each branch clips the + // shared region extract locally, pulls census equity, computes intrinsic + // + extrinsic quality, and appends one joined TractQuality row. + facet CollectTractQuality(tracts: [TractGeo], region_osm_path: String, footprints_path: String, reference_available: Boolean, acs_year: Long) => (records: [TractQuality]) andThen foreach t in $.tracts { + clip = ClipTractOSM(region_osm_path = $.region_osm_path, geometry_wkt = $.t.geometry_wkt) + eq = FetchCensusEquity(geoid = $.t.geoid, acs_year = $.acs_year) + intr = ComputeAttributeQuality(osm_path = clip.osm_path, area_km2 = $.t.area_km2, geoid = $.t.geoid) + extr = ComputeExtrinsicQuality(osm_path = clip.osm_path, geometry_wkt = $.t.geometry_wkt, footprints_path = $.footprints_path, reference_available = $.reference_available) + rec = AssembleTractQuality(geoid = $.t.geoid, geometry_wkt = $.t.geometry_wkt, equity = eq.vars, intrinsic = intr.quality, extrinsic = extr.quality) + yield CollectTractQuality(records = [rec.record]) + } + + // =================================================================== + // Phase 4 - Statistical & spatial analysis across all tracts. + // The four methods are independent, so they run concurrently. + // =================================================================== + event facet SpearmanCorrelation(records: [TractQuality]) => (rho: Double, pvalue: Double) + event facet MoransI(records: [TractQuality]) => (i: Double, pvalue: Double) + event facet GeographicallyWeightedRegression(records: [TractQuality]) => (summary_path: String, mean_r2: Double) + event facet TemporalEvolution(records: [TractQuality], snapshot: String) => (trend: String, gap_change: Double) + event facet CombineStats(spearman_rho: Double, spearman_p: Double, moran_i: Double, moran_p: Double, gwr_summary: String, gwr_r2: Double, temporal_trend: String, temporal_gap_change: Double) => (results: StatResults) + + facet AnalyzeEquity(records: [TractQuality], snapshot: String) => (results: StatResults) andThen { + corr = SpearmanCorrelation(records = $.records) + moran = MoransI(records = $.records) + gwr = GeographicallyWeightedRegression(records = $.records) + temporal = TemporalEvolution(records = $.records, snapshot = $.snapshot) + combined = CombineStats( + spearman_rho = corr.rho, spearman_p = corr.pvalue, + moran_i = moran.i, moran_p = moran.pvalue, + gwr_summary = gwr.summary_path, gwr_r2 = gwr.mean_r2, + temporal_trend = temporal.trend, temporal_gap_change = temporal.gap_change + ) + yield AnalyzeEquity(results = combined.results) + } + + // =================================================================== + // Phase 5 - Actionable reporting: data deserts + map + HTML report. + // =================================================================== + event facet BuildEquityReport(records: [TractQuality], stats: StatResults, title: String) => (report_html: String, map_path: String, data_deserts: [String]) + + // =================================================================== + // Entry point - chains the five phases end to end. + // =================================================================== + workflow DigitalDivideMappingEquity( + region: String = "Oakland, CA", + scale: String = "census_tract", + hypothesis: String = "Lower-income, lower-connectivity tracts show lower OSM building completeness and road density.", + acs_year: Long = 2022, + update: Boolean = false, + benchmark_source: String = "microsoft_open_buildings", + title: String = "Digital Divide and OSM Mapping Equity" + ) => (report_html: String, map_path: String, data_deserts: [String]) andThen { + + // Phase 1 + area = DefineStudyArea(region = $.region, scale = $.scale, hypothesis = $.hypothesis, acs_year = $.acs_year) + + // Phase 2 - region-level acquisition (update-gated), run concurrently + units = ResolveStudyTracts(area = area.area) + osm = FetchRegionOSM(region = $.region, update = $.update) + foot = FetchRegionFootprints(region = $.region, benchmark_source = $.benchmark_source, update = $.update) + + // Phases 2 + 3 - per-tract clip + metrics (fanned out) + perTract = CollectTractQuality( + tracts = units.tracts, + region_osm_path = osm.osm_path, + footprints_path = foot.footprints_path, + reference_available = foot.available, + acs_year = $.acs_year + ) + + // Phase 4 + stats = AnalyzeEquity(records = perTract.records, snapshot = osm.snapshot) + + // Phase 5 + report = BuildEquityReport(records = perTract.records, stats = stats.results, title = $.title) + + yield DigitalDivideMappingEquity( + report_html = report.report_html, + map_path = report.map_path, + data_deserts = report.data_deserts + ) + } + +} diff --git a/examples/osm-equity/handlers/__init__.py b/examples/osm-equity/handlers/__init__.py new file mode 100644 index 00000000..21fd20da --- /dev/null +++ b/examples/osm-equity/handlers/__init__.py @@ -0,0 +1,33 @@ +"""osm.equity mapping-equity handlers - registration aggregator.""" + +from __future__ import annotations + +from .acquisition.acquisition_handlers import register_acquisition_handlers +from .analysis.analysis_handlers import register_analysis_handlers +from .design.design_handlers import register_design_handlers +from .metrics.metrics_handlers import register_metrics_handlers +from .reporting.reporting_handlers import register_reporting_handlers + + +def register_all_handlers(poller) -> None: + """Register every handler with an AgentPoller.""" + register_design_handlers(poller) + register_acquisition_handlers(poller) + register_metrics_handlers(poller) + register_analysis_handlers(poller) + register_reporting_handlers(poller) + + +def register_all_registry_handlers(runner) -> None: + """Register every handler with a RegistryRunner.""" + from .acquisition.acquisition_handlers import register_handlers as reg_acq + from .analysis.analysis_handlers import register_handlers as reg_ana + from .design.design_handlers import register_handlers as reg_des + from .metrics.metrics_handlers import register_handlers as reg_met + from .reporting.reporting_handlers import register_handlers as reg_rep + + reg_des(runner) + reg_acq(runner) + reg_met(runner) + reg_ana(runner) + reg_rep(runner) diff --git a/examples/osm-equity/handlers/acquisition/__init__.py b/examples/osm-equity/handlers/acquisition/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/handlers/acquisition/acquisition_handlers.py b/examples/osm-equity/handlers/acquisition/acquisition_handlers.py new file mode 100644 index 00000000..f4fb8726 --- /dev/null +++ b/examples/osm-equity/handlers/acquisition/acquisition_handlers.py @@ -0,0 +1,84 @@ +"""Phase 2 - data acquisition: ResolveStudyTracts, FetchRegionOSM, +FetchRegionFootprints. Region-level fetches are update-gated (reuse the +cached extract unless `update` is set).""" + +from __future__ import annotations + +import os +from typing import Any + +from handlers.shared import equity_utils as U + +NAMESPACE = "osm.equity" + + +def _truthy(v: Any) -> bool: + if isinstance(v, str): + return v.strip().lower() in {"1", "true", "yes", "on"} + return bool(v) + + +def _log(params, msg, level="success"): + sl = params.get("_step_log") + if sl is not None: + sl.append({"message": msg, "level": level}) + + +def handle_resolve_study_tracts(params: dict[str, Any]) -> dict[str, Any]: + area = params.get("area", {}) + if isinstance(area, str): + import json + + area = json.loads(area) + region = area.get("region", "Oakland, CA") + tracts = U.make_tracts(region) + _log(params, f"Resolved {len(tracts)} census tracts for {region}") + return {"tracts": tracts} + + +def handle_fetch_region_osm(params: dict[str, Any]) -> dict[str, Any]: + region = params.get("region", "Oakland, CA") + update = _truthy(params.get("update", False)) + tracts = U.make_tracts(region) + path, snapshot, cached = U.fetch_region_osm(region, update, tracts) + _log(params, f"Region OSM {'reused (cache)' if cached else 'fetched'} -> {snapshot}") + return {"osm_path": path, "snapshot": snapshot} + + +def handle_fetch_region_footprints(params: dict[str, Any]) -> dict[str, Any]: + region = params.get("region", "Oakland, CA") + source = params.get("benchmark_source", "microsoft_open_buildings") + update = _truthy(params.get("update", False)) + tracts = U.make_tracts(region) + path, available = U.fetch_region_footprints(region, source, update, tracts) + _log( + params, + f"Footprint reference '{source}': {'available' if available else 'NOT available (extrinsic metrics will be N/A)'}", + level="success" if available else "warning", + ) + return {"footprints_path": path, "available": available} + + +_DISPATCH: dict[str, Any] = { + f"{NAMESPACE}.ResolveStudyTracts": handle_resolve_study_tracts, + f"{NAMESPACE}.FetchRegionOSM": handle_fetch_region_osm, + f"{NAMESPACE}.FetchRegionFootprints": handle_fetch_region_footprints, +} + + +def handle(payload: dict) -> dict: + return _DISPATCH[payload["_facet_name"]](payload) + + +def register_handlers(runner) -> None: + for facet_name in _DISPATCH: + runner.register_handler( + facet_name=facet_name, + module_uri=f"file://{os.path.abspath(__file__)}", + entrypoint="handle", + ) + + +def register_acquisition_handlers(poller) -> None: + for facet_name, fn in _DISPATCH.items(): + poller.register(facet_name, fn) diff --git a/examples/osm-equity/handlers/analysis/__init__.py b/examples/osm-equity/handlers/analysis/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/handlers/analysis/analysis_handlers.py b/examples/osm-equity/handlers/analysis/analysis_handlers.py new file mode 100644 index 00000000..83b76da6 --- /dev/null +++ b/examples/osm-equity/handlers/analysis/analysis_handlers.py @@ -0,0 +1,72 @@ +"""Phase 4 - statistics & spatial analysis: SpearmanCorrelation, MoransI, +GeographicallyWeightedRegression, TemporalEvolution, CombineStats.""" + +from __future__ import annotations + +import os +from typing import Any + +from handlers.shared import equity_utils as U + +NAMESPACE = "osm.equity" + + +def handle_spearman(params: dict[str, Any]) -> dict[str, Any]: + rho, p = U.spearman(params.get("records", [])) + return {"rho": rho, "pvalue": p} + + +def handle_morans_i(params: dict[str, Any]) -> dict[str, Any]: + i, p = U.morans_i(params.get("records", [])) + return {"i": i, "pvalue": p} + + +def handle_gwr(params: dict[str, Any]) -> dict[str, Any]: + path, r2 = U.gwr(params.get("records", [])) + return {"summary_path": path, "mean_r2": r2} + + +def handle_temporal(params: dict[str, Any]) -> dict[str, Any]: + trend, gap = U.temporal_evolution(params.get("records", []), params.get("snapshot", "")) + return {"trend": trend, "gap_change": gap} + + +def handle_combine_stats(params: dict[str, Any]) -> dict[str, Any]: + results = { + "spearman_rho": float(params.get("spearman_rho", 0.0)), + "spearman_pvalue": float(params.get("spearman_p", 1.0)), + "morans_i": float(params.get("moran_i", 0.0)), + "morans_pvalue": float(params.get("moran_p", 1.0)), + "gwr_summary_path": params.get("gwr_summary", ""), + "gwr_mean_r2": float(params.get("gwr_r2", 0.0)), + "temporal_trend": params.get("temporal_trend", "n/a"), + "temporal_gap_change": float(params.get("temporal_gap_change", 0.0)), + } + return {"results": results} + + +_DISPATCH: dict[str, Any] = { + f"{NAMESPACE}.SpearmanCorrelation": handle_spearman, + f"{NAMESPACE}.MoransI": handle_morans_i, + f"{NAMESPACE}.GeographicallyWeightedRegression": handle_gwr, + f"{NAMESPACE}.TemporalEvolution": handle_temporal, + f"{NAMESPACE}.CombineStats": handle_combine_stats, +} + + +def handle(payload: dict) -> dict: + return _DISPATCH[payload["_facet_name"]](payload) + + +def register_handlers(runner) -> None: + for facet_name in _DISPATCH: + runner.register_handler( + facet_name=facet_name, + module_uri=f"file://{os.path.abspath(__file__)}", + entrypoint="handle", + ) + + +def register_analysis_handlers(poller) -> None: + for facet_name, fn in _DISPATCH.items(): + poller.register(facet_name, fn) diff --git a/examples/osm-equity/handlers/design/__init__.py b/examples/osm-equity/handlers/design/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/handlers/design/design_handlers.py b/examples/osm-equity/handlers/design/design_handlers.py new file mode 100644 index 00000000..f69b0246 --- /dev/null +++ b/examples/osm-equity/handlers/design/design_handlers.py @@ -0,0 +1,48 @@ +"""Phase 1 - study design: DefineStudyArea.""" + +from __future__ import annotations + +import os +from typing import Any + +NAMESPACE = "osm.equity" + + +def handle_define_study_area(params: dict[str, Any]) -> dict[str, Any]: + area = { + "region": params.get("region", "Oakland, CA"), + "scale": params.get("scale", "census_tract"), + "hypothesis": params.get("hypothesis", ""), + "acs_year": int(params.get("acs_year", 2022)), + } + _log(params, f"Study area defined: {area['region']} ({area['scale']}, ACS {area['acs_year']})") + return {"area": area} + + +def _log(params, msg, level="success"): + sl = params.get("_step_log") + if sl is not None: + sl.append({"message": msg, "level": level}) + + +_DISPATCH: dict[str, Any] = { + f"{NAMESPACE}.DefineStudyArea": handle_define_study_area, +} + + +def handle(payload: dict) -> dict: + return _DISPATCH[payload["_facet_name"]](payload) + + +def register_handlers(runner) -> None: + for facet_name in _DISPATCH: + runner.register_handler( + facet_name=facet_name, + module_uri=f"file://{os.path.abspath(__file__)}", + entrypoint="handle", + ) + + +def register_design_handlers(poller) -> None: + for facet_name, fn in _DISPATCH.items(): + poller.register(facet_name, fn) diff --git a/examples/osm-equity/handlers/metrics/__init__.py b/examples/osm-equity/handlers/metrics/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/handlers/metrics/metrics_handlers.py b/examples/osm-equity/handlers/metrics/metrics_handlers.py new file mode 100644 index 00000000..aa3be38b --- /dev/null +++ b/examples/osm-equity/handlers/metrics/metrics_handlers.py @@ -0,0 +1,88 @@ +"""Phase 3 - per-tract metrics: ClipTractOSM, FetchCensusEquity, +ComputeAttributeQuality, ComputeExtrinsicQuality, AssembleTractQuality.""" + +from __future__ import annotations + +import json +import os +from typing import Any + +from handlers.shared import equity_utils as U + +NAMESPACE = "osm.equity" + + +def _d(v): + return json.loads(v) if isinstance(v, str) else v + + +def _truthy(v: Any) -> bool: + if isinstance(v, str): + return v.strip().lower() in {"1", "true", "yes", "on"} + return bool(v) + + +def handle_clip_tract_osm(params: dict[str, Any]) -> dict[str, Any]: + path = U.clip_tract_osm(params["region_osm_path"], params["geometry_wkt"]) + return {"osm_path": path} + + +def handle_fetch_census_equity(params: dict[str, Any]) -> dict[str, Any]: + vars_ = U.fetch_census_equity(params["geoid"], int(params.get("acs_year", 2022))) + return {"vars": vars_} + + +def handle_compute_attribute_quality(params: dict[str, Any]) -> dict[str, Any]: + q = U.compute_attribute_quality( + params["osm_path"], float(params.get("area_km2", 1.0)), params.get("geoid", "") + ) + return {"quality": q} + + +def handle_compute_extrinsic_quality(params: dict[str, Any]) -> dict[str, Any]: + q = U.compute_extrinsic_quality( + params["osm_path"], + params["geometry_wkt"], + params.get("footprints_path", ""), + _truthy(params.get("reference_available", False)), + params.get("geoid", ""), + ) + return {"quality": q} + + +def handle_assemble_tract_quality(params: dict[str, Any]) -> dict[str, Any]: + record = { + "geoid": params["geoid"], + "geometry_wkt": params["geometry_wkt"], + "equity": _d(params["equity"]), + "intrinsic": _d(params["intrinsic"]), + "extrinsic": _d(params["extrinsic"]), + } + return {"record": record} + + +_DISPATCH: dict[str, Any] = { + f"{NAMESPACE}.ClipTractOSM": handle_clip_tract_osm, + f"{NAMESPACE}.FetchCensusEquity": handle_fetch_census_equity, + f"{NAMESPACE}.ComputeAttributeQuality": handle_compute_attribute_quality, + f"{NAMESPACE}.ComputeExtrinsicQuality": handle_compute_extrinsic_quality, + f"{NAMESPACE}.AssembleTractQuality": handle_assemble_tract_quality, +} + + +def handle(payload: dict) -> dict: + return _DISPATCH[payload["_facet_name"]](payload) + + +def register_handlers(runner) -> None: + for facet_name in _DISPATCH: + runner.register_handler( + facet_name=facet_name, + module_uri=f"file://{os.path.abspath(__file__)}", + entrypoint="handle", + ) + + +def register_metrics_handlers(poller) -> None: + for facet_name, fn in _DISPATCH.items(): + poller.register(facet_name, fn) diff --git a/examples/osm-equity/handlers/reporting/__init__.py b/examples/osm-equity/handlers/reporting/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/handlers/reporting/reporting_handlers.py b/examples/osm-equity/handlers/reporting/reporting_handlers.py new file mode 100644 index 00000000..1595a9ae --- /dev/null +++ b/examples/osm-equity/handlers/reporting/reporting_handlers.py @@ -0,0 +1,55 @@ +"""Phase 5 - actionable reporting: BuildEquityReport.""" + +from __future__ import annotations + +import json +import os +from typing import Any + +from handlers.shared import equity_utils as U + +NAMESPACE = "osm.equity" + + +def handle_build_equity_report(params: dict[str, Any]) -> dict[str, Any]: + records = params.get("records", []) + if isinstance(records, str): + records = json.loads(records) + stats_res = params.get("stats", {}) + if isinstance(stats_res, str): + stats_res = json.loads(stats_res) + title = params.get("title", "Digital Divide and OSM Mapping Equity") + report_html, map_path, deserts = U.build_report(records, stats_res, title) + + sl = params.get("_step_log") + if sl is not None: + sl.append( + { + "message": f"Report built: {len(deserts)} data deserts across {len(records)} tracts", + "level": "success", + } + ) + return {"report_html": report_html, "map_path": map_path, "data_deserts": deserts} + + +_DISPATCH: dict[str, Any] = { + f"{NAMESPACE}.BuildEquityReport": handle_build_equity_report, +} + + +def handle(payload: dict) -> dict: + return _DISPATCH[payload["_facet_name"]](payload) + + +def register_handlers(runner) -> None: + for facet_name in _DISPATCH: + runner.register_handler( + facet_name=facet_name, + module_uri=f"file://{os.path.abspath(__file__)}", + entrypoint="handle", + ) + + +def register_reporting_handlers(poller) -> None: + for facet_name, fn in _DISPATCH.items(): + poller.register(facet_name, fn) diff --git a/examples/osm-equity/handlers/shared/__init__.py b/examples/osm-equity/handlers/shared/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/handlers/shared/equity_utils.py b/examples/osm-equity/handlers/shared/equity_utils.py new file mode 100644 index 00000000..3aa78906 --- /dev/null +++ b/examples/osm-equity/handlers/shared/equity_utils.py @@ -0,0 +1,718 @@ +"""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). +""" + +from __future__ import annotations + +import hashlib +import json +import math +import os +import tempfile +from pathlib import Path +from typing import Any + +import numpy as np +from scipy import stats +from shapely import wkt as shapely_wkt +from shapely.geometry import Point, mapping, shape + +# Amenity vocabulary used to synthesise + score POI diversity. +_AMENITY_TYPES = [ + "cafe", + "restaurant", + "pharmacy", + "school", + "clinic", + "bank", + "library", + "supermarket", + "bar", + "fuel", + "hospital", + "post_office", +] +# Attribute-completeness tags scored per POI (Phase 3.1). +_KEY_TAGS = ["name", "addr:street", "opening_hours"] + +# Rough study-region bounding boxes (lon_min, lat_min, lon_max, lat_max). +# TODO(real): resolve via a geocoder / TIGER PLACE lookup. +_REGION_BBOX = { + "oakland, ca": (-122.34, 37.72, -122.15, 37.85), + "detroit, mi": (-83.29, 42.26, -82.91, 42.45), + "atlanta, ga": (-84.55, 33.65, -84.29, 33.89), + "california": (-124.41, 32.53, -114.13, 42.01), + "los angeles, ca": (-118.67, 33.70, -118.16, 34.34), + "san francisco, ca": (-122.52, 37.70, -122.36, 37.83), +} +_DEFAULT_BBOX = (-122.34, 37.72, -122.15, 37.85) + + +# --------------------------------------------------------------------------- +# Paths / IO +# --------------------------------------------------------------------------- + + +def _slug(region: str) -> str: + return "".join(c if c.isalnum() else "_" for c in region.strip().lower()) + + +def out_dir() -> Path: + base = ( + os.environ.get("FW_CACHE_ROOT") or os.environ.get("FW_OUTPUT_BASE") or tempfile.gettempdir() + ) + d = Path(base) / "osm_equity" + d.mkdir(parents=True, exist_ok=True) + return d + + +def _rng(*parts: Any) -> np.random.Generator: + """A deterministic numpy Generator seeded from the given parts.""" + h = hashlib.sha256("::".join(str(p) for p in parts).encode()).digest() + return np.random.default_rng(int.from_bytes(h[:8], "big")) + + +def _write_geojson(path: Path, features: list[dict]) -> None: + path.write_text(json.dumps({"type": "FeatureCollection", "features": features})) + + +def _read_geojson(path: str) -> list[dict]: + return json.loads(Path(path).read_text()).get("features", []) + + +def grid_size() -> int: + try: + return max(3, int(os.environ.get("FW_EQUITY_GRID", "6"))) + except ValueError: + return 6 + + +# --------------------------------------------------------------------------- +# Phase 1 / 2 - study area, tracts, region OSM, footprints +# --------------------------------------------------------------------------- + + +def region_bbox(region: str) -> tuple[float, float, float, float]: + return _REGION_BBOX.get(region.strip().lower(), _DEFAULT_BBOX) + + +def tract_profile(geoid: str, region: str | None = None) -> dict[str, float]: + """Deterministic per-tract latent profile. + + Income rises smoothly toward the NE corner of the grid (giving Moran's I + a real positive signal) plus a small seeded jitter. Home-internet + subscription tracks income; the two together define the "connectivity" + that drives OSM richness. + """ + n = grid_size() + i, j = tract_ij(geoid) + # smooth spatial gradient in [0, 1] + grad = ((i + j) / (2.0 * (n - 1))) if n > 1 else 0.5 + jitter = float(_rng(geoid, "income").normal(0.0, 0.06)) + wealth = min(1.0, max(0.0, grad + jitter)) + median_income = 28000.0 + wealth * 145000.0 + pct_internet = 0.45 + wealth * 0.5 + float(_rng(geoid, "net").normal(0, 0.03)) + pct_internet = min(0.99, max(0.2, pct_internet)) + return { + "wealth": wealth, + "median_income": round(median_income, 1), + "pct_internet_subscription": round(pct_internet, 4), + } + + +def tract_ij(geoid: str) -> tuple[int, int]: + """Grid coordinates encoded in the last 4 chars of the geoid.""" + tail = geoid[-4:] + return int(tail[:2]), int(tail[2:]) + + +def make_tracts(region: str) -> list[dict[str, Any]]: + """Tile the region bbox into a grid of tract polygons (Phase 2 units). + + TODO(real): replace with a Census TIGER/Line tract fetch for the region. + """ + n = grid_size() + lon0, lat0, lon1, lat1 = region_bbox(region) + dlon = (lon1 - lon0) / n + dlat = (lat1 - lat0) / n + # crude km-per-degree at this latitude for area + midlat = (lat0 + lat1) / 2.0 + km_per_deg_lat = 110.574 + km_per_deg_lon = 111.320 * math.cos(math.radians(midlat)) + cell_km2 = abs(dlon * km_per_deg_lon) * abs(dlat * km_per_deg_lat) + tracts = [] + for i in range(n): + for j in range(n): + x0, y0 = lon0 + i * dlon, lat0 + j * dlat + x1, y1 = x0 + dlon, y0 + dlat + poly_wkt = f"POLYGON (({x0} {y0}, {x1} {y0}, {x1} {y1}, {x0} {y1}, {x0} {y0}))" + geoid = f"06001{i:02d}{j:02d}" + tracts.append( + { + "geoid": geoid, + "name": f"Tract {i}.{j}", + "geometry_wkt": poly_wkt, + "area_km2": round(cell_km2, 4), + } + ) + return tracts + + +def generate_region_osm(region: str, tracts: list[dict]) -> list[dict]: + """Synthesize OSM features for the whole region as GeoJSON. + + Feature density, tag completeness, POI diversity and edit recency all + scale with the tract's connectivity - the digital-divide signal the + study is designed to detect. + + TODO(real): replace with an Overpass/ohsome bbox pull (buildings, highways, + amenities) with full metadata (timestamps) - fetched ONCE for the region. + """ + features: list[dict] = [] + for t in tracts: + geoid = t["geoid"] + prof = tract_profile(geoid, region) + wealth = prof["wealth"] + poly = shapely_wkt.loads(t["geometry_wkt"]) + x0, y0, x1, y1 = poly.bounds + rng = _rng(geoid, "osm") + # counts scale with connectivity + n_buildings = int(15 + wealth * 120 + rng.integers(0, 12)) + n_roads = int(3 + wealth * 12) + n_pois = int(2 + wealth * 30) + # buildings + for _ in range(n_buildings): + px = float(rng.uniform(x0, x1)) + py = float(rng.uniform(y0, y1)) + features.append( + { + "type": "Feature", + "properties": { + "geoid": geoid, + "osm_kind": "building", + "building": "yes", + # wealthier tracts edited more recently + "recency_days": int(rng.uniform(20, 400) * (1.6 - wealth)), + }, + "geometry": mapping(Point(px, py)), + } + ) + # roads as short linestrings + for _ in range(n_roads): + px = float(rng.uniform(x0, x1)) + py = float(rng.uniform(y0, y1)) + qx = px + float(rng.uniform(-0.004, 0.004)) + qy = py + float(rng.uniform(-0.004, 0.004)) + features.append( + { + "type": "Feature", + "properties": { + "geoid": geoid, + "osm_kind": "road", + "highway": "residential", + "recency_days": int(rng.uniform(20, 400) * (1.6 - wealth)), + }, + "geometry": {"type": "LineString", "coordinates": [[px, py], [qx, qy]]}, + } + ) + # POIs - richer, better-tagged, more diverse in connected tracts + n_types = max(1, int(1 + wealth * (len(_AMENITY_TYPES) - 1))) + for k in range(n_pois): + px = float(rng.uniform(x0, x1)) + py = float(rng.uniform(y0, y1)) + amenity = _AMENITY_TYPES[int(rng.integers(0, n_types))] + props = { + "geoid": geoid, + "osm_kind": "poi", + "amenity": amenity, + "recency_days": int(rng.uniform(20, 400) * (1.6 - wealth)), + } + # tag completeness scales with connectivity + for tag in _KEY_TAGS: + if rng.random() < 0.35 + wealth * 0.6: + props[tag] = f"{amenity}-{k}" if tag == "name" else "set" + features.append( + {"type": "Feature", "properties": props, "geometry": mapping(Point(px, py))} + ) + return features + + +def fetch_region_osm(region: str, update: bool, tracts: list[dict]) -> tuple[str, str, bool]: + """Update-gated region OSM fetch. Returns (path, snapshot, was_cached). + + If the extract already exists it is reused unless `update` is set - + exactly the "only fetch if the update flag is set" behaviour, and the + reason region OSM is pulled once rather than per tract. + """ + path = out_dir() / f"osm_{_slug(region)}.geojson" + 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] + snap_path.write_text(snapshot) + return str(path), snapshot, False + + +def fetch_region_footprints( + region: str, source: str, update: bool, tracts: list[dict] +) -> tuple[str, bool]: + """Update-gated authoritative-footprint fetch (extrinsic benchmark). + + Returns (path, available). `available` is false when no reference layer + 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. + """ + # A deployment might have no benchmark for some regions/sources. + available = source.strip().lower() in { + "microsoft_open_buildings", + "google_open_buildings", + "municipal", + } + path = out_dir() / f"footprints_{_slug(region)}.geojson" + if not available: + return "", False + if path.exists() and not update: + return str(path), True + # reference building count per tract = "true" count the OSM sample under-covers + feats = [] + for t in tracts: + geoid = t["geoid"] + poly = shapely_wkt.loads(t["geometry_wkt"]) + x0, y0, x1, y1 = poly.bounds + rng = _rng(geoid, "ref") + # authoritative count is the "complete" total; OSM completeness = osm/ref + true_count = int(60 + rng.integers(0, 60)) + for _ in range(true_count): + px = float(rng.uniform(x0, x1)) + py = float(rng.uniform(y0, y1)) + feats.append( + { + "type": "Feature", + "properties": {"geoid": geoid, "ref": "building"}, + "geometry": mapping(Point(px, py)), + } + ) + _write_geojson(path, feats) + return str(path), True + + +# --------------------------------------------------------------------------- +# Phase 2/3 - census equity, clip, intrinsic + extrinsic quality +# --------------------------------------------------------------------------- + + +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). + """ + prof = tract_profile(geoid, region) + wealth = prof["wealth"] + rng = _rng(geoid, "acs", acs_year) + + def _pct(base, span, lo=0.0, hi=1.0): + return round(min(hi, max(lo, base + span + float(rng.normal(0, 0.02)))), 4) + + return { + "geoid": geoid, + "median_income": prof["median_income"], + # disadvantage indicators fall as wealth rises + "pct_rent_burdened": _pct(0.55, -wealth * 0.35), + "pct_bipoc": _pct(0.75, -wealth * 0.5), + "pct_no_hs_diploma": _pct(0.32, -wealth * 0.28), + "pct_limited_english": _pct(0.24, -wealth * 0.2), + "pct_zero_vehicle": _pct(0.30, -wealth * 0.24), + "pct_internet_subscription": prof["pct_internet_subscription"], + } + + +def clip_tract_osm(region_osm_path: str, geometry_wkt: str) -> str: + """Clip the shared region extract to a tract polygon - LOCAL, no network.""" + poly = shapely_wkt.loads(geometry_wkt) + 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) + tag = hashlib.sha256((region_osm_path + geometry_wkt).encode()).hexdigest()[:10] + path = out_dir() / f"tract_{tag}.geojson" + _write_geojson(path, kept) + return str(path) + + +def _haversine_km(a: tuple[float, float], b: tuple[float, float]) -> float: + lon1, lat1, lon2, lat2 = map(math.radians, [a[0], a[1], b[0], b[1]]) + dlon, dlat = lon2 - lon1, lat2 - lat1 + h = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2 + return 2 * 6371.0 * math.asin(math.sqrt(h)) + + +def compute_attribute_quality(osm_path: str, area_km2: float, geoid: str) -> dict[str, Any]: + """Phase 3.1 intrinsic metrics from the tract's OSM features.""" + feats = _read_geojson(osm_path) + buildings = [f for f in feats if f["properties"].get("osm_kind") == "building"] + roads = [f for f in feats if f["properties"].get("osm_kind") == "road"] + pois = [f for f in feats if f["properties"].get("osm_kind") == "poi"] + + area = max(area_km2, 1e-6) + building_density = len(buildings) / area + + road_km = 0.0 + for r in roads: + coords = r["geometry"]["coordinates"] + road_km += _haversine_km(tuple(coords[0]), tuple(coords[-1])) + road_km_per_km2 = road_km / area + + # attribute completeness: mean fraction of key tags present across POIs + if pois: + completeness = float( + np.mean( + [sum(1 for t in _KEY_TAGS if p["properties"].get(t)) / len(_KEY_TAGS) for p in pois] + ) + ) + else: + completeness = 0.0 + + # POI diversity: Shannon entropy over amenity types, normalised to [0,1] + poi_diversity = 0.0 + if pois: + counts: dict[str, int] = {} + for p in pois: + a = p["properties"].get("amenity", "other") + counts[a] = counts.get(a, 0) + 1 + total = sum(counts.values()) + ent = -sum((c / total) * math.log(c / total) for c in counts.values()) + poi_diversity = ent / math.log(len(_AMENITY_TYPES)) + + recencies = [ + f["properties"].get("recency_days") + for f in feats + if f["properties"].get("recency_days") is not None + ] + median_recency = float(np.median(recencies)) if recencies else 0.0 + + return { + "geoid": geoid, + "building_density": round(building_density, 3), + "road_km_per_km2": round(road_km_per_km2, 4), + "attr_completeness": round(completeness, 4), + "poi_diversity": round(poi_diversity, 4), + "median_recency_days": round(median_recency, 1), + } + + +def compute_extrinsic_quality( + osm_path: str, geometry_wkt: str, footprints_path: str, reference_available: bool, geoid: str +) -> dict[str, Any]: + """Phase 3.2 extrinsic metrics - GATED on reference availability. + + When there is no authoritative footprint layer, has_reference is false + and the numeric metrics are emitted as -1.0 sentinels (explicitly N/A), + never silently zeroed. + """ + if not reference_available or not footprints_path: + return { + "geoid": geoid, + "footprint_completeness": -1.0, + "positional_accuracy_m": -1.0, + "has_reference": False, + } + 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"])) + ] + ref_n = len(ref_buildings) + 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] + 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) + positional = float(np.mean(dists)) + else: + positional = -1.0 + return { + "geoid": geoid, + "footprint_completeness": round(min(completeness, 1.0), 4), + "positional_accuracy_m": round(positional, 2), + "has_reference": True, + } + + +# --------------------------------------------------------------------------- +# Phase 4 - statistics & spatial analysis +# --------------------------------------------------------------------------- + + +def _as_list(records: Any) -> list[dict]: + if isinstance(records, str): + records = json.loads(records) + return list(records or []) + + +def _centroids(records: list[dict]) -> np.ndarray: + pts = [] + for r in records: + c = shapely_wkt.loads(r["geometry_wkt"]).centroid + pts.append([c.x, c.y]) + return np.array(pts) + + +def _income(records: list[dict]) -> np.ndarray: + return np.array([float(r["equity"]["median_income"]) for r in records]) + + +def spearman(records: list[dict]) -> tuple[float, float]: + """Phase 4.1 - Spearman rank correlation: income vs POI diversity.""" + recs = _as_list(records) + inc = _income(recs) + div = np.array([float(r["intrinsic"]["poi_diversity"]) for r in recs]) + if len(recs) < 3: + return 0.0, 1.0 + rho, p = stats.spearmanr(inc, div) + return float(round(rho, 4)), float(round(p, 5)) + + +def _weights(coords: np.ndarray, k: int = 4) -> np.ndarray: + """Row-standardised k-nearest-neighbour spatial weights matrix.""" + n = len(coords) + W = np.zeros((n, n)) + for i in range(n): + d = np.array( + [ + _haversine_km(tuple(coords[i]), tuple(coords[j])) if j != i else np.inf + for j in range(n) + ] + ) + nn = np.argsort(d)[: min(k, n - 1)] + W[i, nn] = 1.0 + rs = W.sum(axis=1, keepdims=True) + rs[rs == 0] = 1.0 + return W / rs + + +def morans_i( + records: list[dict], value: str = "attr_completeness", perms: int = 199 +) -> tuple[float, float]: + """Phase 4.2 - Global Moran's I on a quality metric, with a seeded + permutation p-value. Detects whether 'data deserts' cluster spatially.""" + recs = _as_list(records) + n = len(recs) + if n < 4: + return 0.0, 1.0 + coords = _centroids(recs) + x = np.array([float(r["intrinsic"][value]) for r in recs]) + W = _weights(coords) + z = x - x.mean() + denom = (z**2).sum() + + def _I(zv): + num = (W * np.outer(zv, zv)).sum() + return (n / W.sum()) * (num / denom) if denom else 0.0 + + obs = _I(z) + rng = _rng("morans", n, round(float(x.sum()), 3)) + ge = 1 + for _ in range(perms): + if _I(rng.permutation(z)) >= obs: + ge += 1 + return float(round(obs, 4)), float(round(ge / (perms + 1), 5)) + + +def gwr(records: list[dict]) -> tuple[str, float]: + """Phase 4.3 - a simplified geographically weighted regression of + attribute completeness ~ income, fit locally with a Gaussian kernel over + each tract's neighbourhood. Returns (summary_path, mean_local_r2). + + TODO(real): swap in mgwr.gwr.GWR with an adaptive-bandwidth selector.""" + recs = _as_list(records) + n = len(recs) + coords = _centroids(recs) + x = _income(recs) + y = np.array([float(r["intrinsic"]["attr_completeness"]) for r in recs]) + x = (x - x.mean()) / (x.std() or 1.0) + rows = [] + r2s = [] + # bandwidth = median pairwise distance + dmat = np.array( + [[_haversine_km(tuple(coords[i]), tuple(coords[j])) for j in range(n)] for i in range(n)] + ) + bw = float(np.median(dmat[dmat > 0])) or 1.0 + for i in range(n): + w = np.exp(-(dmat[i] ** 2) / (2 * bw**2)) + W = np.diag(w) + X = np.column_stack([np.ones(n), x]) + try: + beta = np.linalg.solve(X.T @ W @ X, X.T @ W @ y) + except np.linalg.LinAlgError: + continue + pred = X @ beta + ss_res = (w * (y - pred) ** 2).sum() + ss_tot = (w * (y - np.average(y, weights=w)) ** 2).sum() + local_r2 = 1.0 - (ss_res / ss_tot) if ss_tot else 0.0 + r2s.append(local_r2) + rows.append( + { + "geoid": recs[i]["geoid"], + "intercept": round(float(beta[0]), 4), + "income_coef": round(float(beta[1]), 4), + "local_r2": round(float(local_r2), 4), + } + ) + mean_r2 = float(np.mean(r2s)) if r2s else 0.0 + path = out_dir() / "gwr_summary.json" + path.write_text( + json.dumps( + {"bandwidth_km": round(bw, 3), "mean_local_r2": round(mean_r2, 4), "local_fits": rows}, + indent=2, + ) + ) + return str(path), round(mean_r2, 4) + + +def temporal_evolution(records: list[dict], snapshot: str) -> tuple[str, float]: + """Phase 4.4 - is the mapping divide widening? Uses edit recency as the + temporal signal: compares mean feature age in the poorest vs richest + income tercile. A positive gap (poorer areas older) = a persistent divide. + + TODO(real): drive this from the ohsome quality API time series.""" + recs = _as_list(records) + if len(recs) < 3: + return "insufficient_data", 0.0 + inc = _income(recs) + rec_days = np.array([float(r["intrinsic"]["median_recency_days"]) for r in recs]) + order = np.argsort(inc) + k = max(1, len(recs) // 3) + poor = rec_days[order[:k]].mean() + rich = rec_days[order[-k:]].mean() + gap = float(poor - rich) # days older in poor tracts + trend = "widening" if gap > 30 else ("narrowing" if gap < -30 else "stable") + return trend, round(gap, 1) + + +# --------------------------------------------------------------------------- +# Phase 5 - reporting +# --------------------------------------------------------------------------- + + +def data_deserts(records: list[dict]) -> list[str]: + """Bottom-quartile tracts by a composite completeness score = priority + neighbourhoods where missing maps may hinder emergency/transit routing.""" + recs = _as_list(records) + if not recs: + return [] + scored = [] + for r in recs: + iq = r["intrinsic"] + score = ( + 0.5 * iq["attr_completeness"] + + 0.3 * min(iq["building_density"] / 100.0, 1.0) + + 0.2 * iq["poi_diversity"] + ) + scored.append((r["geoid"], score)) + scored.sort(key=lambda s: s[1]) + cut = max(1, len(scored) // 4) + return [g for g, _ in scored[:cut]] + + +def build_report(records: list[dict], stats_res: dict, title: str) -> tuple[str, str, list[str]]: + """Phase 5 - HTML report + a simple GeoJSON choropleth artifact + + the ranked data-desert list, with intervention recommendations.""" + recs = _as_list(records) + deserts = data_deserts(recs) + n_ref = sum(1 for r in recs if r["extrinsic"].get("has_reference")) + + # choropleth artifact: tract polygons coloured by completeness + feats = [] + for r in recs: + feats.append( + { + "type": "Feature", + "properties": { + "geoid": r["geoid"], + "median_income": r["equity"]["median_income"], + "attr_completeness": r["intrinsic"]["attr_completeness"], + "poi_diversity": r["intrinsic"]["poi_diversity"], + "footprint_completeness": r["extrinsic"]["footprint_completeness"], + "has_reference": r["extrinsic"]["has_reference"], + "data_desert": r["geoid"] in deserts, + }, + "geometry": mapping(shapely_wkt.loads(r["geometry_wkt"])), + } + ) + map_path = out_dir() / "equity_choropleth.geojson" + _write_geojson(map_path, feats) + + rho = stats_res.get("spearman_rho", 0.0) + moran = stats_res.get("morans_i", 0.0) + rows = "".join( + f"{r['geoid']}{r['equity']['median_income']:.0f}" + f"{r['intrinsic']['attr_completeness']:.2f}" + f"{r['intrinsic']['poi_diversity']:.2f}" + f"{'N/A' if not r['extrinsic']['has_reference'] else format(r['extrinsic']['footprint_completeness'], '.2f')}" + f"{'YES' if r['geoid'] in deserts else ''}" + for r in sorted(recs, key=lambda x: x["equity"]["median_income"]) + ) + html = f"""{title} + +

{title}

+

Hypothesis under test: lower-income / lower-connectivity tracts show lower OSM completeness.

+

Phase 4 - findings

+ +

Phase 5 - data deserts ({len(deserts)} of {len(recs)} tracts)

+

Priority tracts: {", ".join(deserts) or "none"}

+

Extrinsic reference coverage: {n_ref}/{len(recs)} tracts had an authoritative footprint layer.

+

Recommended interventions

+ +

Per-tract detail

+ +{rows}
GEOIDMedian incomeAttr completenessPOI diversityFootprint compl.Data desert
+""" + report_path = out_dir() / "equity_report.html" + report_path.write_text(html) + return str(report_path), str(map_path), deserts diff --git a/examples/osm-equity/tests/__init__.py b/examples/osm-equity/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/examples/osm-equity/tests/conftest.py b/examples/osm-equity/tests/conftest.py new file mode 100644 index 00000000..f9487a25 --- /dev/null +++ b/examples/osm-equity/tests/conftest.py @@ -0,0 +1,28 @@ +"""Test-level conftest for osm-equity - handler path setup.""" + +import os +import sys + +import pytest + +_EXAMPLE_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) + + +def _ensure_equity_handlers(): + for key in list(sys.modules.keys()): + if key == "handlers" or key.startswith("handlers."): + mod = sys.modules[key] + mod_file = getattr(mod, "__file__", "") or "" + if "osm-equity" not in mod_file: + del sys.modules[key] + if _EXAMPLE_ROOT in sys.path: + sys.path.remove(_EXAMPLE_ROOT) + sys.path.insert(0, _EXAMPLE_ROOT) + + +_ensure_equity_handlers() + + +@pytest.fixture(autouse=True) +def _equity_handlers_on_path(): + _ensure_equity_handlers() diff --git a/examples/osm-equity/tests/test_osm_equity.py b/examples/osm-equity/tests/test_osm_equity.py new file mode 100644 index 00000000..9e1f0215 --- /dev/null +++ b/examples/osm-equity/tests/test_osm_equity.py @@ -0,0 +1,218 @@ +"""Tests for the osm.equity mapping-equity study. + +Covers three levels: + 1. the FFL compiles + validates (relative-scoping model); + 2. the computational utils are correct and deterministic; + 3. an END-TO-END run driving every handler in the workflow's data order + recovers the built-in digital-divide signal (income correlates with + OSM quality) and produces a report + data-desert list. +""" + +from __future__ import annotations + +import os +import subprocess +import sys + +_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") + + +# --------------------------------------------------------------------------- +# 1. FFL compiles & validates +# --------------------------------------------------------------------------- +class TestFFL: + def test_ffl_validates(self): + r = subprocess.run( + [sys.executable, "-m", "facetwork.cli", _FFL, "--check"], + cwd=_REPO_ROOT, + capture_output=True, + text=True, + ) + assert r.returncode == 0, r.stdout + r.stderr + assert "OK" in (r.stdout + r.stderr) + + def test_ffl_emits_json(self, tmp_path): + out = tmp_path / "equity.json" + r = subprocess.run( + [sys.executable, "-m", "facetwork.cli", _FFL, "-o", str(out)], + cwd=_REPO_ROOT, + capture_output=True, + text=True, + ) + assert r.returncode == 0, r.stdout + r.stderr + assert out.exists() + + +# --------------------------------------------------------------------------- +# 2. Utility correctness / determinism +# --------------------------------------------------------------------------- +class TestUtils: + def test_tracts_tile_the_grid(self): + from handlers.shared import equity_utils as U + + tracts = U.make_tracts("Oakland, CA") + n = U.grid_size() + assert len(tracts) == n * n + assert all("geometry_wkt" in t and t["area_km2"] > 0 for t in tracts) + + def test_income_gradient_is_spatial(self): + from handlers.shared import equity_utils as U + + # NE corner richer than SW corner (smooth gradient => Moran signal) + sw = U.tract_profile("06001" + "0000")["median_income"] + ne = U.tract_profile(f"06001{U.grid_size() - 1:02d}{U.grid_size() - 1:02d}")[ + "median_income" + ] + assert ne > sw + + def test_census_equity_determinism(self): + from handlers.shared import equity_utils as U + + a = U.fetch_census_equity("060010203", 2022) + b = U.fetch_census_equity("060010203", 2022) + assert a == b + assert 0 <= a["pct_internet_subscription"] <= 1 + + def test_update_gating_reuses_cache(self, tmp_path, monkeypatch): + from handlers.shared import equity_utils as U + + monkeypatch.setenv("FW_CACHE_ROOT", str(tmp_path)) + tracts = U.make_tracts("Detroit, MI") + p1, s1, cached1 = U.fetch_region_osm("Detroit, MI", update=False, tracts=tracts) + assert cached1 is False # first fetch builds it + p2, s2, cached2 = U.fetch_region_osm("Detroit, MI", update=False, tracts=tracts) + assert cached2 is True and s2 == s1 # reused, no re-fetch + _, s3, cached3 = U.fetch_region_osm("Detroit, MI", update=True, tracts=tracts) + assert cached3 is False # update flag forces a re-fetch + + def test_extrinsic_gating_no_reference(self): + from handlers.shared import equity_utils as U + + q = U.compute_extrinsic_quality("nonexistent", "POINT (0 0)", "", False, "060010000") + assert q["has_reference"] is False + assert q["footprint_completeness"] == -1.0 # explicit N/A, not silently 0 + + def test_footprint_source_unknown_is_unavailable(self, tmp_path, monkeypatch): + from handlers.shared import equity_utils as U + + monkeypatch.setenv("FW_CACHE_ROOT", str(tmp_path)) + tracts = U.make_tracts("Atlanta, GA") + path, avail = U.fetch_region_footprints("Atlanta, GA", "no_such_source", False, tracts) + assert avail is False and path == "" + + +# --------------------------------------------------------------------------- +# 3. End-to-end: drive every handler in the workflow's data order +# --------------------------------------------------------------------------- +class TestEndToEnd: + def _run( + self, + tmp_path, + monkeypatch, + region="Oakland, CA", + update=False, + benchmark="microsoft_open_buildings", + ): + monkeypatch.setenv("FW_CACHE_ROOT", str(tmp_path)) + from handlers.acquisition import acquisition_handlers as ACQ + from handlers.analysis import analysis_handlers as ANA + from handlers.design import design_handlers as DES + from handlers.metrics import metrics_handlers as MET + from handlers.reporting import reporting_handlers as REP + + # Phase 1 + area = DES.handle_define_study_area( + {"region": region, "scale": "census_tract", "hypothesis": "h", "acs_year": 2022} + )["area"] + # Phase 2 + tracts = ACQ.handle_resolve_study_tracts({"area": area})["tracts"] + osm = ACQ.handle_fetch_region_osm({"region": region, "update": update}) + foot = ACQ.handle_fetch_region_footprints( + {"region": region, "benchmark_source": benchmark, "update": update} + ) + # Phases 2+3 fan-out + records = [] + for t in tracts: + clip = MET.handle_clip_tract_osm( + {"region_osm_path": osm["osm_path"], "geometry_wkt": t["geometry_wkt"]} + ) + eq = MET.handle_fetch_census_equity({"geoid": t["geoid"], "acs_year": 2022})["vars"] + intr = MET.handle_compute_attribute_quality( + {"osm_path": clip["osm_path"], "area_km2": t["area_km2"], "geoid": t["geoid"]} + )["quality"] + extr = MET.handle_compute_extrinsic_quality( + { + "osm_path": clip["osm_path"], + "geometry_wkt": t["geometry_wkt"], + "footprints_path": foot["footprints_path"], + "reference_available": foot["available"], + "geoid": t["geoid"], + } + )["quality"] + rec = MET.handle_assemble_tract_quality( + { + "geoid": t["geoid"], + "geometry_wkt": t["geometry_wkt"], + "equity": eq, + "intrinsic": intr, + "extrinsic": extr, + } + )["record"] + records.append(rec) + # Phase 4 + corr = ANA.handle_spearman({"records": records}) + moran = ANA.handle_morans_i({"records": records}) + gwr = ANA.handle_gwr({"records": records}) + temporal = ANA.handle_temporal({"records": records, "snapshot": osm["snapshot"]}) + stats = ANA.handle_combine_stats( + { + "spearman_rho": corr["rho"], + "spearman_p": corr["pvalue"], + "moran_i": moran["i"], + "moran_p": moran["pvalue"], + "gwr_summary": gwr["summary_path"], + "gwr_r2": gwr["mean_r2"], + "temporal_trend": temporal["trend"], + "temporal_gap_change": temporal["gap_change"], + } + )["results"] + # Phase 5 + report = REP.handle_build_equity_report({"records": records, "stats": stats, "title": "T"}) + return records, stats, report + + def test_full_pipeline_runs_and_reports(self, tmp_path, monkeypatch): + records, stats, report = self._run(tmp_path, monkeypatch) + n = __import__("os") # noqa + assert len(records) >= 9 + assert os.path.exists(report["report_html"]) + assert os.path.exists(report["map_path"]) + html = open(report["report_html"]).read() + assert "Mapping Equity" in html or "Digital Divide" in html or "data desert" in html.lower() + assert isinstance(report["data_deserts"], list) and len(report["data_deserts"]) >= 1 + + def test_digital_divide_signal_recovered(self, tmp_path, monkeypatch): + # Income was built to drive POI diversity + completeness, so the study + # must detect a positive, significant relationship and spatial clustering. + records, stats, _ = self._run(tmp_path, monkeypatch) + assert stats["spearman_rho"] > 0.3 + assert stats["spearman_pvalue"] < 0.05 + assert stats["morans_i"] > 0.0 + + def test_data_deserts_are_low_income(self, tmp_path, monkeypatch): + records, _, report = self._run(tmp_path, monkeypatch) + incomes = {r["geoid"]: r["equity"]["median_income"] for r in records} + overall_median = sorted(incomes.values())[len(incomes) // 2] + desert_incomes = [incomes[g] for g in report["data_deserts"]] + assert sum(desert_incomes) / len(desert_incomes) < overall_median + + def test_extrinsic_available_with_real_source(self, tmp_path, monkeypatch): + records, _, _ = self._run(tmp_path, monkeypatch, benchmark="microsoft_open_buildings") + assert all(r["extrinsic"]["has_reference"] for r in records) + assert all(0.0 <= r["extrinsic"]["footprint_completeness"] <= 1.0 for r in records) + + 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)