# ----------------------------------------------------------------------------- # Lou-Ka — Agrégateur de logements à louer (province de Québec) # Auteur : Simon-Pierre Boucher — contact@spboucher.ai # inondation.py : risque d'inondation à l'adresse — BDZI (gouv. du Québec) # # Source : Base de données des zones à risque d'inondation (BDZI, MELCCFP), # donnée ouverte officielle (Données Québec), GeoPackage EPSG:3857. # `build()` éclate les multipolygones de la couche ZOI_s (pleine précision) # en polygones simples indexés R*Tree dans data/inondation.db (+ la couche # Carto_ZI_S : périmètres couverts par une cartographie — sans elle, # « aucune zone » ne veut rien dire), puis `lookup(lat, lng)` fait un test # point-dans-polygone exact (shapely) sur les seuls candidats de la boîte. # # Types BDZI : « Zone de grand courant » (récurrence 0-20 ans, risque # élevé), « Zone de faible courant » (20-100 ans), « Zone de crue # 0-100 ans », variantes « - Pont » et « Autre zone inondable ». # # Usage : python run.py inondation-build (une fois, GPKG requis) # lookup(lat, lng) -> dict (fiche, /api/inondation) # ----------------------------------------------------------------------------- from __future__ import annotations import math import sqlite3 from pathlib import Path DATA = Path(__file__).resolve().parent.parent / "data" DB_PATH = DATA / "inondation.db" GPKG = DATA / "BDZI_GPK.gpkg" # rayon de tolérance : géocodage + emprise du bâtiment NEAR_M = 30.0 # au-delà, on signale quand même une zone toute proche (information utile) WARN_M = 100.0 _SEVERITE = { "Zone de grand courant": ("eleve", "récurrence 0-20 ans"), "Zone de grand courant - Pont": ("eleve", "récurrence 0-20 ans"), "Zone de faible courant": ("modere", "récurrence 20-100 ans"), "Zone de faible courant - Pont": ("modere", "récurrence 20-100 ans"), "Zone de crue 0-100 ans": ("present", "récurrence 0-100 ans"), "Zone de crue 0-100 ans - Pont": ("present", "récurrence 0-100 ans"), "Autre zone inondable": ("present", ""), } _R = 20037508.342789244 def _to_3857(lat: float, lng: float) -> tuple[float, float]: x = lng * _R / 180.0 y = math.log(math.tan((90 + lat) * math.pi / 360.0)) * _R / math.pi return x, y def _gpkg_wkb(blob: bytes) -> bytes: """Retire l'en-tête GeoPackage (magic GP + drapeaux + enveloppe).""" if blob[:2] != b"GP": return blob flags = blob[3] env = (flags >> 1) & 0x07 env_len = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}.get(env, 0) return blob[8 + env_len:] def build(gpkg: Path = GPKG) -> None: """Construit data/inondation.db à partir du GeoPackage BDZI.""" from shapely import wkb as _swkb src = sqlite3.connect(gpkg) con = sqlite3.connect(DB_PATH) con.executescript(""" DROP TABLE IF EXISTS zi; DROP TABLE IF EXISTS zi_rtree; DROP TABLE IF EXISTS couverture; DROP TABLE IF EXISTS couv_rtree; CREATE TABLE zi (id INTEGER PRIMARY KEY, description TEXT, rapport TEXT, date_rapport TEXT, wkb BLOB); CREATE VIRTUAL TABLE zi_rtree USING rtree(id, xmin, xmax, ymin, ymax); CREATE TABLE couverture (id INTEGER PRIMARY KEY, nom TEXT, wkb BLOB); CREATE VIRTUAL TABLE couv_rtree USING rtree(id, xmin, xmax, ymin, ymax); """) nid = 0 for desc, rapport, date_r, blob in src.execute( "SELECT Description, Nm_rapport, Date_rapport, Shape FROM ZOI_s"): geom = _swkb.loads(_gpkg_wkb(blob)) polys = geom.geoms if geom.geom_type == "MultiPolygon" else [geom] for poly in polys: if poly.is_empty: continue nid += 1 con.execute("INSERT INTO zi VALUES (?,?,?,?,?)", (nid, desc, rapport, date_r, poly.wkb)) x0, y0, x1, y1 = poly.bounds con.execute("INSERT INTO zi_rtree VALUES (?,?,?,?,?)", (nid, x0, x1, y0, y1)) if nid % 500 < len(polys): print(f"[inondation] {nid} polygones…", flush=True) cid = 0 for nom, blob in src.execute("SELECT Nom_Carte, Shape FROM Carto_ZI_S"): geom = _swkb.loads(_gpkg_wkb(blob)) polys = geom.geoms if geom.geom_type == "MultiPolygon" else [geom] for poly in polys: if poly.is_empty: continue cid += 1 con.execute("INSERT INTO couverture VALUES (?,?,?)", (cid, nom, poly.wkb)) x0, y0, x1, y1 = poly.bounds con.execute("INSERT INTO couv_rtree VALUES (?,?,?,?,?)", (cid, x0, x1, y0, y1)) con.commit() con.execute("VACUUM") print(f"[inondation] {nid} polygones de zones, {cid} périmètres " f"cartographiés -> {DB_PATH}") con.close() src.close() def lookup(lat: float, lng: float) -> dict | None: """Risque d'inondation BDZI au point (WGS84). None si base absente.""" if not DB_PATH.exists(): return None from shapely import wkb as _swkb from shapely.geometry import Point x, y = _to_3857(lat, lng) # les distances 3857 sont dilatées d'un facteur 1/cos(lat) scale = 1.0 / max(0.2, math.cos(math.radians(lat))) pad = WARN_M * scale pt = Point(x, y) con = sqlite3.connect(f"file:{DB_PATH}?mode=ro", uri=True) zones: list[dict] = [] for zid, desc, rapport, date_r, blob in con.execute( "SELECT z.id, z.description, z.rapport, z.date_rapport, z.wkb " "FROM zi z JOIN zi_rtree r ON z.id = r.id " "WHERE r.xmax >= ? AND r.xmin <= ? AND r.ymax >= ? AND r.ymin <= ?", (x - pad, x + pad, y - pad, y + pad)): poly = _swkb.loads(blob) d = poly.distance(pt) / scale # ~mètres réels if poly.contains(pt): d = 0.0 elif d > WARN_M: continue sev, rec = _SEVERITE.get(desc, ("present", "")) zones.append({"type": desc, "severite": sev, "recurrence": rec, "distance_m": round(d), "date_rapport": (date_r or "")[:10] or None}) couvert = False for (blob,) in con.execute( "SELECT c.wkb FROM couverture c JOIN couv_rtree r ON c.id = r.id " "WHERE r.xmax >= ? AND r.xmin <= ? AND r.ymax >= ? AND r.ymin <= ?", (x, x, y, y)): if _swkb.loads(blob).contains(pt): couvert = True break con.close() zones.sort(key=lambda z: (z["distance_m"], {"eleve": 0, "modere": 1, "present": 2} .get(z["severite"], 3))) dans = [z for z in zones if z["distance_m"] <= NEAR_M] if dans: statut = "en_zone" pire = dans[0]["severite"] elif zones: statut, pire = "a_proximite", zones[0]["severite"] elif couvert: statut, pire = "hors_zone", None else: statut, pire = "non_cartographie", None return {"statut": statut, "severite": pire, "couvert": couvert, "zones": zones[:5]}