1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Geo features cache (flood zone, transit, dissemination-area statistics, CBD distance) for every analytic row.45 python3 scripts/04_geo.py [--synthetic] [--limit N]6→ data/{processed|synthetic}/geo_features.parquet, results/data/geo_summary*.json7Requires data/raw/external/{inondation,transit,quartier}.db (copied from M4M64b:~/apps/immo-ka/data/).8"""9from __future__ import annotations1011import argparse12import json13import time1415import polars as pl1617from agentic_valuation.config import PROCESSED, RESULTS, SYNTHETIC18from agentic_valuation.features.geo import GeoLayers, geo_features_batch192021def main() -> None:22 ap = argparse.ArgumentParser(description=__doc__)23 ap.add_argument("--synthetic", action="store_true")24 ap.add_argument("--limit", type=int, default=None, help="only the first N rows (smoke test)")25 ap.add_argument("--chunk", type=int, default=100_000)26 args = ap.parse_args()27 base = SYNTHETIC if args.synthetic else PROCESSED28 if not GeoLayers.available():29 raise SystemExit("external layers missing: scp M4M64b:~/apps/immo-ka/data/{inondation,transit,quartier}.db data/raw/external/")30 t0 = time.perf_counter()31 layers = GeoLayers.load()32 t_load = time.perf_counter() - t033 analytic = pl.read_parquet(base / "analytic.parquet", columns=["tx_id", "lat", "lng", "x_m", "y_m"])34 if args.limit:35 analytic = analytic.head(args.limit)36 parts = []37 for s0 in range(0, analytic.height, args.chunk):38 parts.append(geo_features_batch(analytic.slice(s0, args.chunk), layers))39 print(f" {min(s0 + args.chunk, analytic.height)}/{analytic.height} rows, {time.perf_counter() - t0:.0f} s", flush=True)40 geo = pl.concat(parts)41 out = base / "geo_features.parquet"42 geo.write_parquet(out, compression="zstd")43 summary = {44 "synthetic": args.synthetic,45 "rows": geo.height,46 "layers_load_seconds": round(t_load, 1),47 "seconds": round(time.perf_counter() - t0, 1),48 "n_flood_polygons": len(layers.flood_geoms),49 "n_coverage_polygons": len(layers.coverage_geoms),50 "n_dissemination_areas": len(layers.da_geoms),51 "n_transit_stops": layers.n_stops,52 "share_flood_mapped": float(geo["flood_mapped"].mean()),53 "share_in_flood_zone_among_mapped": float(geo.filter(pl.col("flood_mapped"))["in_flood_zone"].mean()) if geo["flood_mapped"].any() else None,54 "flood_zone_kinds": {k: int(v) for k, v in geo.group_by("flood_zone_kind").len().iter_rows() if k is not None},55 "share_with_da": float(geo["dauid"].is_not_null().mean()),56 "transit_stops_quantiles": {q: float(geo["transit_stops_within_1km"].quantile(q)) for q in (0.5, 0.75, 0.9)},57 "cbd_counts": {k: int(v) for k, v in geo.group_by("cbd_name").len().iter_rows()},58 "missingness_pct": {c: round(100.0 * geo[c].null_count() / geo.height, 3) for c in geo.columns},59 "output": str(out),60 }61 sp = RESULTS / "data" / ("geo_summary_synthetic.json" if args.synthetic else "geo_summary.json")62 sp.parent.mkdir(parents=True, exist_ok=True)63 sp.write_text(json.dumps(summary, indent=1, ensure_ascii=False))64 print(json.dumps({k: summary[k] for k in ("rows", "layers_load_seconds", "seconds", "share_flood_mapped", "share_in_flood_zone_among_mapped", "share_with_da", "cbd_counts")}, ensure_ascii=False, indent=1))656667if __name__ == "__main__":68 main()69