Aislo/E07_StabArea/E07_StabArea_Engine.py

"""E07 광역 자연사면 위험 — 무한사면 안전율 · SINMAP 안정 지수 (요구 345 · 351 · PLAN 536-2).

순수 계산(storage 모름) — D 예정노선 재해 안전성도 같은 엔진을 씀(한 벌).
근거 = 지식DB `original/예정노선/02` §1.1 · §2 ① · `original/안정성평가/04`(붕적토 범위).

  지형   DEM → 고른 격자(cell_m) · 경사 β(Horn 3×3) · D8 흐름(웅덩이 채움 · 수신 셀 = B04 배수유역
         함수 재사용) · 단위 집수면적 a = 상류 셀 수 × cell²/cell(등고선 폭 = cell)
  포화비 m  시나리오 둘
         sinmap(기본 · SINMAP 표준) = 습윤지수 W = min(x·a / sinβ, 1) · x = R/T ∈ [1/(T/R)_hi, 1/(T/R)_lo]
         saturated(최악 · 강우 극한) = 지표 포화 m = 1
  안전율 Fs = (c + c_r + [γt z(1−m) + (γsat−γw) m z] cos²β tanφ)
            / ([γt z(1−m) + γsat m z] sinβ cosβ)
  안정 지수 SI(SINMAP · Pack 1998) — c ∈ [c, c_hi] · tanφ ∈ [tanφ, tanφ_hi] · x ∈ [x_lo, x_hi] 고른 분포
         하한 Fs_min = (c, tanφ, x_hi 가장 젖음) ≥ 1 → SI = Fs_min · 상한 Fs_max = (c_hi, tanφ_hi, x_lo) < 1
         → SI = 0 · 그 사이 → SI = P(Fs > 1)(c · tanφ 닫힌 꼴 · x 는 중점 적분 X_STEPS 마디)
  등급  1 안정(SI ≥ 1.5) · 2 보통 안정(≥ 1.25) · 3 준안정(≥ 1.0) · 4 위험 하한(≥ 0.5)
        · 5 위험 상한(> 0) · 6 불안정(SI = 0) · 0 = 자료 없음

격자 꼴 = 행 0 이 남쪽(y0) · 셀 (iy, ix) 중심 = (x0 + ix·dx, y0 + iy·dx) — E04 래스터와 같음.
"""

from __future__ import annotations

import math
from dataclasses import dataclass
from types import SimpleNamespace
from typing import Any

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.spatial import cKDTree

from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
    compute_receivers,
    condition_surface,
)

GAMMA_W = 9.81  # kN/m³
SI_EDGES = (1.5, 1.25, 1.0, 0.5)  # 등급 1~4 아래 경계 · 5 = (0, 0.5) · 6 = 0
CLASS_KEYS = ("", "stable", "moderate", "quasi", "lower", "upper", "defended")
RISK_CLASS = 4  # 이 등급부터 위험(SI < 1)
FLAT_SLOPE_RAD = 1e-4  # 이보다 평평하면 활동 없음(Fs = ∞)
X_STEPS = 20  # SI 확률의 x = R/T 중점 적분 마디


@dataclass(frozen=True)
class AreaParams:
    """한 번 해석 조건 — 단위 kN/m³ · kPa · ° · m · mm/h · m/s."""

    gamma_t: float
    gamma_sat: float
    c: float  # 점착력 아래(결정 Fs · SI 범위 아래)
    c_hi: float  # SI 범위 위
    phi: float
    phi_hi: float
    c_root: float = 0.0  # 뿌리 점착(나지 = 0)
    soil_depth_m: float = 1.5
    scenario: str = "sinmap"  # sinmap = 습윤지수(T/R 범위) · saturated = 지표 포화 m = 1(최악)
    tr_lo: float = 2000.0  # T/R 아래(m · 젖은 쪽)
    tr_hi: float = 3000.0  # T/R 위(m · 마른 쪽)


# ---------------------------------------------------------------------------
# 지형
# ---------------------------------------------------------------------------


