""" build_nature_zones.py ===================== Download park, forest, and nature area polygons from the Overpass API (OpenStreetMap) and compute per-IRIS nature land coverage for the target area covering the full département Nord (59) — from Dunkerque (coast) to Maubeuge (Belgian/Walloon border), including Lille, Douai, Valenciennes, Cambrai, Avesnes-sur-Helpe and all 1 345 IRIS. Two-pass Overpass strategy: Pass 1 — ways (out geom): geometry returned directly per closed way. Pass 2 — relations (out body;>;out geom qt): member ways assembled into outer rings. OSM tags collected (polygon areas only — closed ways and multipolygon relations): leisure = park | nature_reserve | garden | recreation_ground landuse = forest | grass | meadow | orchard | allotments | village_green natural = wood | scrub | grassland | heath | wetland Coverage: département Nord (59), bbox 50.02–51.10°N / 2.30–4.30°E. Outputs: data/2026/iris_nature.json {code_iris: {nature_pct, nature_m2_approx, iris_m2_approx, categories}} nature_zones table in MySQL (raw polygon records — osm_id, category, name, bbox, coordinates JSON) Usage: python build_nature_zones.py python build_nature_zones.py --force """ from __future__ import annotations import argparse import json import logging import math import os import urllib.parse import urllib.request from datetime import datetime, timezone import numpy as np from shapely.geometry import MultiPolygon, Polygon, shape as shapely_shape from shapely.ops import unary_union from shapely.strtree import STRtree from domain.core.mysql_db import get_connection, reset_table logging.basicConfig( level=logging.INFO, format="%(asctime)s [%(levelname)s] %(name)s — %(message)s", datefmt="%H:%M:%S", ) logger = logging.getLogger("build_nature_zones") # --------------------------------------------------------------------------- # Constants # --------------------------------------------------------------------------- # Overpass bbox (min_lat,min_lon,max_lat,max_lon) — full département Nord (59) # Dunkerque (W coast) → Maubeuge (E border) / Bray-Dunes (N) → Avesnes (S) TARGET_BBOX = "50.02,2.30,51.10,4.30" OVERPASS_URL = "https://overpass-api.de/api/interpreter" _LEISURE_REGEX = "^(park|nature_reserve|garden|recreation_ground)$" _LANDUSE_REGEX = "^(forest|grass|meadow|orchard|allotments|village_green)$" _NATURAL_REGEX = "^(wood|scrub|grassland|heath|wetland)$" _Q_WAYS = """\ [out:json][timeout:500][bbox:{bbox}]; ( way["leisure"~"{leisure}"]; way["landuse"~"{landuse}"]; way["natural"~"{natural}"]; ); out geom; """ _Q_RELATIONS = """\ [out:json][timeout:500][bbox:{bbox}]; ( relation["leisure"~"{leisure}"]; relation["landuse"~"{landuse}"]; relation["natural"~"{natural}"]; ); out body; >; out geom qt; """ TAG_CATEGORY: dict[tuple[str, str], str] = { ("leisure", "park"): "parc", ("leisure", "nature_reserve"): "reserve_naturelle", ("leisure", "garden"): "jardin", ("leisure", "recreation_ground"): "terrain_recreation", ("landuse", "forest"): "foret", ("landuse", "grass"): "espace_vert", ("landuse", "meadow"): "prairie", ("landuse", "orchard"): "verger", ("landuse", "allotments"): "jardins_familiaux", ("landuse", "village_green"): "espace_vert", ("natural", "wood"): "bois", ("natural", "scrub"): "friche", ("natural", "grassland"): "prairie", ("natural", "heath"): "lande", ("natural", "wetland"): "zone_humide", } # Area conversion at 50.65°N: 1 deg² ≈ this many m² # (ratio = nature_m2/iris_m2 is accurate regardless; absolute m² are approximate) _M2_PER_DEG2 = 111_194.0 * (111_194.0 * math.cos(math.radians(50.65))) _IRIS_GEOJSON = os.path.join("data", "2026", "iris_nord.geojson") _OUTPUT_JSON = os.path.join("data", "2026", "iris_nature.json") _DDL = """ CREATE TABLE nature_zones ( id INT NOT NULL AUTO_INCREMENT, osm_id BIGINT NOT NULL, osm_type VARCHAR(10) NOT NULL, category VARCHAR(30) NOT NULL, name TEXT, geom_type VARCHAR(20) NOT NULL, coordinates MEDIUMTEXT NOT NULL, bbox_min_lng DOUBLE NOT NULL, bbox_min_lat DOUBLE NOT NULL, bbox_max_lng DOUBLE NOT NULL, bbox_max_lat DOUBLE NOT NULL, fetched_at VARCHAR(40) NOT NULL, PRIMARY KEY (id), UNIQUE KEY uq_osm_id_type (osm_id, osm_type) ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 """ _DDL_MERGED = """ CREATE TABLE nature_zones_merged ( id INT NOT NULL AUTO_INCREMENT, geom_type VARCHAR(20) NOT NULL, coordinates MEDIUMTEXT NOT NULL, bbox_min_lng DOUBLE NOT NULL, bbox_min_lat DOUBLE NOT NULL, bbox_max_lng DOUBLE NOT NULL, bbox_max_lat DOUBLE NOT NULL, built_at VARCHAR(40) NOT NULL, PRIMARY KEY (id) ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 """ _IDX_NATURE_BBOX = ( "CREATE INDEX idx_nature_bbox ON nature_zones" "(bbox_min_lng, bbox_max_lng, bbox_min_lat, bbox_max_lat)" ) _IDX_NATURE_CAT = "CREATE INDEX idx_nature_cat ON nature_zones (category)" _IDX_MERGED_BBOX = ( "CREATE INDEX idx_nature_merged_bbox ON nature_zones_merged" "(bbox_min_lng, bbox_max_lng, bbox_min_lat, bbox_max_lat)" ) # --------------------------------------------------------------------------- # Overpass fetching # --------------------------------------------------------------------------- def _overpass(query: str) -> list[dict]: body = urllib.parse.urlencode({"data": query}).encode("utf-8") req = urllib.request.Request( OVERPASS_URL, data=body, headers={ "Content-Type": "application/x-www-form-urlencoded", "User-Agent": "agent-immobilier/1.0 (nature zone builder)", }, ) with urllib.request.urlopen(req, timeout=500) as resp: return json.loads(resp.read()).get("elements", []) def fetch_ways() -> list[dict]: logger.info("Overpass — Pass 1 (ways) bbox %s …", TARGET_BBOX) q = _Q_WAYS.format( bbox=TARGET_BBOX, leisure=_LEISURE_REGEX, landuse=_LANDUSE_REGEX, natural=_NATURAL_REGEX, ) els = _overpass(q) ways = [e for e in els if e.get("type") == "way"] logger.info(" %d ways reçus", len(ways)) return ways def fetch_relations() -> tuple[list[dict], dict[int, list[dict]]]: """ Returns (relations, way_geom_by_id): relations — list of relation elements (with members) way_geom_by_id — {way_id: [{lat, lon}, ...]} for all member ways returned """ logger.info("Overpass — Pass 2 (relations) bbox %s …", TARGET_BBOX) q = _Q_RELATIONS.format( bbox=TARGET_BBOX, leisure=_LEISURE_REGEX, landuse=_LANDUSE_REGEX, natural=_NATURAL_REGEX, ) els = _overpass(q) relations: list[dict] = [] way_geom: dict[int, list[dict]] = {} for el in els: t = el.get("type") if t == "relation": relations.append(el) elif t == "way" and el.get("geometry"): way_geom[el["id"]] = el["geometry"] logger.info(" %d relations, %d member ways reçus", len(relations), len(way_geom)) return relations, way_geom # --------------------------------------------------------------------------- # Polygon building helpers # --------------------------------------------------------------------------- def _classify(tags: dict) -> str | None: for key in ("leisure", "landuse", "natural"): val = tags.get(key, "") cat = TAG_CATEGORY.get((key, val)) if cat: return cat return None def _coords_from_geom(geom_nodes: list[dict]) -> list[tuple[float, float]]: """Convert Overpass {lat,lon} list → [(lon, lat), ...] for Shapely (x=lon, y=lat).""" return [(n["lon"], n["lat"]) for n in geom_nodes if "lon" in n and "lat" in n] def _make_polygon(coords: list[tuple[float, float]]) -> Polygon | None: """Build and validate a Shapely Polygon; try buffer(0) for invalid rings.""" if len(coords) < 3: return None # Ensure ring is closed if coords[0] != coords[-1]: coords = coords + [coords[0]] try: poly = Polygon(coords) if not poly.is_valid: poly = poly.buffer(0) if poly.is_empty or not poly.is_valid: return None return poly except Exception: return None def _assemble_outer_ring( member_way_ids: list[int], way_geom: dict[int, list[dict]], ) -> list[tuple[float, float]] | None: """ Chain member way coordinate lists into a single outer ring. OSM relation outer rings are built from one or more ways whose endpoints must connect. This assembles them greedily by matching endpoints. Returns the ordered (lon, lat) ring, or None if assembly fails. """ segs: list[list[tuple[float, float]]] = [] for wid in member_way_ids: raw = way_geom.get(wid) if raw: segs.append(_coords_from_geom(raw)) if not segs: return None if len(segs) == 1: return segs[0] assembled = list(segs.pop(0)) max_tries = len(segs) * len(segs) + 1 tries = 0 while segs and tries < max_tries: tries += 1 end = assembled[-1] matched = False for i, seg in enumerate(segs): if not seg: continue if seg[0] == end: assembled.extend(seg[1:]) segs.pop(i) matched = True break if seg[-1] == end: assembled.extend(reversed(seg[:-1])) segs.pop(i) matched = True break if not matched: # Try reversing assembled to see if we can connect to the start assembled = list(reversed(assembled)) end = assembled[-1] for i, seg in enumerate(segs): if not seg: continue if seg[0] == end: assembled.extend(seg[1:]) segs.pop(i) matched = True break if seg[-1] == end: assembled.extend(reversed(seg[:-1])) segs.pop(i) matched = True break if not matched: break # give up — disconnected ring return assembled if not segs else None # --------------------------------------------------------------------------- # Feature parsing # --------------------------------------------------------------------------- def parse_ways(ways: list[dict]) -> list[dict]: """ Convert way elements into feature dicts with a Shapely polygon. Returns list of {osm_id, osm_type, category, name, polygon}. """ features = [] skipped = 0 for el in ways: tags = el.get("tags") or {} cat = _classify(tags) if cat is None: skipped += 1 continue raw_geom = el.get("geometry") or [] coords = _coords_from_geom(raw_geom) poly = _make_polygon(coords) if poly is None or poly.area == 0: skipped += 1 continue features.append({ "osm_id": el["id"], "osm_type": "way", "category": cat, "name": (tags.get("name") or "").strip() or None, "polygon": poly, }) by_cat: dict[str, int] = {} for f in features: by_cat[f["category"]] = by_cat.get(f["category"], 0) + 1 logger.info( "Ways → %d polygones valides | %d ignorés | %s", len(features), skipped, " ".join(f"{k}:{v}" for k, v in sorted(by_cat.items())), ) return features def parse_relations( relations: list[dict], way_geom: dict[int, list[dict]], ) -> list[dict]: """ Convert relation elements into feature dicts with a Shapely polygon. Only outer rings are assembled; inner rings (holes) are ignored (conservative: slightly over-estimates coverage for parks with internal roads/buildings). """ features = [] skipped = 0 for rel in relations: tags = rel.get("tags") or {} cat = _classify(tags) if cat is None: skipped += 1 continue outer_ids = [ m["ref"] for m in (rel.get("members") or []) if m.get("type") == "way" and m.get("role") in ("outer", "") ] if not outer_ids: skipped += 1 continue coords = _assemble_outer_ring(outer_ids, way_geom) if coords is None: # Fallback: concatenate all outer way coords (may self-intersect; buffer(0) fixes) coords = [] for wid in outer_ids: raw = way_geom.get(wid) if raw: coords.extend(_coords_from_geom(raw)) poly = _make_polygon(coords) if poly is None or poly.area == 0: skipped += 1 continue features.append({ "osm_id": rel["id"], "osm_type": "relation", "category": cat, "name": (tags.get("name") or "").strip() or None, "polygon": poly, }) by_cat: dict[str, int] = {} for f in features: by_cat[f["category"]] = by_cat.get(f["category"], 0) + 1 logger.info( "Relations → %d polygones valides | %d ignorés | %s", len(features), skipped, " ".join(f"{k}:{v}" for k, v in sorted(by_cat.items())), ) return features # --------------------------------------------------------------------------- # IRIS loading # --------------------------------------------------------------------------- def load_iris(path: str) -> tuple[list[str], list[str], list[object]]: """Load IRIS polygons from iris_nord.geojson. Returns (codes, noms, polygons).""" logger.info("Chargement IRIS : %s", path) with open(path, encoding="utf-8") as f: gj = json.load(f) codes, noms, polys = [], [], [] for feat in gj["features"]: try: p = feat.get("properties") or {} code = str(p.get("code_iris", "")).zfill(9) nom = str(p.get("nom_iris") or "") poly = shapely_shape(feat["geometry"]) if not poly.is_valid: poly = poly.buffer(0) codes.append(code) noms.append(nom) polys.append(poly) except Exception: pass logger.info(" %d polygones IRIS chargés", len(polys)) return codes, noms, polys # --------------------------------------------------------------------------- # Per-IRIS coverage computation # --------------------------------------------------------------------------- def compute_iris_coverage( features: list[dict], iris_codes: list[str], iris_noms: list[str], iris_polys: list[object], ) -> dict[str, dict]: """ Compute per-IRIS nature coverage. For each IRIS, collect all intersecting nature polygons, take their union (avoids double-counting overlapping OSM features), then compute intersection with the IRIS polygon. Returns {code_iris: {nature_pct, nature_m2_approx, iris_m2_approx, categories}}. """ iris_arr = np.array(iris_polys, dtype=object) iris_tree = STRtree(iris_arr) # Collect per-IRIS: {idx: [(polygon, category), ...]} iris_candidates: dict[int, list[tuple]] = {} total = len(features) logger.info("Indexation des candidats IRIS (%d polygones nature) …", total) for i, feat in enumerate(features): poly = feat["polygon"] cat = feat["category"] for idx in iris_tree.query(poly): if iris_arr[idx].intersects(poly): if idx not in iris_candidates: iris_candidates[idx] = [] iris_candidates[idx].append((poly, cat)) if (i + 1) % 500 == 0: logger.info(" %d / %d polygones indexés", i + 1, total) logger.info("Union + intersection pour %d IRIS …", len(iris_candidates)) result: dict[str, dict] = {} for idx, pairs in iris_candidates.items(): iris_poly = iris_arr[idx] i_area = iris_poly.area if i_area <= 0: continue nat_polys = [p for p, _ in pairs] cats = {c for _, c in pairs} # Union all nature polygons touching this IRIS (removes double-counting) try: nat_union = unary_union(nat_polys) inter = iris_poly.intersection(nat_union) except Exception: continue if inter.is_empty: continue n_area = inter.area pct = round(min(100.0, 100.0 * n_area / i_area), 2) code = iris_codes[idx] result[code] = { "nature_pct": pct, "nature_m2_approx": round(n_area * _M2_PER_DEG2), "iris_m2_approx": round(i_area * _M2_PER_DEG2), "categories": sorted(cats), } # IRIS with zero nature coverage → include with pct=0 for completeness for idx, code in enumerate(iris_codes): if code not in result: i_area = iris_arr[idx].area result[code] = { "nature_pct": 0.0, "nature_m2_approx": 0, "iris_m2_approx": round(i_area * _M2_PER_DEG2), "categories": [], } return result # --------------------------------------------------------------------------- # Database save # --------------------------------------------------------------------------- def save_db(features: list[dict], conn) -> None: fetched_at = datetime.now(timezone.utc).isoformat() reset_table(conn, "nature_zones", _DDL, indexes=[_IDX_NATURE_BBOX, _IDX_NATURE_CAT]) rows = [] for feat in features: poly = feat["polygon"] bounds = poly.bounds # (min_lng, min_lat, max_lng, max_lat) # Exterior ring coordinates as JSON [[lon,lat], ...] if isinstance(poly, Polygon): ext_coords = list(poly.exterior.coords) elif isinstance(poly, MultiPolygon): ext_coords = list(max(poly.geoms, key=lambda p: p.area).exterior.coords) else: ext_coords = [] coords_json = json.dumps([[[round(c[0], 7), round(c[1], 7)] for c in ext_coords]], separators=(",", ":")) rows.append(( feat["osm_id"], feat["osm_type"], feat["category"], feat.get("name"), "Polygon", coords_json, bounds[0], bounds[1], bounds[2], bounds[3], fetched_at, )) with conn.cursor() as cur: cur.executemany( """INSERT IGNORE INTO nature_zones (osm_id, osm_type, category, name, geom_type, coordinates, bbox_min_lng, bbox_min_lat, bbox_max_lng, bbox_max_lat, fetched_at) VALUES (%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s)""", rows, ) conn.commit() logger.info("Table nature_zones : %d polygones insérés → MySQL", len(rows)) # --------------------------------------------------------------------------- # Pre-merge build (fast-path for /nature-polygons API) # --------------------------------------------------------------------------- def build_nature_merged(conn) -> None: """Pre-compute buffer+union+simplify and store in nature_zones_merged table. This runs the full Shapely pipeline once at build time so that the /nature-polygons API endpoint can serve results with a pure SQL bbox query instead of recomputing Shapely geometry on every request. """ logger.info("Pré-calcul nature_zones_merged …") with conn.cursor() as cur: cur.execute( """SELECT geom_type, coordinates, (bbox_max_lng - bbox_min_lng) * (bbox_max_lat - bbox_min_lat) AS bbox_area FROM nature_zones ORDER BY bbox_area DESC""" ) rows = cur.fetchall() if not rows: logger.warning("nature_zones vide — pré-calcul ignoré") return # Drop bottom 30% by bbox proxy area (same filter as the API) areas = sorted(r["bbox_area"] for r in rows) threshold = areas[int(len(areas) * 0.30)] rows = [r for r in rows if r["bbox_area"] >= threshold] logger.info(" %d polygones retenus (top 70%% par surface)", len(rows)) polys = [] for r in rows: coords = json.loads(r["coordinates"]) if coords and not isinstance(coords[0][0], list): coords = [coords] try: poly = shapely_shape({"type": r["geom_type"], "coordinates": coords}) if not poly.is_valid: poly = poly.buffer(0) if poly.is_valid and not poly.is_empty: polys.append(poly) except Exception: continue if not polys: logger.warning("Aucun polygone valide — pré-calcul ignoré") return logger.info(" buffer+union+simplify sur %d polygones …", len(polys)) MERGE_DIST = 0.0006 SHRINK_DIST = 0.0003 SIMPLIFY_TOL = 0.0001 MIN_AREA = 5e-7 try: merged = unary_union([p.buffer(MERGE_DIST) for p in polys]).buffer(-SHRINK_DIST) final = merged.simplify(SIMPLIFY_TOL, preserve_topology=True) except Exception: try: final = unary_union(polys) except Exception: logger.error("Fusion échouée — pré-calcul ignoré") return def _collect(geom): t = geom.geom_type if t in ("Polygon", "MultiPolygon"): if geom.area >= MIN_AREA: yield geom elif t == "GeometryCollection": for g in geom.geoms: yield from _collect(g) def _to_coords(geom): if geom.geom_type == "Polygon": rings = [list(geom.exterior.coords)] + [list(r.coords) for r in geom.interiors] return "Polygon", [[[round(x, 6), round(y, 6)] for x, y in ring] for ring in rings] polys_out = [] for p in geom.geoms: rings = [list(p.exterior.coords)] + [list(r.coords) for r in p.interiors] polys_out.append([[[round(x, 6), round(y, 6)] for x, y in ring] for ring in rings]) return "MultiPolygon", polys_out built_at = datetime.now(timezone.utc).isoformat() reset_table(conn, "nature_zones_merged", _DDL_MERGED, indexes=[_IDX_MERGED_BBOX]) count = 0 with conn.cursor() as cur: for geom in _collect(final): gtype, coords = _to_coords(geom) bounds = geom.bounds # (min_lng, min_lat, max_lng, max_lat) cur.execute( """INSERT INTO nature_zones_merged (geom_type, coordinates, bbox_min_lng, bbox_min_lat, bbox_max_lng, bbox_max_lat, built_at) VALUES (%s,%s,%s,%s,%s,%s,%s)""", (gtype, json.dumps(coords, separators=(",", ":")), bounds[0], bounds[1], bounds[2], bounds[3], built_at), ) count += 1 conn.commit() logger.info("Table nature_zones_merged : %d polygones → MySQL", count) # --------------------------------------------------------------------------- # Main # --------------------------------------------------------------------------- def main() -> None: parser = argparse.ArgumentParser( description="Calcule la couverture naturelle/verte par IRIS (Overpass API)" ) parser.add_argument( "--force", action="store_true", help="Reconstruire même si iris_nature.json est déjà présent", ) args = parser.parse_args() if os.path.exists(_OUTPUT_JSON) and not args.force: import json as _j with open(_OUTPUT_JSON, encoding="utf-8") as _f: existing = _j.load(_f) print(f"iris_nature.json déjà présent ({len(existing)} IRIS) — utilisez --force pour reconstruire.") return print() print("=" * 60) print(" Nature Zones — Département Nord (59) — couverture complète") print("=" * 60) print(f" Bbox : {TARGET_BBOX} (50.02–51.10°N / 2.30–4.30°E)") print(f" Source : Overpass API (OpenStreetMap)") print(f" IRIS : {_IRIS_GEOJSON}") print(f" Output : {_OUTPUT_JSON}") print() # -- Fetch ways = fetch_ways() relations, rel_way_geom = fetch_relations() # -- Parse way_features = parse_ways(ways) rel_features = parse_relations(relations, rel_way_geom) all_features = way_features + rel_features logger.info("Total polygones nature : %d", len(all_features)) # -- Load IRIS iris_codes, iris_noms, iris_polys = load_iris(_IRIS_GEOJSON) # -- Compute coverage coverage = compute_iris_coverage(all_features, iris_codes, iris_noms, iris_polys) # -- Summary stats covered = sum(1 for v in coverage.values() if v["nature_pct"] > 0) pcts = [v["nature_pct"] for v in coverage.values() if v["nature_pct"] > 0] mean_pct = sum(pcts) / len(pcts) if pcts else 0.0 high_coverage = sorted( [(code, v) for code, v in coverage.items() if v["nature_pct"] >= 20], key=lambda x: -x[1]["nature_pct"], )[:10] print() print(f" IRIS couverts (>0%) : {covered} / {len(coverage)}") print(f" Couverture moyenne : {mean_pct:.1f}% (IRIS avec nature seulement)") if high_coverage: print() print(" Top 10 IRIS les plus verts :") for code, v in high_coverage: cats = "+".join(v["categories"][:3]) print(f" {code} {v['nature_pct']:>5.1f}% {round(v['nature_m2_approx']/10000, 1):>7.1f} ha [{cats}]") # -- Save JSON os.makedirs(os.path.dirname(_OUTPUT_JSON), exist_ok=True) with open(_OUTPUT_JSON, "w", encoding="utf-8") as f: json.dump(coverage, f, separators=(",", ":"), ensure_ascii=False) logger.info("iris_nature.json → %s (%d IRIS)", _OUTPUT_JSON, len(coverage)) # -- Save DB + pre-merged table conn = get_connection() try: save_db(all_features, conn) build_nature_merged(conn) finally: conn.close() print() print("=" * 60) print(f" Polygones insérés : {len(all_features):,}") print(f" IRIS enrichis : {covered:,} / {len(coverage):,}") print(f" iris_nature.json : {_OUTPUT_JSON}") print(f" Base : MySQL (MYSQL_DATABASE env var)") print(f" nature_zones_merged : pré-calculé (/nature-polygons rapide)") print() print(" Prochaines étapes :") print(" python build_pression_iris.py # intégrer greenery_pct dans BPS") print("=" * 60) if __name__ == "__main__": main()