spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Figures du chapitre 38 — SIG et risques géospatiaux (cas 1 Gatineau et cas 4 portefeuille, DONNÉES SYNTHÉTIQUES, graine 20261031).3Cartes geopandas/Matplotlib sans fond de tuiles (contextily absent) ; coordonnées projetées NAD83 / MTM zone 9 (cas 1) et4NAD83 / Québec Lambert (cas 4), exprimées en kilomètres ; toutes les géométries « zones », « corridor » et « pôles » sont FICTIVES.5Sorties (assets/figures/ch38/) :6  ch38-fig-carte-cas1.pdf      600 ventes, zones d'exposition fictives, corridor fictif et ses tampons, pôles fictifs, sujet7  ch38-fig-choroplethe.pdf     médiane du prix au m² par cellule hexagonale selon trois classifications + résidus moyens (divergent)8  ch38-fig-faux-positifs.pdf   distribution des t de polygones aléatoires et part de « significatifs » selon le rayon et la spécification9  ch38-fig-tampon.pdf          coefficient de l'indicatrice de corridor selon le rayon du tampon (IC à 95 %)10  ch38-fig-cas4.pdf            48 immeubles du cas 4 par RMR et part de la valeur en zone inondable11Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch38_figures.py"""12import sys, pathlib, warnings13ROOT = pathlib.Path(__file__).resolve().parents[3]14sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))15import numpy as np16import matplotlib.pyplot as plt17from matplotlib.colors import ListedColormap, BoundaryNorm18from matplotlib.lines import Line2D19from matplotlib.patches import Patch20from scipy.stats import norm21from traite.style import setup, palette, savefig, LARGEUR, DEMI, fmt_milliers22from traite import geo23import ch38_sig as c382425SECTEURS = ["Hull", "Gatineau", "Aylmer", "Buckingham", "Masson-Angers"]262728def _km(ax, xy, marge=1.5):29    ax.set_aspect("equal"); ax.set_xlim(xy[:, 0].min() - marge, xy[:, 0].max() + marge); ax.set_ylim(xy[:, 1].min() - marge, xy[:, 1].max() + marge)30    ax.grid(True, lw=0.4)313233def _poly_km(ax, geom, **kw):34    """Trace un (multi)polygone shapely en km (coordonnées MTM divisées par 1000)."""35    polys = getattr(geom, "geoms", [geom])36    for p in polys:37        ext = np.asarray(p.exterior.coords) / 100038        ax.fill(ext[:, 0], ext[:, 1], **kw)39        for hole in p.interiors:40            hh = np.asarray(hole.coords) / 1000; ax.fill(hh[:, 0], hh[:, 1], color=palette.surface, lw=0)414243def fig_carte_cas1(r):44    b, gf, v = r["base"], r["geom"], r["ventes"]; xy = b["xy_m"] / 1000; cs = np.array([b["xs"], b["ys"]]) / 100045    fig, ax = plt.subplots(figsize=(LARGEUR, 4.3))46    zA, zB = gf["zones"].geometry.iloc[0], gf["zones"].geometry.iloc[1]47    _poly_km(ax, zB, color=palette.seq[1], alpha=0.55, lw=0, zorder=1); _poly_km(ax, zA, color=palette.seq[5], alpha=0.6, lw=0, zorder=1)48    riv = np.asarray(gf["cours_eau"].coords) / 1000; ax.plot(riv[:, 0], riv[:, 1], color=palette.seq[8], lw=1.1, zorder=2)49    t300 = gf["tampons"].set_index("rayon_m").geometry.loc[c38.RAYON_CORRIDOR_REF_M]50    _poly_km(ax, t300, color=palette.cat[1], alpha=0.3, lw=0, zorder=1)51    cor = np.asarray(gf["corridor"].coords) / 1000; ax.plot(cor[:, 0], cor[:, 1], color=palette.cat[1], lw=1.4, zorder=2)52    hors = (v["classe_expo"] == "hors zone").values53    ax.scatter(xy[hors, 0], xy[hors, 1], s=6, color=palette.ink2, alpha=0.45, lw=0, zorder=3)54    mA = (v["classe_expo"] == "A").values; mB = (v["classe_expo"] == "B").values55    ax.scatter(xy[mB, 0], xy[mB, 1], s=14, color=palette.seq[4], lw=0, zorder=4); ax.scatter(xy[mA, 0], xy[mA, 1], s=16, color=palette.seq[9], lw=0, zorder=4)56    pc = np.column_stack([gf["poles"].geometry.x.values, gf["poles"].geometry.y.values]) / 100057    ax.scatter(pc[:, 0], pc[:, 1], marker="^", s=18 + gf["poles"]["masse"].values * 1.2, color=palette.cat[2], edgecolor=palette.ink, lw=0.6, zorder=5)58    for (x, y), nom in zip(pc, gf["poles"]["pole"]):59        ax.annotate(nom, (x, y), xytext=(5, 4), textcoords="offset points", fontsize=7, color=palette.ink2)60    ax.scatter(cs[0], cs[1], marker="*", s=130, color=palette.ink, zorder=6)61    for s in SECTEURS:62        m = (v["secteur"] == s).values; ax.text(np.median(xy[m, 0]), np.median(xy[m, 1]) + 3.2, s, fontsize=7, color=palette.ink2, ha="center",63                                                 bbox=dict(boxstyle="round,pad=0.15", fc=palette.surface, ec="none", alpha=0.8))64    _km(ax, xy); ax.set_xlabel("Est, NAD83 / MTM zone 9 (km) — coordonnées synthétiques"); ax.set_ylabel("Nord (km)")65    leg = [Patch(color=palette.seq[5], alpha=0.6, label="Zone d'exposition fictive, classe A (≤ 350 m)"), Patch(color=palette.seq[1], alpha=0.55, label="Classe B (350–900 m)"),66           Line2D([], [], color=palette.cat[1], lw=1.4, label="Corridor de transport fictif et tampon de 300 m"),67           Line2D([], [], marker="^", color="none", markerfacecolor=palette.cat[2], markeredgecolor=palette.ink, markersize=7, label="Pôles fictifs (taille selon la masse)"),68           Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[9], markersize=5, label="Vente en classe A"),69           Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[4], markersize=5, label="Vente en classe B"),70           Line2D([], [], marker="o", color="none", markerfacecolor=palette.ink2, markersize=4, alpha=0.6, label="Autre vente"),71           Line2D([], [], marker="*", color="none", markerfacecolor=palette.ink, markersize=10, label="Sujet")]72    ax.legend(handles=leg, fontsize=6.8, loc="upper left", ncol=2, framealpha=0.9, frameon=True, edgecolor="none", facecolor=palette.surface)73    savefig(fig, "ch38/ch38-fig-carte-cas1.pdf")747576def _carte_classes(ax, cells, z, bornes, cmap_colors, titre, fmt=lambda x: fmt_milliers(x) + " $"):77    cmap = ListedColormap(cmap_colors); normc = BoundaryNorm(bornes, cmap.N)78    for geom, val in zip(cells.geometry, z):79        ext = np.asarray(geom.exterior.coords) / 100080        ax.fill(ext[:, 0], ext[:, 1], color=cmap(normc(val)), lw=0.3, ec=palette.surface)81    ax.set_title(titre, fontsize=8.5)82    handles = [Patch(color=cmap(i), label=f"{fmt(bornes[i])} – {fmt(bornes[i + 1])}") for i in range(len(bornes) - 1)]83    ax.legend(handles=handles, fontsize=5.8, loc="upper left", frameon=True, framealpha=0.9, edgecolor="none", facecolor=palette.surface, handlelength=1.2, borderpad=0.3, labelspacing=0.25)848586def fig_choroplethe(r):87    g = r["grille"]; ok = g["cartographiees"]; occ = g["occupees"]; b = r["base"]; xy = b["xy_m"] / 100088    z = ok["prix_m2_med"].values; k = c38.K_CLASSES89    seq = [palette.seq[i] for i in np.linspace(0, 9, k).round().astype(int)]90    fig, axes = plt.subplots(2, 2, figsize=(LARGEUR, 4.0), sharex=True, sharey=True)91    titres = {"intervalles égaux": "(a) Intervalles égaux", "quantiles": "(b) Quantiles", "Jenks": "(c) Fisher–Jenks"}92    for ax, (nom, cl) in zip(axes.flat[:3], g["classif"].items()):93        for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry:94            ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3)95        _carte_classes(ax, ok, z, cl["bornes"], seq, f"{titres[nom]}, GVF = {cl['gvf']:.3f}".replace(".", ","))96        ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0)97    ax = axes.flat[3]; res = ok["residu_moy"].values; vmax = float(np.abs(res).max()); bornes = np.linspace(-vmax, vmax, 7)98    div = [palette.div["neg"], "#7fa8e0", "#c5d4ea", "#ead0cf", "#eb9a99", palette.div["pos"]]99    for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry:100        ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3)101    _carte_classes(ax, ok, res, bornes, div, "(d) Résidu moyen, modèle de référence", fmt=lambda x: f"{x:+.3f}".replace(".", ","))102    ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0)103    for ax in axes.flat:104        _km(ax, xy, marge=1.0)105    for ax in axes[1]:106        ax.set_xlabel("Est, MTM 9 (km)")107    for ax in axes[:, 0]:108        ax.set_ylabel("Nord (km)")109    fig.text(0.5, 0.0, f"Cellules hexagonales de {g['aire_cellule_km2']:.2f} km² ; hachures : cellules de moins de {c38.MIN_VENTES_CELLULE} ventes (non cartographiées) ; points : ventes", ha="center", fontsize=7, color=palette.ink2)110    fig.subplots_adjust(hspace=0.28, wspace=0.06, bottom=0.12, top=0.95)111    savefig(fig, "ch38/ch38-fig-choroplethe.pdf")112113114def fig_faux_positifs(r):115    fp = r["faux_positifs"]; tab = fp["tableau"]; T = fp["t"]116    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.1))117    ax = axes[0]; bins = np.linspace(-8, 8, 49); x = np.linspace(-8, 8, 400)118    for i, (spec, col) in enumerate([("référence", palette.cat[0]), ("sans secteur", palette.cat[1])]):119        t = np.array(T[spec][1.0])[:, 0]120        ax.hist(t, bins=bins, density=True, histtype="stepfilled", alpha=0.45, color=col, lw=0, label=f"{spec} (écart-type {t.std():.2f})")121    ax.plot(x, norm.pdf(x), color=palette.ink, lw=1.0, ls="--", label="Loi normale centrée réduite")122    for s in (-1.96, 1.96):123        ax.axvline(s, color=palette.ink2, lw=0.7, ls=":")124    ax.set_xlabel("Statistique $t$ de l'indicatrice d'un disque aléatoire de 1 km"); ax.set_ylabel("Densité"); ax.set_title("(a) Distribution des $t$ (346 disques)", fontsize=8.5)125    ax.legend(fontsize=6.8, loc="upper right")126    ax = axes[1]127    for i, spec in enumerate(["référence", "sans secteur", "sans secteur ni gradient"]):128        w = tab[(tab["spec"] == spec) & (tab["erreur_type"] == "classique")].sort_values("rayon_km")129        ax.plot(w["rayon_km"], w["part_5pct"] * 100, marker="o", color=palette.cat[i], label=spec)130    ax.axhline(5, color=palette.ink2, lw=0.8, ls="--"); ax.text(2.0, 6.0, "seuil nominal 5 %", fontsize=7, color=palette.ink2, ha="right")131    ax.set_xlabel("Rayon du disque aléatoire (km)"); ax.set_ylabel("Part des coefficients « significatifs » à 5 % (%)"); ax.set_ylim(0, 85)132    ax.set_xticks([0.5, 1.0, 2.0]); ax.set_title("(b) Faux positifs selon la spécification", fontsize=8.5); ax.legend(fontsize=6.8, loc="upper left")133    fig.tight_layout(); savefig(fig, "ch38/ch38-fig-faux-positifs.pdf")134135136def fig_tampon(r):137    tab = r["tampon"]138    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))139    for i, (spec, dx) in enumerate([("référence", -0.012), ("sans secteur", 0.012)]):140        w = tab[tab["spec"] == spec].sort_values("rayon_m"); x = np.log10(w["rayon_m"]) + dx141        ax.errorbar(x, w["coef"] * 100, yerr=1.96 * w["se"] * 100, fmt="o-", color=palette.cat[i], capsize=2, lw=1.2, ms=4, label=spec)142    ax.axhline(0, color=palette.ink, lw=0.8)143    w = tab[tab["spec"] == "référence"].sort_values("rayon_m")144    for xr, n in zip(w["rayon_m"], w["n_dans_tampon"]):145        ax.text(np.log10(xr), 13.5, f"n = {n}", fontsize=6.5, ha="center", color=palette.ink2)146    ax.set_xticks(np.log10(w["rayon_m"])); ax.set_xticklabels([str(int(v)) for v in w["rayon_m"]])147    ax.set_xlabel("Rayon du tampon autour du corridor fictif (m, échelle logarithmique)"); ax.set_ylabel("Coefficient de l'indicatrice (%) et IC à 95 %")148    ax.set_ylim(-20, 16); ax.legend(fontsize=7.5, loc="lower right")149    savefig(fig, "ch38/ch38-fig-tampon.pdf")150151152def fig_cas4(r):153    c4 = r["cas4"]; g = c4["gdf"]; par = c4["par_rmr"]154    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.6), gridspec_kw={"width_ratios": [1.55, 1]})155    ax = axes[0]; x = g.geometry.x.values / 1000; y = g.geometry.y.values / 1000; val = g["valeur_2026"].values156    s = 12 + 110 * val / val.max()157    ax.scatter(x, y, s=s, color=palette.cat[0], alpha=0.55, lw=0.4, edgecolor=palette.ink, zorder=3)158    m = g["zone_inondable"].values == 1159    ax.scatter(x[m], y[m], s=s[m] + 60, facecolor="none", edgecolor=palette.cat[5], lw=1.6, zorder=4)160    m2 = g["contamination_sol"].values == 1161    ax.scatter(x[m2], y[m2], marker="s", s=s[m2] + 70, facecolor="none", edgecolor=palette.cat[3], lw=1.3, zorder=4)162    for rmr, w in g.groupby("rmr"):163        cx, cy = w.geometry.x.mean() / 1000, w.geometry.y.mean() / 1000164        dx, dy, ha = (0, 14, "center") if rmr not in ("Montréal", "Sherbrooke") else ((0, -20, "center") if rmr == "Montréal" else (14, -4, "left"))165        ax.annotate(f"{rmr} ({len(w)})", (cx, cy), xytext=(dx, dy), textcoords="offset points", fontsize=7, ha=ha, color=palette.ink2)166    ax.set_aspect("equal"); ax.grid(True, lw=0.4); ax.set_xlabel("Est, NAD83 / Québec Lambert (km)"); ax.set_ylabel("Nord (km)")167    ax.set_title("(a) 48 immeubles, taille selon la valeur 2026", fontsize=8.5)168    leg = [Line2D([], [], marker="o", color="none", markerfacecolor=palette.cat[0], alpha=0.6, markersize=7, label="Immeuble"),169           Line2D([], [], marker="o", color="none", markerfacecolor="none", markeredgecolor=palette.cat[5], markeredgewidth=1.6, markersize=9, label="Zone inondable (3)"),170           Line2D([], [], marker="s", color="none", markerfacecolor="none", markeredgecolor=palette.cat[3], markeredgewidth=1.3, markersize=9, label="Contamination du sol (2)")]171    ax.legend(handles=leg, fontsize=6.8, loc="upper left")172    ax.set_ylim(y.min() - 55, y.max() + 60); ax.set_xlim(x.min() - 60, x.max() + 90)173    ax = axes[1]; ordre = par.index[::-1]; yy = np.arange(len(ordre))174    ax.barh(yy, par.loc[ordre, "part_inond_rmr"] * 100, color=palette.cat[5], alpha=0.85, height=0.6)175    for yi, rmr in zip(yy, ordre):176        p = par.loc[rmr, "part_inond_rmr"] * 100; ax.text(p + 0.8, yi, f"{p:.1f} %".replace(".", ","), va="center", fontsize=7)177    ax.set_yticks(yy); ax.set_yticklabels(ordre, fontsize=7.5); ax.set_xlim(0, 36); ax.set_xlabel("Part de la valeur de la RMR en zone inondable (%)")178    ax.set_title(f"(b) Ensemble : {c4['total']['part_zone_inondable'] * 100:.1f} % de la valeur".replace(".", ","), fontsize=8.5)179    fig.tight_layout(); savefig(fig, "ch38/ch38-fig-cas4.pdf")180181182def main():183    warnings.filterwarnings("ignore"); setup(); r = c38.tout()184    fig_carte_cas1(r); fig_choroplethe(r); fig_faux_positifs(r); fig_tampon(r); fig_cas4(r)185    print("figures ch38 écrites dans", ROOT / "assets/figures/ch38")186187188if __name__ == "__main__":189    main()190