def resample_dem(dem: dict[str, Any], cell_m: float) -> tuple[np.ndarray, float, float]:
    """DTM npz(x · y 오름 축 · z (ny, nx) · valid_mask) → 고른 격자 (z, x0, y0) · 밖 = NaN."""
    x = np.asarray(dem["x"], float)
    y = np.asarray(dem["y"], float)
    z = np.where(np.asarray(dem["valid_mask"], bool), np.asarray(dem["z"], float), np.nan)
    nx = int(math.floor((x[-1] - x[0]) / cell_m)) + 1
    ny = int(math.floor((y[-1] - y[0]) / cell_m)) + 1
    gx = x[0] + np.arange(nx) * cell_m
    gy = y[0] + np.arange(ny) * cell_m
    yy, xx = np.meshgrid(gy, gx, indexing="ij")
    f = RegularGridInterpolator((y, x), z, bounds_error=False, fill_value=np.nan)
    return f(np.column_stack([yy.ravel(), xx.ravel()])).reshape(ny, nx), float(x[0]), float(y[0])


def slope_rad(z: np.ndarray, cell_m: float) -> np.ndarray:
    """Horn 3×3 경사(rad) — 이웃에 NaN 이 있거나 격자 끝이면 NaN."""
    p = np.pad(z, 1, constant_values=np.nan)

    def at(dy: int, dx: int) -> np.ndarray:
        return p[1 + dy : 1 + dy + z.shape[0], 1 + dx : 1 + dx + z.shape[1]]

    gx = ((at(-1, 1) + 2 * at(0, 1) + at(1, 1)) - (at(-1, -1) + 2 * at(0, -1) + at(1, -1))) / (
        8 * cell_m
    )
    gy = ((at(1, -1) + 2 * at(1, 0) + at(1, 1)) - (at(-1, -1) + 2 * at(-1, 0) + at(-1, 1))) / (
        8 * cell_m
    )
    return np.arctan(np.hypot(gx, gy))


def catchment_area(z: np.ndarray, cell_m: float) -> np.ndarray:
    """D8 단위 집수면적 a(m) = 자기 포함 상류 셀 수 × cell. 밖 = NaN.

    ponytail: 지형 격자 끝 밖 상류는 모름(복사 DTM = 노선 띠) — 띠 끝 셀은 a 가 작게 나옴 ·
    넓은 DEM(D 예정노선 등)이면 그대로 맞음.
    """
    filled, valid = condition_surface(z.astype(np.float32))
    spec = SimpleNamespace(n_rows=z.shape[0], n_cols=z.shape[1], cell_m=cell_m)
    receiver, _ = compute_receivers(spec, filled, valid)
    flat = filled.ravel()
    order = np.flatnonzero(valid.ravel())
    order = order[np.argsort(-flat[order], kind="stable")].tolist()  # 높은 데서 낮은 데로
    acc = valid.ravel().astype(float).tolist()
    recv = receiver.tolist()
    for i in order:
        r = recv[i]
        if r != i:
            acc[r] += acc[i]
    out = np.asarray(acc).reshape(z.shape) * cell_m
    return np.where(valid, out, np.nan)


# ---------------------------------------------------------------------------
# 무한사면 · SINMAP
# ---------------------------------------------------------------------------


def wetness(beta: np.ndarray, a: np.ndarray, x: float) -> np.ndarray:
    """습윤지수 W = min(x·a / sinβ, 1) · x = R/T(1/m)."""
    with np.errstate(divide="ignore", invalid="ignore"):
        m = x * a / np.sin(beta)
    return np.clip(np.nan_to_num(m, nan=1.0, posinf=1.0), 0.0, 1.0)


def x_range(p: AreaParams) -> tuple[float, float]:
    """x = R/T 범위(마른 쪽, 젖은 쪽) — T/R 위 · 아래의 역수."""
    return 1.0 / max(p.tr_hi, p.tr_lo), 1.0 / p.tr_lo


def wet_m(beta: np.ndarray, a: np.ndarray, p: AreaParams, x: float | None = None) -> np.ndarray:
    """시나리오 포화비 — saturated = 1 · sinmap = W(x · 없으면 가장 젖은 x_hi)."""
    if p.scenario == "saturated":
        return np.ones_like(beta)
    return wetness(beta, a, x_range(p)[1] if x is None else x)


