269 lines
12 KiB
Python
269 lines
12 KiB
Python
"""536-2 — E07 광역 자연사면 위험 `E07_StabArea_Engine` · `_Router`.
|
||
|
||
무한사면 해석해(건조 · 포화 · 점착) · SINMAP 안정 지수(닫힌 꼴 = 표본 · 등급 경계) ·
|
||
D8 집수면적 견본(기운 평면 · V 골짜기) · 습윤지수 · SI 하한 / 상한 해석해 · 시나리오 둘(표준 · 최악) ·
|
||
가짜 복사본으로 띠 · 측점표 · 저장본 · 지문.
|
||
"""
|
||
|
||
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)
|
||
SAT = dict(scenario="saturated")
|
||
|
||
|
||
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, **SAT
|
||
)
|
||
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, **SAT)
|
||
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_습윤지수():
|
||
beta = np.radians(np.array([30.0, 30.0, 0.0]))
|
||
a = np.array([100.0, 1e6, 50.0])
|
||
w = engine.wetness(beta, a, 1 / 2000)
|
||
assert w[0] == pytest.approx(100 / (2000 * 0.5)) # R a / (T sinβ)
|
||
assert w[1] == 1.0 and w[2] == 1.0 # 넘침 · 평평 = 포화
|
||
p = engine.AreaParams(**SOIL, tr_lo=2000, tr_hi=3000)
|
||
assert engine.x_range(p) == (1 / 3000, 1 / 2000)
|
||
assert engine.wet_m(beta, a, p).tolist() == w.tolist() # 기본 = 가장 젖은 쪽
|
||
assert engine.wet_m(beta, a, engine.AreaParams(**SOIL, **SAT)).tolist() == [1.0] * 3
|
||
|
||
|
||
SINMAP = dict(gamma_t=18.0, gamma_sat=18.0, c=2.0, c_hi=6.0, phi=30.0, phi_hi=40.0)
|
||
|
||
|
||
def _sinmap_fs(beta, a, c, tanphi, x, z=1.5, g=18.0):
|
||
"""Pack 1998 식 — FS = [C + cosθ(1 − W r) tanφ] / sinθ · C = c / (h ρs g) · h = z cosθ."""
|
||
w = min(x * a / math.sin(beta), 1.0)
|
||
h = z * math.cos(beta)
|
||
return (c / (h * g) + math.cos(beta) * (1 - w * 9.81 / g) * tanphi) / math.sin(beta)
|
||
|
||
|
||
def test_SI_하한_상한_해석해():
|
||
"""γt = γsat 이면 엔진 = SINMAP 원식 · 하한 = (c, φ 아래 · T/R 아래) · 상한 = (위 · T/R 위)."""
|
||
p = engine.AreaParams(**SINMAP, tr_lo=1000, tr_hi=3000)
|
||
t1, t2 = math.tan(math.radians(30)), math.tan(math.radians(40))
|
||
beta = np.radians(np.array([15.0, 33.0, 38.0, 45.0, 70.0]))
|
||
area = np.array([20.0, 400.0, 300.0, 800.0, 5000.0])
|
||
si = engine.stability_index(beta, area, p)
|
||
rng = np.random.default_rng(7)
|
||
n = 300_000
|
||
c, t, x = rng.uniform(2, 6, n), rng.uniform(t1, t2, n), rng.uniform(1 / 3000, 1 / 1000, n)
|
||
seen = set()
|
||
for k, (b, ak) in enumerate(zip(beta, area)):
|
||
lo = _sinmap_fs(b, ak, 2.0, t1, 1 / 1000)
|
||
hi = _sinmap_fs(b, ak, 6.0, t2, 1 / 3000)
|
||
fs = engine.infinite_slope_fs(
|
||
beta[k : k + 1], engine.wet_m(beta[k : k + 1], area[k : k + 1], p), p
|
||
)
|
||
assert fs[0] == pytest.approx(lo, rel=1e-12) # 결정 Fs = 하한
|
||
if lo >= 1:
|
||
assert si[k] == pytest.approx(lo, rel=1e-12)
|
||
seen.add("lo")
|
||
elif hi < 1:
|
||
assert si[k] == 0.0
|
||
seen.add("zero")
|
||
else:
|
||
w = np.minimum(x * ak / math.sin(b), 1)
|
||
h = 1.5 * math.cos(b)
|
||
mc = (c / (h * 18) + math.cos(b) * (1 - w * 9.81 / 18) * t) / math.sin(b)
|
||
assert si[k] == pytest.approx((mc > 1).mean(), abs=3e-3) # P(Fs > 1) · T/R 도 범위
|
||
seen.add("prob")
|
||
assert seen == {"lo", "zero", "prob"}
|
||
|
||
|
||
def test_시나리오_둘_최악이_더_나쁨():
|
||
beta = np.radians(np.linspace(5, 60, 56))
|
||
area = np.full_like(beta, 50.0)
|
||
std = engine.stability_index(beta, area, engine.AreaParams(**SINMAP))
|
||
worst = engine.stability_index(beta, area, engine.AreaParams(**SINMAP, **SAT))
|
||
assert np.all(worst <= std + 1e-12) # 완전 포화 = 아래 경계
|
||
assert (engine.classify(worst) >= engine.RISK_CLASS).sum() > (
|
||
engine.classify(std) >= engine.RISK_CLASS
|
||
).sum()
|
||
# 다 젖는 집수(x·a / sinβ ≥ 1) 면 두 시나리오가 같음
|
||
big = np.full_like(beta, 1e7)
|
||
assert np.allclose(
|
||
engine.stability_index(beta, big, engine.AreaParams(**SINMAP)), worst, rtol=1e-12
|
||
)
|
||
|
||
|
||
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_측점_등급_문턱_경계():
|
||
"""위험 셀 비율 < 5 % 낮음 · 5 % 이상 20 % 미만 보통 · 20 % 이상 높음 (경계 포함은 윗등급)."""
|
||
cases = [(0, 100, "low"), (4, 100, "low"), (5, 100, "mid"), (19, 100, "mid")]
|
||
cases += [(20, 100, "high"), (100, 100, "high"), (1, 20, "mid"), (1, 21, "low"), (1, 5, "high")]
|
||
for risk, cells, want in cases:
|
||
assert engine.station_level(risk, cells) == want, (risk, cells)
|
||
# 사용자가 바꾼 문턱
|
||
assert engine.station_level(9, 100, 10, 30) == "low"
|
||
assert engine.station_level(30, 100, 10, 30) == "high"
|
||
|
||
|
||
def test_측점_영향_구역과_비율():
|
||
"""측점 0 · 10 · 30 → 영향 구역 = 이웃과 중간 누가로 나눔 · 띠 안 셀만 · 자료 없음 제외."""
|
||
chain = np.arange(0.0, 40.0, 1.0)[None, :] # 1 × 40 셀 · 셀 i 누가 = i
|
||
cls = np.ones(chain.shape, np.uint8)
|
||
cls[0, [0, 1, 2, 3, 4]] = 4 # 측점 0 영역(0~5 · 정확히 중간 5 는 앞 측점) 6 셀 중 5 → 83 %
|
||
cls[0, [6, 7]] = 5 # 측점 10 영역(6~20) 셀 14(19 는 띠 밖) · 위험 6 · 7 · 20 → 21 %
|
||
cls[0, [20, 21]] = 6 # 20 은 측점 10 몫 · 측점 30 영역(21~39) 은 21 하나
|
||
cls[0, 39] = 0 # 자료 없음 → 셀 수에서 빠짐(18 셀) 1/18 → 5.6 %
|
||
inside = np.ones(chain.shape, bool)
|
||
inside[0, 19] = False # 띠 밖 → 빠짐
|
||
out = {k: np.ones(chain.shape) for k in ("fs", "si", "slope")} | {"cls": cls}
|
||
rows = engine.station_table([30.0, 0.0, 10.0], chain, inside, out, 2.0)
|
||
assert [r["chainage_m"] for r in rows] == [0.0, 10.0, 30.0]
|
||
assert [r["cells"] for r in rows] == [6, 14, 18]
|
||
assert [r["level"] for r in rows] == ["high", "high", "mid"]
|
||
assert rows[0]["risk_ratio"] == round(5 / 6, 4) and rows[1]["risk_ratio"] == round(3 / 14, 4)
|
||
assert rows[2]["risk_ratio"] == round(1 / 18, 4) and rows[2]["risk_area_m2"] == 4.0
|
||
# 문턱을 올리면 등급이 내려감
|
||
again = engine.station_table([0.0, 10.0, 30.0], chain, inside, out, 2.0, 15.0, 50.0)
|
||
assert [r["level"] for r in again] == ["high", "mid", "low"]
|
||
|
||
|
||
def test_측점_문턱_설정_검증():
|
||
base = router_mod.AreaRequest.model_fields
|
||
assert (base["level_lo_pct"].default, base["level_hi_pct"].default) == (5.0, 20.0)
|
||
ok = dict(gamma_t=18.0, gamma_sat=20.0, c=0.0, c_hi=0.0, phi=30.0, phi_hi=30.0)
|
||
router_mod.AreaRequest(**ok, level_lo_pct=10, level_hi_pct=10)
|
||
with pytest.raises(ValueError):
|
||
router_mod.AreaRequest(**ok, level_lo_pct=30, level_hi_pct=10)
|
||
|
||
|
||
def test_길_띠_측점표_저장본(project):
|
||
base = router_mod.defaults(project)
|
||
assert base["scenario"] == "sinmap" and (base["tr_lo"], base["tr_hi"]) == (2000.0, 3000.0)
|
||
settings = {**base, "band_m": 30.0, "cell_m": 2.0, "scenario": "saturated", "c_root": 0.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 all(r["level"] == "high" and r["risk_ratio"] > 0.99 for r in rows) # 완전 포화 30° 비탈
|
||
assert summary["station_levels"] == {"low": 0, "mid": 0, "high": 6}
|
||
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
|