diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py index 10b87ca10..8da3b6deb 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py @@ -1,37 +1,46 @@ -"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 요구 443). +"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 요구 443 · 486). 0차곡 = 뚜렷한 유로가 없는 곡두의 집수지형(`노선_배치.md:31` · 교본 `3. 임도망의 배치 절차.md:156`) → 사업 DEM(평활 안 한 `models_plan/dtm_sheet_<방법>.npz`)을 2 m 로 솎아 D8 흐름누적 · 위 집수 넓이가 `min_m2` 이상 · `max_m2` 미만인 칸 = 0차곡(그보다 크면 1차 계류) → 세선화 · 선 뽑기(B05 `_mask_to_polylines`) · 경위도 선 + 가닥 · 길이 합. 임계값 = 지식DB 수치 없음 → 정함(브레인) 0.1 ha · 1 ha · 바꿈 = 요청 `min_m2` · `max_m2`. -파생 저장본(열쇠 = 격자 지문 · 임계 · 칸). +파생 저장본(열쇠 = 격자 지문 · 임계 · 칸 · 쓴 도엽 누적 지문). +흐름누적 = 공용 지형 창고 도엽마다(도엽 + 여유 200 m · 창고에 있는 이웃 도엽만) 한 번 셈 → +`storage/_shared/terrain/<도엽>/zero_acc__<칸>m.npz`(열쇠 = 쓴 도엽 지형 열쇠) · +임계 바꿈 = 저장 누적에서 칸 고르기 · 선 뽑기만 · 사업 DEM 이 창고 격자가 아니면(옛 받기) 통째 셈. 계산 자리 = ③ 서버 단독. """ from __future__ import annotations +import json import math from pathlib import Path from typing import Any import numpy as np +from common_util import common_util_terrain_store as terrain_store +from common_util.common_util_atomic import atomic_write_npz from common_util.common_util_derived_store import ( derived_key, file_stamp, read_derived, write_derived, ) +from common_util.common_util_json import atomic_write_json from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD ZERO_ORDER_KIND = "plan_zero_order" -ZERO_ORDER_VERSION = 1 +ZERO_ORDER_VERSION = 2 ZERO_ORDER_MIN_M2 = 1_000.0 # 집수 시작(0.1 ha) ZERO_ORDER_MAX_M2 = 10_000.0 # 이 넓이부터 1차 계류(1 ha) ZERO_ORDER_CELL_M = 2.0 # 흐름누적 칸(1 m 격자는 솎음) ZERO_ORDER_MIN_LEN_M = 10.0 # 세선화 갈림목 토막 거름 +ZERO_ORDER_MARGIN_M = 200.0 # 도엽 누적 여유(통째 셈과 0차곡 칸 차 0 · 분석 486 3-1) +ACC_CAP = 65535 # 저장 누적 상한(칸 수 · 2 m 면 26 ha · 임계 1 ha 보다 넉넉) SIMPLIFY_M = 1.0 _D8 = ((-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)) @@ -103,6 +112,28 @@ def zero_order_mask( return valid & (area >= min_m2) & (area < max_m2) +def mask_lines( + acc: np.ndarray, + z: np.ndarray, + valid: np.ndarray, + x: np.ndarray, + y: np.ndarray, + cell_m: float, + min_m2: float = ZERO_ORDER_MIN_M2, + max_m2: float = ZERO_ORDER_MAX_M2, +) -> list[list[tuple[float, float]]]: + """흐름누적 → 0차곡 칸 → 선(작업 좌표계 m 점 목록 · `ZERO_ORDER_MIN_LEN_M` 미만 거름).""" + from B05_Profile.B05_Profile_Engine_Skeleton import _mask_to_polylines + + mask = zero_order_mask(acc, valid, cell_m, min_m2, max_m2) + lines = [[(p[0], p[1]) for p in line["polyline"]] for line in _mask_to_polylines(mask, x, y, z)] + return [ + line + for line in lines + if sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= ZERO_ORDER_MIN_LEN_M + ] + + def zero_order_lines( z: np.ndarray, valid: np.ndarray, @@ -112,16 +143,108 @@ def zero_order_lines( min_m2: float = ZERO_ORDER_MIN_M2, max_m2: float = ZERO_ORDER_MAX_M2, ) -> list[list[tuple[float, float]]]: - """격자(행 = 남 → 북) → 0차곡 선(작업 좌표계 m 점 목록 · `ZERO_ORDER_MIN_LEN_M` 미만 거름).""" - from B05_Profile.B05_Profile_Engine_Skeleton import _mask_to_polylines + """격자(행 = 남 → 북) 통째 흐름누적 → 0차곡 선.""" + return mask_lines(d8_accumulation(z, valid), z, valid, x, y, cell_m, min_m2, max_m2) - mask = zero_order_mask(d8_accumulation(z, valid), valid, cell_m, min_m2, max_m2) - lines = [[(p[0], p[1]) for p in line["polyline"]] for line in _mask_to_polylines(mask, x, y, z)] - return [ - line - for line in lines - if sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= ZERO_ORDER_MIN_LEN_M - ] + +def _read_json(path: Path) -> Any: + try: + return json.loads(path.read_text(encoding="utf-8")) + except (OSError, ValueError): + return None + + +def _acc_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: + folder = terrain_store.TERRAIN_ROOT / str(sheet_no) + stem = f"zero_acc_{int(epsg)}_{cell_m:g}m" + return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock" + + +def _dem_key(sheet_no: str, epsg: int, cell_m: float) -> dict[str, Any] | None: + """창고 지형 열쇠(셈 시간 뺌) — 안 셌거나 원천이 바뀌었으면 None.""" + stored = terrain_store._stored_ok(sheet_no, epsg, cell_m) + return None if stored is None else {k: v for k, v in stored.items() if k != "seconds"} + + +def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL_M) -> Path | None: + """도엽 하나의 흐름누적(창고) — 도엽 + 여유 200 m 지형(창고에 있는 도엽만 · 없는 이웃 = + 출구)으로 셈 · 도엽 칸만 저장(uint16 · `ACC_CAP` 로 자름) · 열쇠 = 쓴 도엽 지형 열쇠 · + 지형이 없으면 None.""" + own = _dem_key(sheet_no, epsg, cell_m) + dem_npz = terrain_store._names(sheet_no, epsg, cell_m)[0] + if own is None or own.get("empty") or not dem_npz.is_file(): + return None + with np.load(dem_npz) as data: + ix0, iy0 = int(data["ix0"]), int(data["iy0"]) + rows, cols = data["z"].shape + pad = math.ceil(ZERO_ORDER_MARGIN_M / cell_m) + box = ( + (ix0 - pad) * cell_m, + (iy0 - pad) * cell_m, + (ix0 + cols + pad) * cell_m, + (iy0 + rows + pad) * cell_m, + ) + used = {} + for no in terrain_store.sheets_for_box(box, epsg): + dem = _dem_key(no, epsg, cell_m) + if dem is not None and not dem.get("empty"): + used[no] = dem + key = {"version": ZERO_ORDER_VERSION, "margin_m": ZERO_ORDER_MARGIN_M, "dem": used} + npz, key_path, lock_path = _acc_names(sheet_no, epsg, cell_m) + if npz.is_file() and _read_json(key_path) == key: + return npz + with terrain_store._Lock(lock_path): + if npz.is_file() and _read_json(key_path) == key: # 기다리는 동안 남이 셌을 수 있음 + return npz + grid = terrain_store.load_terrain(epsg, sheets=sorted(used), box=box, cell_m=cell_m) + acc = d8_accumulation(grid["z"], grid["valid_mask"])[pad : pad + rows, pad : pad + cols] + key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게 + atomic_write_npz( + npz, + acc=np.minimum(acc, ACC_CAP).astype(np.uint16), + ix0=np.int64(ix0), + iy0=np.int64(iy0), + cell_m=np.float64(cell_m), + ) + atomic_write_json(key_path, key) + return npz + + +def store_accumulation_paths( + x: np.ndarray, y: np.ndarray, epsg: int, cell_m: float +) -> list[Path] | None: + """사업 DEM 격자(창고 격자 위)를 덮는 도엽 누적 — 창고 격자 아님 · 지형 안 센 도엽 = None.""" + fx, fy = x[0] / cell_m - 0.5, y[0] / cell_m - 0.5 + if cell_m != ZERO_ORDER_CELL_M or abs(fx - round(fx)) > 1e-6 or abs(fy - round(fy)) > 1e-6: + return None + half = cell_m / 2 + box = (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half) + paths: list[Path] = [] + for no in terrain_store.sheets_for_box(box, epsg): + dem = _dem_key(no, epsg, cell_m) + if dem is None: + return None + path = None if dem.get("empty") else sheet_accumulation(no, epsg, cell_m) + if path is not None: + paths.append(path) + return paths + + +def paste_accumulation( + paths: list[Path], x: np.ndarray, y: np.ndarray, cell_m: float +) -> np.ndarray: + """도엽 누적을 사업 격자(y 오름차순)에 이어 붙임 — 없는 칸 0.""" + ix0, iy0 = round(x[0] / cell_m - 0.5), round(y[0] / cell_m - 0.5) + acc = np.zeros((len(y), len(x))) + for path in paths: + with np.load(path) as data: + px, py, part = int(data["ix0"]), int(data["iy0"]), data["acc"] + r0, c0 = max(py, iy0), max(px, ix0) + r1, c1 = min(py + part.shape[0], iy0 + len(y)), min(px + part.shape[1], ix0 + len(x)) + if r0 < r1 and c0 < c1: + src = part[r0 - py : r1 - py, c0 - px : c1 - px] + np.copyto(acc[r0 - iy0 : r1 - iy0, c0 - ix0 : c1 - ix0], src, where=src > 0) + return acc def zero_order_layer( @@ -131,7 +254,8 @@ def zero_order_layer( min_m2: float = ZERO_ORDER_MIN_M2, max_m2: float = ZERO_ORDER_MAX_M2, ) -> dict[str, Any]: - """0차곡 겹 — 지문(격자 파일 · 임계 · 칸) 같으면 저장본 · 경위도 선 + 가닥 · 길이 합 m.""" + """0차곡 겹 — 지문(격자 파일 · 임계 · 칸 · 도엽 누적) 같으면 저장본 · 경위도 선 + 가닥 · + 길이 합 m.""" from shapely.geometry import LineString, mapping from shapely.ops import transform @@ -141,6 +265,12 @@ def zero_order_layer( stamp = file_stamp(npz) if stamp is None: raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.") + with np.load(npz) as dem: + x, y = np.asarray(dem["x"], float), np.asarray(dem["y"], float) + cell = float(dem["resolution"][0]) if "resolution" in dem.files else float(x[1] - x[0]) + step = max(1, round(ZERO_ORDER_CELL_M / cell)) + x, y, cell_m = x[::step], y[::step], cell * step + paths = store_accumulation_paths(x, y, epsg, cell_m) key = derived_key( ZERO_ORDER_KIND, ZERO_ORDER_VERSION, @@ -152,23 +282,25 @@ def zero_order_layer( ZERO_ORDER_CELL_M, ZERO_ORDER_MIN_LEN_M, SIMPLIFY_M, + None if paths is None else [f"{p.parent.name}:{file_stamp(p)}" for p in paths], ) saved = read_derived(project_root, ZERO_ORDER_KIND, key) if saved is not None: return saved["value"] with np.load(npz) as dem: - x, y = np.asarray(dem["x"], float), np.asarray(dem["y"], float) - cell = float(dem["resolution"][0]) if "resolution" in dem.files else float(x[1] - x[0]) - step = max(1, round(ZERO_ORDER_CELL_M / cell)) z = np.asarray(dem["z"], float)[::step, ::step] valid = np.asarray(dem["valid_mask"], bool)[::step, ::step] - lines = zero_order_lines(z, valid, x[::step], y[::step], cell * step, min_m2, max_m2) + if paths is None: + lines = zero_order_lines(z, valid, x, y, cell_m, min_m2, max_m2) + else: + acc = paste_accumulation(paths, x, y, cell_m) + lines = mask_lines(acc, z, valid, x, y, cell_m, min_m2, max_m2) _, to_ll = _transforms(epsg) shapes = [LineString(line).simplify(SIMPLIFY_M) for line in lines] value = { "min_m2": float(min_m2), "max_m2": float(max_m2), - "cell_m": cell * step, + "cell_m": cell_m, "count": len(shapes), "length_m": round(sum(s.length for s in shapes), 1), "collection": { diff --git a/resources/tester/test_608_plan_zero_order_store.py b/resources/tester/test_608_plan_zero_order_store.py new file mode 100644 index 000000000..c399e4b0a --- /dev/null +++ b/resources/tester/test_608_plan_zero_order_store.py @@ -0,0 +1,104 @@ +"""D 0차곡 흐름누적 = 도엽별 공용 창고(PLAN 608-2) — 도엽 + 여유 200 m 조각 = 통째 셈 · +처음만 셈 · 이웃 도엽 지형이 새로 생기면 다시 · 임계 바꿈은 누적 그대로 · 창고 격자 아닌 DEM = 통째. + +./venv/Scripts/python.exe -m pytest resources/tester/test_608_plan_zero_order_store.py -q +""" + +import numpy as np +import pytest + +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 D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zo + +EPSG = 5187 +WEST, EAST = "36901074", "36901075" +CELL = 20.0 # 시험은 굵은 칸(여유 200 m = 10 칸) + + +def _height(x, y): + """기운 면 + 골 · 마루(파장 240 m) — 웅덩이 없음 · 작은 집수.""" + return 0.08 * x + 6.0 * np.sin(y / 38.0) + 2.0 * np.sin(x / 55.0) + + +def _put_sheet(no: str) -> None: + ix0, iy0, own = store._owned(no, EPSG, CELL) + rows, cols = own.shape + gx, gy = np.meshgrid((np.arange(cols) + ix0 + 0.5) * CELL, (np.arange(rows) + iy0 + 0.5) * CELL) + z = _height(gx - 300000.0, gy).astype(np.float32) + z[~own] = np.nan + npz, key_path, _ = store._names(no, EPSG, CELL) + atomic_write_npz(npz, z=z, ix0=np.int64(ix0), iy0=np.int64(iy0), cell_m=np.float64(CELL)) + atomic_write_json(key_path, {**store._key(EPSG, CELL, {}), "source": "gpkg", "seconds": 1}) + + +@pytest.fixture +def shared(tmp_path, monkeypatch): + monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path / "terrain") + calls = [] + real = zo.d8_accumulation + monkeypatch.setattr(zo, "d8_accumulation", lambda *a: calls.append(1) or real(*a)) + return calls + + +def _acc(no: str): + with np.load(zo.sheet_accumulation(no, EPSG, CELL)) as data: + return int(data["ix0"]), int(data["iy0"]), np.asarray(data["acc"], float) + + +def test_tiles_equal_whole_and_saved(shared) -> None: + calls = shared + _put_sheet(WEST) + assert zo.sheet_accumulation(EAST, EPSG, CELL) is None # 지형 안 센 도엽 + lone = _acc(WEST)[2] + assert len(calls) == 1 + _acc(WEST) + assert len(calls) == 1 # 두 번째 = 저장본 + _put_sheet(EAST) # 이웃 지형이 생김 = 다시 셈 + pieces = [_acc(WEST), _acc(EAST)] + assert len(calls) == 3 + + whole = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL) + ix0, iy0 = round(whole["x"][0] / CELL - 0.5), round(whole["y"][0] / CELL - 0.5) + acc = zo.d8_accumulation(whole["z"], whole["valid_mask"]) + tiled = np.zeros_like(acc) + for px, py, part in pieces: + dst = tiled[py - iy0 : py - iy0 + part.shape[0], px - ix0 : px - ix0 + part.shape[1]] + np.copyto(dst, part, where=part > 0) + valid = whole["valid_mask"] + cells = CELL * CELL + for low, high in ((2, 8), (1, 10)): # 0차곡 칸 = 집수 칸 수 [low, high) + a = zo.zero_order_mask(acc, valid, CELL, low * cells, high * cells) + b = zo.zero_order_mask(tiled, valid, CELL, low * cells, high * cells) + assert a.any() and np.array_equal(a, b) + # 이웃 없이 셌던 서쪽 도엽은 동쪽 경계에서 달랐음(이웃이 들어와 고쳐짐) + assert not np.array_equal(lone, pieces[0][2]) + + +def test_layer_uses_store_and_threshold_reuses(shared, tmp_path, monkeypatch) -> None: + """사업 DEM 이 창고 격자 위면 도엽 누적을 이어 붙임 · 임계 바꿈 = 누적 다시 안 셈.""" + calls = shared + monkeypatch.setattr(zo, "ZERO_ORDER_CELL_M", CELL) + _put_sheet(WEST) + _put_sheet(EAST) + grid = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL) + x, y, z = grid["x"][5:-5], grid["y"][5:-5], grid["z"][5:-5, 5:-5] + models = tmp_path / "models" + models.mkdir() + np.savez( + models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz", + x=x, + y=y, + z=np.nan_to_num(z, nan=0.0), + valid_mask=np.isfinite(z), + resolution=np.array([CELL], np.float32), + ) + cells = CELL * CELL + first = zo.zero_order_layer(tmp_path, models, EPSG, 2 * cells, 8 * cells) + assert first["count"] >= 1 and first["cell_m"] == CELL + assert len(calls) == 2 # 도엽 둘 한 번씩 + zo.zero_order_layer(tmp_path, models, EPSG, 1 * cells, 10 * cells) + assert len(calls) == 2 # 임계 바꿈 = 저장 누적에서 칸 고르기만 + # 창고 격자에서 반 칸 어긋난 DEM = 통째 셈 + assert zo.store_accumulation_paths(x + CELL / 2, y, EPSG, CELL) is None