#!/usr/bin/env python3 """Chapitre 38 — SIG et risques géospatiaux (cas 1 Gatineau et cas 4 portefeuille, DONNÉES SYNTHÉTIQUES, graine 20261031). Reproduit tous les chiffres cités dans manuscript/chapters/ch38.tex : 1. systèmes de référence : projection des 600 ventes du cas 1 en NAD83 / MTM zone 9 (EPSG 32189) ; comparaison avec la projection équirectangulaire locale du chapitre 19, le haversine et la géodésique WGS84 ; facteur d'échelle ; 2. géométries FICTIVES déclarées : « cours d'eau fictif » et ses deux classes d'exposition (tampons), « corridor de transport fictif » et ses tampons, trois « pôles » fictifs pour l'accessibilité ; 3. opérations SIG : jointure spatiale (point dans polygone), tampons, distances (euclidienne, Manhattan, au corridor), agrégation par cellule hexagonale, pondération surfacique, interpolation IDW ; 4. accessibilité gravitaire (Hansen) ; 5. intégration au modèle hédonique de référence (traite.cas1.hedonique) : variables d'exposition, de corridor et d'accessibilité — l'effet attendu est NUL (le générateur ne contient pas ces variables) ; 6. faux positifs spatiaux : polygones aléatoires et part de coefficients « significatifs » selon la spécification ; 7. sensibilité au rayon du tampon et à la méthode de classification des cartes choroplèthes ; 8. cas 4 : 48 immeubles par RMR, part de la valeur en zone inondable, coefficient « inondable » des effets régionaux. Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch38_sig.py""" from __future__ import annotations import sys, pathlib, warnings import numpy as np, pandas as pd sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1])) from traite import cas1, cas4, geo from traite import hedonique as h ROOT = pathlib.Path(__file__).resolve().parents[3] GRAINE = 20261031 EPSG_WGS84, EPSG_MTM9, EPSG_LAMBERT = 4326, 32189, 32198 LON0_MTM9 = -76.5 # méridien central de la zone MTM 9 RAYON_HEX_M = 1200.0 # rayon circonscrit des cellules hexagonales (m) ; aire ≈ 3,74 km² MIN_VENTES_CELLULE = 5 RAYONS_CORRIDOR_M = [100, 150, 200, 300, 400, 500, 750, 1000] RAYON_CORRIDOR_REF_M = 300 N_POLYGONES_ALEATOIRES = 500 RAYONS_POLYGONES_KM = [0.5, 1.0, 2.0] K_CLASSES = 5 # géométries FICTIVES (coordonnées MTM zone 9, mètres) — aucune correspondance avec un cours d'eau ou une route réels COURS_EAU_FICTIF = [(351000, 5028200), (358000, 5029000), (363500, 5029400), (369000, 5031200), (376000, 5033500), (384000, 5036800), (393000, 5039500)] LARGEUR_CLASSE_A_M, LARGEUR_CLASSE_B_M = 350.0, 900.0 # classe A (intensité élevée) : ≤ 350 m ; classe B (modérée) : 350–900 m CORRIDOR_FICTIF = [(359500, 5025500), (362500, 5030000), (365500, 5034500), (369500, 5040500), (373000, 5046000)] MASSES_POLES = {"Pôle central": 100.0, "Pôle est": 40.0, "Pôle ouest": 30.0} FRICTION = 0.25 # λ (par km) de l'indice gravitaire exponentiel _CACHE: dict = {} def _gpd(): import geopandas as gpd return gpd # ----------------------------------------------------------------------------- 1. projections def charger_base() -> dict: """600 ventes du cas 1 et sujet en GeoDataFrame WGS84 → NAD83 / MTM zone 9 ; comparaisons de distances.""" from shapely.geometry import Point from pyproj import Geod, CRS gpd = _gpd() d = cas1.charger(); v = d["ventes"].copy(); s = d["sujet"] g = gpd.GeoDataFrame(v, geometry=gpd.points_from_xy(v["lon"], v["lat"]), crs=EPSG_WGS84).to_crs(EPSG_MTM9) gs = gpd.GeoSeries([Point(s["lon"], s["lat"])], crs=EPSG_WGS84).to_crs(EPSG_MTM9) xs, ys = float(gs.x.iloc[0]), float(gs.y.iloc[0]) xy_m = np.column_stack([g.geometry.x.values, g.geometry.y.values]) # distances sujet → ventes : MTM 9, équirectangulaire locale (ch. 19), haversine, géodésique WGS84 d_mtm = np.sqrt(((xy_m - [xs, ys]) ** 2).sum(1)) / 1000 lat0, lon0 = float(v["lat"].mean()), float(v["lon"].mean()) xy_eq = h.coordonnees_km(v["lat"].values, v["lon"].values, lat0, lon0); cs_eq = h.coordonnees_km([s["lat"]], [s["lon"]], lat0, lon0)[0] d_eq = np.sqrt(((xy_eq - cs_eq) ** 2).sum(1)) d_hav = geo.haversine_km(s["lat"], s["lon"], v["lat"].values, v["lon"].values) geod = Geod(ellps="WGS84"); _, _, d_geod = geod.inv(np.full(len(v), s["lon"]), np.full(len(v), s["lat"]), v["lon"].values, v["lat"].values) d_geod = d_geod / 1000 # distance en degrés « naïve » (erreur classique) d_deg = np.sqrt((v["lat"].values - s["lat"]) ** 2 + (v["lon"].values - s["lon"]) ** 2) k_sujet = float(geo.facteur_echelle_tm(s["lat"], s["lon"], LON0_MTM9)) crs = CRS.from_epsg(EPSG_MTM9) comp = pd.DataFrame({"mtm9": d_mtm, "equirect": d_eq, "haversine": d_hav, "geodesique": d_geod}) ecarts = {c: {"max_m": float(np.abs(comp[c] - comp["geodesique"]).max() * 1000), "max_rel_pct": float((np.abs(comp[c] / comp["geodesique"] - 1)).max() * 100)} for c in ["mtm9", "equirect", "haversine"]} # distance de Manhattan (réseau orthogonal idéalisé) vs euclidienne man = (np.abs(xy_m - [xs, ys]).sum(1)) / 1000; ratio = man / np.where(d_mtm > 0, d_mtm, np.nan) return {"ventes": g, "sujet": s, "xs": xs, "ys": ys, "xy_m": xy_m, "comparaison_distances": comp, "ecarts": ecarts, "k_sujet": k_sujet, "crs_nom": crs.name, "d_deg_vs_km": pd.DataFrame({"deg": d_deg, "km": d_geod}), "manhattan": {"moyenne_ratio": float(np.nanmean(ratio)), "max_ratio": float(np.nanmax(ratio)), "mediane_ratio": float(np.nanmedian(ratio))}, "etendue_km": ((xy_m[:, 0].max() - xy_m[:, 0].min()) / 1000, (xy_m[:, 1].max() - xy_m[:, 1].min()) / 1000)} # ----------------------------------------------------------------------------- 2. géométries fictives def geometries_fictives(b: dict) -> dict: from shapely.geometry import LineString, Point gpd = _gpd() riv = LineString(COURS_EAU_FICTIF); corr = LineString(CORRIDOR_FICTIF) zA = riv.buffer(LARGEUR_CLASSE_A_M); zB = riv.buffer(LARGEUR_CLASSE_B_M).difference(zA) zones = gpd.GeoDataFrame({"classe": ["A", "B"], "intensite": ["élevée (fictive)", "modérée (fictive)"]}, geometry=[zA, zB], crs=EPSG_MTM9) tampons = gpd.GeoDataFrame({"rayon_m": RAYONS_CORRIDOR_M}, geometry=[corr.buffer(r) for r in RAYONS_CORRIDOR_M], crs=EPSG_MTM9) v = b["ventes"] # pôles fictifs : le pôle central est placé à la vente de plus petite distance au centre du générateur ; les deux autres # aux médianes des coordonnées des secteurs Gatineau (est) et Aylmer (ouest) i0 = int(v["dist_centre_km"].values.argmin()); pc = (float(v.geometry.x.iloc[i0]), float(v.geometry.y.iloc[i0])) pe = (float(v.loc[v["secteur"] == "Gatineau"].geometry.x.median()), float(v.loc[v["secteur"] == "Gatineau"].geometry.y.median())) po = (float(v.loc[v["secteur"] == "Aylmer"].geometry.x.median()), float(v.loc[v["secteur"] == "Aylmer"].geometry.y.median())) poles = gpd.GeoDataFrame({"pole": list(MASSES_POLES), "masse": list(MASSES_POLES.values())}, geometry=[Point(pc), Point(pe), Point(po)], crs=EPSG_MTM9) return {"cours_eau": riv, "corridor": corr, "zones": zones, "tampons": tampons, "poles": poles, "aires_km2": {"A": zA.area / 1e6, "B": zB.area / 1e6, "corridor_300": corr.buffer(RAYON_CORRIDOR_REF_M).area / 1e6}, "longueurs_km": {"cours_eau": riv.length / 1000, "corridor": corr.length / 1000}} # ----------------------------------------------------------------------------- 3. opérations SIG def jointure_et_distances(b: dict, gf: dict) -> "pd.DataFrame": """Jointure spatiale (point dans polygone) avec les classes d'exposition, distances au corridor et au cours d'eau, indicatrices de tampon, accessibilité gravitaire.""" gpd = _gpd(); v = b["ventes"].copy() j = gpd.sjoin(v[["geometry"]], gf["zones"][["classe", "geometry"]], how="left", predicate="within") j = j[~j.index.duplicated(keep="first")] v["classe_expo"] = j["classe"].reindex(v.index).fillna("hors zone") v["zone_A"] = (v["classe_expo"] == "A").astype(float); v["zone_B"] = (v["classe_expo"] == "B").astype(float) v["dist_cours_eau_m"] = v.geometry.distance(gf["cours_eau"]); v["dist_corridor_m"] = v.geometry.distance(gf["corridor"]) for r in RAYONS_CORRIDOR_M: v[f"corr_{r}"] = (v["dist_corridor_m"] <= r).astype(float) v["ln_dist_corridor"] = np.log(1.0 + v["dist_corridor_m"] / 100.0) # ln(1 + d/100 m) coords = np.column_stack([v.geometry.x.values, v.geometry.y.values]) / 1000 pc = np.column_stack([gf["poles"].geometry.x.values, gf["poles"].geometry.y.values]) / 1000 v["acces"] = geo.accessibilite_gravitaire(coords, pc, gf["poles"]["masse"].values, FRICTION) v["ln_acces"] = np.log(v["acces"]); v["dist_pole_min_km"] = np.sqrt(((coords[:, None, :] - pc[None, :, :]) ** 2).sum(-1)).min(1) # contrôle : le lancer de rayon de traite.geo donne le même résultat que shapely pour la zone A (polygone extérieur sans trou) ext = np.asarray(gf["zones"].geometry.iloc[0].exterior.coords)[:-1] v["zone_A_rayon"] = geo.point_dans_polygone(np.column_stack([v.geometry.x.values, v.geometry.y.values]), ext).astype(float) return v def tableau_secteurs(v) -> pd.DataFrame: g = v.groupby("secteur").agg(n=("id_vente", "count"), zone_A=("zone_A", "sum"), zone_B=("zone_B", "sum"), corridor_300=(f"corr_{RAYON_CORRIDOR_REF_M}", "sum"), dist_corr_med_km=("dist_corridor_m", lambda x: x.median() / 1000), acces_moy=("acces", "mean"), prix_m2_med=("prix_vente", lambda x: np.nan)) for sct, w in v.groupby("secteur"): g.loc[sct, "prix_m2_med"] = float((w["prix_vente"] / w["superficie_habitable_m2"]).median()) g = g.reindex(cas1.SECTEURS) tot = pd.Series({"n": len(v), "zone_A": v["zone_A"].sum(), "zone_B": v["zone_B"].sum(), "corridor_300": v[f"corr_{RAYON_CORRIDOR_REF_M}"].sum(), "dist_corr_med_km": v["dist_corridor_m"].median() / 1000, "acces_moy": v["acces"].mean(), "prix_m2_med": float((v["prix_vente"] / v["superficie_habitable_m2"]).median())}, name="Total") return pd.concat([g, tot.to_frame().T]) def grille_et_agregation(b: dict, v) -> dict: """Agrégation des ventes par cellule hexagonale (médiane du prix au m², effectif, résidu moyen du modèle de référence), puis pondération surfacique vers une grille carrée grossière comparée au comptage exact.""" from shapely.geometry import Polygon, box gpd = _gpd() xy = b["xy_m"]; xmin, ymin, xmax, ymax = xy[:, 0].min(), xy[:, 1].min(), xy[:, 0].max(), xy[:, 1].max() cells = geo.grille_hexagonale(xmin - RAYON_HEX_M, ymin - RAYON_HEX_M, xmax, ymax, RAYON_HEX_M) hexes = gpd.GeoDataFrame({"cellule": range(len(cells))}, geometry=[Polygon(c) for c in cells], crs=EPSG_MTM9) m = cas1.hedonique(b["ventes"].drop(columns="geometry")) pts = v[["geometry", "prix_vente", "superficie_habitable_m2"]].copy(); pts["prix_m2"] = pts["prix_vente"] / pts["superficie_habitable_m2"] pts["residu"] = m["residus"] j = gpd.sjoin(pts, hexes, how="inner", predicate="within") agg = j.groupby("cellule").agg(n=("prix_m2", "count"), prix_m2_med=("prix_m2", "median"), residu_moy=("residu", "mean")) hexes = hexes.join(agg, on="cellule"); hexes["n"] = hexes["n"].fillna(0).astype(int) occ = hexes[hexes["n"] > 0].copy(); ok = occ[occ["n"] >= MIN_VENTES_CELLULE].copy() # classification de la médiane du prix au m² (cellules avec n ≥ 5) z = ok["prix_m2_med"].values bornes = {"intervalles égaux": geo.bornes_intervalles_egaux(z, K_CLASSES), "quantiles": geo.bornes_quantiles(z, K_CLASSES), "Jenks": geo.bornes_jenks(z, K_CLASSES)} classif = {k: {"bornes": bv, "gvf": geo.gvf(z, bv), "effectifs": np.bincount(geo.classer(z, bv), minlength=K_CLASSES)} for k, bv in bornes.items()} # pondération surfacique : des hexagones (n) vers une grille carrée de 4 km ; comparaison au comptage exact des points pas = 4000.0; cols = np.arange(xmin - 0.5 * pas, xmax + pas, pas); rows = np.arange(ymin - 0.5 * pas, ymax + pas, pas) # origine décalée : aucun point sur une frontière carres = gpd.GeoDataFrame(geometry=[box(x, y, x + pas, y + pas) for x in cols for y in rows], crs=EPSG_MTM9) carres["exact"] = [int(v.geometry.within(c).sum()) for c in carres.geometry] carres = carres[carres["exact"] > 0].copy(); est = [] for c in carres.geometry: tot = 0.0 for _, hx in occ.iterrows(): inter = hx.geometry.intersection(c).area if inter > 0: tot += hx["n"] * inter / hx.geometry.area est.append(tot) carres["pond_surf"] = est; carres["erreur"] = carres["pond_surf"] - carres["exact"] assert int(hexes["n"].sum()) == len(v), "toutes les ventes doivent tomber dans une cellule" return {"hexes": hexes, "occupees": occ, "cartographiees": ok, "classif": classif, "n_cellules": len(hexes), "n_occupees": len(occ), "n_cartographiees": len(ok), "aire_cellule_km2": float(hexes.geometry.iloc[0].area / 1e6), "carres": carres, "pond_surf": {"eam": float(carres["erreur"].abs().mean()), "max": float(carres["erreur"].abs().max()), "total_exact": int(carres["exact"].sum()), "total_est": float(carres["pond_surf"].sum()), "n_carres": len(carres)}, "modele": m} def interpolation_idw(b: dict, grille: dict) -> dict: """IDW du résidu du modèle sans secteur au sujet (comparaison avec le krigeage du chapitre 19) et validation croisée.""" m = grille["modele"]; Xr = m["X"].drop(columns=[c for c in m["X"].columns if c.startswith("sect_")]) r = h.mco(Xr.values, m["y"]); res = r["residus"]; xy = b["xy_m"] / 1000; cs = np.array([[b["xs"], b["ys"]]]) / 1000 out = {} for p in [1.0, 2.0, 3.0]: au_sujet = float(geo.idw(xy, res, cs, puissance=p, voisins=40)[0]) loo = np.array([geo.idw(np.delete(xy, i, 0), np.delete(res, i), xy[i:i + 1], puissance=p, voisins=40)[0] for i in range(len(res))]) out[p] = {"sujet": au_sujet, "eqm_loo": float(np.sqrt(np.mean((res - loo) ** 2)))} out["eqm_sans_correction"] = float(np.sqrt(np.mean(res ** 2))) return out # ----------------------------------------------------------------------------- 5. modèle hédonique augmenté def _X_ref(v): return cas1.matrice(v.drop(columns="geometry") if "geometry" in v else v) def modeles_augmentes(v) -> dict: """Modèle de référence + variables géospatiales (indicatrices d'exposition, corridor ≤ 300 m, ln distance au corridor, ln accessibilité), pour trois spécifications de la localisation : référence (secteurs + gradient), sans secteur, sans secteur ni gradient.""" X = _X_ref(v); y = np.log(v["prix_vente"].values.astype(float)) extra = {"zone_A": v["zone_A"].values, "zone_B": v["zone_B"].values, f"corr_{RAYON_CORRIDOR_REF_M}": v[f"corr_{RAYON_CORRIDOR_REF_M}"].values, "ln_dist_corridor": v["ln_dist_corridor"].values, "ln_acces": v["ln_acces"].values} specs = {"référence": X, "sans secteur": X.drop(columns=[c for c in X.columns if c.startswith("sect_")]), "sans secteur ni gradient": X.drop(columns=[c for c in X.columns if c.startswith("sect_")] + ["dist_centre"])} from scipy.stats import t as student out = {} for nom, Xs in specs.items(): Xa = Xs.copy() for k, col in extra.items(): Xa[k] = col r = h.mco(Xa.values, y); n, k = r["n"], r["k"] coef = pd.DataFrame({"coef": r["beta"], "se": r["se_hc1"]}, index=Xa.columns); coef["t"] = coef["coef"] / coef["se"] coef["p"] = 2 * student.sf(np.abs(coef["t"]), n - k); coef["ic_bas"] = coef["coef"] - 1.96 * coef["se"]; coef["ic_haut"] = coef["coef"] + 1.96 * coef["se"] coef["effet_pct"] = (np.exp(coef["coef"]) - 1) * 100 base = h.mco(Xs.values, y) # test F conjoint des cinq variables géospatiales (homoscédastique, à titre indicatif) q = len(extra); rss_r = base["residus"] @ base["residus"]; rss_u = r["residus"] @ r["residus"] F = ((rss_r - rss_u) / q) / (rss_u / (n - k)) from scipy.stats import f as fdist out[nom] = {"coef": coef.loc[list(extra)], "r2": r["r2"], "r2_base": base["r2"], "sigma": float(np.sqrt(r["sigma2"])), "F": float(F), "p_F": float(fdist.sf(F, q, n - k)), "n": n, "k": k, "mde": {k2: geo.effet_minimal_detectable(float(coef.loc[k2, "se"])) for k2 in extra}} return out def sensibilite_tampon(v) -> pd.DataFrame: """Coefficient de l'indicatrice « à moins de r m du corridor » selon le rayon r, modèle de référence et modèle sans secteur.""" X = _X_ref(v); y = np.log(v["prix_vente"].values.astype(float)); Xr = X.drop(columns=[c for c in X.columns if c.startswith("sect_")]) rows = [] for r in RAYONS_CORRIDOR_M: col = v[f"corr_{r}"].values; n_in = int(col.sum()) for nom, Xs in [("référence", X), ("sans secteur", Xr)]: Xa = Xs.copy(); Xa["corr"] = col; m = h.mco(Xa.values, y); i = list(Xa.columns).index("corr") rows.append({"rayon_m": r, "n_dans_tampon": n_in, "spec": nom, "coef": m["beta"][i], "se": m["se_hc1"][i], "t": m["beta"][i] / m["se_hc1"][i]}) return pd.DataFrame(rows) # ----------------------------------------------------------------------------- 6. faux positifs spatiaux def faux_positifs(b: dict, v, n_poly: int = N_POLYGONES_ALEATOIRES, graine: int = GRAINE) -> dict: """Pour chaque rayon, n_poly disques centrés sur une vente tirée au hasard ; indicatrice « dans le disque » ajoutée au modèle ; part des |t| > 1,96 (et > 2,576) selon la spécification. Sous le générateur (aucun effet), la part attendue est 5 % (1 %) si la spécification est correcte ; les résidus spatialement corrélés des spécifications sans secteur gonflent cette part.""" rng = np.random.default_rng(graine); X = _X_ref(v); y = np.log(v["prix_vente"].values.astype(float)); xy = b["xy_m"] specs = {"référence": X, "sans secteur": X.drop(columns=[c for c in X.columns if c.startswith("sect_")]), "sans secteur ni gradient": X.drop(columns=[c for c in X.columns if c.startswith("sect_")] + ["dist_centre"])} rows = []; tstats = {s: {r: [] for r in RAYONS_POLYGONES_KM} for s in specs} for r in RAYONS_POLYGONES_KM: centres = rng.integers(0, len(v), size=n_poly); n_ok = 0 for c in centres: d = np.sqrt(((xy - xy[c]) ** 2).sum(1)) / 1000; ind = (d <= r).astype(float); n_in = int(ind.sum()) if n_in < 10 or n_in > len(v) - 10: continue n_ok += 1 for s, Xs in specs.items(): Xa = np.column_stack([Xs.values, ind]); m = h.mco(Xa, y) tstats[s][r].append((float(m["beta"][-1] / m["se"][-1]), float(m["beta"][-1] / m["se_hc1"][-1]))) # (t classique, t HC1) for s in specs: tt = np.array(tstats[s][r]); seuil_bonf = _z(0.05 / len(tt)) for j, et in enumerate(["classique", "HC1"]): t = tt[:, j] rows.append({"rayon_km": r, "spec": s, "erreur_type": et, "n_polygones": len(t), "part_5pct": float(np.mean(np.abs(t) > 1.96)), "part_1pct": float(np.mean(np.abs(t) > 2.576)), "part_bonferroni": float(np.mean(np.abs(t) > seuil_bonf)), "ecart_type_t": float(t.std()), "t_max_abs": float(np.abs(t).max())}) return {"tableau": pd.DataFrame(rows), "t": tstats} def _z(alpha_bilateral: float) -> float: from scipy.stats import norm return float(norm.ppf(1 - alpha_bilateral / 2)) # ----------------------------------------------------------------------------- 8. cas 4 def cas4_geo() -> dict: """48 immeubles du cas 4 : projection en NAD83 / Québec Lambert (EPSG 32198), part de la valeur en zone inondable par RMR, coefficient « inondable » des effets régionaux, distances entre centroïdes de RMR.""" from pyproj import Geod gpd = _gpd(); d = cas4.charger(); pf = d["portefeuille"]; som = cas4.sommaire(pf); eff = cas4.effets_regionaux(pf) g = gpd.GeoDataFrame(pf, geometry=gpd.points_from_xy(pf["lon"], pf["lat"]), crs=EPSG_WGS84).to_crs(EPSG_LAMBERT) par = pf.groupby("rmr").agg(n=("id_immeuble", "count"), valeur=("valeur_2026", "sum"), n_inond=("zone_inondable", "sum"), n_sol=("contamination_sol", "sum"), val_inond=("valeur_2026", lambda x: x[pf.loc[x.index, "zone_inondable"] == 1].sum()), lat=("lat", "mean"), lon=("lon", "mean")) par["part_inond_rmr"] = par["val_inond"] / par["valeur"]; par["part_inond_total"] = par["val_inond"] / pf["valeur_2026"].sum() par = par.sort_values("valeur", ascending=False) geod = Geod(ellps="WGS84"); cen = par[["lat", "lon"]] dist = pd.DataFrame(index=cen.index, columns=cen.index, dtype=float) for a in cen.index: for bb in cen.index: dist.loc[a, bb] = geod.inv(cen.loc[a, "lon"], cen.loc[a, "lat"], cen.loc[bb, "lon"], cen.loc[bb, "lat"])[2] / 1000 inond = pf[pf["zone_inondable"] == 1][["id_immeuble", "rmr", "nb_logements", "valeur_2026", "taux_global", "etat_1a5"]] # taux global moyen des immeubles inondables vs non, à état et RMR comparables : le coefficient du modèle le mesure return {"gdf": g, "par_rmr": par, "total": som["total"], "effets": eff, "coef_inondable": float(eff.loc["inondable", "coef"]), "se_inondable": float(eff.loc["inondable", "se_hc1"]), "effet_pct_inondable": float(eff.loc["inondable", "effet_pct"] * 100), "coef_sol": float(eff.loc["sol", "coef"]), "effet_pct_sol": float(eff.loc["sol", "effet_pct"] * 100), "distances_km": dist, "inondables": inond, "n_inond": int(pf["zone_inondable"].sum()), "n_sol": int(pf["contamination_sol"].sum()), "etendue_km": (float((g.geometry.x.max() - g.geometry.x.min()) / 1000), float((g.geometry.y.max() - g.geometry.y.min()) / 1000))} # ----------------------------------------------------------------------------- capitalisation du risque (exemple numérique du texte) def exemple_capitalisation(P25=(0.70, 0.20, 0.07), dommage_pct=(0.15, 0.08, 0.03), R=0.055, prix=637953.0) -> pd.DataFrame: """Décote capitalisée d'une perte annuelle espérée : ΔV/V = −p_1 · D / R, avec p_1 déduit de la probabilité sur 25 ans (bornes de classe).""" rows = [] for p, D in zip(P25, dommage_pct): p1 = geo.probabilite_annuelle(p); rows.append({"P25": p, "p_annuelle": p1, "dommage_pct": D, "perte_annuelle_pct": p1 * D, "decote_pct": p1 * D / R, "decote_dollars": p1 * D / R * prix}) return pd.DataFrame(rows) # ----------------------------------------------------------------------------- orchestration def tout() -> dict: if "tout" in _CACHE: return _CACHE["tout"] warnings.filterwarnings("ignore") b = charger_base(); gf = geometries_fictives(b); v = jointure_et_distances(b, gf) grille = grille_et_agregation(b, v); idw_ = interpolation_idw(b, grille) res = {"base": b, "geom": gf, "ventes": v, "secteurs": tableau_secteurs(v), "grille": grille, "idw": idw_, "modeles": modeles_augmentes(v), "tampon": sensibilite_tampon(v), "faux_positifs": faux_positifs(b, v), "cas4": cas4_geo(), "capitalisation": exemple_capitalisation()} _CACHE["tout"] = res return res def main(): r = tout(); b = r["base"]; pd.set_option("display.width", 200); pd.set_option("display.max_columns", 30) print("=== 1. Projection (NAD83 / MTM zone 9, EPSG 32189) :", b["crs_nom"]) print(f"sujet : x = {b['xs']:.1f} m, y = {b['ys']:.1f} m ; facteur d'échelle k = {b['k_sujet']:.6f} ; étendue {b['etendue_km'][0]:.1f} × {b['etendue_km'][1]:.1f} km") print("écarts max à la géodésique WGS84 (sujet → 600 ventes) :", {k: (round(x["max_m"], 1), round(x["max_rel_pct"], 4)) for k, x in b["ecarts"].items()}) print("Manhattan / euclidienne :", {k: round(x, 3) for k, x in b["manhattan"].items()}) print("\n=== 2. Géométries fictives :", {k: round(x, 2) for k, x in r["geom"]["aires_km2"].items()}, {k: round(x, 1) for k, x in r["geom"]["longueurs_km"].items()}) print("\n=== 3. Jointure spatiale et distances par secteur\n", r["secteurs"].round(3).to_string()) v = r["ventes"]; print("contrôle lancer de rayon = shapely :", bool((v["zone_A"] == v["zone_A_rayon"]).all())) g = r["grille"]; print(f"\n=== 4. Grille hexagonale : {g['n_cellules']} cellules ({g['aire_cellule_km2']:.2f} km²), {g['n_occupees']} occupées, {g['n_cartographiees']} avec n ≥ {MIN_VENTES_CELLULE}") for k, c in g["classif"].items(): print(f" {k:20s} bornes = {np.round(c['bornes'], 0)} ; effectifs = {c['effectifs']} ; GVF = {c['gvf']:.3f}") print(" pondération surfacique :", {k: (round(x, 2) if isinstance(x, float) else x) for k, x in g["pond_surf"].items()}) print("\n=== IDW du résidu sans secteur au sujet :", {k: (x if not isinstance(x, dict) else {kk: round(vv, 4) for kk, vv in x.items()}) for k, x in r["idw"].items()}) print("\n=== 5. Modèles augmentés") for nom, m in r["modeles"].items(): print(f"-- {nom} : R² = {m['r2']:.4f} (base {m['r2_base']:.4f}), σ = {m['sigma']:.4f}, F({len(m['coef'])}, {m['n'] - m['k']}) = {m['F']:.3f}, p = {m['p_F']:.3f}") print(m["coef"].round(4).to_string()); print(" EMD (80 %) :", {k: round(x, 4) for k, x in m["mde"].items()}) print("\n=== 6. Sensibilité au rayon du tampon\n", r["tampon"].round(4).to_string()) print("\n=== 7. Faux positifs spatiaux\n", r["faux_positifs"]["tableau"].round(3).to_string()) c4 = r["cas4"]; print("\n=== 8. Cas 4 par RMR\n", c4["par_rmr"].round(4).to_string()) print(f"part de la valeur en zone inondable : {c4['total']['part_zone_inondable']*100:.2f} % ; {c4['n_inond']} immeubles ; coef inondable = {c4['coef_inondable']:.4f} ({c4['se_inondable']:.4f}) → {c4['effet_pct_inondable']:.2f} % ; sol {c4['effet_pct_sol']:.2f} %") print("distances géodésiques entre centroïdes de RMR (km)\n", c4["distances_km"].round(0).to_string()) print("\n=== Capitalisation du risque (exemple)\n", r["capitalisation"].round(4).to_string()) if __name__ == "__main__": main()