"""E07 광역 자연사면 위험 — 무한사면 안전율 · SINMAP 안정 지수 (요구 345 · 351 · PLAN 536-2). 순수 계산(storage 모름) — D 예정노선 재해 안전성도 같은 엔진을 씀(한 벌). 근거 = 지식DB `original/예정노선/02` §1.1 · §2 ① · `original/안정성평가/04`(붕적토 범위). 지형 DEM → 고른 격자(cell_m) · 경사 β(Horn 3×3) · D∞ 흐름(Tarboton 1997 · SINMAP 원본 · 웅덩이 채움 = B04 배수유역 함수 재사용 · 내리막 없는 평지만 D8) · 단위 집수면적 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.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) LEVEL_KEYS = ("low", "mid", "high") # 측점 등급 — 영향 구역 안 위험 셀 비율 기준 LEVEL_LO_PCT = 5.0 # 낮음 < 이 % ≤ 보통 LEVEL_HI_PCT = 20.0 # 보통 < 이 % ≤ 높음(이 % 이상 = 높음) 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 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)) #: D∞ 삼각 면 여덟(Tarboton 1997 표 1) — (e1 = 곁 이웃 행·열, e2 = 대각 이웃 행·열) _FACETS = ( ((0, 1), (-1, 1)), ((-1, 0), (-1, 1)), ((-1, 0), (-1, -1)), ((0, -1), (-1, -1)), ((0, -1), (1, -1)), ((1, 0), (1, -1)), ((1, 0), (1, 1)), ((0, 1), (1, 1)), ) def dinf_receivers(z: np.ndarray, cell_m: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """D∞ — 셀마다 가장 가파른 면의 받는 셀 둘(평면 번호 · 없으면 −1) · e1 몫 = 1 − r/(π/4). 면 기울기 s1 = (e0 − e1)/d · s2 = (e1 − e2)/d · r = atan(s2/s1) 를 [0, π/4] 로 자름(끝이면 한 셀로). NaN(밖) 이웃이 낀 면은 안 봄.""" ny, nx = z.shape p = np.pad(z.astype(float), 1, constant_values=np.nan) iy, ix = np.indices((ny, nx)) best = np.zeros((ny, nx)) r1 = np.full((ny, nx), -1, np.int64) r2 = np.full((ny, nx), -1, np.int64) f1 = np.ones((ny, nx)) quarter = math.pi / 4 with np.errstate(invalid="ignore"): for (a, b), (c, d) in _FACETS: e1 = p[1 + a : 1 + a + ny, 1 + b : 1 + b + nx] e2 = p[1 + c : 1 + c + ny, 1 + d : 1 + d + nx] s1 = (z - e1) / cell_m s2 = (e1 - e2) / cell_m r = np.arctan2(s2, s1) s = np.where( r < 0, s1, np.where(r > quarter, (z - e2) / (cell_m * math.sqrt(2)), np.hypot(s1, s2)), ) up = np.isfinite(s) & (s > best) best[up] = s[up] r1[up] = ((iy + a) * nx + ix + b)[up] r2[up] = ((iy + c) * nx + ix + d)[up] f1[up] = 1.0 - np.clip(r, 0.0, quarter)[up] / quarter return r1, r2, f1 def catchment_area(z: np.ndarray, cell_m: float) -> np.ndarray: """D∞ 단위 집수면적 a(m) = 자기 포함 상류 셀 몫 합 × cell. 밖 = NaN · 내리막 없는 셀(채운 평지) = D8. ponytail: 지형 격자 끝 밖 상류는 모름 — E07 은 도엽 + 여유 400 m 로 셈(`E07_StabArea_Sheet` · 포화 근처라 차 0.1 %). """ 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) r1, r2, f1 = dinf_receivers(np.where(valid, filled, np.nan), cell_m) d8 = r1 < 0 r1[d8] = receiver.reshape(z.shape)[d8] r2[d8] = -1 f1[d8] = 1.0 flat = filled.ravel() order = np.flatnonzero(valid.ravel()) order = order[np.argsort(-flat[order], kind="stable")].tolist() # 높은 데서 낮은 데로 acc = valid.ravel().astype(float).tolist() one, two, share = r1.ravel().tolist(), r2.ravel().tolist(), f1.ravel().tolist() for i in order: j, f = one[i], share[i] if j != i: acc[j] += acc[i] * f if f < 1.0: acc[two[i]] += acc[i] * (1.0 - f) 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) return stability(beta, a, p) def stability( beta: np.ndarray, a: np.ndarray, p: AreaParams, si: np.ndarray | None = None ) -> dict[str, np.ndarray]: """경사 · 집수면적 → {slope · a · m · fs · si · cls} — si 를 주면(같은 물성 미리 셈) 그대로.""" m = wet_m(beta, a, p) # 가장 젖은 쪽(Fs = 하한과 같은 조건) fs = infinite_slope_fs(beta, m, p) si = stability_index(beta, a, p) if si is None else np.array(si, float) 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_level( risk_cells: int, cells: int, lo_pct: float = LEVEL_LO_PCT, hi_pct: float = LEVEL_HI_PCT ) -> str: """측점 등급 — 위험(SI < 1) 셀 비율 < lo_pct 낮음 · lo_pct 이상 hi_pct 미만 보통 · hi_pct 이상 높음. (정수 곱 비교 — 문턱 경계에서 소수 오차 없음)""" if risk_cells * 100 >= hi_pct * cells: return "high" return "mid" if risk_cells * 100 >= lo_pct * cells else "low" def station_table( chainages: list[float], cell_chain: np.ndarray, inside: np.ndarray, out: dict[str, np.ndarray], cell_m: float, lo_pct: float = LEVEL_LO_PCT, hi_pct: float = LEVEL_HI_PCT, ) -> 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] risk = int((cls >= RISK_CLASS).sum()) 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(risk * cell_m**2, 1), "risk_ratio": round(risk / n, 4), "level": station_level(risk, n, lo_pct, hi_pct), "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) }