329 lines
14 KiB
Python
329 lines
14 KiB
Python
"""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)
|
||
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 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_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)
|
||
}
|