#!/usr/bin/env python3
"""Figures du chapitre 38 — SIG et risques géospatiaux (cas 1 Gatineau et cas 4 portefeuille, DONNÉES SYNTHÉTIQUES, graine 20261031).
Cartes geopandas/Matplotlib sans fond de tuiles (contextily absent) ; coordonnées projetées NAD83 / MTM zone 9 (cas 1) et
NAD83 / Québec Lambert (cas 4), exprimées en kilomètres ; toutes les géométries « zones », « corridor » et « pôles » sont FICTIVES.
Sorties (assets/figures/ch38/) :
  ch38-fig-carte-cas1.pdf      600 ventes, zones d'exposition fictives, corridor fictif et ses tampons, pôles fictifs, sujet
  ch38-fig-choroplethe.pdf     médiane du prix au m² par cellule hexagonale selon trois classifications + résidus moyens (divergent)
  ch38-fig-faux-positifs.pdf   distribution des t de polygones aléatoires et part de « significatifs » selon le rayon et la spécification
  ch38-fig-tampon.pdf          coefficient de l'indicatrice de corridor selon le rayon du tampon (IC à 95 %)
  ch38-fig-cas4.pdf            48 immeubles du cas 4 par RMR et part de la valeur en zone inondable
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/figures/ch38_figures.py"""
import sys, pathlib, warnings
ROOT = pathlib.Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT / "analysis/python")); sys.path.insert(0, str(ROOT / "analysis/python/cases"))
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
from scipy.stats import norm
from traite.style import setup, palette, savefig, LARGEUR, DEMI, fmt_milliers
from traite import geo
import ch38_sig as c38

SECTEURS = ["Hull", "Gatineau", "Aylmer", "Buckingham", "Masson-Angers"]


def _km(ax, xy, marge=1.5):
    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)
    ax.grid(True, lw=0.4)


def _poly_km(ax, geom, **kw):
    """Trace un (multi)polygone shapely en km (coordonnées MTM divisées par 1000)."""
    polys = getattr(geom, "geoms", [geom])
    for p in polys:
        ext = np.asarray(p.exterior.coords) / 1000
        ax.fill(ext[:, 0], ext[:, 1], **kw)
        for hole in p.interiors:
            hh = np.asarray(hole.coords) / 1000; ax.fill(hh[:, 0], hh[:, 1], color=palette.surface, lw=0)


def fig_carte_cas1(r):
    b, gf, v = r["base"], r["geom"], r["ventes"]; xy = b["xy_m"] / 1000; cs = np.array([b["xs"], b["ys"]]) / 1000
    fig, ax = plt.subplots(figsize=(LARGEUR, 4.3))
    zA, zB = gf["zones"].geometry.iloc[0], gf["zones"].geometry.iloc[1]
    _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)
    riv = np.asarray(gf["cours_eau"].coords) / 1000; ax.plot(riv[:, 0], riv[:, 1], color=palette.seq[8], lw=1.1, zorder=2)
    t300 = gf["tampons"].set_index("rayon_m").geometry.loc[c38.RAYON_CORRIDOR_REF_M]
    _poly_km(ax, t300, color=palette.cat[1], alpha=0.3, lw=0, zorder=1)
    cor = np.asarray(gf["corridor"].coords) / 1000; ax.plot(cor[:, 0], cor[:, 1], color=palette.cat[1], lw=1.4, zorder=2)
    hors = (v["classe_expo"] == "hors zone").values
    ax.scatter(xy[hors, 0], xy[hors, 1], s=6, color=palette.ink2, alpha=0.45, lw=0, zorder=3)
    mA = (v["classe_expo"] == "A").values; mB = (v["classe_expo"] == "B").values
    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)
    pc = np.column_stack([gf["poles"].geometry.x.values, gf["poles"].geometry.y.values]) / 1000
    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)
    for (x, y), nom in zip(pc, gf["poles"]["pole"]):
        ax.annotate(nom, (x, y), xytext=(5, 4), textcoords="offset points", fontsize=7, color=palette.ink2)
    ax.scatter(cs[0], cs[1], marker="*", s=130, color=palette.ink, zorder=6)
    for s in SECTEURS:
        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",
                                                 bbox=dict(boxstyle="round,pad=0.15", fc=palette.surface, ec="none", alpha=0.8))
    _km(ax, xy); ax.set_xlabel("Est, NAD83 / MTM zone 9 (km) — coordonnées synthétiques"); ax.set_ylabel("Nord (km)")
    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)"),
           Line2D([], [], color=palette.cat[1], lw=1.4, label="Corridor de transport fictif et tampon de 300 m"),
           Line2D([], [], marker="^", color="none", markerfacecolor=palette.cat[2], markeredgecolor=palette.ink, markersize=7, label="Pôles fictifs (taille selon la masse)"),
           Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[9], markersize=5, label="Vente en classe A"),
           Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[4], markersize=5, label="Vente en classe B"),
           Line2D([], [], marker="o", color="none", markerfacecolor=palette.ink2, markersize=4, alpha=0.6, label="Autre vente"),
           Line2D([], [], marker="*", color="none", markerfacecolor=palette.ink, markersize=10, label="Sujet")]
    ax.legend(handles=leg, fontsize=6.8, loc="upper left", ncol=2, framealpha=0.9, frameon=True, edgecolor="none", facecolor=palette.surface)
    savefig(fig, "ch38/ch38-fig-carte-cas1.pdf")


