spb/traite-evaluation-qc
Public
1#!/usr/bin/env python32"""Chapitre 38 — SIG et risques géospatiaux (cas 1 Gatineau et cas 4 portefeuille, DONNÉES SYNTHÉTIQUES, graine 20261031).34Reproduit tous les chiffres cités dans manuscript/chapters/ch38.tex :5 1. systèmes de référence : projection des 600 ventes du cas 1 en NAD83 / MTM zone 9 (EPSG 32189) ; comparaison avec la6 projection équirectangulaire locale du chapitre 19, le haversine et la géodésique WGS84 ; facteur d'échelle ;7 2. géométries FICTIVES déclarées : « cours d'eau fictif » et ses deux classes d'exposition (tampons), « corridor de8 transport fictif » et ses tampons, trois « pôles » fictifs pour l'accessibilité ;9 3. opérations SIG : jointure spatiale (point dans polygone), tampons, distances (euclidienne, Manhattan, au corridor),10 agrégation par cellule hexagonale, pondération surfacique, interpolation IDW ;11 4. accessibilité gravitaire (Hansen) ;12 5. intégration au modèle hédonique de référence (traite.cas1.hedonique) : variables d'exposition, de corridor et13 d'accessibilité — l'effet attendu est NUL (le générateur ne contient pas ces variables) ;14 6. faux positifs spatiaux : polygones aléatoires et part de coefficients « significatifs » selon la spécification ;15 7. sensibilité au rayon du tampon et à la méthode de classification des cartes choroplèthes ;16 8. cas 4 : 48 immeubles par RMR, part de la valeur en zone inondable, coefficient « inondable » des effets régionaux.17Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch38_sig.py"""18from __future__ import annotations19import sys, pathlib, warnings20import numpy as np, pandas as pd21sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))22from traite import cas1, cas4, geo23from traite import hedonique as h2425ROOT = pathlib.Path(__file__).resolve().parents[3]26GRAINE = 2026103127EPSG_WGS84, EPSG_MTM9, EPSG_LAMBERT = 4326, 32189, 3219828LON0_MTM9 = -76.5 # méridien central de la zone MTM 929RAYON_HEX_M = 1200.0 # rayon circonscrit des cellules hexagonales (m) ; aire ≈ 3,74 km²30MIN_VENTES_CELLULE = 531RAYONS_CORRIDOR_M = [100, 150, 200, 300, 400, 500, 750, 1000]32RAYON_CORRIDOR_REF_M = 30033N_POLYGONES_ALEATOIRES = 50034RAYONS_POLYGONES_KM = [0.5, 1.0, 2.0]35K_CLASSES = 536# géométries FICTIVES (coordonnées MTM zone 9, mètres) — aucune correspondance avec un cours d'eau ou une route réels37COURS_EAU_FICTIF = [(351000, 5028200), (358000, 5029000), (363500, 5029400), (369000, 5031200), (376000, 5033500), (384000, 5036800), (393000, 5039500)]38LARGEUR_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 m39CORRIDOR_FICTIF = [(359500, 5025500), (362500, 5030000), (365500, 5034500), (369500, 5040500), (373000, 5046000)]40MASSES_POLES = {"Pôle central": 100.0, "Pôle est": 40.0, "Pôle ouest": 30.0}41FRICTION = 0.25 # λ (par km) de l'indice gravitaire exponentiel4243_CACHE: dict = {}444546def _gpd():47 import geopandas as gpd48 return gpd495051# ----------------------------------------------------------------------------- 1. projections5253def charger_base() -> dict:54 """600 ventes du cas 1 et sujet en GeoDataFrame WGS84 → NAD83 / MTM zone 9 ; comparaisons de distances."""55 from shapely.geometry import Point56 from pyproj import Geod, CRS57 gpd = _gpd()58 d = cas1.charger(); v = d["ventes"].copy(); s = d["sujet"]59 g = gpd.GeoDataFrame(v, geometry=gpd.points_from_xy(v["lon"], v["lat"]), crs=EPSG_WGS84).to_crs(EPSG_MTM9)60 gs = gpd.GeoSeries([Point(s["lon"], s["lat"])], crs=EPSG_WGS84).to_crs(EPSG_MTM9)61 xs, ys = float(gs.x.iloc[0]), float(gs.y.iloc[0])62 xy_m = np.column_stack([g.geometry.x.values, g.geometry.y.values])63 # distances sujet → ventes : MTM 9, équirectangulaire locale (ch. 19), haversine, géodésique WGS8464 d_mtm = np.sqrt(((xy_m - [xs, ys]) ** 2).sum(1)) / 100065 lat0, lon0 = float(v["lat"].mean()), float(v["lon"].mean())66 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]67 d_eq = np.sqrt(((xy_eq - cs_eq) ** 2).sum(1))68 d_hav = geo.haversine_km(s["lat"], s["lon"], v["lat"].values, v["lon"].values)69 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)70 d_geod = d_geod / 100071 # distance en degrés « naïve » (erreur classique)72 d_deg = np.sqrt((v["lat"].values - s["lat"]) ** 2 + (v["lon"].values - s["lon"]) ** 2)73 k_sujet = float(geo.facteur_echelle_tm(s["lat"], s["lon"], LON0_MTM9))74 crs = CRS.from_epsg(EPSG_MTM9)75 comp = pd.DataFrame({"mtm9": d_mtm, "equirect": d_eq, "haversine": d_hav, "geodesique": d_geod})76 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)}77 for c in ["mtm9", "equirect", "haversine"]}78 # distance de Manhattan (réseau orthogonal idéalisé) vs euclidienne79 man = (np.abs(xy_m - [xs, ys]).sum(1)) / 1000; ratio = man / np.where(d_mtm > 0, d_mtm, np.nan)80 return {"ventes": g, "sujet": s, "xs": xs, "ys": ys, "xy_m": xy_m, "comparaison_distances": comp, "ecarts": ecarts,81 "k_sujet": k_sujet, "crs_nom": crs.name, "d_deg_vs_km": pd.DataFrame({"deg": d_deg, "km": d_geod}),82 "manhattan": {"moyenne_ratio": float(np.nanmean(ratio)), "max_ratio": float(np.nanmax(ratio)), "mediane_ratio": float(np.nanmedian(ratio))},83 "etendue_km": ((xy_m[:, 0].max() - xy_m[:, 0].min()) / 1000, (xy_m[:, 1].max() - xy_m[:, 1].min()) / 1000)}848586# ----------------------------------------------------------------------------- 2. géométries fictives8788def geometries_fictives(b: dict) -> dict:89 from shapely.geometry import LineString, Point90 gpd = _gpd()91 riv = LineString(COURS_EAU_FICTIF); corr = LineString(CORRIDOR_FICTIF)92 zA = riv.buffer(LARGEUR_CLASSE_A_M); zB = riv.buffer(LARGEUR_CLASSE_B_M).difference(zA)93 zones = gpd.GeoDataFrame({"classe": ["A", "B"], "intensite": ["élevée (fictive)", "modérée (fictive)"]}, geometry=[zA, zB], crs=EPSG_MTM9)94 tampons = gpd.GeoDataFrame({"rayon_m": RAYONS_CORRIDOR_M}, geometry=[corr.buffer(r) for r in RAYONS_CORRIDOR_M], crs=EPSG_MTM9)95 v = b["ventes"]96 # pôles fictifs : le pôle central est placé à la vente de plus petite distance au centre du générateur ; les deux autres97 # aux médianes des coordonnées des secteurs Gatineau (est) et Aylmer (ouest)98 i0 = int(v["dist_centre_km"].values.argmin()); pc = (float(v.geometry.x.iloc[i0]), float(v.geometry.y.iloc[i0]))99 pe = (float(v.loc[v["secteur"] == "Gatineau"].geometry.x.median()), float(v.loc[v["secteur"] == "Gatineau"].geometry.y.median()))100 po = (float(v.loc[v["secteur"] == "Aylmer"].geometry.x.median()), float(v.loc[v["secteur"] == "Aylmer"].geometry.y.median()))101 poles = gpd.GeoDataFrame({"pole": list(MASSES_POLES), "masse": list(MASSES_POLES.values())}, geometry=[Point(pc), Point(pe), Point(po)], crs=EPSG_MTM9)102 return {"cours_eau": riv, "corridor": corr, "zones": zones, "tampons": tampons, "poles": poles,103 "aires_km2": {"A": zA.area / 1e6, "B": zB.area / 1e6, "corridor_300": corr.buffer(RAYON_CORRIDOR_REF_M).area / 1e6},104 "longueurs_km": {"cours_eau": riv.length / 1000, "corridor": corr.length / 1000}}105106107# ----------------------------------------------------------------------------- 3. opérations SIG108109def jointure_et_distances(b: dict, gf: dict) -> "pd.DataFrame":110 """Jointure spatiale (point dans polygone) avec les classes d'exposition, distances au corridor et au cours d'eau,111 indicatrices de tampon, accessibilité gravitaire."""112 gpd = _gpd(); v = b["ventes"].copy()113 j = gpd.sjoin(v[["geometry"]], gf["zones"][["classe", "geometry"]], how="left", predicate="within")114 j = j[~j.index.duplicated(keep="first")]115 v["classe_expo"] = j["classe"].reindex(v.index).fillna("hors zone")116 v["zone_A"] = (v["classe_expo"] == "A").astype(float); v["zone_B"] = (v["classe_expo"] == "B").astype(float)117 v["dist_cours_eau_m"] = v.geometry.distance(gf["cours_eau"]); v["dist_corridor_m"] = v.geometry.distance(gf["corridor"])118 for r in RAYONS_CORRIDOR_M:119 v[f"corr_{r}"] = (v["dist_corridor_m"] <= r).astype(float)120 v["ln_dist_corridor"] = np.log(1.0 + v["dist_corridor_m"] / 100.0) # ln(1 + d/100 m)121 coords = np.column_stack([v.geometry.x.values, v.geometry.y.values]) / 1000122 pc = np.column_stack([gf["poles"].geometry.x.values, gf["poles"].geometry.y.values]) / 1000123 v["acces"] = geo.accessibilite_gravitaire(coords, pc, gf["poles"]["masse"].values, FRICTION)124 v["ln_acces"] = np.log(v["acces"]); v["dist_pole_min_km"] = np.sqrt(((coords[:, None, :] - pc[None, :, :]) ** 2).sum(-1)).min(1)125 # 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)126 ext = np.asarray(gf["zones"].geometry.iloc[0].exterior.coords)[:-1]127 v["zone_A_rayon"] = geo.point_dans_polygone(np.column_stack([v.geometry.x.values, v.geometry.y.values]), ext).astype(float)128 return v129130131def tableau_secteurs(v) -> pd.DataFrame:132 g = v.groupby("secteur").agg(n=("id_vente", "count"), zone_A=("zone_A", "sum"), zone_B=("zone_B", "sum"),133 corridor_300=(f"corr_{RAYON_CORRIDOR_REF_M}", "sum"), dist_corr_med_km=("dist_corridor_m", lambda x: x.median() / 1000),134 acces_moy=("acces", "mean"), prix_m2_med=("prix_vente", lambda x: np.nan))135 for sct, w in v.groupby("secteur"):136 g.loc[sct, "prix_m2_med"] = float((w["prix_vente"] / w["superficie_habitable_m2"]).median())137 g = g.reindex(cas1.SECTEURS)138 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(),139 "dist_corr_med_km": v["dist_corridor_m"].median() / 1000, "acces_moy": v["acces"].mean(),140 "prix_m2_med": float((v["prix_vente"] / v["superficie_habitable_m2"]).median())}, name="Total")141 return pd.concat([g, tot.to_frame().T])142143144def grille_et_agregation(b: dict, v) -> dict:145 """Agrégation des ventes par cellule hexagonale (médiane du prix au m², effectif, résidu moyen du modèle de référence),146 puis pondération surfacique vers une grille carrée grossière comparée au comptage exact."""147 from shapely.geometry import Polygon, box148 gpd = _gpd()149 xy = b["xy_m"]; xmin, ymin, xmax, ymax = xy[:, 0].min(), xy[:, 1].min(), xy[:, 0].max(), xy[:, 1].max()150 cells = geo.grille_hexagonale(xmin - RAYON_HEX_M, ymin - RAYON_HEX_M, xmax, ymax, RAYON_HEX_M)151 hexes = gpd.GeoDataFrame({"cellule": range(len(cells))}, geometry=[Polygon(c) for c in cells], crs=EPSG_MTM9)152 m = cas1.hedonique(b["ventes"].drop(columns="geometry"))153 pts = v[["geometry", "prix_vente", "superficie_habitable_m2"]].copy(); pts["prix_m2"] = pts["prix_vente"] / pts["superficie_habitable_m2"]154 pts["residu"] = m["residus"]155 j = gpd.sjoin(pts, hexes, how="inner", predicate="within")156 agg = j.groupby("cellule").agg(n=("prix_m2", "count"), prix_m2_med=("prix_m2", "median"), residu_moy=("residu", "mean"))157 hexes = hexes.join(agg, on="cellule"); hexes["n"] = hexes["n"].fillna(0).astype(int)158 occ = hexes[hexes["n"] > 0].copy(); ok = occ[occ["n"] >= MIN_VENTES_CELLULE].copy()159 # classification de la médiane du prix au m² (cellules avec n ≥ 5)160 z = ok["prix_m2_med"].values161 bornes = {"intervalles égaux": geo.bornes_intervalles_egaux(z, K_CLASSES), "quantiles": geo.bornes_quantiles(z, K_CLASSES), "Jenks": geo.bornes_jenks(z, K_CLASSES)}162 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()}163 # pondération surfacique : des hexagones (n) vers une grille carrée de 4 km ; comparaison au comptage exact des points164 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ère165 carres = gpd.GeoDataFrame(geometry=[box(x, y, x + pas, y + pas) for x in cols for y in rows], crs=EPSG_MTM9)166 carres["exact"] = [int(v.geometry.within(c).sum()) for c in carres.geometry]167 carres = carres[carres["exact"] > 0].copy(); est = []168 for c in carres.geometry:169 tot = 0.0170 for _, hx in occ.iterrows():171 inter = hx.geometry.intersection(c).area172 if inter > 0:173 tot += hx["n"] * inter / hx.geometry.area174 est.append(tot)175 carres["pond_surf"] = est; carres["erreur"] = carres["pond_surf"] - carres["exact"]176 assert int(hexes["n"].sum()) == len(v), "toutes les ventes doivent tomber dans une cellule"177 return {"hexes": hexes, "occupees": occ, "cartographiees": ok, "classif": classif, "n_cellules": len(hexes), "n_occupees": len(occ), "n_cartographiees": len(ok),178 "aire_cellule_km2": float(hexes.geometry.iloc[0].area / 1e6), "carres": carres,179 "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)},180 "modele": m}181182183def interpolation_idw(b: dict, grille: dict) -> dict:184 """IDW du résidu du modèle sans secteur au sujet (comparaison avec le krigeage du chapitre 19) et validation croisée."""185 m = grille["modele"]; Xr = m["X"].drop(columns=[c for c in m["X"].columns if c.startswith("sect_")])186 r = h.mco(Xr.values, m["y"]); res = r["residus"]; xy = b["xy_m"] / 1000; cs = np.array([[b["xs"], b["ys"]]]) / 1000187 out = {}188 for p in [1.0, 2.0, 3.0]:189 au_sujet = float(geo.idw(xy, res, cs, puissance=p, voisins=40)[0])190 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))])191 out[p] = {"sujet": au_sujet, "eqm_loo": float(np.sqrt(np.mean((res - loo) ** 2)))}192 out["eqm_sans_correction"] = float(np.sqrt(np.mean(res ** 2)))193 return out194195196# ----------------------------------------------------------------------------- 5. modèle hédonique augmenté197198def _X_ref(v):199 return cas1.matrice(v.drop(columns="geometry") if "geometry" in v else v)200201202def modeles_augmentes(v) -> dict:203 """Modèle de référence + variables géospatiales (indicatrices d'exposition, corridor ≤ 300 m, ln distance au corridor, ln accessibilité),204 pour trois spécifications de la localisation : référence (secteurs + gradient), sans secteur, sans secteur ni gradient."""205 X = _X_ref(v); y = np.log(v["prix_vente"].values.astype(float))206 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,207 "ln_dist_corridor": v["ln_dist_corridor"].values, "ln_acces": v["ln_acces"].values}208 specs = {"référence": X, "sans secteur": X.drop(columns=[c for c in X.columns if c.startswith("sect_")]),209 "sans secteur ni gradient": X.drop(columns=[c for c in X.columns if c.startswith("sect_")] + ["dist_centre"])}210 from scipy.stats import t as student211 out = {}212 for nom, Xs in specs.items():213 Xa = Xs.copy()214 for k, col in extra.items():215 Xa[k] = col216 r = h.mco(Xa.values, y); n, k = r["n"], r["k"]217 coef = pd.DataFrame({"coef": r["beta"], "se": r["se_hc1"]}, index=Xa.columns); coef["t"] = coef["coef"] / coef["se"]218 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"]219 coef["effet_pct"] = (np.exp(coef["coef"]) - 1) * 100220 base = h.mco(Xs.values, y)221 # test F conjoint des cinq variables géospatiales (homoscédastique, à titre indicatif)222 q = len(extra); rss_r = base["residus"] @ base["residus"]; rss_u = r["residus"] @ r["residus"]223 F = ((rss_r - rss_u) / q) / (rss_u / (n - k))224 from scipy.stats import f as fdist225 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,226 "mde": {k2: geo.effet_minimal_detectable(float(coef.loc[k2, "se"])) for k2 in extra}}227 return out228229230def sensibilite_tampon(v) -> pd.DataFrame:231 """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."""232 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_")])233 rows = []234 for r in RAYONS_CORRIDOR_M:235 col = v[f"corr_{r}"].values; n_in = int(col.sum())236 for nom, Xs in [("référence", X), ("sans secteur", Xr)]:237 Xa = Xs.copy(); Xa["corr"] = col; m = h.mco(Xa.values, y); i = list(Xa.columns).index("corr")238 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]})239 return pd.DataFrame(rows)240241242# ----------------------------------------------------------------------------- 6. faux positifs spatiaux243244def faux_positifs(b: dict, v, n_poly: int = N_POLYGONES_ALEATOIRES, graine: int = GRAINE) -> dict:245 """Pour chaque rayon, n_poly disques centrés sur une vente tirée au hasard ; indicatrice « dans le disque » ajoutée au modèle ;246 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 %)247 si la spécification est correcte ; les résidus spatialement corrélés des spécifications sans secteur gonflent cette part."""248 rng = np.random.default_rng(graine); X = _X_ref(v); y = np.log(v["prix_vente"].values.astype(float)); xy = b["xy_m"]249 specs = {"référence": X, "sans secteur": X.drop(columns=[c for c in X.columns if c.startswith("sect_")]),250 "sans secteur ni gradient": X.drop(columns=[c for c in X.columns if c.startswith("sect_")] + ["dist_centre"])}251 rows = []; tstats = {s: {r: [] for r in RAYONS_POLYGONES_KM} for s in specs}252 for r in RAYONS_POLYGONES_KM:253 centres = rng.integers(0, len(v), size=n_poly); n_ok = 0254 for c in centres:255 d = np.sqrt(((xy - xy[c]) ** 2).sum(1)) / 1000; ind = (d <= r).astype(float); n_in = int(ind.sum())256 if n_in < 10 or n_in > len(v) - 10:257 continue258 n_ok += 1259 for s, Xs in specs.items():260 Xa = np.column_stack([Xs.values, ind]); m = h.mco(Xa, y)261 tstats[s][r].append((float(m["beta"][-1] / m["se"][-1]), float(m["beta"][-1] / m["se_hc1"][-1]))) # (t classique, t HC1)262 for s in specs:263 tt = np.array(tstats[s][r]); seuil_bonf = _z(0.05 / len(tt))264 for j, et in enumerate(["classique", "HC1"]):265 t = tt[:, j]266 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)),267 "part_bonferroni": float(np.mean(np.abs(t) > seuil_bonf)), "ecart_type_t": float(t.std()), "t_max_abs": float(np.abs(t).max())})268 return {"tableau": pd.DataFrame(rows), "t": tstats}269270271def _z(alpha_bilateral: float) -> float:272 from scipy.stats import norm273 return float(norm.ppf(1 - alpha_bilateral / 2))274275276# ----------------------------------------------------------------------------- 8. cas 4277278def cas4_geo() -> dict:279 """48 immeubles du cas 4 : projection en NAD83 / Québec Lambert (EPSG 32198), part de la valeur en zone inondable par RMR,280 coefficient « inondable » des effets régionaux, distances entre centroïdes de RMR."""281 from pyproj import Geod282 gpd = _gpd(); d = cas4.charger(); pf = d["portefeuille"]; som = cas4.sommaire(pf); eff = cas4.effets_regionaux(pf)283 g = gpd.GeoDataFrame(pf, geometry=gpd.points_from_xy(pf["lon"], pf["lat"]), crs=EPSG_WGS84).to_crs(EPSG_LAMBERT)284 par = pf.groupby("rmr").agg(n=("id_immeuble", "count"), valeur=("valeur_2026", "sum"), n_inond=("zone_inondable", "sum"), n_sol=("contamination_sol", "sum"),285 val_inond=("valeur_2026", lambda x: x[pf.loc[x.index, "zone_inondable"] == 1].sum()), lat=("lat", "mean"), lon=("lon", "mean"))286 par["part_inond_rmr"] = par["val_inond"] / par["valeur"]; par["part_inond_total"] = par["val_inond"] / pf["valeur_2026"].sum()287 par = par.sort_values("valeur", ascending=False)288 geod = Geod(ellps="WGS84"); cen = par[["lat", "lon"]]289 dist = pd.DataFrame(index=cen.index, columns=cen.index, dtype=float)290 for a in cen.index:291 for bb in cen.index:292 dist.loc[a, bb] = geod.inv(cen.loc[a, "lon"], cen.loc[a, "lat"], cen.loc[bb, "lon"], cen.loc[bb, "lat"])[2] / 1000293 inond = pf[pf["zone_inondable"] == 1][["id_immeuble", "rmr", "nb_logements", "valeur_2026", "taux_global", "etat_1a5"]]294 # taux global moyen des immeubles inondables vs non, à état et RMR comparables : le coefficient du modèle le mesure295 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"]),296 "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),297 "distances_km": dist, "inondables": inond, "n_inond": int(pf["zone_inondable"].sum()), "n_sol": int(pf["contamination_sol"].sum()),298 "etendue_km": (float((g.geometry.x.max() - g.geometry.x.min()) / 1000), float((g.geometry.y.max() - g.geometry.y.min()) / 1000))}299300301# ----------------------------------------------------------------------------- capitalisation du risque (exemple numérique du texte)302303def 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:304 """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)."""305 rows = []306 for p, D in zip(P25, dommage_pct):307 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})308 return pd.DataFrame(rows)309310311# ----------------------------------------------------------------------------- orchestration312313def tout() -> dict:314 if "tout" in _CACHE:315 return _CACHE["tout"]316 warnings.filterwarnings("ignore")317 b = charger_base(); gf = geometries_fictives(b); v = jointure_et_distances(b, gf)318 grille = grille_et_agregation(b, v); idw_ = interpolation_idw(b, grille)319 res = {"base": b, "geom": gf, "ventes": v, "secteurs": tableau_secteurs(v), "grille": grille, "idw": idw_, "modeles": modeles_augmentes(v),320 "tampon": sensibilite_tampon(v), "faux_positifs": faux_positifs(b, v), "cas4": cas4_geo(), "capitalisation": exemple_capitalisation()}321 _CACHE["tout"] = res322 return res323324325def main():326 r = tout(); b = r["base"]; pd.set_option("display.width", 200); pd.set_option("display.max_columns", 30)327 print("=== 1. Projection (NAD83 / MTM zone 9, EPSG 32189) :", b["crs_nom"])328 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")329 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()})330 print("Manhattan / euclidienne :", {k: round(x, 3) for k, x in b["manhattan"].items()})331 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()})332 print("\n=== 3. Jointure spatiale et distances par secteur\n", r["secteurs"].round(3).to_string())333 v = r["ventes"]; print("contrôle lancer de rayon = shapely :", bool((v["zone_A"] == v["zone_A_rayon"]).all()))334 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}")335 for k, c in g["classif"].items():336 print(f" {k:20s} bornes = {np.round(c['bornes'], 0)} ; effectifs = {c['effectifs']} ; GVF = {c['gvf']:.3f}")337 print(" pondération surfacique :", {k: (round(x, 2) if isinstance(x, float) else x) for k, x in g["pond_surf"].items()})338 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()})339 print("\n=== 5. Modèles augmentés")340 for nom, m in r["modeles"].items():341 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}")342 print(m["coef"].round(4).to_string()); print(" EMD (80 %) :", {k: round(x, 4) for k, x in m["mde"].items()})343 print("\n=== 6. Sensibilité au rayon du tampon\n", r["tampon"].round(4).to_string())344 print("\n=== 7. Faux positifs spatiaux\n", r["faux_positifs"]["tableau"].round(3).to_string())345 c4 = r["cas4"]; print("\n=== 8. Cas 4 par RMR\n", c4["par_rmr"].round(4).to_string())346 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} %")347 print("distances géodésiques entre centroïdes de RMR (km)\n", c4["distances_km"].round(0).to_string())348 print("\n=== Capitalisation du risque (exemple)\n", r["capitalisation"].round(4).to_string())349350351# ----------------------------------------------------------------------------- exercices 6 et 7 (solutions/ch38.tex)352353def exercice_6(n_poly: int = 1000, rayon_km: float = 3.0) -> dict:354 """Disques aléatoires de 3 km (1 000 tirages) pour les trois spécifications, puis comparaison classique / HC1 / HC3 à 0,5 km (référence)."""355 r = tout(); b, v = r["base"], r["ventes"]356 global RAYONS_POLYGONES_KM357 sauve = RAYONS_POLYGONES_KM; RAYONS_POLYGONES_KM = [rayon_km]358 try:359 fp = faux_positifs(b, v, n_poly=n_poly)["tableau"]360 finally:361 RAYONS_POLYGONES_KM = sauve362 rng = np.random.default_rng(GRAINE); X = _X_ref(v); y = np.log(v["prix_vente"].values.astype(float)); xy = b["xy_m"]363 ts = {"classique": [], "HC1": [], "HC3": []}364 for cc in rng.integers(0, len(v), size=N_POLYGONES_ALEATOIRES):365 d = np.sqrt(((xy - xy[cc]) ** 2).sum(1)) / 1000; ind = (d <= 0.5).astype(float); n_in = int(ind.sum())366 if n_in < 10 or n_in > len(v) - 10:367 continue368 Xa = np.column_stack([X.values, ind]); m = h.mco(Xa, y); e = m["residus"]369 XtXi = np.linalg.pinv(Xa.T @ Xa); lev = np.einsum("ij,jk,ik->i", Xa, XtXi, Xa)370 V3 = XtXi @ ((Xa * ((e / (1 - lev)) ** 2)[:, None]).T @ Xa) @ XtXi371 ts["classique"].append(m["beta"][-1] / m["se"][-1]); ts["HC1"].append(m["beta"][-1] / m["se_hc1"][-1]); ts["HC3"].append(m["beta"][-1] / np.sqrt(V3[-1, -1]))372 comp = {k: {"n": len(t), "part_5pct": float(np.mean(np.abs(t) > 1.96)), "ecart_type": float(np.std(t))} for k, t in ts.items()}373 return {"tableau_3km": fp, "comparaison_05km": comp}374375376def exercice_7() -> pd.DataFrame:377 """Médiane du prix au m² par cellule pour trois rayons hexagonaux et une grille carrée : cellules, bornes de Jenks, GVF, part intercellulaire de la variance."""378 from shapely.geometry import Polygon, box379 gpd = _gpd(); r = tout(); b, v = r["base"], r["ventes"]; xy = b["xy_m"]380 pts = v[["geometry", "prix_vente", "superficie_habitable_m2"]].copy(); pts["prix_m2"] = pts["prix_vente"] / pts["superficie_habitable_m2"]381 xmin, ymin, xmax, ymax = xy[:, 0].min(), xy[:, 1].min(), xy[:, 0].max(), xy[:, 1].max(); tot = ((pts["prix_m2"] - pts["prix_m2"].mean()) ** 2).sum()382 def evaluer(cells, nom):383 g = gpd.GeoDataFrame({"cellule": range(len(cells))}, geometry=cells, crs=EPSG_MTM9)384 j = gpd.sjoin(pts, g, how="inner", predicate="within"); assert len(j) == len(pts)385 agg = j.groupby("cellule").agg(n=("prix_m2", "count"), med=("prix_m2", "median"), moy=("prix_m2", "mean")); ok = agg[agg["n"] >= MIN_VENTES_CELLULE]386 z = ok["med"].values; bj = geo.bornes_jenks(z, K_CLASSES); inter = float((agg["n"] * (agg["moy"] - pts["prix_m2"].mean()) ** 2).sum())387 return {"grille": nom, "cellules": len(g), "occupees": len(agg), "n_min5": len(ok), "bornes_jenks": np.round(bj, 0), "gvf": geo.gvf(z, bj), "part_inter": inter / tot}388 rows = []389 for R in (800.0, 1200.0, 2000.0):390 rows.append(evaluer([Polygon(cc) for cc in geo.grille_hexagonale(xmin - R, ymin - R, xmax, ymax, R)], f"hexagones R = {R:.0f} m"))391 pas = 2000.0392 rows.append(evaluer([box(x, y_, x + pas, y_ + pas) for x in np.arange(xmin - 0.5 * pas, xmax + pas, pas) for y_ in np.arange(ymin - 0.5 * pas, ymax + pas, pas)], "carrés de 2 km"))393 return pd.DataFrame(rows)394395396if __name__ == "__main__":397 main()398