From fc4f3ef1ae5ae625441532114b3c8796a428fdce Mon Sep 17 00:00:00 2001 From: umsangdon Date: Fri, 30 Oct 2026 22:13:07 +0900 Subject: [PATCH] =?UTF-8?q?feat(E07=20=EA=B4=91=EC=97=AD):=20=EC=A7=80?= =?UTF-8?q?=ED=98=95=20=3D=20=EA=B3=B5=EC=9A=A9=20=EC=A7=80=ED=98=95=20?= =?UTF-8?q?=EC=B0=BD=EA=B3=A0=20=C2=B7=20=EB=8F=84=EC=97=BD=EB=A7=88?= =?UTF-8?q?=EB=8B=A4=20=EA=B2=BD=EC=82=AC=20=C2=B7=20=EC=A7=91=EC=88=98?= =?UTF-8?q?=EB=A9=B4=EC=A0=81=20=C2=B7=20=EA=B8=B0=EB=B3=B8=20=EB=AC=BC?= =?UTF-8?q?=EC=84=B1=20SI=20=EB=AF=B8=EB=A6=AC=20=C2=B7=20=EB=AC=BC?= =?UTF-8?q?=EC=84=B1=20=EB=B0=94=EA=BF=88=20=3D=20=EC=B9=B8=20=EA=B3=84?= =?UTF-8?q?=EC=82=B0=EB=A7=8C=20(PLAN=20610-1=20=C2=B7=20=EC=9A=94?= =?UTF-8?q?=EA=B5=AC=20487=20=C2=B7=20488)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 새 E07_StabArea_Sheet — 도엽 지형 + 여유 400 m 로 경사 · D∞ 집수면적 → 창고 도엽 폴더 e07__<격자>m.npz · 기본 물성(마스터 E02 붕적토) SI = e07si_ · 열쇠 = 판 + 여유 + 이웃 도엽 지형 열쇠 - 라우터 — LAS DTM 복사본 안 씀 · 노선 띠 상자 → 도엽 이어 붙임 · 좌표계 = 동결 DB 지형 모델 crs_epsg · 물성이 기본과 같으면 미리 셈 SI · 저장 상자 = 띠 안 지형 있는 칸 - 엔진 — stability() 로 나눔 · 안 쓰는 resample_dem 지움 - E07 안내 문구 — 지형 = 전국 등고선 공용 지형 - 시험 — 536_2 가짜 복사본 DTM → 가짜 창고(LAS 안 씀) · 새 610_1(도엽 이음 = 통째 · 위쪽 유역 안 잘림 · 두 번째 · 물성 바꿈 = 흐름 안 셈 · 지형 열쇠 바뀌면 다시) Co-Authored-By: Claude Opus 5.5 --- E07_StabArea/E07_StabArea_Engine.py | 28 ++-- E07_StabArea/E07_StabArea_Router.py | 91 +++++++---- E07_StabArea/E07_StabArea_Sheet.py | 183 +++++++++++++++++++++++ resources/tester/test_536_2_e07_area.py | 21 +-- resources/tester/test_610_1_e07_sheet.py | 138 +++++++++++++++++ ui_template/ui_template_locale_e07.ts | 4 +- 6 files changed, 400 insertions(+), 65 deletions(-) create mode 100644 E07_StabArea/E07_StabArea_Sheet.py create mode 100644 resources/tester/test_610_1_e07_sheet.py diff --git a/E07_StabArea/E07_StabArea_Engine.py b/E07_StabArea/E07_StabArea_Engine.py index e6ce7ee40..f116dc49b 100644 --- a/E07_StabArea/E07_StabArea_Engine.py +++ b/E07_StabArea/E07_StabArea_Engine.py @@ -28,7 +28,6 @@ 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 ( @@ -69,20 +68,6 @@ class AreaParams: # --------------------------------------------------------------------------- -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) @@ -148,8 +133,8 @@ def dinf_receivers(z: np.ndarray, cell_m: float) -> tuple[np.ndarray, np.ndarray def catchment_area(z: np.ndarray, cell_m: float) -> np.ndarray: """D∞ 단위 집수면적 a(m) = 자기 포함 상류 셀 몫 합 × cell. 밖 = NaN · 내리막 없는 셀(채운 평지) = D8. - ponytail: 지형 격자 끝 밖 상류는 모름(복사 DTM = 노선 띠) — 띠 끝 셀은 a 가 작게 나옴 · - 넓은 DEM(D 예정노선 등)이면 그대로 맞음. + 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) @@ -284,9 +269,16 @@ def analyze_grid(z: np.ndarray, cell_m: float, p: AreaParams) -> dict[str, np.nd """고른 격자 지형 → {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) + 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 diff --git a/E07_StabArea/E07_StabArea_Router.py b/E07_StabArea/E07_StabArea_Router.py index 8cedd5840..ef54f4af2 100644 --- a/E07_StabArea/E07_StabArea_Router.py +++ b/E07_StabArea/E07_StabArea_Router.py @@ -5,8 +5,9 @@ - `GET …/stab/area/progress` — 진행률 - `GET …/stab/area` — 기본값(E02 붕적토) · 물성 목록 · 노선 · 측점 · 저장본(격자 · 측점표 · 요약) -계산 자리 = ③ 서버 단독(엔진 `E07_StabArea_Engine` · 순수). 지형 = 복사본 원지형 DTM -(노선면 합성 안 함 — 노선 위 · 아래 자연사면). +계산 자리 = ③ 서버 단독(엔진 `E07_StabArea_Engine` · 순수). 지형 = 공용 지형 창고(전국 등고선 2 m · +LAS 안 씀 · PLAN 610-1) — 도엽마다 경사 · 집수면적 · 기본 물성 SI 를 미리(`E07_StabArea_Sheet`) · +물성이 다르면 저장된 경사 · 집수면적으로 칸 계산만. 파생 저장본 = `E07_StabArea/_derived/area.json`(지문이 같으면 안 셈). 권한 = main.py `protected_with_company`(지정받은 회사 = GET 만). @@ -43,11 +44,16 @@ from config.config_db import get_db_pool from E01_StabRegister import E01_StabRegister_Source as source from E02_StabGround import E02_StabGround_Model as ground_model from E07_StabArea import E07_StabArea_Engine as engine +from E07_StabArea import E07_StabArea_Sheet as sheet logger = logging.getLogger(__name__) router = APIRouter(prefix="/api/projects", tags=["E07 Stab Area"]) -ENGINE_FILES = ("E07_StabArea_Engine.py", "E07_StabArea_Router.py") # 엔진 판 = 이 파일들 해시 +ENGINE_FILES = ( + "E07_StabArea_Engine.py", + "E07_StabArea_Router.py", + "E07_StabArea_Sheet.py", +) # 엔진 판 = 이 파일들 해시 MATERIAL_KEYS = ("gamma_t", "gamma_sat", "c", "phi") RESULT = Path("E07_StabArea") / "_derived" / "area.json" DEFAULT_MATERIAL = "colluvium" # 표층 붕괴에 가장 취약(지식DB 안정성평가/04 붕적토) @@ -106,7 +112,31 @@ def _ground(project_root: Path) -> dict[str, Any]: def defaults(project_root: Path) -> dict[str, Any]: """기본 조건 — E02 붕적토(없으면 첫 물성) · 범위 위 = + 5 kPa · + 3°.""" - materials = _ground(project_root).get("materials") or [] + return _defaults_for(_ground(project_root).get("materials") or []) + + +PARAM_KEYS = ( + "gamma_t", + "gamma_sat", + "c", + "c_hi", + "phi", + "phi_hi", + "c_root", + "soil_depth_m", + "scenario", + "tr_lo", + "tr_hi", +) + + +def default_params() -> engine.AreaParams: + """공용 창고에 미리 셀 기본 물성 — 마스터 E02 붕적토 + 기본 조건.""" + base = _defaults_for(ground_model.pick(ground_model.load_defaults()).get("materials") or []) + return engine.AreaParams(**{k: base[k] for k in PARAM_KEYS}) + + +def _defaults_for(materials: list[dict[str, Any]]) -> dict[str, Any]: mat = next((m for m in materials if m.get("id") == DEFAULT_MATERIAL), None) or materials[0] return { **{ @@ -167,6 +197,12 @@ def _route_line(root: Path) -> np.ndarray: return np.array([[p["model_x"], p["model_y"], p["chainage_m"]] for p in points], float) +def _epsg(root: Path) -> int: + """작업 좌표계 — 동결 DB 지형 모델 crs_epsg(B 지형 PRJ 와 같은 값) · 없으면 중부원점.""" + models = source.db_snapshot(root).get("surface_models") or [] + return int(next((m["crs_epsg"] for m in models if m.get("crs_epsg")), 5186)) + + def _b64(values: np.ndarray, dtype: str) -> str: return base64.b64encode(np.ascontiguousarray(values, dtype=dtype).tobytes()).decode("ascii") @@ -186,29 +222,22 @@ def run( """해석 한 번 → 저장본 쓰고 요약 돌려줌.""" started = time.perf_counter() cell = float(settings["cell_m"]) - report(5, "지형 격자") - with np.load(source.terrain_grid_path(root)) as npz: - dem = {k: npz[k] for k in ("x", "y", "z", "valid_mask")} - z, x0, y0 = engine.resample_dem(dem, cell) - report(25, "경사 · 집수면적 · 안전율") - keys = ( - "gamma_t", - "gamma_sat", - "c", - "c_hi", - "phi", - "phi_hi", - "c_root", - "soil_depth_m", - "scenario", - "tr_lo", - "tr_hi", + line = _route_line(root) + pad = float(settings["band_m"]) + cell + box = (line[:, 0].min() - pad, line[:, 1].min() - pad) + box = (*box, line[:, 0].max() + pad, line[:, 1].max() + pad) + base = default_params() + grid = sheet.load_area( + _epsg(root), box, cell, base, lambda pct, stage: report(2 + pct * 68 // 100, stage) ) - out = engine.analyze_grid(z, cell, engine.AreaParams(**{k: settings[k] for k in keys})) - report(75, "노선 띠 · 측점") + x0, y0 = grid["x0"], grid["y0"] + report(72, "안전율 · 안정 지수") + params = engine.AreaParams(**{k: settings[k] for k in PARAM_KEYS}) + out = engine.stability(grid["slope"], grid["a"], params, grid["si"] if params == base else None) + report(80, "노선 띠 · 측점") stations = sorted(source.stations(root), key=lambda s: float(s["chainage_m"])) inside, chain = engine.band_cells( - z.shape, x0, y0, cell, _route_line(root), float(settings["band_m"]) + grid["slope"].shape, x0, y0, cell, line, float(settings["band_m"]) ) rows = engine.station_table( [s["chainage_m"] for s in stations], @@ -222,14 +251,14 @@ def run( labels = {round(float(s["chainage_m"]), 3): s.get("label") for s in stations} for row in rows: row["label"] = labels.get(round(row["chainage_m"], 3)) - # 띠 상자로 잘라 저장 - iy, ix = np.nonzero(inside) - if not len(iy): - raise ValueError("노선 띠가 지형 격자와 겹치지 않습니다.") - win = (slice(iy.min(), iy.max() + 1), slice(ix.min(), ix.max() + 1)) - keep = inside[win] & (out["cls"][win] > 0) - areas = engine.class_areas(out["cls"], inside, cell) + # 띠 안 지형 있는 상자로 잘라 저장 have = inside & (out["cls"] > 0) + iy, ix = np.nonzero(have) + if not len(iy): + raise ValueError("노선 띠에 지형이 없습니다.") + win = (slice(iy.min(), iy.max() + 1), slice(ix.min(), ix.max() + 1)) + keep = have[win] + areas = engine.class_areas(out["cls"], inside, cell) summary = { "areas_m2": areas, "band_area_m2": round(float(inside.sum()) * cell**2, 1), diff --git a/E07_StabArea/E07_StabArea_Sheet.py b/E07_StabArea/E07_StabArea_Sheet.py new file mode 100644 index 000000000..41e72e541 --- /dev/null +++ b/E07_StabArea/E07_StabArea_Sheet.py @@ -0,0 +1,183 @@ +"""E07 광역 — 도엽마다 경사 · 집수면적 · 기본 물성 안정 지수를 공용 창고에 미리 +(PLAN 610-1 · 요구 487 · 488). + +지형 원천 = 공용 지형 창고(`common_util_terrain_store` · 전국 등고선 2 m) — LAS DTM 복사본 안 씀. +자리 = 창고 도엽 폴더 `storage/_shared/terrain/<도엽>/` + e07__<격자>m.npz 경사 β(rad) · 단위 집수면적 a(m) — 도엽 지형 + 여유 400 m 로 셈 → + 도엽 안쪽 칸만 저장 → 이어 붙여도 위쪽 유역 안 잘림(분석 486 3-1) + e07si__<격자>m.npz 기본 물성(마스터 E02 붕적토) 안정 지수 — 물성이 다르면 부르는 쪽이 + 저장된 β · a 로 칸 계산만(흐름 다시 안 셈) +열쇠 = 판 + 여유 + 쓴 도엽 지형 열쇠(전국 파일 지문 · 방식 판) · SI 는 + 물성. +격자 = 창고와 같음(칸 중심 (k + ½)·cell · 행 0 = 남쪽). +""" + +from __future__ import annotations + +import json +import math +from dataclasses import asdict +from pathlib import Path +from typing import Any, Callable + +import numpy as np + +from common_util import common_util_terrain_store as store +from common_util.common_util_atomic import atomic_write_npz +from common_util.common_util_json import atomic_write_json +from E07_StabArea import E07_StabArea_Engine as engine + +SHEET_VER = 1 +MARGIN_M = 400.0 + +Report = Callable[[int, str], None] +Box = tuple[float, float, float, float] + + +def _paths(sheet_no: str, epsg: int, cell_m: float) -> dict[str, Path]: + folder = store.TERRAIN_ROOT / str(sheet_no) + stem = f"{int(epsg)}_{cell_m:g}m" + return { + "flow": folder / f"e07_{stem}.npz", + "key": folder / f"e07_{stem}.key.json", + "si": folder / f"e07si_{stem}.npz", + "lock": folder / f".e07_{stem}.lock", + } + + +def _dem_stamp(sheet_no: str, epsg: int, cell_m: float) -> Any: + """도엽 지형 열쇠(셈 시간 뺌) — 지형이 다시 셈되면 달라짐.""" + path = store.TERRAIN_ROOT / str(sheet_no) / f"dem_{int(epsg)}_{cell_m:g}m.key.json" + try: + key = json.loads(path.read_text(encoding="utf-8")) + except (OSError, ValueError): + return None + key.pop("seconds", None) + return key + + +def _read_key(path: Path) -> Any: + try: + return json.loads(path.read_text(encoding="utf-8")) + except (OSError, ValueError): + return None + + +def ensure_sheet( + sheet_no: str, + epsg: int, + cell_m: float, + defaults: engine.AreaParams, + report: Report | None = None, +) -> Path | None: + """도엽 하나의 E07 흐름(β · a) + 기본 물성 SI — 있으면 그대로 · 없거나 지형이 바뀌었으면 셈. + + 빈 도엽(바다) = None · 지형 받기 실패 = RuntimeError(저장 안 함 · 다음에 다시).""" + + def say(pct: int, stage: str) -> None: + if report is not None: + report(pct, stage) + + sheet_no = str(sheet_no) + paths = _paths(sheet_no, epsg, cell_m) + say(0, f"지형 · {sheet_no}") + dem = store.ensure_sheet_terrain(sheet_no, epsg, cell_m, lambda p, s: say(p * 30 // 100, s)) + if dem is None: + return None + with np.load(dem) as data: + ix0, iy0 = int(data["ix0"]), int(data["iy0"]) + core_ok = np.isfinite(data["z"]) + ny, nx = core_ok.shape + pad = math.ceil(MARGIN_M / cell_m) + box = ( + (ix0 - pad) * cell_m, + (iy0 - pad) * cell_m, + (ix0 + nx + pad) * cell_m, + (iy0 + ny + pad) * cell_m, + ) + around = store.sheets_for_box(box, epsg) + stamp: dict[str, Any] = {} + for i, no in enumerate(around): # 이웃 지형(처음이면 셈 · 있으면 열쇠만 봄) + say(30 + 40 * i // len(around), f"이웃 지형 {i + 1}/{len(around)} · {no}") + store.ensure_sheet_terrain(no, epsg, cell_m) + stamp[no] = _dem_stamp(no, epsg, cell_m) + key = {"ver": SHEET_VER, "margin_m": MARGIN_M, "dem": stamp} + if not (_read_key(paths["key"]) == key and paths["flow"].is_file()): + with store._Lock(paths["lock"]): + if not (_read_key(paths["key"]) == key and paths["flow"].is_file()): + say(70, f"경사 · 집수면적 · {sheet_no}") + region = store.load_terrain(epsg, sheets=around, box=box, cell_m=cell_m) + if region["failed"]: + raise RuntimeError(f"도엽 {sheet_no} 이웃 지형 실패: {region['failed']}") + z = region["z"] + inner = (slice(pad, pad + ny), slice(pad, pad + nx)) + beta = engine.slope_rad(z, cell_m)[inner] + area = engine.catchment_area(z, cell_m)[inner] + beta[~core_ok] = np.nan + area[~core_ok] = np.nan + paths["key"].unlink(missing_ok=True) + atomic_write_npz( + paths["flow"], + slope=beta.astype(np.float32), + a=area.astype(np.float32), + ix0=np.int64(ix0), + iy0=np.int64(iy0), + ) + atomic_write_json(paths["key"], key) + si_key = json.dumps({"flow": key, "params": asdict(defaults)}, sort_keys=True) + if _si_key(paths["si"]) != si_key: + say(90, f"기본 물성 안정 지수 · {sheet_no}") + with np.load(paths["flow"]) as data: + si = engine.stability_index( + data["slope"].astype(float), data["a"].astype(float), defaults + ) + atomic_write_npz(paths["si"], si=si.astype(np.float32), key=np.array(si_key)) + say(100, f"완료 · {sheet_no}") + return paths["flow"] + + +def _si_key(path: Path) -> str | None: + try: + with np.load(path) as data: + return str(data["key"]) + except (OSError, KeyError, ValueError): + return None + + +def load_area( + epsg: int, + box: Box, + cell_m: float, + defaults: engine.AreaParams, + report: Report | None = None, +) -> dict[str, Any]: + """TM 상자 → 이어 붙인 {x0 · y0(첫 칸 중심) · slope · a · si(기본 물성)} · 없는 칸 = NaN.""" + sheets = store.sheets_for_box(box, epsg) + ix0, iy0 = math.ceil(box[0] / cell_m - 0.5), math.ceil(box[1] / cell_m - 0.5) + ix1, iy1 = math.floor(box[2] / cell_m - 0.5) + 1, math.floor(box[3] / cell_m - 0.5) + 1 + shape = (max(iy1 - iy0, 0), max(ix1 - ix0, 0)) + out = {k: np.full(shape, np.nan) for k in ("slope", "a", "si")} + total = max(len(sheets), 1) + for i, no in enumerate(sheets): + + def inner(pct: int, stage: str, i: int = i) -> None: + if report is not None: + report((i * 100 + pct) // total, f"{stage} ({i + 1}/{total})") + + path = ensure_sheet(no, epsg, cell_m, defaults, inner) + if path is None: + continue + with np.load(path) as data: + px, py = int(data["ix0"]), int(data["iy0"]) + pieces = {"slope": data["slope"], "a": data["a"]} + with np.load(_paths(no, epsg, cell_m)["si"]) as data: + pieces["si"] = data["si"] + h, w = pieces["slope"].shape + r0, c0, r1, c1 = max(py, iy0), max(px, ix0), min(py + h, iy1), min(px + w, ix1) + if r0 >= r1 or c0 >= c1: + continue + for k, v in pieces.items(): + src = v[r0 - py : r1 - py, c0 - px : c1 - px] + np.copyto( + out[k][r0 - iy0 : r1 - iy0, c0 - ix0 : c1 - ix0], src, where=~np.isnan(src) + ) # SI ∞ = 평지 + return {"x0": (ix0 + 0.5) * cell_m, "y0": (iy0 + 0.5) * cell_m, "sheets": sheets, **out} diff --git a/resources/tester/test_536_2_e07_area.py b/resources/tester/test_536_2_e07_area.py index 947be7abb..6a439b117 100644 --- a/resources/tester/test_536_2_e07_area.py +++ b/resources/tester/test_536_2_e07_area.py @@ -11,6 +11,7 @@ import math import numpy as np import pytest +from test_610_1_e07_sheet import fake_store from E07_StabArea import E07_StabArea_Engine as engine from E07_StabArea import E07_StabArea_Router as router_mod @@ -178,19 +179,12 @@ def _decode(text: str, dtype: str, grid: dict) -> np.ndarray: @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, - ) + + def slope30(x, y): # 동쪽 끝(x > 180) = 지형 없음 + inside = (x >= 0) & (x <= 180) & (y >= 0) & (y <= 100) + return np.where(inside, x * math.tan(math.radians(30)), np.nan) + + fake_store(monkeypatch, tmp_path, slope30) # 공용 지형 창고(PLAN 610-1) 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)) @@ -200,7 +194,6 @@ def project(tmp_path, monkeypatch): 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 diff --git a/resources/tester/test_610_1_e07_sheet.py b/resources/tester/test_610_1_e07_sheet.py new file mode 100644 index 000000000..f2d535f0d --- /dev/null +++ b/resources/tester/test_610_1_e07_sheet.py @@ -0,0 +1,138 @@ +"""PLAN 610-1 — E07 광역 = 공용 지형 창고 · 도엽마다 경사 · 집수면적 · 기본 물성 SI 미리. + +가짜 창고(정사각 도엽 · 함수 지형)로 — 도엽 이어 붙임 = 통째 셈 · 위쪽 유역 안 잘림 · +두 번째 = 흐름 안 셈 · 물성 바꿈 = 흐름 안 셈 · 지형 열쇠가 바뀌면 다시 셈. +""" + +from __future__ import annotations + +import json +import math +from dataclasses import replace + +import numpy as np +import pytest + +from common_util import common_util_terrain_store as store +from E07_StabArea import E07_StabArea_Engine as engine +from E07_StabArea import E07_StabArea_Sheet as sheet + +PARAMS = engine.AreaParams(gamma_t=18.0, gamma_sat=20.0, c=2.0, c_hi=7.0, phi=30.0, phi_hi=33.0) + + +def _grid(box, cell): + ix0, iy0 = math.ceil(box[0] / cell - 0.5), math.ceil(box[1] / cell - 0.5) + ix1, iy1 = math.floor(box[2] / cell - 0.5) + 1, math.floor(box[3] / cell - 0.5) + 1 + return ix0, iy0, (np.arange(ix0, ix1) + 0.5) * cell, (np.arange(iy0, iy1) + 0.5) * cell + + +def fake_store(monkeypatch, tmp_path, zfun, side_m=100.0): + """정사각 도엽(side_m · 번호 'i_j') · 지형 = zfun(x, y)(NaN = 없음) · 셈 횟수를 돌려줌.""" + monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path / "terrain") + calls = {"dem": 0} + + def sheets_for_box(box, epsg): + i0, i1 = math.floor(box[0] / side_m), math.ceil(box[2] / side_m) + j0, j1 = math.floor(box[1] / side_m), math.ceil(box[3] / side_m) + return [f"{i}_{j}" for j in range(j0, j1) for i in range(i0, i1)] + + def ensure_sheet_terrain(no, epsg, cell_m=2.0, on_progress=None): + i, j = (int(v) for v in no.split("_")) + box = (i * side_m, j * side_m, (i + 1) * side_m, (j + 1) * side_m) + ix0, iy0, xs, ys = _grid(box, cell_m) + z = zfun(*np.meshgrid(xs, ys)).astype(np.float32) + if not np.isfinite(z).any(): + return None + path = store.TERRAIN_ROOT / no / f"dem_{epsg}_{cell_m:g}m.npz" + if not path.is_file(): + calls["dem"] += 1 + path.parent.mkdir(parents=True, exist_ok=True) + np.savez(path, z=z, ix0=np.int64(ix0), iy0=np.int64(iy0)) + (path.parent / f"dem_{epsg}_{cell_m:g}m.key.json").write_text(json.dumps({"algo": 1})) + return path + + def load_terrain(epsg, *, sheets=None, box=None, cell_m=2.0, on_progress=None): + _, _, xs, ys = _grid(box, cell_m) + z = zfun(*np.meshgrid(xs, ys)).astype(np.float32) + return {"x": xs, "y": ys, "z": z, "valid_mask": np.isfinite(z), "failed": {}} + + monkeypatch.setattr(store, "sheets_for_box", sheets_for_box) + monkeypatch.setattr(store, "ensure_sheet_terrain", ensure_sheet_terrain) + monkeypatch.setattr(store, "load_terrain", load_terrain) + return calls + + +def _hills(x, y): + """300 × 300 m 언덕(밖 = 없음).""" + z = 0.3 * x + 0.1 * y + 12 * np.sin(x / 37) * np.cos(y / 29) + return np.where((x > 0) & (x < 300) & (y > 0) & (y < 300), z, np.nan) + + +@pytest.fixture() +def counted(monkeypatch): + seen = {"n": 0} + real = engine.catchment_area + + def counting(z, cell): + seen["n"] += 1 + return real(z, cell) + + monkeypatch.setattr(engine, "catchment_area", counting) + return seen + + +def test_도엽_이어붙임_통째와_같음(monkeypatch, tmp_path): + fake_store(monkeypatch, tmp_path, _hills) + got = sheet.load_area(5186, (0.0, 0.0, 300.0, 300.0), 2.0, PARAMS) + _, _, xs, ys = _grid((0.0, 0.0, 300.0, 300.0), 2.0) + z = _hills(*np.meshgrid(xs, ys)) + assert (got["x0"], got["y0"]) == (1.0, 1.0) + np.testing.assert_allclose(got["slope"], engine.slope_rad(z, 2.0), atol=1e-5) # float32 저장 + np.testing.assert_allclose(got["a"], engine.catchment_area(z, 2.0), rtol=1e-5) + np.testing.assert_allclose( + got["si"], engine.stability_index(got["slope"], got["a"], PARAMS), rtol=1e-6 + ) + + +def test_위쪽_유역_안_잘림(monkeypatch, tmp_path): + """노선 띠만 지형(옛 LAS 복사본)이면 띠 위쪽 끝 a 가 작게 나옴 — 창고(여유 400 m)는 안 잘림.""" + fake_store(monkeypatch, tmp_path, _hills) + band = (100.0, 100.0, 200.0, 200.0) + got = sheet.load_area(5186, band, 2.0, PARAMS) + _, _, xs, ys = _grid(band, 2.0) + cut = engine.catchment_area(_hills(*np.meshgrid(xs, ys)), 2.0) + # 물은 서쪽으로 — 띠 동쪽 끝(위쪽) 칸은 띠만이면 a 가 거의 자기 칸(2 m) + east = np.s_[:, -3] + assert np.median(cut[east]) < 10 and np.median(got["a"][east]) > 5 * np.median(cut[east]) + + +def test_두번째_물성바꿈은_흐름_안_셈(monkeypatch, tmp_path, counted): + calls = fake_store(monkeypatch, tmp_path, _hills) + box = (50.0, 50.0, 150.0, 150.0) + first = sheet.load_area(5186, box, 2.0, PARAMS) + assert counted["n"] == 4 and calls["dem"] == 9 # 도엽 넷 · 지형 = 언덕 아홉 장 + again = sheet.load_area(5186, box, 2.0, PARAMS) + assert counted["n"] == 4 + np.testing.assert_array_equal(first["a"], again["a"]) + # 마스터 기본 물성이 바뀜 → SI 만 다시(흐름 그대로) + other = replace(PARAMS, c=4.0) + moved = sheet.load_area(5186, box, 2.0, other) + assert counted["n"] == 4 + np.testing.assert_allclose( + moved["si"], engine.stability_index(moved["slope"], moved["a"], other), rtol=1e-6 + ) + # 미리 셈 SI 를 쓰나 칸 계산하나 결과 같음 + a = engine.stability(first["slope"], first["a"], PARAMS, first["si"]) + b = engine.stability(first["slope"], first["a"], PARAMS) + np.testing.assert_array_equal(a["cls"], b["cls"]) + # 이웃 지형 열쇠가 바뀌면 그 도엽 흐름 다시 + key = store.TERRAIN_ROOT / "2_2" / "dem_5186_2m.key.json" + key.write_text(json.dumps({"algo": 2})) + sheet.load_area(5186, box, 2.0, PARAMS) + assert counted["n"] == 8 + + +def test_빈_도엽(monkeypatch, tmp_path): + fake_store(monkeypatch, tmp_path, lambda x, y: np.full_like(x, np.nan)) + got = sheet.load_area(5186, (0.0, 0.0, 100.0, 100.0), 2.0, PARAMS) + assert np.isnan(got["a"]).all() and got["slope"].shape == (50, 50) diff --git a/ui_template/ui_template_locale_e07.ts b/ui_template/ui_template_locale_e07.ts index 577d387d2..ef3e6b10b 100644 --- a/ui_template/ui_template_locale_e07.ts +++ b/ui_template/ui_template_locale_e07.ts @@ -11,8 +11,8 @@ export const ui_locales_e07 = { /* --- E07 광역 위험 지도 --- */ E07_NoProject: ["프로젝트를 먼저 고르세요.", "Pick a project first."], E07_Hint: [ - "노선 띠 안 자연사면의 표층 붕괴 위험 — 무한사면 안전율과 SINMAP 안정 지수(SI)로 등급을 매깁니다. 지형 = 복사한 원지형.", - "Shallow-landslide hazard of natural slopes in the route band — infinite-slope Fs and SINMAP stability index (SI). Terrain = copied original DTM.", + "노선 띠 안 자연사면의 표층 붕괴 위험 — 무한사면 안전율과 SINMAP 안정 지수(SI)로 등급을 매깁니다. 지형 = 전국 등고선 공용 지형(2 m).", + "Shallow-landslide hazard of natural slopes in the route band — infinite-slope Fs and SINMAP stability index (SI). Terrain = shared national-contour terrain (2 m).", ], E07_ReadOnly: ["보기만 허용된 프로젝트입니다.", "This project is view-only."], E07_Card_Soil: ["토층 물성", "Soil layer"],