def _carte_classes(ax, cells, z, bornes, cmap_colors, titre, fmt=lambda x: fmt_milliers(x) + " $"):
    cmap = ListedColormap(cmap_colors); normc = BoundaryNorm(bornes, cmap.N)
    for geom, val in zip(cells.geometry, z):
        ext = np.asarray(geom.exterior.coords) / 1000
        ax.fill(ext[:, 0], ext[:, 1], color=cmap(normc(val)), lw=0.3, ec=palette.surface)
    ax.set_title(titre, fontsize=8.5)
    handles = [Patch(color=cmap(i), label=f"{fmt(bornes[i])} – {fmt(bornes[i + 1])}") for i in range(len(bornes) - 1)]
    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)


def fig_choroplethe(r):
    g = r["grille"]; ok = g["cartographiees"]; occ = g["occupees"]; b = r["base"]; xy = b["xy_m"] / 1000
    z = ok["prix_m2_med"].values; k = c38.K_CLASSES
    seq = [palette.seq[i] for i in np.linspace(0, 9, k).round().astype(int)]
    fig, axes = plt.subplots(2, 2, figsize=(LARGEUR, 4.0), sharex=True, sharey=True)
    titres = {"intervalles égaux": "(a) Intervalles égaux", "quantiles": "(b) Quantiles", "Jenks": "(c) Fisher–Jenks"}
    for ax, (nom, cl) in zip(axes.flat[:3], g["classif"].items()):
        for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry:
            ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3)
        _carte_classes(ax, ok, z, cl["bornes"], seq, f"{titres[nom]}, GVF = {cl['gvf']:.3f}".replace(".", ","))
        ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0)
    ax = axes.flat[3]; res = ok["residu_moy"].values; vmax = float(np.abs(res).max()); bornes = np.linspace(-vmax, vmax, 7)
    div = [palette.div["neg"], "#7fa8e0", "#c5d4ea", "#ead0cf", "#eb9a99", palette.div["pos"]]
    for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry:
        ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3)
    _carte_classes(ax, ok, res, bornes, div, "(d) Résidu moyen, modèle de référence", fmt=lambda x: f"{x:+.3f}".replace(".", ","))
    ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0)
    for ax in axes.flat:
        _km(ax, xy, marge=1.0)
    for ax in axes[1]:
        ax.set_xlabel("Est, MTM 9 (km)")
    for ax in axes[:, 0]:
        ax.set_ylabel("Nord (km)")
    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)
    fig.subplots_adjust(hspace=0.28, wspace=0.06, bottom=0.12, top=0.95)
    savefig(fig, "ch38/ch38-fig-choroplethe.pdf")


def fig_faux_positifs(r):
    fp = r["faux_positifs"]; tab = fp["tableau"]; T = fp["t"]
    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.1))
    ax = axes[0]; bins = np.linspace(-8, 8, 49); x = np.linspace(-8, 8, 400)
    for i, (spec, col) in enumerate([("référence", palette.cat[0]), ("sans secteur", palette.cat[1])]):
        t = np.array(T[spec][1.0])[:, 0]
        ax.hist(t, bins=bins, density=True, histtype="stepfilled", alpha=0.45, color=col, lw=0, label=f"{spec} (écart-type {t.std():.2f})")
    ax.plot(x, norm.pdf(x), color=palette.ink, lw=1.0, ls="--", label="Loi normale centrée réduite")
    for s in (-1.96, 1.96):
        ax.axvline(s, color=palette.ink2, lw=0.7, ls=":")
    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)
    ax.legend(fontsize=6.8, loc="upper right")
    ax = axes[1]
    for i, spec in enumerate(["référence", "sans secteur", "sans secteur ni gradient"]):
        w = tab[(tab["spec"] == spec) & (tab["erreur_type"] == "classique")].sort_values("rayon_km")
        ax.plot(w["rayon_km"], w["part_5pct"] * 100, marker="o", color=palette.cat[i], label=spec)
    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")
    ax.set_xlabel("Rayon du disque aléatoire (km)"); ax.set_ylabel("Part des coefficients « significatifs » à 5 % (%)"); ax.set_ylim(0, 85)
    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")
    fig.tight_layout(); savefig(fig, "ch38/ch38-fig-faux-positifs.pdf")