def _terms(beta: np.ndarray, m: np.ndarray, p: AreaParams) -> tuple[np.ndarray, np.ndarray]:
    """Fs = (c + c_r)·A + tanφ·B 의 A · B(평평한 셀 = ∞)."""
    z = p.soil_depth_m
    cos, sin = np.cos(beta), np.sin(beta)
    normal = (p.gamma_t * z * (1 - m) + (p.gamma_sat - GAMMA_W) * m * z) * cos**2
    shear = (p.gamma_t * z * (1 - m) + p.gamma_sat * m * z) * sin * cos
    flat = beta < FLAT_SLOPE_RAD
    with np.errstate(divide="ignore", invalid="ignore"):
        a = np.where(flat, np.inf, 1.0 / shear)
        b = np.where(flat, np.inf, normal / shear)
    return a, b


def infinite_slope_fs(beta: np.ndarray, m: np.ndarray, p: AreaParams) -> np.ndarray:
    """결정 안전율(c · φ = 범위 아래 = 보수) — 지식DB 02 §1.1 식."""
    a, b = _terms(beta, m, p)
    with np.errstate(invalid="ignore"):
        return np.where(
            np.isinf(a), np.inf, (p.c + p.c_root) * a + math.tan(math.radians(p.phi)) * b
        )


def _cdf_sum_uniform(
    s: np.ndarray, a1: np.ndarray, a2: np.ndarray, b1: np.ndarray, b2: np.ndarray
) -> np.ndarray:
    """P(X + Y ≤ s) · X ~ U[a1, a2] · Y ~ U[b1, b2] 독립 — 폭 0 이면 그 변수는 상수."""
    wa, wb = a2 - a1, b2 - b1
    eps = 1e-12

    def g(u: np.ndarray) -> np.ndarray:
        return np.maximum(u, 0.0) ** 2 / 2

    with np.errstate(divide="ignore", invalid="ignore"):
        both = (g(s - a1 - b1) - g(s - a2 - b1) - g(s - a1 - b2) + g(s - a2 - b2)) / (wa * wb)
        only_a = (s - a1 - b1) / wa
        only_b = (s - a1 - b1) / wb
    step = (a1 + b1 <= s).astype(float)
    out = np.where(
        (wa > eps) & (wb > eps), both, np.where(wa > eps, only_a, np.where(wb > eps, only_b, step))
    )
    return np.clip(np.nan_to_num(out, nan=0.0), 0.0, 1.0)


