spb/traite-evaluation-qc

1"""Outils géospatiaux génériques du traité (chapitre 38 ; réutilisables aux chapitres 43 et 48).2Fonctions pures NumPy, sans dépendance à geopandas : distances (haversine, équirectangulaire locale), facteur d'échelle3d'une projection transverse de Mercator, point dans polygone (lancer de rayon), interpolation par l'inverse de la4distance, grille hexagonale, classification des cartes choroplèthes (intervalles égaux, quantiles, Fisher–Jenks, GVF),5accessibilité gravitaire, conversion d'une probabilité sur 25 ans en probabilité annuelle, effet minimal détectable.6Tests : analysis/python/tests/test_ch38.py. Ajouté le 2026-10-31 ; n'altère aucune fonction existante."""7from __future__ import annotations8import numpy as np910R_TERRE_KM = 6371.0088          # rayon moyen de la Terre (IUGG, rayon volumétrique), km11MTM_K0 = 0.9999                 # facteur d'échelle au méridien central des zones MTM du Québec (EPSG 32187-32189)12UTM_K0 = 0.9996131415def haversine_km(lat1, lon1, lat2, lon2, rayon: float = R_TERRE_KM) -> np.ndarray:16    """Distance orthodromique (formule du haversine) sur une sphère de rayon `rayon`, en km. Entrées en degrés."""17    p1, p2 = np.radians(np.asarray(lat1, float)), np.radians(np.asarray(lat2, float))18    dphi = p2 - p1; dlam = np.radians(np.asarray(lon2, float) - np.asarray(lon1, float))19    a = np.sin(dphi / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dlam / 2) ** 220    return 2 * rayon * np.arcsin(np.sqrt(np.clip(a, 0.0, 1.0)))212223def equirectangulaire_km(lat, lon, lat0: float, lon0: float) -> np.ndarray:24    """Projection équirectangulaire locale autour de (lat0, lon0) : x = R cos(lat0) Δλ, y = R Δφ (km). Même rayon que25    haversine_km, de sorte que les deux mesures coïncident au premier ordre autour du point de référence."""26    lat = np.radians(np.asarray(lat, float)); lon = np.radians(np.asarray(lon, float))27    x = R_TERRE_KM * np.cos(np.radians(lat0)) * (lon - np.radians(lon0)); y = R_TERRE_KM * (lat - np.radians(lat0))28    return np.column_stack([x, y])293031def facteur_echelle_tm(lat, lon, lon0: float, k0: float = MTM_K0) -> np.ndarray:32    """Facteur d'échelle ponctuel d'une projection transverse de Mercator (approximation sphérique au second ordre) :33    k ≈ k0 [1 + ½ (Δλ cos φ)²], Δλ en radians. Une distance au sol de 1 km est représentée par k km sur la carte."""34    dl = np.radians(np.asarray(lon, float) - lon0); phi = np.radians(np.asarray(lat, float))35    return k0 * (1.0 + 0.5 * (dl * np.cos(phi)) ** 2)363738def point_dans_polygone(points: np.ndarray, polygone: np.ndarray) -> np.ndarray:39    """Test point-dans-polygone par lancer de rayon (parité du nombre de croisements d'une demi-droite horizontale).40    points : (m, 2) ; polygone : (n, 2) sommets ordonnés (fermeture implicite). Retourne un tableau booléen (m,)."""41    P = np.asarray(points, float); V = np.asarray(polygone, float)42    x, y = P[:, 0], P[:, 1]; n = len(V); dedans = np.zeros(len(P), bool)43    for i in range(n):44        x1, y1 = V[i]; x2, y2 = V[(i + 1) % n]45        croise = (y1 > y) != (y2 > y)46        with np.errstate(divide="ignore", invalid="ignore"):47            xi = x1 + (y - y1) * (x2 - x1) / (y2 - y1)48        dedans ^= croise & (x < xi)49    return dedans505152def idw(coords_obs: np.ndarray, z: np.ndarray, coords_cibles: np.ndarray, puissance: float = 2.0, voisins: int | None = None) -> np.ndarray:53    """Interpolation par l'inverse de la distance : ẑ(s0) = Σ w_i z_i / Σ w_i, w_i = d_i^(−p) (valeur exacte si d_i = 0)."""54    co = np.asarray(coords_obs, float); z = np.asarray(z, float); ct = np.asarray(coords_cibles, float)55    d = np.sqrt(((ct[:, None, :] - co[None, :, :]) ** 2).sum(-1)); out = np.zeros(len(ct))56    for i in range(len(ct)):57        idx = np.argsort(d[i])[:voisins] if voisins else np.arange(len(co))58        di = d[i, idx]59        if di.min() == 0:60            out[i] = z[idx][di.argmin()]; continue61        w = di ** (-puissance); out[i] = (w @ z[idx]) / w.sum()62    return out636465def grille_hexagonale(xmin: float, ymin: float, xmax: float, ymax: float, rayon: float) -> list[np.ndarray]:66    """Pavage hexagonal (hexagones « pointe en haut » de rayon circonscrit `rayon`) couvrant la boîte englobante.67    Retourne la liste des sommets (6, 2) de chaque cellule ; aire d'une cellule = (3√3/2) rayon²."""68    w = np.sqrt(3) * rayon; hstep = 1.5 * rayon; cells = []69    ang = np.radians(np.arange(30, 390, 60))70    j = 0; y = ymin71    while y <= ymax + rayon:72        off = (w / 2) if (j % 2) else 0.0; x = xmin + off73        while x <= xmax + w:74            cells.append(np.column_stack([x + rayon * np.cos(ang), y + rayon * np.sin(ang)])); x += w75        y += hstep; j += 176    return cells777879def bornes_intervalles_egaux(z, k: int) -> np.ndarray:80    z = np.asarray(z, float); return np.linspace(z.min(), z.max(), k + 1)818283def bornes_quantiles(z, k: int) -> np.ndarray:84    z = np.asarray(z, float); return np.quantile(z, np.linspace(0, 1, k + 1))858687def bornes_jenks(z, k: int) -> np.ndarray:88    """Classification de Fisher–Jenks (programmation dynamique exacte) : minimise la somme intra-classe des carrés des écarts.89    Retourne k + 1 bornes (min, ruptures, max)."""90    z = np.sort(np.asarray(z, float)); n = len(z)91    if k >= n:92        return np.concatenate([[z[0]], z[1:], [z[-1]]])93    cs = np.concatenate([[0.0], np.cumsum(z)]); cs2 = np.concatenate([[0.0], np.cumsum(z ** 2)])94    def ssd(i, j):  # somme des carrés des écarts de z[i:j]95        m = j - i; s = cs[j] - cs[i]; return (cs2[j] - cs2[i]) - s * s / m96    INF = np.inf97    cout = np.full((k + 1, n + 1), INF); arg = np.zeros((k + 1, n + 1), int); cout[0, 0] = 0.098    for c in range(1, k + 1):99        for j in range(c, n + 1):100            best, bi = INF, -1101            for i in range(c - 1, j):102                v = cout[c - 1, i] + ssd(i, j)103                if v < best:104                    best, bi = v, i105            cout[c, j] = best; arg[c, j] = bi106    bornes = [z[-1]]; j = n107    for c in range(k, 0, -1):108        i = arg[c, j]; bornes.append(z[i]); j = i109    return np.array(bornes[::-1])110111112def gvf(z, bornes) -> float:113    """Goodness of variance fit = 1 − SDCM/SDAM (somme des carrés intra-classe / somme totale des carrés)."""114    z = np.asarray(z, float); b = np.asarray(bornes, float)115    cl = np.clip(np.searchsorted(b, z, side="right") - 1, 0, len(b) - 2)116    sdcm = sum(((z[cl == c] - z[cl == c].mean()) ** 2).sum() for c in range(len(b) - 1) if (cl == c).any())117    return 1.0 - sdcm / ((z - z.mean()) ** 2).sum()118119120def classer(z, bornes) -> np.ndarray:121    """Indice de classe (0 … k−1) de chaque valeur selon des bornes croissantes ; la dernière borne est incluse."""122    z = np.asarray(z, float); b = np.asarray(bornes, float)123    return np.clip(np.searchsorted(b, z, side="right") - 1, 0, len(b) - 2)124125126def accessibilite_gravitaire(coords: np.ndarray, coords_poles: np.ndarray, masses, friction: float, forme: str = "exp") -> np.ndarray:127    """Indice de Hansen A_i = Σ_k E_k f(d_ik), f(d) = exp(−λ d) (forme « exp ») ou d^(−λ) (forme « puissance », d en km, d ≥ 0,1)."""128    c = np.asarray(coords, float); p = np.asarray(coords_poles, float); E = np.asarray(masses, float)129    d = np.sqrt(((c[:, None, :] - p[None, :, :]) ** 2).sum(-1))130    f = np.exp(-friction * d) if forme == "exp" else np.maximum(d, 0.1) ** (-friction)131    return f @ E132133134def probabilite_annuelle(p_horizon: float, horizon: int = 25) -> float:135    """Probabilité annuelle d'au moins une occurrence équivalente à une probabilité p sur `horizon` années (événements indépendants) :136    p_1 = 1 − (1 − p)^(1/horizon)."""137    return 1.0 - (1.0 - p_horizon) ** (1.0 / horizon)138139140def effet_minimal_detectable(se: float, alpha: float = 0.05, puissance: float = 0.80) -> float:141    """Effet minimal détectable (bilatéral) : (z_{1−α/2} + z_{puissance}) · se."""142    from scipy.stats import norm143    return float((norm.ppf(1 - alpha / 2) + norm.ppf(puissance)) * se)144