def fig_tampon(r):
    tab = r["tampon"]
    fig, ax = plt.subplots(figsize=(LARGEUR, 3.0))
    for i, (spec, dx) in enumerate([("référence", -0.012), ("sans secteur", 0.012)]):
        w = tab[tab["spec"] == spec].sort_values("rayon_m"); x = np.log10(w["rayon_m"]) + dx
        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)
    ax.axhline(0, color=palette.ink, lw=0.8)
    w = tab[tab["spec"] == "référence"].sort_values("rayon_m")
    for xr, n in zip(w["rayon_m"], w["n_dans_tampon"]):
        ax.text(np.log10(xr), 13.5, f"n = {n}", fontsize=6.5, ha="center", color=palette.ink2)
    ax.set_xticks(np.log10(w["rayon_m"])); ax.set_xticklabels([str(int(v)) for v in w["rayon_m"]])
    ax.set_xlabel("Rayon du tampon autour du corridor fictif (m, échelle logarithmique)"); ax.set_ylabel("Coefficient de l'indicatrice (%) et IC à 95 %")
    ax.set_ylim(-20, 16); ax.legend(fontsize=7.5, loc="lower right")
    savefig(fig, "ch38/ch38-fig-tampon.pdf")


def fig_cas4(r):
    c4 = r["cas4"]; g = c4["gdf"]; par = c4["par_rmr"]
    fig, axes = plt.subplots(1, 2, figsize=(LARGEUR, 3.6), gridspec_kw={"width_ratios": [1.55, 1]})
    ax = axes[0]; x = g.geometry.x.values / 1000; y = g.geometry.y.values / 1000; val = g["valeur_2026"].values
    s = 12 + 110 * val / val.max()
    ax.scatter(x, y, s=s, color=palette.cat[0], alpha=0.55, lw=0.4, edgecolor=palette.ink, zorder=3)
    m = g["zone_inondable"].values == 1
    ax.scatter(x[m], y[m], s=s[m] + 60, facecolor="none", edgecolor=palette.cat[5], lw=1.6, zorder=4)
    m2 = g["contamination_sol"].values == 1
    ax.scatter(x[m2], y[m2], marker="s", s=s[m2] + 70, facecolor="none", edgecolor=palette.cat[3], lw=1.3, zorder=4)
    for rmr, w in g.groupby("rmr"):
        cx, cy = w.geometry.x.mean() / 1000, w.geometry.y.mean() / 1000
        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"))
        ax.annotate(f"{rmr} ({len(w)})", (cx, cy), xytext=(dx, dy), textcoords="offset points", fontsize=7, ha=ha, color=palette.ink2)
    ax.set_aspect("equal"); ax.grid(True, lw=0.4); ax.set_xlabel("Est, NAD83 / Québec Lambert (km)"); ax.set_ylabel("Nord (km)")
    ax.set_title("(a) 48 immeubles, taille selon la valeur 2026", fontsize=8.5)
    leg = [Line2D([], [], marker="o", color="none", markerfacecolor=palette.cat[0], alpha=0.6, markersize=7, label="Immeuble"),
           Line2D([], [], marker="o", color="none", markerfacecolor="none", markeredgecolor=palette.cat[5], markeredgewidth=1.6, markersize=9, label="Zone inondable (3)"),
           Line2D([], [], marker="s", color="none", markerfacecolor="none", markeredgecolor=palette.cat[3], markeredgewidth=1.3, markersize=9, label="Contamination du sol (2)")]
    ax.legend(handles=leg, fontsize=6.8, loc="upper left")
    ax.set_ylim(y.min() - 55, y.max() + 60); ax.set_xlim(x.min() - 60, x.max() + 90)
    ax = axes[1]; ordre = par.index[::-1]; yy = np.arange(len(ordre))
    ax.barh(yy, par.loc[ordre, "part_inond_rmr"] * 100, color=palette.cat[5], alpha=0.85, height=0.6)
    for yi, rmr in zip(yy, ordre):
        p = par.loc[rmr, "part_inond_rmr"] * 100; ax.text(p + 0.8, yi, f"{p:.1f} %".replace(".", ","), va="center", fontsize=7)
    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 (%)")
    ax.set_title(f"(b) Ensemble : {c4['total']['part_zone_inondable'] * 100:.1f} % de la valeur".replace(".", ","), fontsize=8.5)
    fig.tight_layout(); savefig(fig, "ch38/ch38-fig-cas4.pdf")


def main():
    warnings.filterwarnings("ignore"); setup(); r = c38.tout()
    fig_carte_cas1(r); fig_choroplethe(r); fig_faux_positifs(r); fig_tampon(r); fig_cas4(r)
    print("figures ch38 écrites dans", ROOT / "assets/figures/ch38")


if __name__ == "__main__":
    main()