eomsangdonandClaude Opus 5.5
0fd6437628
- 엔진 E07_StabArea_Engine — DTM 고른 격자 · Horn 경사 · D8 집수면적(B04 웅덩이 채움 · 수신 셀 재사용) · 무한사면 Fs(포화비 m) · SINMAP 안정 지수(c · tanφ 고른 분포 닫힌 꼴) · 6 등급 · 노선 띠 · 측점표 - 길 E07_StabArea_Router — POST 해석(뒤에서 · 진행률 · 지문 같으면 cached) · GET 기본값(E02 붕적토 · 범위 위 +5 kPa · +3°) · 파생 저장본 _derived/area.json(격자 base64) - E05 보고서 s7 광역요약 — 조건 · 등급별 면적 · 위험 측점 줄 · 양식 · 자리표 덧붙임 · 해석 지문에 E07 저장본 - 시험 test_536_2 새로 · test_534_3 두 파일은 시트 7 · PDF 쪽 +1 로 고침(광역요약 시트 더함) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
284 lines
12 KiB
Python
284 lines
12 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 = 1(지표 포화 · 기본) | min(R·a / (T·sinβ), 1) · T = K·z(정상류 · SINMAP · SHALSTAB)
|
||
안전율 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] 고른 분포
|
||
Fs_min ≥ 1 → SI = Fs_min · Fs_max < 1 → SI = 0
|
||
· 그 사이 → SI = P(Fs > 1)(두 고른 분포 합 · 닫힌 꼴)
|
||
등급 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 = ∞)
|
||
|
||
|
||
@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 = "saturated" # saturated = 지표 포화 m = 1 · rain = 정상류 강우
|
||
rain_mm_h: float = 30.0
|
||
k_m_s: float = 1e-4 # 포화 투수계수
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 지형
|
||
# ---------------------------------------------------------------------------
|
||
|
||
|
||
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, p: AreaParams) -> np.ndarray:
|
||
"""포화비 m(0~1) — 지표 포화 = 1 · 강우 = min(R·a/(T sinβ), 1)."""
|
||
if p.scenario == "saturated":
|
||
return np.ones_like(beta)
|
||
recharge = p.rain_mm_h / 1000.0 / 3600.0 # m/s
|
||
transmissivity = p.k_m_s * p.soil_depth_m # m²/s
|
||
with np.errstate(divide="ignore", invalid="ignore"):
|
||
m = recharge * a / (transmissivity * np.sin(beta))
|
||
return np.clip(np.nan_to_num(m, nan=1.0, posinf=1.0), 0.0, 1.0)
|
||
|
||
|
||
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, m: np.ndarray, p: AreaParams) -> np.ndarray:
|
||
"""SINMAP 안정 지수 — Fs_min ≥ 1 → Fs_min · Fs_max < 1 → 0 · 그 사이 → P(Fs > 1)."""
|
||
a, b = _terms(beta, m, p)
|
||
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
|
||
flat = np.isinf(a)
|
||
a = np.where(flat, 0.0, a)
|
||
b = np.where(flat, 0.0, b)
|
||
lo, hi = c1 * a + t1 * b, c2 * a + t2 * b
|
||
prob = 1.0 - _cdf_sum_uniform(np.ones_like(a), c1 * a, c2 * a, t1 * b, t2 * b)
|
||
si = np.where(lo >= 1.0, lo, np.where(hi < 1.0, 0.0, prob))
|
||
return np.where(flat, 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 == "rain" else np.full(z.shape, np.nan)
|
||
m = wetness(beta, a, p)
|
||
fs = infinite_slope_fs(beta, m, p)
|
||
si = stability_index(beta, m, 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)
|
||
}
|