def stability_index(beta: np.ndarray, area: np.ndarray, p: AreaParams) -> np.ndarray:
    """SINMAP 안정 지수 — 하한 Fs_min ≥ 1 → Fs_min · 상한 Fs_max < 1 → 0 · 그 사이 → P(Fs > 1)."""
    t1, t2 = math.tan(math.radians(p.phi)), math.tan(math.radians(max(p.phi_hi, p.phi)))
    c1, c2 = p.c + p.c_root, max(p.c_hi, p.c) + p.c_root
    x_lo, x_hi = x_range(p)

    def ab(m: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        a, b = _terms(beta, m, p)
        return np.where(np.isinf(a), 0.0, a), np.where(np.isinf(b), 0.0, b)

    a, b = ab(wet_m(beta, area, p, x_hi))
    lo = c1 * a + t1 * b
    a, b = ab(wet_m(beta, area, p, x_lo))
    hi = c2 * a + t2 * b
    xs = (
        [x_hi]
        if p.scenario == "saturated"
        else x_lo + (x_hi - x_lo) * (np.arange(X_STEPS) + 0.5) / X_STEPS
    )
    prob = np.zeros_like(beta, dtype=float)
    one = np.ones_like(prob)
    for x in xs:
        a, b = ab(wet_m(beta, area, p, x))
        prob += 1.0 - _cdf_sum_uniform(one, c1 * a, c2 * a, t1 * b, t2 * b)
    prob /= len(xs)
    si = np.where(lo >= 1.0, lo, np.where(hi < 1.0, 0.0, prob))
    return np.where(beta < FLAT_SLOPE_RAD, np.inf, si)


def classify(si: np.ndarray) -> np.ndarray:
    """SI → 등급 1~6(0 = 자료 없음 · NaN)."""
    cls = np.full(si.shape, 6, dtype=np.uint8)
    cls[si > 0] = 5
    for k, edge in enumerate(reversed(SI_EDGES)):  # 0.5 → 4 · 1.0 → 3 · 1.25 → 2 · 1.5 → 1
        cls[si >= edge] = 4 - k
    cls[np.isnan(si)] = 0
    return cls


def analyze_grid(z: np.ndarray, cell_m: float, p: AreaParams) -> dict[str, np.ndarray]:
    """고른 격자 지형 → {slope(rad) · a · m · fs · si · cls} · 밖 = NaN / 0."""
    beta = slope_rad(z, cell_m)
    a = catchment_area(z, cell_m) if p.scenario != "saturated" else np.full(z.shape, np.nan)
    m = wet_m(beta, a, p)  # 가장 젖은 쪽(Fs = 하한과 같은 조건)
    fs = infinite_slope_fs(beta, m, p)
    si = stability_index(beta, a, p)
    nodata = np.isnan(beta)
    fs[nodata] = np.nan
    si[nodata] = np.nan
    m = np.where(nodata, np.nan, m)
    return {"slope": beta, "a": a, "m": m, "fs": fs, "si": si, "cls": classify(si)}


# ---------------------------------------------------------------------------
# 노선 띠 · 측점
# ---------------------------------------------------------------------------


def dense_line(line: np.ndarray, step_m: float = 0.5) -> np.ndarray:
    """중심선 [x, y, 누가] 를 step 간격으로 촘촘히(꺾인 점 사이 셀도 바른 누가)."""
    out = [line[:1]]
    for a, b in zip(line[:-1], line[1:]):
        n = max(1, int(math.ceil(math.hypot(*(b[:2] - a[:2])) / step_m)))
        out.append(a + (b - a) * np.linspace(0, 1, n + 1)[1:, None])
    return np.vstack(out)


def band_cells(
    shape: tuple[int, int], x0: float, y0: float, cell_m: float, line: np.ndarray, half_m: float
) -> tuple[np.ndarray, np.ndarray]:
    """노선 띠(중심선 ± half_m) 셀 → (띠 안, 셀의 가장 가까운 중심선 누가)."""
    dense = dense_line(np.asarray(line, float))
    yy, xx = np.meshgrid(
        y0 + np.arange(shape[0]) * cell_m, x0 + np.arange(shape[1]) * cell_m, indexing="ij"
    )
    dist, at = cKDTree(dense[:, :2]).query(np.column_stack([xx.ravel(), yy.ravel()]))
    return (dist <= half_m).reshape(shape), dense[at, 2].reshape(shape)


def station_table(
    chainages: list[float],
    cell_chain: np.ndarray,
    inside: np.ndarray,
    out: dict[str, np.ndarray],
    cell_m: float,
) -> list[dict[str, Any]]:
    """측점마다 맡은 띠 셀(이웃 측점과 중간 누가로 나눔) —
    최소 Fs · SI · 가장 나쁜 등급 · 위험 면적."""
    chs = np.asarray(sorted(chainages), float)
    if not len(chs):
        return []
    mids = (chs[:-1] + chs[1:]) / 2
    owner = np.searchsorted(mids, cell_chain)
    have = inside & (out["cls"] > 0)
    rows = []
    for k, ch in enumerate(chs):
        sel = have & (owner == k)
        n = int(sel.sum())
        if not n:
            rows.append({"chainage_m": float(ch), "cells": 0})
            continue
        cls = out["cls"][sel]
        fs = out["fs"][sel]
        si = out["si"][sel]
        rows.append(
            {
                "chainage_m": float(ch),
                "cells": n,
                "fs_min": round(float(np.min(fs)), 3),
                "si_min": round(float(np.min(si)), 3),
                "worst_class": int(cls.max()),
                "risk_area_m2": round(float((cls >= RISK_CLASS).sum()) * cell_m**2, 1),
                "slope_max_deg": round(float(np.degrees(np.nanmax(out["slope"][sel]))), 1),
            }
        )
    return rows


def class_areas(cls: np.ndarray, inside: np.ndarray, cell_m: float) -> dict[str, float]:
    """띠 안 등급별 면적 m²."""
    return {
        CLASS_KEYS[k]: round(float(((cls == k) & inside).sum()) * cell_m**2, 1) for k in range(1, 7)
    }