Aislo/resources/tester/test_536_2_e07_area.py

"""536-2 — E07 광역 자연사면 위험 `E07_StabArea_Engine` · `_Router`.

무한사면 해석해(건조 · 포화 · 점착) · SINMAP 안정 지수(닫힌 꼴 = 표본 · 등급 경계) ·
D8 집수면적 견본(기운 평면 · V 골짜기) · 강우 포화비 · 가짜 복사본으로 띠 · 측점표 · 저장본 · 지문.
"""

import base64
import json
import math

import numpy as np
import pytest

from E07_StabArea import E07_StabArea_Engine as engine
from E07_StabArea import E07_StabArea_Router as router_mod

SOIL = dict(gamma_t=18.0, gamma_sat=20.0, c=0.0, c_hi=0.0, phi=35.0, phi_hi=35.0)


def _fs(beta_deg: float, m: float, **kw) -> float:
    p = engine.AreaParams(**{**SOIL, **kw})
    return float(engine.infinite_slope_fs(np.array([math.radians(beta_deg)]), np.array([m]), p)[0])


def test_무한사면_해석해_건조_포화_점착():
    t = math.tan(math.radians(35)) / math.tan(math.radians(30))
    assert _fs(30, 0.0) == pytest.approx(t, rel=1e-12)  # 건조 c = 0 → tanφ / tanβ
    assert _fs(30, 1.0) == pytest.approx((20 - 9.81) / 20 * t, rel=1e-12)  # 포화 → γ' / γsat 배
    # 손계산 — c 5 · 뿌리 2 · z 1.5 · β 35° · m 0.5
    b, z = math.radians(35), 1.5
    w = 18 * z * 0.5 + 20 * z * 0.5
    n = (18 * z * 0.5 + (20 - 9.81) * z * 0.5) * math.cos(b) ** 2
    hand = (7 + n * math.tan(math.radians(35))) / (w * math.sin(b) * math.cos(b))
    assert _fs(35, 0.5, c=5.0, c_root=2.0) == pytest.approx(hand, rel=1e-12)
    assert math.isinf(_fs(0, 1.0))  # 평평 = 활동 없음


def test_SINMAP_안정지수_닫힌꼴_표본():
    p = engine.AreaParams(gamma_t=17.5, gamma_sat=18.5, c=0.0, c_hi=8.0, phi=25.0, phi_hi=38.0)
    beta = np.radians(np.array([10.0, 30.0, 35.0, 40.0, 60.0]))
    m = np.ones_like(beta)
    si = engine.stability_index(beta, m, p)
    a, b = engine._terms(beta, m, p)
    rng = np.random.default_rng(1)
    c = rng.uniform(0, 8, 400_000)
    t = rng.uniform(math.tan(math.radians(25)), math.tan(math.radians(38)), 400_000)
    for k in range(len(beta)):
        lo = 0 * a[k] + math.tan(math.radians(25)) * b[k]
        fs = c * a[k] + t * b[k]
        if lo >= 1:
            assert si[k] == pytest.approx(lo)  # Fs_min ≥ 1 → SI = Fs_min
        elif fs.max() < 1:
            assert si[k] == 0.0
        else:
            assert si[k] == pytest.approx((fs > 1).mean(), abs=3e-3)  # P(Fs > 1)
    assert si[0] >= 1 and 0 < si[1] < 1 and si[4] == 0  # 세 갈래 모두 지남
    # 폭 0(범위 없음) — 결정 Fs 와 같거나 0
    q = engine.AreaParams(**SOIL)
    fs = engine.infinite_slope_fs(beta, m, q)
    si0 = engine.stability_index(beta, m, q)
    assert np.allclose(si0, np.where(fs >= 1, fs, 0.0))


def test_SINMAP_등급_경계():
    si = np.array([np.inf, 1.5, 1.4999, 1.25, 1.2499, 1.0, 0.9999, 0.5, 0.4999, 1e-9, 0.0, np.nan])
    assert engine.classify(si).tolist() == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 0]


def test_D8_집수면적_기운평면_V골짜기():
    cell = 2.0
    ny, nx = 5, 7
    xx = np.arange(nx)[None, :] * cell
    plane = np.repeat(100.0 - 0.5 * xx, ny, axis=0)  # 동쪽으로 기움 → 행마다 동쪽으로 흐름
    a = engine.catchment_area(plane, cell)
    assert np.allclose(a, (np.arange(nx)[None, :] + 1) * cell)
    yy = np.arange(9)[:, None] * cell
    valley = 100.0 + np.abs(yy - 8.0) * 1.0 - 0.1 * np.arange(nx)[None, :] * cell
    a = engine.catchment_area(valley, cell)
    assert a[4, -1] == pytest.approx(9 * nx * cell)  # 전부가 골 출구로
    assert np.all(a[4, :-1] < a[4, 1:])  # 골 따라 늘어남


