diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_Terrain.py b/D03_PlanAnalysis/D03_PlanAnalysis_Terrain.py index 6827181ef..b01159626 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_Terrain.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_Terrain.py @@ -1,5 +1,5 @@ """D 예정노선 지형 받기 — 사업 전체 한 벌(도엽 · 위성 · DEM) (요구 345 · 389 · 463 · PLAN 530-1 · -561-2 · 595). +561-2 · 595 · 608-1). 지형은 D 프로젝트 storage 안 `B04_PreProcess/` 배치 그대로(지형 코드 주인 = B04 · 고칠 곳 0): @@ -10,7 +10,9 @@ - 병합 = 지난번 병합 도엽 ∪ 이번 도엽(옛 도엽 유지 · 열쇠 같으면 건너뜀) - 위성 = 범위 도엽 도곽 ∪ 옛 위성 범위(덮으면 건너뜀) - 지적도 = 지형 일에서 뺌 — 겹을 켤 때 GET 이 받음(`D03_PlanAnalysis_Layers.ensure_cadastral`) -- DEM = 한 벌 `models_plan/`(연도는 D03 나눔 뒤에 정해짐) · 받은 범위 = `plan_terrain.json` +- DEM = 공용 지형 창고(`common_util_terrain_store` · 도엽마다 처음 한 번 셈) 도엽을 받은 범위로 + 이어 붙여 한 벌 `models_plan/`(뒤 소비자가 읽는 꼴 그대로 · 열쇠 = 이어 붙인 값 지문) · + 도엽 zip = 보조 겹(하천 · 도로 등)만 · 받은 범위 = `plan_terrain.json` (상자 · 여유 · 좌표계 · 위치 칸 가운데) - 작업 좌표계 = 지점 epsg · `B03_FileInput/input/prj/plan.prj` 한 장(B04 · B05 좌표계 창구가 읽음) - 받기 전 범위 보기 = `plan_extent`(내려받기 없음) · 받기는 사용자가 누를 때만 @@ -20,6 +22,7 @@ from __future__ import annotations +import hashlib import json import logging import math @@ -39,9 +42,15 @@ from B04_PreProcess.B04_PreProcess_Engine_SheetStore import ( get_project_map_sheets_dir, get_sheet_path, ) -from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import build_sheet_surface_model +from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import ( + _key_path, + _method_model_entry, + _write_method_model, +) from B04_PreProcess.B04_PreProcess_Engine_VWorld import download_vworld_satellite_map +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec from common_util.common_util_json import atomic_write_json +from common_util.common_util_terrain_store import load_terrain from config.config_system_terrain import ( PLAN_TERRAIN_GRID_M, SHEET_SURFACE_DEFAULT_METHOD, @@ -296,35 +305,27 @@ def ensure_plan_terrain( before = {s for s in wanted if get_sheet_path(store, s) is not None} report(10, f"2/5 수치지형도 도엽 {len(wanted) - len(before)}매 받는 중") got = ensure_sheets(store, wanted) - if not got["available"]: - raise ValueError(f"수치지형도 도엽을 하나도 받지 못했습니다: {got['failed']}") + if not got["available"]: # 도엽 = 보조 겹만 · 지형은 창고가 셈 + logger.warning("D 수치지형도 도엽 받기 실패: %s", got["failed"]) report(45, "3/5 도엽 레이어 병합 중") merged = _merge_sheets(processed, store, list(got["available"]), keep_old=have is not None) report(60, "4/5 위성 지도 받는 중") ring = [(p[1], p[0]) for p in box_ring(box, epsg)] satellite = _satellite(processed, prj, ring, epsg, keep_old=have is not None) - report(65, f"5/5 사업 지형면 만드는 중 — 칸 {cells_of(box) / 1e4:.0f}만 · 수 분 걸릴 수 있음") + report(65, f"5/5 사업 지형면 — 칸 {cells_of(box) / 1e4:.0f}만") models_dir = models_dir_of(project_root) models_dir.mkdir(parents=True, exist_ok=True) - # B 절취 = 넘긴 꼭짓점 상자 + `SHEET_SURFACE_MARGIN_M` → 그만큼 안쪽으로 넘김(여유 ≥ 그 값) - inner = SHEET_SURFACE_MARGIN_M - outline = np.array( - [[box[0] + inner, box[1] + inner], [box[2] - inner, box[3] - inner]], dtype=np.float64 - ) - models = build_sheet_surface_model( - project_root, - processed, - models_dir, - outline, - f"EPSG:{int(epsg)}", - [SHEET_SURFACE_DEFAULT_METHOD], - PLAN_TERRAIN_GRID_M, - ) + + def store_report(percent: int, message: str) -> None: + report(65 + percent * 34 // 100, f"5/5 {message}") + + grid = load_terrain(epsg, box=box, cell_m=PLAN_TERRAIN_GRID_M, on_progress=store_report) + models = _write_store_model(project_root, models_dir, grid) if not models: raise ValueError( - "사업 지형면을 만들지 못했습니다 — 범위 안 도엽 등고선 높이 확인" - f"(병합 도엽 {len(merged)}매)" + "사업 지형면을 만들지 못했습니다 — 범위 안 등고선 확인" + f"(도엽 {len(grid['sheets'])}매 · 실패 {grid['failed'] or '없음'})" ) meta: dict[str, Any] = {"box": list(box), "margin_m": margin, "epsg": int(epsg)} kept = None if reset else received_meta(project_root).get("center") @@ -342,4 +343,40 @@ def ensure_plan_terrain( "satellite": satellite, "box_m": list(box), "models": models, + "terrain": {"computed": grid["computed"], "failed": grid["failed"]}, } + + +def _write_store_model( + project_root: Path, models_dir: Path, grid: dict[str, Any] +) -> list[dict[str, Any]]: + """창고 격자 → `dtm_sheet_<방법>.npz` · 미리보기 · 평활본(B 꼴) — 값 지문 같으면 그대로.""" + z, x, y = grid["z"], grid["x"], grid["y"] + if not z.size or not grid["valid_mask"].any(): + return [] + cell = grid["cell_m"] + method = SHEET_SURFACE_DEFAULT_METHOD + key = { + "source": "terrain_store", + "origin": [float(x[0]), float(y[0])], + "shape": list(z.shape), + "z": hashlib.sha1(np.ascontiguousarray(z).tobytes()).hexdigest(), + } + entry = _method_model_entry(project_root, models_dir, method, cell) + key_path = _key_path(models_dir, method) + files = [entry["model_file_path"], *(layer["file_path"] for layer in entry["layers"])] + if _read_json(key_path) == key and all((project_root / f).is_file() for f in files): + return [entry] + key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게 + spec = GridSpec( + x_min=float(x[0]) - cell / 2, + y_max=float(y[-1]) + cell / 2, + cell_m=cell, + n_rows=z.shape[0], + n_cols=z.shape[1], + ) + model = _write_method_model(project_root, models_dir, spec, z[::-1, :], method, cell) + if model is None: + return [] + atomic_write_json(key_path, key) + return [model] 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_530_1_plan_terrain.py b/resources/tester/test_530_1_plan_terrain.py index 74af4f6a1..68839f677 100644 --- a/resources/tester/test_530_1_plan_terrain.py +++ b/resources/tester/test_530_1_plan_terrain.py @@ -1,6 +1,6 @@ """D 예정노선 지형 받기(PLAN 530-1 · 561-2 · 595) — 사업 전체 한 벌 · 범위 보기 = 내려받기 없음 · 옛 도엽 유지 · 새 도엽만 · 범위는 줄지 않음 · 빠짐 없음 · 지적도 뺌 · 새로(reset) = 옛 범위 버림 · -받은 범위 좌표계 · 칸 상한 · 거리. +받은 범위 좌표계 · 칸 상한 · 거리 · DEM = 공용 지형 창고 이어 붙임(608-1 · 값 같으면 안 다시 씀). ./venv/Scripts/python.exe -m pytest resources/tester/test_530_1_plan_terrain.py -q """ @@ -8,6 +8,7 @@ import json from pathlib import Path +import numpy as np import pytest from pyproj import Transformer @@ -46,15 +47,22 @@ def fake(tmp_path, monkeypatch): (output_dir / "vworld_satellite_meta.json").write_text(json.dumps(meta), "utf-8") return output_dir / "vworld_satellite.png" - def dem(root, processed, models_dir, route_xy, crs, methods=None, grid_m=None): - seen["dem"].append((models_dir.name, route_xy.tolist(), crs)) - seen["grid_m"] = grid_m - return [{"model_file_path": f"B04_PreProcess/{models_dir.name}/dtm.npz"}] + def load_terrain(epsg, *, sheets=None, box=None, cell_m=2.0, on_progress=None): + seen["dem"].append(list(box)) + seen["grid_m"] = cell_m + on_progress(50, "지형 계산") + x = (np.arange(10) + np.floor(box[0] / cell_m) + 0.5) * cell_m # 작은 격자 · 원점 칸 배수 + y = (np.arange(8) + np.floor(box[1] / cell_m) + 0.5) * cell_m + z = seen["z"] + 0.1 * np.arange(80, dtype=np.float32).reshape(8, 10) + grid = {"x": x, "y": y, "z": z, "valid_mask": np.isfinite(z), "cell_m": cell_m} + return grid | {"epsg": epsg, "sheets": ["1"], "computed": [], "failed": {}} + + seen["z"] = np.float32(100.0) monkeypatch.setattr(store_mod, "download_sheet", download_sheet) monkeypatch.setattr(terrain, "parse_and_merge_sheets", merge) monkeypatch.setattr(terrain, "download_vworld_satellite_map", satellite) - monkeypatch.setattr(terrain, "build_sheet_surface_model", dem) + monkeypatch.setattr(terrain, "load_terrain", load_terrain) return tmp_path, seen @@ -63,10 +71,8 @@ def _zips(root): def _dem_box(seen) -> list[float]: - """가짜 DEM 이 받은 꼭짓점 상자 + B 절취 여유 = 실제 DEM 상자.""" - (x0, y0), (x1, y1) = seen["dem"][-1][1] - m = SHEET_SURFACE_MARGIN_M - return [x0 - m, y0 - m, x1 + m, y1 + m] + """창고에 물은 상자 = 받은 범위 그대로.""" + return seen["dem"][-1] def test_one_set_and_extend_keeps_old(fake) -> None: @@ -75,7 +81,13 @@ def test_one_set_and_extend_keeps_old(fake) -> None: old = _zips(root) assert sorted(seen["download"]) == old == sorted(first["sheets"]["downloaded"]) assert first["satellite"] == "done" and len(first["models"]) == 1 - assert seen["dem"][-1][0] == "models_plan" # 연도 없는 한 벌 + model = first["models"][0]["model_file_path"] + assert model == "B04_PreProcess/models_plan/dtm_sheet_laplace.npz" # 연도 없는 한 벌 + with np.load(root / model) as dem: # 뒤 소비자 꼴 — y 오름차순 · 칸 2 m + assert dem["z"].shape == (8, 10) and dem["y"][1] > dem["y"][0] + assert float(dem["resolution"][0]) == 2.0 and dem["valid_mask"].all() + assert dem["z"][1, 0] == pytest.approx(101.0) # 창고 값 그대로(뒤집힘 없음) + assert (root / model.replace(".npz", "_smooth.npz")).is_file() m = terrain.DEFAULT_MARGIN_M box = [min(A[0], B[0]) - m, min(A[1], B[1]) - m, max(A[0], B[0]) + m, max(A[1], B[1]) + m] assert _dem_box(seen) == pytest.approx(box) == first["box_m"] @@ -103,6 +115,16 @@ def test_one_set_and_extend_keeps_old(fake) -> None: assert seen["download"] == [] and third["satellite"] == "skip" assert len(seen["merge"]) == merges and third["box_m"] == grown + # 창고 값 같음 = DEM 안 다시 씀 · 값 바뀜(창고 다시 셈) = 다시 씀 + npz = root / third["models"][0]["model_file_path"] + stamp = npz.stat().st_mtime_ns + terrain.ensure_plan_terrain(root, [A, B], EPSG) + assert npz.stat().st_mtime_ns == stamp + seen["z"] = np.float32(200.0) + terrain.ensure_plan_terrain(root, [A, B], EPSG) + with np.load(npz) as dem: + assert float(dem["z"][0, 0]) == 200.0 + def test_reset_drops_old_box_sheets_satellite(fake) -> None: """595-2 — 새로 = 옛 받은 범위 · 병합 도엽 · 위성 범위 버림 · 좌표계 · 가운데 기록 · @@ -174,7 +196,8 @@ def test_box_sheets_leave_no_gap(fake) -> None: def test_empty_surface_raises_and_keeps_no_box(fake, monkeypatch) -> None: root, _ = fake - monkeypatch.setattr(terrain, "build_sheet_surface_model", lambda *a, **k: []) + seen = fake[1] + seen["z"] = np.float32(np.nan) # 창고에 등고선 없음(빈 도엽) with pytest.raises(ValueError, match="사업 지형면을 만들지 못했습니다"): terrain.ensure_plan_terrain(root, [A, B], EPSG) assert terrain.received_box(root) is None 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