spb/lou-ka

1# -----------------------------------------------------------------------------2# Lou-Ka — Agrégateur de logements à louer (province de Québec)3# Auteur : Simon-Pierre Boucher — contact@spboucher.ai4# inondation.py : risque d'inondation à l'adresse — BDZI (gouv. du Québec)5#6#   Source : Base de données des zones à risque d'inondation (BDZI, MELCCFP),7#   donnée ouverte officielle (Données Québec), GeoPackage EPSG:3857.8#   `build()` éclate les multipolygones de la couche ZOI_s (pleine précision)9#   en polygones simples indexés R*Tree dans data/inondation.db (+ la couche10#   Carto_ZI_S : périmètres couverts par une cartographie — sans elle,11#   « aucune zone » ne veut rien dire), puis `lookup(lat, lng)` fait un test12#   point-dans-polygone exact (shapely) sur les seuls candidats de la boîte.13#14#   Types BDZI : « Zone de grand courant » (récurrence 0-20 ans, risque15#   élevé), « Zone de faible courant » (20-100 ans), « Zone de crue16#   0-100 ans », variantes « - Pont » et « Autre zone inondable ».17#18#   Usage : python run.py inondation-build   (une fois, GPKG requis)19#           lookup(lat, lng) -> dict (fiche, /api/inondation)20# -----------------------------------------------------------------------------21from __future__ import annotations2223import math24import sqlite325from pathlib import Path2627DATA = Path(__file__).resolve().parent.parent / "data"28DB_PATH = DATA / "inondation.db"29GPKG = DATA / "BDZI_GPK.gpkg"3031# rayon de tolérance : géocodage + emprise du bâtiment32NEAR_M = 30.033# au-delà, on signale quand même une zone toute proche (information utile)34WARN_M = 100.03536_SEVERITE = {37    "Zone de grand courant": ("eleve", "récurrence 0-20 ans"),38    "Zone de grand courant - Pont": ("eleve", "récurrence 0-20 ans"),39    "Zone de faible courant": ("modere", "récurrence 20-100 ans"),40    "Zone de faible courant - Pont": ("modere", "récurrence 20-100 ans"),41    "Zone de crue 0-100 ans": ("present", "récurrence 0-100 ans"),42    "Zone de crue 0-100 ans - Pont": ("present", "récurrence 0-100 ans"),43    "Autre zone inondable": ("present", ""),44}4546_R = 20037508.342789244474849def _to_3857(lat: float, lng: float) -> tuple[float, float]:50    x = lng * _R / 180.051    y = math.log(math.tan((90 + lat) * math.pi / 360.0)) * _R / math.pi52    return x, y535455def _gpkg_wkb(blob: bytes) -> bytes:56    """Retire l'en-tête GeoPackage (magic GP + drapeaux + enveloppe)."""57    if blob[:2] != b"GP":58        return blob59    flags = blob[3]60    env = (flags >> 1) & 0x0761    env_len = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}.get(env, 0)62    return blob[8 + env_len:]636465def build(gpkg: Path = GPKG) -> None:66    """Construit data/inondation.db à partir du GeoPackage BDZI."""67    from shapely import wkb as _swkb6869    src = sqlite3.connect(gpkg)70    con = sqlite3.connect(DB_PATH)71    con.executescript("""72        DROP TABLE IF EXISTS zi; DROP TABLE IF EXISTS zi_rtree;73        DROP TABLE IF EXISTS couverture; DROP TABLE IF EXISTS couv_rtree;74        CREATE TABLE zi (id INTEGER PRIMARY KEY, description TEXT,75                         rapport TEXT, date_rapport TEXT, wkb BLOB);76        CREATE VIRTUAL TABLE zi_rtree USING rtree(id, xmin, xmax, ymin, ymax);77        CREATE TABLE couverture (id INTEGER PRIMARY KEY, nom TEXT, wkb BLOB);78        CREATE VIRTUAL TABLE couv_rtree USING rtree(id, xmin, xmax, ymin, ymax);79    """)80    nid = 081    for desc, rapport, date_r, blob in src.execute(82            "SELECT Description, Nm_rapport, Date_rapport, Shape FROM ZOI_s"):83        geom = _swkb.loads(_gpkg_wkb(blob))84        polys = geom.geoms if geom.geom_type == "MultiPolygon" else [geom]85        for poly in polys:86            if poly.is_empty:87                continue88            nid += 189            con.execute("INSERT INTO zi VALUES (?,?,?,?,?)",90                        (nid, desc, rapport, date_r, poly.wkb))91            x0, y0, x1, y1 = poly.bounds92            con.execute("INSERT INTO zi_rtree VALUES (?,?,?,?,?)",93                        (nid, x0, x1, y0, y1))94        if nid % 500 < len(polys):95            print(f"[inondation] {nid} polygones…", flush=True)96    cid = 097    for nom, blob in src.execute("SELECT Nom_Carte, Shape FROM Carto_ZI_S"):98        geom = _swkb.loads(_gpkg_wkb(blob))99        polys = geom.geoms if geom.geom_type == "MultiPolygon" else [geom]100        for poly in polys:101            if poly.is_empty:102                continue103            cid += 1104            con.execute("INSERT INTO couverture VALUES (?,?,?)",105                        (cid, nom, poly.wkb))106            x0, y0, x1, y1 = poly.bounds107            con.execute("INSERT INTO couv_rtree VALUES (?,?,?,?,?)",108                        (cid, x0, x1, y0, y1))109    con.commit()110    con.execute("VACUUM")111    print(f"[inondation] {nid} polygones de zones, {cid} périmètres "112          f"cartographiés -> {DB_PATH}")113    con.close()114    src.close()115116117def lookup(lat: float, lng: float) -> dict | None:118    """Risque d'inondation BDZI au point (WGS84). None si base absente."""119    if not DB_PATH.exists():120        return None121    from shapely import wkb as _swkb122    from shapely.geometry import Point123124    x, y = _to_3857(lat, lng)125    # les distances 3857 sont dilatées d'un facteur 1/cos(lat)126    scale = 1.0 / max(0.2, math.cos(math.radians(lat)))127    pad = WARN_M * scale128    pt = Point(x, y)129130    con = sqlite3.connect(f"file:{DB_PATH}?mode=ro", uri=True)131    zones: list[dict] = []132    for zid, desc, rapport, date_r, blob in con.execute(133            "SELECT z.id, z.description, z.rapport, z.date_rapport, z.wkb "134            "FROM zi z JOIN zi_rtree r ON z.id = r.id "135            "WHERE r.xmax >= ? AND r.xmin <= ? AND r.ymax >= ? AND r.ymin <= ?",136            (x - pad, x + pad, y - pad, y + pad)):137        poly = _swkb.loads(blob)138        d = poly.distance(pt) / scale        # ~mètres réels139        if poly.contains(pt):140            d = 0.0141        elif d > WARN_M:142            continue143        sev, rec = _SEVERITE.get(desc, ("present", ""))144        zones.append({"type": desc, "severite": sev, "recurrence": rec,145                      "distance_m": round(d),146                      "date_rapport": (date_r or "")[:10] or None})147    couvert = False148    for (blob,) in con.execute(149            "SELECT c.wkb FROM couverture c JOIN couv_rtree r ON c.id = r.id "150            "WHERE r.xmax >= ? AND r.xmin <= ? AND r.ymax >= ? AND r.ymin <= ?",151            (x, x, y, y)):152        if _swkb.loads(blob).contains(pt):153            couvert = True154            break155    con.close()156157    zones.sort(key=lambda z: (z["distance_m"],158                              {"eleve": 0, "modere": 1, "present": 2}159                              .get(z["severite"], 3)))160    dans = [z for z in zones if z["distance_m"] <= NEAR_M]161    if dans:162        statut = "en_zone"163        pire = dans[0]["severite"]164    elif zones:165        statut, pire = "a_proximite", zones[0]["severite"]166    elif couvert:167        statut, pire = "hors_zone", None168    else:169        statut, pire = "non_cartographie", None170    return {"statut": statut, "severite": pire, "couvert": couvert,171            "zones": zones[:5]}172