- 새 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>
139 lines
5.9 KiB
Python
139 lines
5.9 KiB
Python
"""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)
|