feat(E07 광역): 지형 = 공용 지형 창고 · 도엽마다 경사 · 집수면적 · 기본 물성 SI 미리 · 물성 바꿈 = 칸 계산만 (PLAN 610-1 · 요구 487 · 488)

- 새 E07_StabArea_Sheet — 도엽 지형 + 여유 400 m 로 경사 · D∞ 집수면적 → 창고 도엽 폴더 e07_<epsg>_<격자>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 <noreply@anthropic.com>
2026-10-30 22:13:07 +09:00
parent 5b7fc964a7
commit fc4f3ef1ae
6 changed files with 400 additions and 65 deletions
@@ -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
@@ -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)