def test_강우_포화비():
    p = engine.AreaParams(**SOIL, scenario="rain", rain_mm_h=36.0, k_m_s=1e-4, soil_depth_m=2.0)
    beta = np.radians(np.array([30.0, 30.0]))
    a = np.array([10.0, 1e6])
    m = engine.wetness(beta, a, p)
    assert m[0] == pytest.approx(1e-5 * 10 / (2e-4 * 0.5))  # R a / (T sinβ)
    assert m[1] == 1.0
    assert engine.wetness(beta, a, engine.AreaParams(**SOIL)).tolist() == [1.0, 1.0]


def _decode(text: str, dtype: str, grid: dict) -> np.ndarray:
    return np.frombuffer(base64.b64decode(text), dtype=dtype).reshape(grid["ny"], grid["nx"])


@pytest.fixture()
def project(tmp_path, monkeypatch):
    """가짜 E 복사본 — 서쪽이 낮은 30° 비탈(동서 200 m · 남북 100 m) · 노선 = 남북 중심선."""
    x = np.arange(0.0, 200.5, 0.5)
    y = np.arange(0.0, 100.5, 0.5)
    z = np.repeat((x * math.tan(math.radians(30)))[None, :], len(y), axis=0)
    valid = np.ones_like(z, bool)
    valid[:, x > 180] = False  # 동쪽 끝 = 지형 없음
    dem = tmp_path / "dtm.npz"
    np.savez(
        dem,
        x=x.astype(np.float32),
        y=y.astype(np.float32),
        z=z.astype(np.float32),
        valid_mask=valid,
    )
    stations = [
        {"chainage_m": float(c), "label": f"No.{i}", "center_x": 100.0, "center_y": float(c)}
        for i, c in enumerate(range(0, 101, 20))
    ]
    route = [
        {"sequence_num": i, "model_x": 100.0, "model_y": float(c), "chainage_m": float(c)}
        for i, c in enumerate((0, 100))
    ]
    src = router_mod.source
    monkeypatch.setattr(src, "terrain_grid_path", lambda _r: dem)
    monkeypatch.setattr(src, "stations", lambda _r: stations)
    monkeypatch.setattr(src, "db_snapshot", lambda _r: {"route_points": route})
    return tmp_path


def test_길_띠_측점표_저장본(project):
    settings = {**router_mod.defaults(project), "band_m": 30.0, "cell_m": 2.0}
    assert settings["material_id"] == "colluvium" and settings["c_hi"] == settings["c"] + 5
    steps = []
    summary = router_mod.run(project, project, settings, "fp", lambda p, m: steps.append(p))
    doc = router_mod.load(project)
    assert steps == sorted(steps) and doc["fingerprint"] == "fp"
    g = doc["grid"]
    cls = _decode(doc["cls"], "u1", g)
    fs = _decode(doc["fs"], "<u2", g)
    # 붕적토 포화 30° — 손계산 Fs
    b = math.radians(30)
    n = (18.5 - 9.81) * 1.5 * math.cos(b) ** 2
    hand = (5 + n * math.tan(math.radians(27))) / (18.5 * 1.5 * math.sin(b) * math.cos(b))
    inner = fs[cls > 0]
    assert np.median(inner) / 1000 == pytest.approx(hand, abs=2e-3)
    assert summary["fs_min"] == pytest.approx(hand, abs=5e-3)
    # 띠 = 중심선 ± 30 m(격자 상자도 띠만큼) · 등급 면적 합 = 덮인 면적
    assert g["nx"] * g["dx"] <= 62 + 1e-9 and g["ny"] * g["dx"] <= 102 + 1e-9
    assert sum(summary["areas_m2"].values()) == pytest.approx(summary["covered_area_m2"])
    assert summary["covered_area_m2"] == pytest.approx((cls > 0).sum() * 4.0)
    rows = doc["stations"]
    assert [r["chainage_m"] for r in rows] == [0, 20, 40, 60, 80, 100]
    assert all(r["label"] == f"No.{i}" and r["worst_class"] >= 4 for i, r in enumerate(rows))
    assert summary["risk_stations"] == [0, 20, 40, 60, 80, 100]
    assert sum(r["cells"] for r in rows) == (cls > 0).sum()
    # 지문 = 설정에 따라 바뀜
    assert router_mod.fingerprint("s", settings) != router_mod.fingerprint(
        "s", {**settings, "c": 6}
    )
    json.dumps(doc)  # 그대로 JSON