1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Geospatial context (evidence id G1) from static external layers copied from the Immo-Ka node3(``data/raw/external/``):45* ``inondation.db`` — Quebec flood-zone polygons (``zi``: WKB in EPSG:3857, kinds "grand courant" 0–20 yr,6 "faible courant" 20–100 yr, …) and the mapped coverage areas (``couverture``). ``in_flood_zone`` is True/False7 only where a flood map exists; None (with a note) outside the mapped coverage.8* ``transit.db`` — public-transit stops (bus/metro/REM/train, WGS84): count within 1 km.9* ``quartier.db`` — Statistics Canada 2021 dissemination-area polygons (GeoJSON rings, WGS84) and profile10 statistics (population, density, median age, median income, share of renters).11* distance to the nearest major CBD (Montréal, Québec, Gatineau, Sherbrooke, Trois-Rivières, Saguenay).1213These layers are **time-invariant** in this study (one snapshot each, census 2021) and are given to every arm;14they contain no transaction information, so they cannot leak the subject's price. Batch computation15(``geo_features_batch``) uses shapely STRtrees and a KD-tree so that all 632k rows are processed in minutes;16``geo_context`` serves a single subject from the same trees.17"""18from __future__ import annotations1920import json21import sqlite322from dataclasses import dataclass23from pathlib import Path24from typing import Any2526import numpy as np27import polars as pl2829from ..config import RAW30from ..data.analytic_schema import haversine_km, project_xy31from ..schemas import GeoContext3233EXTERNAL = RAW / "external"34CBD: dict[str, tuple[float, float]] = {35 "Montréal": (45.5017, -73.5673),36 "Québec": (46.8139, -71.2080),37 "Gatineau": (45.4765, -75.7013),38 "Sherbrooke": (45.4042, -71.8929),39 "Trois-Rivières": (46.3432, -72.5477),40 "Saguenay": (48.4284, -71.0590),41}42TRANSIT_RADIUS_M = 1000.04344GEO_COLUMNS: dict[str, pl.DataType] = {45 "tx_id": pl.String,46 "flood_mapped": pl.Boolean,47 "in_flood_zone": pl.Boolean,48 "flood_zone_kind": pl.String,49 "transit_stops_within_1km": pl.Int32,50 "dauid": pl.String,51 "da_population": pl.Int64,52 "da_population_density": pl.Float64,53 "da_median_age": pl.Float64,54 "da_median_income": pl.Float64,55 "da_pct_renters": pl.Float64,56 "cbd_name": pl.String,57 "distance_to_cbd_km": pl.Float64,58}596061def _flood_kind(description: str | None) -> str | None:62 if not description:63 return None64 d = description.lower()65 if "grand courant" in d:66 return "grand_courant_0_20"67 if "faible courant" in d:68 return "faible_courant_20_100"69 if "0-100" in d:70 return "crue_0_100"71 return "autre"727374@dataclass75class GeoLayers:76 """Loaded external layers (shapely trees). Build once with ``GeoLayers.load()``."""7778 flood_tree: Any79 flood_geoms: list[Any]80 flood_kinds: list[str | None]81 coverage_tree: Any82 coverage_geoms: list[Any]83 da_tree: Any84 da_geoms: list[Any]85 da_ids: list[str]86 da_stats: dict[str, dict[str, Any]]87 transit_tree: Any88 n_stops: int89 to_3857: Any9091 @classmethod92 def available(cls, root: Path = EXTERNAL) -> bool:93 return all((root / f).exists() for f in ("inondation.db", "transit.db", "quartier.db"))9495 @classmethod96 def load(cls, root: Path = EXTERNAL) -> "GeoLayers":97 from pyproj import Transformer98 from scipy.spatial import cKDTree99 from shapely import STRtree, wkb100 from shapely.geometry import MultiPolygon, Polygon101102 to_3857 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)103 con = sqlite3.connect(root / "inondation.db")104 geoms, kinds = [], []105 for _id, desc, blob in con.execute("SELECT id, description, wkb FROM zi"):106 g = wkb.loads(blob)107 if g.is_empty:108 continue109 geoms.append(g)110 kinds.append(_flood_kind(desc))111 cov = [wkb.loads(b) for (b,) in con.execute("SELECT wkb FROM couverture")]112 con.close()113114 con = sqlite3.connect(root / "quartier.db")115 da_geoms, da_ids = [], []116 for dauid, poly in con.execute("SELECT dauid, poly FROM da_poly"):117 rings = json.loads(poly)118 try:119 if rings and isinstance(rings[0][0][0], (int, float)): # Polygon: [ring][pt][lng,lat]120 g = Polygon(rings[0], rings[1:])121 else: # MultiPolygon: [poly][ring][pt]122 g = MultiPolygon([Polygon(p[0], p[1:]) for p in rings])123 except (IndexError, TypeError, ValueError):124 continue125 if not g.is_valid:126 g = g.buffer(0)127 da_geoms.append(g)128 da_ids.append(str(dauid))129 stats = {}130 for dauid, pop, dens, age, inc, renters in con.execute("SELECT dauid, population, densite, age_median, revenu_median, pct_locataires FROM da_stats"):131 stats[str(dauid)] = {"population": pop, "density": dens, "median_age": age, "median_income": inc, "pct_renters": renters}132 con.close()133134 con = sqlite3.connect(root / "transit.db")135 stops = np.array(con.execute("SELECT lat, lng FROM arrets WHERE lat IS NOT NULL AND lng IS NOT NULL").fetchall(), dtype=float)136 con.close()137 sx, sy = project_xy(stops[:, 0], stops[:, 1])138 ok = np.isfinite(sx) & np.isfinite(sy)139 return cls(140 flood_tree=STRtree(geoms), flood_geoms=geoms, flood_kinds=kinds,141 coverage_tree=STRtree(cov), coverage_geoms=cov,142 da_tree=STRtree(da_geoms), da_geoms=da_geoms, da_ids=da_ids, da_stats=stats,143 transit_tree=cKDTree(np.column_stack([sx[ok], sy[ok]])), n_stops=int(ok.sum()), to_3857=to_3857,144 )145146147def nearest_cbd(lat: np.ndarray, lng: np.ndarray) -> tuple[np.ndarray, np.ndarray]:148 names = list(CBD)149 d = np.column_stack([haversine_km(lat, lng, *CBD[n]) for n in names])150 j = d.argmin(axis=1)151 return np.array(names, dtype=object)[j], d[np.arange(len(lat)), j]152153154def geo_features_batch(df: pl.DataFrame, layers: GeoLayers) -> pl.DataFrame:155 """Vectorised geo features for every row of ``df`` (needs tx_id, lat, lng, x_m, y_m) → ``GEO_COLUMNS``."""156 from shapely import points157158 n = df.height159 lat = df["lat"].to_numpy().astype(float)160 lng = df["lng"].to_numpy().astype(float)161 x = df["x_m"].to_numpy().astype(float)162 y = df["y_m"].to_numpy().astype(float)163164 # flood (EPSG:3857)165 mx, my = layers.to_3857.transform(lng, lat)166 pts = points(np.column_stack([mx, my]))167 mapped = np.zeros(n, dtype=bool)168 qi, _ = layers.coverage_tree.query(pts, predicate="within")169 mapped[np.unique(qi)] = True170 in_zone = np.zeros(n, dtype=bool)171 kind = np.full(n, None, dtype=object)172 qi, gi = layers.flood_tree.query(pts, predicate="within")173 in_zone[np.unique(qi)] = True174 # most severe kind when a point lies in several polygons175 sev = {"grand_courant_0_20": 3, "crue_0_100": 2, "faible_courant_20_100": 1, "autre": 0, None: -1}176 for a, b in zip(qi, gi):177 k = layers.flood_kinds[b]178 if kind[a] is None or sev.get(k, -1) > sev.get(kind[a], -1):179 kind[a] = k180 mapped |= in_zone181182 # transit183 counts = layers.transit_tree.query_ball_point(np.column_stack([x, y]), r=TRANSIT_RADIUS_M, return_length=True)184185 # dissemination areas (WGS84)186 wpts = points(np.column_stack([lng, lat]))187 dauid = np.full(n, None, dtype=object)188 qi, gi = layers.da_tree.query(wpts, predicate="within")189 for a, b in zip(qi, gi):190 if dauid[a] is None:191 dauid[a] = layers.da_ids[b]192 stats = [layers.da_stats.get(d) if d is not None else None for d in dauid]193194 cbd_name, cbd_km = nearest_cbd(lat, lng)195 out = pl.DataFrame(196 {197 "tx_id": df["tx_id"],198 "flood_mapped": mapped,199 "in_flood_zone": np.where(mapped, in_zone, None).tolist(),200 "flood_zone_kind": kind.tolist(),201 "transit_stops_within_1km": np.asarray(counts, dtype=np.int32),202 "dauid": dauid.tolist(),203 "da_population": [None if s is None else s["population"] for s in stats],204 "da_population_density": [None if s is None else s["density"] for s in stats],205 "da_median_age": [None if s is None else s["median_age"] for s in stats],206 "da_median_income": [None if s is None else s["median_income"] for s in stats],207 "da_pct_renters": [None if s is None else s["pct_renters"] for s in stats],208 "cbd_name": cbd_name.tolist(),209 "distance_to_cbd_km": np.round(cbd_km, 3),210 },211 schema=GEO_COLUMNS,212 )213 return out214215216def geo_context(subject: dict[str, Any] | pl.DataFrame, layers: GeoLayers | None, cached: pl.DataFrame | None = None) -> GeoContext | None:217 """GeoContext for one subject. Uses the cached ``geo_features.parquet`` row when provided, else the layers.218 Returns None when neither is available."""219 from .index import subject_to_dict220221 s = subject_to_dict(subject)222 row: dict[str, Any] | None = None223 if cached is not None and "tx_id" in s:224 hit = cached.filter(pl.col("tx_id") == s["tx_id"])225 if hit.height:226 row = hit.row(0, named=True)227 if row is None:228 if layers is None:229 return None230 one = pl.DataFrame({"tx_id": [s.get("tx_id", "subject")], "lat": [float(s["lat"])], "lng": [float(s["lng"])]})231 if s.get("x_m") is not None:232 one = one.with_columns(pl.lit(float(s["x_m"])).alias("x_m"), pl.lit(float(s["y_m"])).alias("y_m"))233 else:234 xx, yy = project_xy(one["lat"].to_numpy(), one["lng"].to_numpy())235 one = one.with_columns(pl.Series("x_m", xx), pl.Series("y_m", yy))236 row = geo_features_batch(one, layers).row(0, named=True)237 notes: list[str] = []238 if not row["flood_mapped"]:239 notes.append("no provincial flood map covers this location; flood exposure unknown")240 elif row["in_flood_zone"]:241 notes.append(f"flood zone kind: {row['flood_zone_kind']}")242 if row["dauid"] is None:243 notes.append("no 2021 dissemination area matched (statistics unavailable)")244 notes.append(f"CBD reference: {row['cbd_name']}")245 return GeoContext(246 in_flood_zone=row["in_flood_zone"],247 transit_stops_within_1km=row["transit_stops_within_1km"],248 da_population_density=row["da_population_density"],249 da_median_income=row["da_median_income"],250 distance_to_cbd_km=row["distance_to_cbd_km"],251 notes=notes,252 )253