# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca """Geospatial context (evidence id G1) from static external layers copied from the Immo-Ka node (``data/raw/external/``): * ``inondation.db`` — Quebec flood-zone polygons (``zi``: WKB in EPSG:3857, kinds "grand courant" 0–20 yr, "faible courant" 20–100 yr, …) and the mapped coverage areas (``couverture``). ``in_flood_zone`` is True/False only where a flood map exists; None (with a note) outside the mapped coverage. * ``transit.db`` — public-transit stops (bus/metro/REM/train, WGS84): count within 1 km. * ``quartier.db`` — Statistics Canada 2021 dissemination-area polygons (GeoJSON rings, WGS84) and profile statistics (population, density, median age, median income, share of renters). * distance to the nearest major CBD (Montréal, Québec, Gatineau, Sherbrooke, Trois-Rivières, Saguenay). These layers are **time-invariant** in this study (one snapshot each, census 2021) and are given to every arm; they contain no transaction information, so they cannot leak the subject's price. Batch computation (``geo_features_batch``) uses shapely STRtrees and a KD-tree so that all 632k rows are processed in minutes; ``geo_context`` serves a single subject from the same trees. """ from __future__ import annotations import json import sqlite3 from dataclasses import dataclass from pathlib import Path from typing import Any import numpy as np import polars as pl from ..config import RAW from ..data.analytic_schema import haversine_km, project_xy from ..schemas import GeoContext EXTERNAL = RAW / "external" CBD: dict[str, tuple[float, float]] = { "Montréal": (45.5017, -73.5673), "Québec": (46.8139, -71.2080), "Gatineau": (45.4765, -75.7013), "Sherbrooke": (45.4042, -71.8929), "Trois-Rivières": (46.3432, -72.5477), "Saguenay": (48.4284, -71.0590), } TRANSIT_RADIUS_M = 1000.0 GEO_COLUMNS: dict[str, pl.DataType] = { "tx_id": pl.String, "flood_mapped": pl.Boolean, "in_flood_zone": pl.Boolean, "flood_zone_kind": pl.String, "transit_stops_within_1km": pl.Int32, "dauid": pl.String, "da_population": pl.Int64, "da_population_density": pl.Float64, "da_median_age": pl.Float64, "da_median_income": pl.Float64, "da_pct_renters": pl.Float64, "cbd_name": pl.String, "distance_to_cbd_km": pl.Float64, } def _flood_kind(description: str | None) -> str | None: if not description: return None d = description.lower() if "grand courant" in d: return "grand_courant_0_20" if "faible courant" in d: return "faible_courant_20_100" if "0-100" in d: return "crue_0_100" return "autre" @dataclass class GeoLayers: """Loaded external layers (shapely trees). Build once with ``GeoLayers.load()``.""" flood_tree: Any flood_geoms: list[Any] flood_kinds: list[str | None] coverage_tree: Any coverage_geoms: list[Any] da_tree: Any da_geoms: list[Any] da_ids: list[str] da_stats: dict[str, dict[str, Any]] transit_tree: Any n_stops: int to_3857: Any @classmethod def available(cls, root: Path = EXTERNAL) -> bool: return all((root / f).exists() for f in ("inondation.db", "transit.db", "quartier.db")) @classmethod def load(cls, root: Path = EXTERNAL) -> "GeoLayers": from pyproj import Transformer from scipy.spatial import cKDTree from shapely import STRtree, wkb from shapely.geometry import MultiPolygon, Polygon to_3857 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) con = sqlite3.connect(root / "inondation.db") geoms, kinds = [], [] for _id, desc, blob in con.execute("SELECT id, description, wkb FROM zi"): g = wkb.loads(blob) if g.is_empty: continue geoms.append(g) kinds.append(_flood_kind(desc)) cov = [wkb.loads(b) for (b,) in con.execute("SELECT wkb FROM couverture")] con.close() con = sqlite3.connect(root / "quartier.db") da_geoms, da_ids = [], [] for dauid, poly in con.execute("SELECT dauid, poly FROM da_poly"): rings = json.loads(poly) try: if rings and isinstance(rings[0][0][0], (int, float)): # Polygon: [ring][pt][lng,lat] g = Polygon(rings[0], rings[1:]) else: # MultiPolygon: [poly][ring][pt] g = MultiPolygon([Polygon(p[0], p[1:]) for p in rings]) except (IndexError, TypeError, ValueError): continue if not g.is_valid: g = g.buffer(0) da_geoms.append(g) da_ids.append(str(dauid)) stats = {} for dauid, pop, dens, age, inc, renters in con.execute("SELECT dauid, population, densite, age_median, revenu_median, pct_locataires FROM da_stats"): stats[str(dauid)] = {"population": pop, "density": dens, "median_age": age, "median_income": inc, "pct_renters": renters} con.close() con = sqlite3.connect(root / "transit.db") stops = np.array(con.execute("SELECT lat, lng FROM arrets WHERE lat IS NOT NULL AND lng IS NOT NULL").fetchall(), dtype=float) con.close() sx, sy = project_xy(stops[:, 0], stops[:, 1]) ok = np.isfinite(sx) & np.isfinite(sy) return cls( flood_tree=STRtree(geoms), flood_geoms=geoms, flood_kinds=kinds, coverage_tree=STRtree(cov), coverage_geoms=cov, da_tree=STRtree(da_geoms), da_geoms=da_geoms, da_ids=da_ids, da_stats=stats, transit_tree=cKDTree(np.column_stack([sx[ok], sy[ok]])), n_stops=int(ok.sum()), to_3857=to_3857, ) def nearest_cbd(lat: np.ndarray, lng: np.ndarray) -> tuple[np.ndarray, np.ndarray]: names = list(CBD) d = np.column_stack([haversine_km(lat, lng, *CBD[n]) for n in names]) j = d.argmin(axis=1) return np.array(names, dtype=object)[j], d[np.arange(len(lat)), j] def geo_features_batch(df: pl.DataFrame, layers: GeoLayers) -> pl.DataFrame: """Vectorised geo features for every row of ``df`` (needs tx_id, lat, lng, x_m, y_m) → ``GEO_COLUMNS``.""" from shapely import points n = df.height lat = df["lat"].to_numpy().astype(float) lng = df["lng"].to_numpy().astype(float) x = df["x_m"].to_numpy().astype(float) y = df["y_m"].to_numpy().astype(float) # flood (EPSG:3857) mx, my = layers.to_3857.transform(lng, lat) pts = points(np.column_stack([mx, my])) mapped = np.zeros(n, dtype=bool) qi, _ = layers.coverage_tree.query(pts, predicate="within") mapped[np.unique(qi)] = True in_zone = np.zeros(n, dtype=bool) kind = np.full(n, None, dtype=object) qi, gi = layers.flood_tree.query(pts, predicate="within") in_zone[np.unique(qi)] = True # most severe kind when a point lies in several polygons sev = {"grand_courant_0_20": 3, "crue_0_100": 2, "faible_courant_20_100": 1, "autre": 0, None: -1} for a, b in zip(qi, gi): k = layers.flood_kinds[b] if kind[a] is None or sev.get(k, -1) > sev.get(kind[a], -1): kind[a] = k mapped |= in_zone # transit counts = layers.transit_tree.query_ball_point(np.column_stack([x, y]), r=TRANSIT_RADIUS_M, return_length=True) # dissemination areas (WGS84) wpts = points(np.column_stack([lng, lat])) dauid = np.full(n, None, dtype=object) qi, gi = layers.da_tree.query(wpts, predicate="within") for a, b in zip(qi, gi): if dauid[a] is None: dauid[a] = layers.da_ids[b] stats = [layers.da_stats.get(d) if d is not None else None for d in dauid] cbd_name, cbd_km = nearest_cbd(lat, lng) out = pl.DataFrame( { "tx_id": df["tx_id"], "flood_mapped": mapped, "in_flood_zone": np.where(mapped, in_zone, None).tolist(), "flood_zone_kind": kind.tolist(), "transit_stops_within_1km": np.asarray(counts, dtype=np.int32), "dauid": dauid.tolist(), "da_population": [None if s is None else s["population"] for s in stats], "da_population_density": [None if s is None else s["density"] for s in stats], "da_median_age": [None if s is None else s["median_age"] for s in stats], "da_median_income": [None if s is None else s["median_income"] for s in stats], "da_pct_renters": [None if s is None else s["pct_renters"] for s in stats], "cbd_name": cbd_name.tolist(), "distance_to_cbd_km": np.round(cbd_km, 3), }, schema=GEO_COLUMNS, ) return out def geo_context(subject: dict[str, Any] | pl.DataFrame, layers: GeoLayers | None, cached: pl.DataFrame | None = None) -> GeoContext | None: """GeoContext for one subject. Uses the cached ``geo_features.parquet`` row when provided, else the layers. Returns None when neither is available.""" from .index import subject_to_dict s = subject_to_dict(subject) row: dict[str, Any] | None = None if cached is not None and "tx_id" in s: hit = cached.filter(pl.col("tx_id") == s["tx_id"]) if hit.height: row = hit.row(0, named=True) if row is None: if layers is None: return None one = pl.DataFrame({"tx_id": [s.get("tx_id", "subject")], "lat": [float(s["lat"])], "lng": [float(s["lng"])]}) if s.get("x_m") is not None: one = one.with_columns(pl.lit(float(s["x_m"])).alias("x_m"), pl.lit(float(s["y_m"])).alias("y_m")) else: xx, yy = project_xy(one["lat"].to_numpy(), one["lng"].to_numpy()) one = one.with_columns(pl.Series("x_m", xx), pl.Series("y_m", yy)) row = geo_features_batch(one, layers).row(0, named=True) notes: list[str] = [] if not row["flood_mapped"]: notes.append("no provincial flood map covers this location; flood exposure unknown") elif row["in_flood_zone"]: notes.append(f"flood zone kind: {row['flood_zone_kind']}") if row["dauid"] is None: notes.append("no 2021 dissemination area matched (statistics unavailable)") notes.append(f"CBD reference: {row['cbd_name']}") return GeoContext( in_flood_zone=row["in_flood_zone"], transit_stops_within_1km=row["transit_stops_within_1km"], da_population_density=row["da_population_density"], da_median_income=row["da_median_income"], distance_to_cbd_km=row["distance_to_cbd_km"], notes=notes, )