""" geo_utils.py ============ Utilitaires géographiques bas niveau : - Conversion GML 3.1.1 → GeoJSON (WFS 1.1.0 / EPSG:4326) - Helpers SQLite pour la table flood_zones - Point-in-polygon (ray casting) """ from __future__ import annotations import json import logging import xml.etree.ElementTree as ET from domain.db import get_db logger = logging.getLogger(__name__) # --------------------------------------------------------------------------- # GML 3.1.1 → GeoJSON # WFS 1.1.0 + EPSG:4326 retourne (lat, lon) ; _parse_pos_list détecte et # permute en (lon, lat) tel que requis par GeoJSON. # --------------------------------------------------------------------------- def _parse_pos_list(text: str) -> list: nums = list(map(float, text.split())) pairs = [[nums[i], nums[i + 1]] for i in range(0, len(nums) - 1, 2)] # France : lon ∈ [-5, 8], lat ∈ [42, 51]. Si la 1ʳᵉ valeur > 10 → c'est une latitude. if pairs and pairs[0][0] > 10: return [[p[1], p[0]] for p in pairs] return pairs def _parse_ring(ring_el, GML: str) -> list: pl = ring_el.find(f"{{{GML}}}posList") if pl is not None and pl.text: return _parse_pos_list(pl.text.strip()) co = ring_el.find(f"{{{GML}}}coordinates") if co is not None and co.text: return [list(map(float, p.split(","))) for p in co.text.strip().split()] return [] def _parse_polygon_coords(p_el, GML: str) -> list: rings = [] ext = p_el.find(f"{{{GML}}}exterior/{{{GML}}}LinearRing") if ext is None: ext = p_el.find(f"{{{GML}}}outerBoundaryIs/{{{GML}}}LinearRing") if ext is not None: rings.append(_parse_ring(ext, GML)) for tag in (f"{{{GML}}}interior", f"{{{GML}}}innerBoundaryIs"): for intr in p_el.findall(f"{tag}/{{{GML}}}LinearRing"): rings.append(_parse_ring(intr, GML)) return rings def _parse_gml_feature(feat_el, GML: str) -> dict | None: ms = feat_el.find(f".//{{{GML}}}MultiSurface") if ms is not None: polys = [] for sm in ms.findall(f"{{{GML}}}surfaceMember"): p = sm.find(f"{{{GML}}}Polygon") if p is not None: c = _parse_polygon_coords(p, GML) if c: polys.append(c) if polys: return {"type": "MultiPolygon", "coordinates": polys} mp = feat_el.find(f".//{{{GML}}}MultiPolygon") if mp is not None: polys = [] for pm in mp.findall(f"{{{GML}}}polygonMember"): p = pm.find(f"{{{GML}}}Polygon") if p is not None: c = _parse_polygon_coords(p, GML) if c: polys.append(c) if polys: return {"type": "MultiPolygon", "coordinates": polys} p = feat_el.find(f".//{{{GML}}}Polygon") if p is not None: c = _parse_polygon_coords(p, GML) if c: return {"type": "Polygon", "coordinates": c} return None def _gml_to_geojson(gml_bytes: bytes) -> dict: GML = "http://www.opengis.net/gml" root = ET.fromstring(gml_bytes) features = [] for tag in (f"{{{GML}}}featureMember", f"{{{GML}}}featureMembers"): for member in root.iter(tag): for feat in member: geom = _parse_gml_feature(feat, GML) if geom: features.append({"type": "Feature", "geometry": geom, "properties": {}}) logger.info("GML parse — %d features extraits", len(features)) if features: geom0 = features[0]["geometry"] sample = ( geom0["coordinates"][0][0] if geom0["type"] == "Polygon" else geom0["coordinates"][0][0][0] ) logger.debug("GML — premier coord (doit être [lon,lat] ~[1-5, 50-51]) : %s", sample) return {"type": "FeatureCollection", "features": features} def _wfs_params(typename: str, bbox: str, max_features: int = 1) -> dict: """Paramètres GetFeature WFS 1.1.0 (sans outputFormat → GML par défaut).""" return { "SERVICE": "WFS", "VERSION": "1.1.0", "REQUEST": "GetFeature", "typeName": typename, "BBOX": bbox, "maxFeatures": max_features, } # --------------------------------------------------------------------------- # Flood zones — helpers MySQL (table flood_zones, construite par build_flood_zones.py) # --------------------------------------------------------------------------- def _flood_table_exists(conn) -> bool: return conn.execute( "SELECT 1 FROM information_schema.tables" " WHERE table_schema=DATABASE() AND table_name='flood_zones'" ).fetchone() is not None def _ray_cast(lng: float, lat: float, ring: list) -> bool: """Test point-dans-polygone (ray casting) pour un anneau GeoJSON [lng, lat].""" inside = False j = len(ring) - 1 for i in range(len(ring)): xi, yi = ring[i][0], ring[i][1] xj, yj = ring[j][0], ring[j][1] if (yi > lat) != (yj > lat) and lng < (xj - xi) * (lat - yi) / (yj - yi) + xi: inside = not inside j = i return inside def _point_in_feature(lng: float, lat: float, geom_type: str, coords: list) -> bool: """True si (lng, lat) est à l'intérieur d'un Polygon ou MultiPolygon (trous exclus).""" if geom_type == "Polygon": return _ray_cast(lng, lat, coords[0]) and not any( _ray_cast(lng, lat, h) for h in coords[1:] ) if geom_type == "MultiPolygon": for poly in coords: if _ray_cast(lng, lat, poly[0]) and not any( _ray_cast(lng, lat, h) for h in poly[1:] ): return True return False def flood_polygons_from_db( scenario: str, min_lng: float, min_lat: float, max_lng: float, max_lat: float, ) -> dict: """GeoJSON FeatureCollection depuis flood_zones, filtré par scenario + bbox.""" with get_db() as conn: if not _flood_table_exists(conn): logger.warning( "Table flood_zones absente | scénario=%s — exécutez build_flood_zones.py", scenario, ) return {"type": "FeatureCollection", "features": []} rows = conn.execute( """SELECT geom_type, coordinates FROM flood_zones WHERE scenario = ? AND bbox_max_lng > ? AND bbox_min_lng < ? AND bbox_max_lat > ? AND bbox_min_lat < ?""", (scenario, min_lng, max_lng, min_lat, max_lat), ).fetchall() features = [ { "type": "Feature", "geometry": { "type": r["geom_type"], "coordinates": json.loads(r["coordinates"]), }, "properties": {}, } for r in rows ] logger.info("flood_zones DB | scénario=%s → %d features", scenario, len(features)) return {"type": "FeatureCollection", "features": features} def flood_zones_from_db(lng: float, lat: float) -> dict[str, bool | None]: """ Test d'appartenance à une zone TRI débordement pour les 3 scénarios. Retourne True / False / None (pas de données pour ce scénario). """ with get_db() as conn: if not _flood_table_exists(conn): logger.debug("flood_zones absente — retourne None pour tous les scénarios") return {"frequent": None, "moyen": None, "rare": None} result: dict[str, bool | None] = {} for scenario in ("frequent", "moyen", "rare"): rows = conn.execute( "SELECT geom_type, coordinates FROM flood_zones WHERE scenario = ?", (scenario,), ).fetchall() if not rows: result[scenario] = None continue result[scenario] = any( _point_in_feature(lng, lat, r["geom_type"], json.loads(r["coordinates"])) for r in rows ) logger.debug("flood_zones_from_db | lng=%s lat=%s | résultat=%s", lng, lat, result) return result def simplify_coords(coords: list, precision: int = 4) -> list: """Round coordinate values to reduce GeoJSON payload size (~30%).""" if not coords: return coords if isinstance(coords[0], list): return [simplify_coords(c, precision) for c in coords] return [round(coords[0], precision), round(coords[1], precision)]