diff --git a/common_util/common_util_terrain_store.py b/common_util/common_util_terrain_store.py new file mode 100644 index 000000000..dad1b3263 --- /dev/null +++ b/common_util/common_util_terrain_store.py @@ -0,0 +1,417 @@ +"""공용 지형 창고 — 도엽(1:5,000)마다 2 m 지형을 처음 쓸 때 한 번 셈 · 영구 저장 +(PLAN 607-2 · 607-3). + +자리 = `storage/_shared/terrain/<도엽>/dem__<격자>m.npz` + `.key.json` +(모든 작업 폴더가 보는 정션 · 시놀로지로 데스크탑에도). +원천 = 전국 등고선 파일(EPSG:5179) · 그 도엽에 등고선이 없으면(섬 등) +그때만 도엽 zip 을 받아 같은 창고 폴더에(B04 받기 함수 · `prune_sheets` 안 부름). +도엽 + 여유 200 m 등고선으로 셈 → 도엽 안쪽 칸(칸 중심 위경도가 그 도엽)만 저장 → 이어 붙이면 +통째 셈과 같음. 격자 원점 = 좌표 격자 배수(칸 모서리 = k·cell_m). +열쇠 = 계산 방식 판 + 여유 + 전국 파일 지문(크기 · 수정 시각) + 쓴 도엽 zip 해시. +""" + +from __future__ import annotations + +import hashlib +import json +import logging +import math +import os +import time +from pathlib import Path +from typing import Any, Callable + +import numpy as np +from pyproj import Transformer + +from B04_PreProcess.B04_PreProcess_Engine_MapSheet import ( + SHEET5K_SIZE_DEG, + latlon_to_sheet5k, + neighbors_3x3, + sheet5k_to_bounds, +) +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ELEVATION_KEYS, GridSpec +from common_util.common_util_atomic import atomic_write_npz +from common_util.common_util_json import atomic_write_json +from common_util.common_util_terrain_surface import ( + TERRAIN_ALGO, + build_contour_surface, + segments_from_lines, +) +from config.config_system import PROJECT_ROOT, STORAGE_BASE_DIR + +logger = logging.getLogger(__name__) + +TERRAIN_ROOT = Path(STORAGE_BASE_DIR) / "_shared" / "terrain" +NATIONAL_GPKG = Path(PROJECT_ROOT) / "resources" / "data_global_contours" / "national_contours.gpkg" +TERRAIN_MARGIN_M = 200.0 +# ponytail: 잠금 파일이 이보다 오래되면 죽은 셈으로 보고 깸 · 도엽 하나 셈 1~5 분 +LOCK_STALE_S = 15 * 60 + +Progress = Callable[[int, str], None] | None +Box = tuple[float, float, float, float] +Segments = tuple[np.ndarray, np.ndarray, np.ndarray] +_NO_SEGMENTS: Segments = (np.zeros((0, 2)), np.zeros((0, 2)), np.zeros(0)) + + +def _names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: + folder = TERRAIN_ROOT / str(sheet_no) + stem = f"dem_{int(epsg)}_{cell_m:g}m" + return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock" + + +def _gpkg_stamp() -> str | None: + try: + stat = NATIONAL_GPKG.stat() + except OSError: + return None + return f"{stat.st_size}:{stat.st_mtime_ns}" + + +def _sha256(path: Path) -> str | None: + try: + return hashlib.sha256(path.read_bytes()).hexdigest() + except OSError: + return None + + +def _key(epsg: int, cell_m: float, zips: dict[str, str | None]) -> dict[str, Any]: + return { + "algo": TERRAIN_ALGO, + "margin_m": TERRAIN_MARGIN_M, + "epsg": int(epsg), + "cell_m": float(cell_m), + "gpkg": _gpkg_stamp(), + "zips": zips, + } + + +def _stored_ok(sheet_no: str, epsg: int, cell_m: float) -> dict[str, Any] | None: + """저장 열쇠가 지금 원천과 같으면 그 열쇠 — 아니면 None(다시 셈).""" + npz, key_path, _ = _names(sheet_no, epsg, cell_m) + try: + stored = json.loads(key_path.read_text(encoding="utf-8")) + except (OSError, ValueError): + return None + if not stored.get("empty") and not npz.is_file(): + return None + zips = {s: _sha256(TERRAIN_ROOT / s / f"{s}.zip") for s in stored.get("zips") or {}} + now = _key(epsg, cell_m, zips) + if now["gpkg"] is None: # 전국 파일 없는 PC — 있는 저장본을 믿음(시놀로지 덮어쓰기 핑퐁 막음) + now["gpkg"] = stored.get("gpkg") + stored_key = {k: stored.get(k) for k in now} + return stored if stored_key == now else None + + +class _Lock: + """도엽 · 좌표계 · 격자 하나를 한 셈만 — 잠금 파일(서버 여럿 · 스레드 여럿 공통).""" + + def __init__(self, path: Path) -> None: + self.path = path + + def __enter__(self) -> "_Lock": + self.path.parent.mkdir(parents=True, exist_ok=True) + while True: + try: + fd = os.open(self.path, os.O_CREAT | os.O_EXCL | os.O_WRONLY) + os.write(fd, f"{os.getpid()} {time.time():.0f}".encode()) + os.close(fd) + return self + except FileExistsError: + try: + if time.time() - self.path.stat().st_mtime > LOCK_STALE_S: + self.path.unlink(missing_ok=True) + continue + except OSError: + continue + time.sleep(0.3) + + def __exit__(self, *exc: Any) -> None: + self.path.unlink(missing_ok=True) + + +def sheets_for_box(box: Box, epsg: int) -> list[str]: + """TM 상자(x0, y0, x1, y1)에 닿는 도엽 번호(북서부터).""" + x0, y0, x1, y1 = box + t = np.linspace(0.0, 1.0, 21) + xs = np.concatenate([x0 + (x1 - x0) * t, np.full(21, x1), x1 - (x1 - x0) * t, np.full(21, x0)]) + ys = np.concatenate([np.full(21, y0), y0 + (y1 - y0) * t, np.full(21, y1), y1 - (y1 - y0) * t]) + lon, lat = Transformer.from_crs(f"EPSG:{int(epsg)}", "EPSG:4326", always_xy=True).transform( + xs, ys + ) + step = SHEET5K_SIZE_DEG + sheets: list[str] = [] + for i in range(math.floor(lat.max() / step), math.floor(lat.min() / step) - 1, -1): + for j in range(math.floor(lon.min() / step), math.floor(lon.max() / step) + 1): + no = latlon_to_sheet5k((i + 0.5) * step, (j + 0.5) * step) + if no not in sheets: + sheets.append(no) + return sheets + + +def _owned(sheet_no: str, epsg: int, cell_m: float) -> tuple[int, int, np.ndarray]: + """도엽 안쪽 칸 — (첫 칸 열 번호 ix0 · 첫 칸 행 번호 iy0 · 소유 표 (ny, nx) y 오름차순).""" + lon0, lat0, lon1, lat1 = sheet5k_to_bounds(sheet_no) + to_tm = Transformer.from_crs("EPSG:4326", f"EPSG:{int(epsg)}", always_xy=True) + t = np.linspace(0.0, 1.0, 21) + lon = np.concatenate([lon0 + (lon1 - lon0) * t, np.full(21, lon1), lon0 + (lon1 - lon0) * t]) + lat = np.concatenate([np.full(21, lat0), lat0 + (lat1 - lat0) * t, np.full(21, lat1)]) + lon = np.concatenate([lon, np.full(21, lon0)]) + lat = np.concatenate([lat, lat0 + (lat1 - lat0) * t]) + xs, ys = to_tm.transform(lon, lat) + ix0, ix1 = math.floor(xs.min() / cell_m) - 1, math.ceil(xs.max() / cell_m) + 1 + iy0, iy1 = math.floor(ys.min() / cell_m) - 1, math.ceil(ys.max() / cell_m) + 1 + cx = (np.arange(ix0, ix1) + 0.5) * cell_m + cy = (np.arange(iy0, iy1) + 0.5) * cell_m + gx, gy = np.meshgrid(cx, cy) + clon, clat = Transformer.from_crs(f"EPSG:{int(epsg)}", "EPSG:4326", always_xy=True).transform( + gx, gy + ) + own = (clon >= lon0) & (clon < lon1) & (clat >= lat0) & (clat < lat1) + return ix0, iy0, own + + +def _gpkg_segments(epsg: int, box: Box) -> Segments: + """전국 등고선에서 TM 상자에 닿는 선분 — 끝 a · b (M,2) 사업 좌표 · 높이 z (M,).""" + if not NATIONAL_GPKG.is_file(): + return _NO_SEGMENTS + import pyogrio + import shapely + + x0, y0, x1, y1 = box + to_5179 = Transformer.from_crs(f"EPSG:{int(epsg)}", "EPSG:5179", always_xy=True) + t = np.linspace(0.0, 1.0, 21) + bx, by = to_5179.transform( + np.concatenate([x0 + (x1 - x0) * t, x0 + (x1 - x0) * t, np.full(21, x0), np.full(21, x1)]), + np.concatenate([np.full(21, y0), np.full(21, y1), y0 + (y1 - y0) * t, y0 + (y1 - y0) * t]), + ) + pad = 20.0 + lo_x, lo_y, hi_x, hi_y = bx.min() - pad, by.min() - pad, bx.max() + pad, by.max() + pad + frame = pyogrio.read_dataframe( + NATIONAL_GPKG, bbox=(lo_x, lo_y, hi_x, hi_y), columns=["CTRLN_HG"] + ) + if frame.empty: + return _NO_SEGMENTS + parts, owner = shapely.get_parts(frame.geometry.values, return_index=True) + levels = frame["CTRLN_HG"].to_numpy(dtype=np.float64)[owner] + coords, part = shapely.get_coordinates(parts, return_index=True) + # 상자 밖 선분은 버림(격자 밖 선분은 결과에 안 닿음 · 꼭짓점 수백만 변환을 줄임) + same = part[:-1] == part[1:] + a, b = coords[:-1], coords[1:] + hit = ( + same + & (np.maximum(a[:, 0], b[:, 0]) >= lo_x) + & (np.minimum(a[:, 0], b[:, 0]) <= hi_x) + & (np.maximum(a[:, 1], b[:, 1]) >= lo_y) + & (np.minimum(a[:, 1], b[:, 1]) <= hi_y) + ) + seg = np.flatnonzero(hit) + if not len(seg): + return _NO_SEGMENTS + to_tm = Transformer.from_crs("EPSG:5179", f"EPSG:{int(epsg)}", always_xy=True) + ax, ay = to_tm.transform(a[seg, 0], a[seg, 1]) + bx2, by2 = to_tm.transform(b[seg, 0], b[seg, 1]) + return np.column_stack([ax, ay]), np.column_stack([bx2, by2]), levels[part[seg]] + + +def _zip_lines(sheet_nos: list[str], epsg: int) -> list[tuple[np.ndarray, float]]: + """창고 폴더의 도엽 zip 등고선 → (꼭짓점 · 높이) 목록.""" + from B04_PreProcess.B04_PreProcess_Engine_SheetParser import SHEET_LAYERS, _read_layer_from_zip + from B04_PreProcess.B04_PreProcess_Engine_SheetStore import _load_index + from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import iter_linestrings + + prefix, code = SHEET_LAYERS["등고선"] + lines: list[tuple[np.ndarray, float]] = [] + for no in sheet_nos: + folder = TERRAIN_ROOT / no + entry = _load_index(folder)["sheets"].get(no) + if not entry: + continue + frame = _read_layer_from_zip(folder / entry["file"], prefix, code, entry["crs"]) + if frame is None: + continue + frame = frame.to_crs(f"EPSG:{int(epsg)}") + column = next((k for k in ELEVATION_KEYS if k in frame.columns), None) + if column is None: + continue + for geometry, level in zip(frame.geometry, frame[column]): + if geometry is None or level is None: + continue + for line in iter_linestrings(geometry): + lines.append((np.asarray(line.coords)[:, :2], float(level))) + return lines + + +def _fetch_zips(sheet_no: str, report: Callable[[int, str], None]) -> tuple[dict, str | None]: + """빠진 곳 — 도엽 + 이웃 8 장 zip 을 창고 폴더마다 받음. ({도엽: 해시}, 실패 까닭).""" + from B04_PreProcess.B04_PreProcess_Engine_SheetStore import ensure_sheets + + zips: dict[str, str | None] = {} + neighbors = neighbors_3x3(sheet_no) + for i, no in enumerate([sheet_no] + [n for n in neighbors if n != sheet_no]): + report(10 + 40 * i // 9, f"도엽 받기 {i + 1}/9 · {no}") + result = ensure_sheets(TERRAIN_ROOT / no, [no]) + if no in result["available"]: + zips[no] = _sha256(TERRAIN_ROOT / no / f"{no}.zip") + elif no == sheet_no: + reason = result["failed"].get(no, "") + # 검색 결과 없음 = 바다 등 진짜 빈 도엽 · 그 밖(로그인 · 망)은 다음에 다시 + return zips, None if reason.startswith("LookupError") else reason + return zips, None + + +def ensure_sheet_terrain( + sheet_no: str, epsg: int, cell_m: float = 2.0, on_progress: Progress = None +) -> Path | None: + """도엽 하나의 창고 지형 — 있으면 그대로 · 없거나 원천이 바뀌었으면 셈. 빈 도엽이면 None. + + 받기 실패(로그인 · 망)는 RuntimeError — 저장 안 함(다음에 다시). + """ + sheet_no = str(sheet_no) + + def report(pct: int, stage: str) -> None: + if on_progress is not None: + on_progress(pct, stage) + + npz, key_path, lock_path = _names(sheet_no, epsg, cell_m) + stored = _stored_ok(sheet_no, epsg, cell_m) + if stored is not None: + return None if stored.get("empty") else npz + with _Lock(lock_path): + stored = _stored_ok(sheet_no, epsg, cell_m) # 기다리는 동안 남이 셌을 수 있음 + if stored is not None: + return None if stored.get("empty") else npz + started = time.monotonic() + report(2, f"전국 등고선 읽기 · {sheet_no}") + ix0, iy0, own = _owned(sheet_no, epsg, cell_m) + pad = math.ceil(TERRAIN_MARGIN_M / cell_m) + rows, cols = own.shape + spec = GridSpec( + x_min=(ix0 - pad) * cell_m, + y_max=(iy0 + rows + pad) * cell_m, + cell_m=cell_m, + n_rows=rows + 2 * pad, + n_cols=cols + 2 * pad, + ) + region = ( + spec.x_min, + spec.y_max - spec.n_rows * cell_m, + spec.x_min + spec.n_cols * cell_m, + spec.y_max, + ) + a, b, z = _gpkg_segments(epsg, region) + x_lo, y_lo = ix0 * cell_m, iy0 * cell_m + x_hi, y_hi = (ix0 + cols) * cell_m, (iy0 + rows) * cell_m + has_own = bool( + np.any( + (np.maximum(a[:, 0], b[:, 0]) >= x_lo) + & (np.minimum(a[:, 0], b[:, 0]) <= x_hi) + & (np.maximum(a[:, 1], b[:, 1]) >= y_lo) + & (np.minimum(a[:, 1], b[:, 1]) <= y_hi) + ) + ) + zips: dict[str, str | None] = {} + if not has_own: + zips, failure = _fetch_zips(sheet_no, report) + if failure: + raise RuntimeError(f"도엽 {sheet_no} 받기 실패: {failure}") + za, zb, zz = segments_from_lines(_zip_lines(list(zips), epsg)) + a, b, z = np.vstack([a, za]), np.vstack([b, zb]), np.concatenate([z, zz]) + report(55, f"지형 계산 · {sheet_no}") + surface = build_contour_surface(spec, a, b, z) if len(z) else None + key = _key(epsg, cell_m, zips) + key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게 + if surface is None: + npz.unlink(missing_ok=True) + atomic_write_json(key_path, {**key, "empty": True}) + report(100, f"빈 도엽 · {sheet_no}") + return None + core = surface[pad : pad + rows, pad : pad + cols][::-1, :].astype(np.float32) + core[~own] = np.nan + report(90, f"저장 · {sheet_no}") + atomic_write_npz( + npz, z=core, ix0=np.int64(ix0), iy0=np.int64(iy0), cell_m=np.float64(cell_m) + ) + atomic_write_json( + key_path, + { + **key, + "source": "zip" if zips else "gpkg", + "seconds": round(time.monotonic() - started, 1), + }, + ) + report(100, f"완료 · {sheet_no}") + logger.info( + "공용 지형 %s (EPSG:%s · %gm) %.1fs", sheet_no, epsg, cell_m, time.monotonic() - started + ) + return npz + + +def load_terrain( + epsg: int, + *, + sheets: list[str] | None = None, + box: Box | None = None, + cell_m: float = 2.0, + on_progress: Progress = None, +) -> dict[str, Any]: + """도엽 목록 또는 TM 상자 → 이어 붙인 격자(x · y 오름차순 · z (ny, nx) · valid_mask) + 메타.""" + if sheets is None: + if box is None: + raise ValueError("sheets 또는 box 가 필요함") + sheets = sheets_for_box(box, epsg) + sheets = [str(s) for s in sheets] + pieces: list[tuple[int, int, np.ndarray]] = [] + computed: list[str] = [] + failed: dict[str, str] = {} + total = max(len(sheets), 1) + for i, no in enumerate(sheets): + fresh = _stored_ok(no, epsg, cell_m) is None + + def inner(pct: int, stage: str, i: int = i) -> None: + if on_progress is not None: + on_progress(int((i * 100 + pct) / total), f"{stage} ({i + 1}/{total})") + + try: + path = ensure_sheet_terrain(no, epsg, cell_m, inner if on_progress else None) + except Exception as exc: # noqa: BLE001 — 한 장이 죽어도 나머지는 이어 붙임 + failed[no] = f"{type(exc).__name__}: {exc}" + continue + if fresh: + computed.append(no) + if path is None: + continue + with np.load(path) as data: + pieces.append((int(data["ix0"]), int(data["iy0"]), data["z"])) + if box is not None: + 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 + elif pieces: + ix0 = min(p[0] for p in pieces) + iy0 = min(p[1] for p in pieces) + ix1 = max(p[0] + p[2].shape[1] for p in pieces) + iy1 = max(p[1] + p[2].shape[0] for p in pieces) + else: + ix0 = iy0 = ix1 = iy1 = 0 + z = np.full((max(iy1 - iy0, 0), max(ix1 - ix0, 0)), np.nan, dtype=np.float32) + for px, py, pz in pieces: + r0, c0 = max(py, iy0), max(px, ix0) + r1, c1 = min(py + pz.shape[0], iy1), min(px + pz.shape[1], ix1) + if r0 >= r1 or c0 >= c1: + continue + src = pz[r0 - py : r1 - py, c0 - px : c1 - px] + dst = z[r0 - iy0 : r1 - iy0, c0 - ix0 : c1 - ix0] + np.copyto(dst, src, where=np.isfinite(src)) + if on_progress is not None: + on_progress(100, "지형 준비 끝") + return { + "x": (np.arange(ix0, ix1) + 0.5) * cell_m, + "y": (np.arange(iy0, iy1) + 0.5) * cell_m, + "z": z, + "valid_mask": np.isfinite(z), + "cell_m": float(cell_m), + "epsg": int(epsg), + "sheets": sheets, + "computed": computed, + "failed": failed, + } diff --git a/common_util/common_util_terrain_surface.py b/common_util/common_util_terrain_surface.py new file mode 100644 index 000000000..18ac58cb9 --- /dev/null +++ b/common_util/common_util_terrain_surface.py @@ -0,0 +1,180 @@ +"""등고선 → 지형 격자(라플라스) — 공용 지형 창고가 쓰는 흠 고친 판(PLAN 607-1 · 요구 486~491). + +B04 `B04_PreProcess_Engine_SheetSurface.py` 라플라스 길을 본뜸 · B 파일은 안 고침 +(B 는 옛 판 그대로). +고친 흠 둘: + ① 굽기 — rasterio `all_touched=False` 는 선 방향 · 조각에 따라 다른 칸을 칠함(같은 선 뒤집어도 + 최대 9 m) → 선분 끝을 정렬한 뒤 큰 축 칸 중심마다 한 칸 + 꼭짓점 칸(방향 무관 · 굵기 1 칸). + ② 폐합 안쪽 채움 — 봉우리 · 웅덩이 방향 · 기울기를 같은 높이 선 전체가 아니라 그 고리 둘레로만 + → 범위가 바뀌어도 같은 고리는 같은 값(조각 = 통째). +보간 자체(거리 비례 + 라플라스 400 회)는 B04 `build_laplace` 를 그대로 부름. +""" + +from __future__ import annotations + +import logging +from typing import Any + +import numpy as np + +from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import build_laplace +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec + +logger = logging.getLogger(__name__) + +# 계산 방식 판 — 바뀌면 올림 · 창고 열쇠가 어긋나 저장본을 다시 셈. +TERRAIN_ALGO = 1 + +# 채움 기울기를 재는 고리 밖 띠(칸) — B04 와 같은 값. +_BAND_CELLS = 8 + + +def segments_from_lines(lines: list[tuple[np.ndarray, float]]) -> tuple[np.ndarray, ...]: + """(꼭짓점 (N,2) · 높이) 목록 → 선분 끝 a · b (M,2) · 높이 z (M,).""" + starts, ends, levels = [], [], [] + for coords, level in lines: + coords = np.asarray(coords, dtype=np.float64)[:, :2] + if len(coords) < 2: + continue + starts.append(coords[:-1]) + ends.append(coords[1:]) + levels.append(np.full(len(coords) - 1, float(level))) + if not starts: + empty = np.zeros((0, 2)) + return empty, empty, np.zeros(0) + return np.vstack(starts), np.vstack(ends), np.concatenate(levels) + + +def burn_segments(spec: GridSpec, a: np.ndarray, b: np.ndarray, z: np.ndarray) -> np.ndarray: + """선분을 격자에 굽는다 — 칸 = 그 위를 지나는 선의 높이 · 그 밖 NaN · 겹치면 낮은 쪽. + + 방향 무관: 선분 끝을 (x, y) 사전 차례로 맞춘 뒤 큰 축의 칸 중심마다 한 칸을 칠함 · + 꼭짓점 칸도 칠해 이음매 빈칸을 막음 · 선분마다 따로라 조각 · 이어 붙임에도 같음. + """ + burned = np.full((spec.n_rows, spec.n_cols), np.inf, dtype=np.float64) + if len(z): + swap = (a[:, 0] > b[:, 0]) | ((a[:, 0] == b[:, 0]) & (a[:, 1] > b[:, 1])) + p = np.where(swap[:, None], b, a) + q = np.where(swap[:, None], a, b) + # 연속 칸 좌표 — 칸 c 의 중심 = c + 0.5 + u0 = (p[:, 0] - spec.x_min) / spec.cell_m + u1 = (q[:, 0] - spec.x_min) / spec.cell_m + v0 = (spec.y_max - p[:, 1]) / spec.cell_m + v1 = (spec.y_max - q[:, 1]) / spec.cell_m + x_major = np.abs(u1 - u0) >= np.abs(v1 - v0) + rows, cols, levels = [], [], [] + for major, (s0, s1, t0, t1) in ((True, (u0, u1, v0, v1)), (False, (v0, v1, u0, u1))): + pick = x_major if major else ~x_major + s0, s1, t0, t1, zz = s0[pick], s1[pick], t0[pick], t1[pick], z[pick] + lo = np.ceil(np.minimum(s0, s1) - 0.5).astype(np.int64) + hi = np.floor(np.maximum(s0, s1) - 0.5).astype(np.int64) + count = np.clip(hi - lo + 1, 0, None) + seg = np.repeat(np.arange(len(lo)), count) + step = np.arange(int(count.sum())) - np.repeat(np.cumsum(count) - count, count) + major_idx = lo[seg] + step + span = s1[seg] - s0[seg] + t = (major_idx + 0.5 - s0[seg]) / np.where(span == 0, 1.0, span) + minor_idx = np.floor(t0[seg] + t * (t1[seg] - t0[seg])).astype(np.int64) + r, c = (minor_idx, major_idx) if major else (major_idx, minor_idx) + rows.append(r) + cols.append(c) + levels.append(zz[seg]) + for ends in (a, b): # 꼭짓점 칸 + rows.append(np.floor((spec.y_max - ends[:, 1]) / spec.cell_m).astype(np.int64)) + cols.append(np.floor((ends[:, 0] - spec.x_min) / spec.cell_m).astype(np.int64)) + levels.append(z) + r, c, zz = np.concatenate(rows), np.concatenate(cols), np.concatenate(levels) + inside = (r >= 0) & (r < spec.n_rows) & (c >= 0) & (c < spec.n_cols) + np.minimum.at(burned, (r[inside], c[inside]), zz[inside]) + burned[np.isinf(burned)] = np.nan + return burned.astype(np.float32) + + +def _window(component: np.ndarray, pad: int) -> tuple[slice, slice]: + rows = np.flatnonzero(component.any(axis=1)) + cols = np.flatnonzero(component.any(axis=0)) + return ( + slice(max(rows[0] - pad, 0), rows[-1] + pad + 1), + slice(max(cols[0] - pad, 0), cols[-1] + pad + 1), + ) + + +def resolve_enclosed_interiors( + burned: np.ndarray, present: list[float], surface: np.ndarray, cell_m: float, interval_m: float +) -> int: + """폐합 등고선 안쪽(봉우리 · 웅덩이)을 그 고리 바깥 사면 경사로 연장한다(in-place). + + B04 `_resolve_enclosed_interiors` 와 같은 규칙 · 다른 점 = 방향(바깥 평균)과 기울기 띠를 + 그 칸 묶음을 감싼 고리 칸 바깥으로만 잼(같은 높이 다른 선은 안 섞음). 채운 곳 수를 돌려줌. + """ + from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label + + constrained = np.isfinite(burned) + limit = max(interval_m - 0.5, 0.5) + pad = _BAND_CELLS + 2 + eight = np.ones((3, 3), dtype=bool) + resolved = 0 + for level in present: + mask = burned == level + interior = binary_fill_holes(mask) & ~mask + if not interior.any(): + continue + components, count = label(interior) + for component_id in range(1, count + 1): + whole = components == component_id + if (constrained & whole).any(): + continue # 안에 다른 제약이 있으면 마루가 아니다(보통의 감싸는 링) + win = _window(whole, pad) + component = whole[win] + line = mask[win] + local = surface[win] + # 그 묶음을 감싼 고리 칸 → 고리 바로 바깥 + loop = binary_dilation(component, structure=eight) & line + ring = binary_dilation(loop, structure=eight) & ~component & ~line + ring &= np.isfinite(local) + if not ring.any(): + continue + direction = 1.0 if float(local[ring].mean()) < level else -1.0 + inner = distance_transform_edt(component, sampling=cell_m) + outer_distance = distance_transform_edt(~(loop | component), sampling=cell_m) + band = (outer_distance > 0) & (outer_distance <= _BAND_CELLS * cell_m) + band &= np.isfinite(local) & ~component & ~line + if band.any(): + slope = float( + np.mean(np.abs(level - local[band].astype(np.float64)) / outer_distance[band]) + ) + else: + slope = 0.0 + if slope > 1e-3: + offset = np.minimum(slope * inner[component], limit) + else: + peak = float(inner.max()) + offset = (limit / 2.0) * (inner[component] / peak) if peak > 0 else 0.0 + local[component] = level + direction * offset + resolved += 1 + return resolved + + +def build_contour_surface( + spec: GridSpec, a: np.ndarray, b: np.ndarray, z: np.ndarray +) -> np.ndarray | None: + """선분 → (R, C) float32 표고(북쪽 첫 행 · 없음 NaN). 등고 두 단 미만이면 None.""" + burned = burn_segments(spec, a, b, z) + present = sorted(np.unique(burned[np.isfinite(burned)]).tolist()) + if len(present) < 2: + return None + interval_m = float(np.diff(np.array(present)).min()) + surface = build_laplace(spec, burned, None, spec.cell_m) + resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m) + return surface + + +def surface_grid(spec: GridSpec, surface: np.ndarray) -> dict[str, Any]: + """DtmGridSampler 꼴 — x · y 오름차순 · z (ny, nx) · valid_mask.""" + z = surface[::-1, :].astype(np.float32) + return { + "x": spec.cell_centers_x(), + "y": spec.cell_centers_y()[::-1], + "z": z, + "valid_mask": np.isfinite(z), + } diff --git a/resources/tester/test_607_terrain_store.py b/resources/tester/test_607_terrain_store.py new file mode 100644 index 000000000..a2973d9c9 --- /dev/null +++ b/resources/tester/test_607_terrain_store.py @@ -0,0 +1,222 @@ +"""공용 지형 창고 · 흠 고친 지형 계산(PLAN 607-1~3) — 굽기 방향 무관 · 고리 둘레 채움 · +조각 = 통째 · 처음만 셈 · 같이 두 번 = 한 번 · 원천 바뀌면 다시 · 빠진 곳 = 도엽 zip · +받기 실패는 저장 안 함. + +./venv/Scripts/python.exe -m pytest resources/tester/test_607_terrain_store.py -q +""" + +import threading + +import numpy as np +import pytest +from pyproj import Transformer + +from B04_PreProcess.B04_PreProcess_Engine_MapSheet import sheet5k_to_bounds +from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import _resolve_enclosed_interiors +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec +from common_util import common_util_terrain_store as store +from common_util import common_util_terrain_surface as surf + +EPSG = 5187 +WEST, EAST = "36901074", "36901075" + + +def _circle(cx, cy, r, n=48): + t = np.linspace(0.0, 2 * np.pi, n + 1) + return np.column_stack([cx + r * np.cos(t), cy + r * np.sin(t)]) + + +def _field_lines(x0, y0, hill): + """기운 면(대각 등고선 · 50 m 마다 5 m) + 봉우리 고리 둘 — (꼭짓점, 높이) 목록.""" + normal, along = np.array([0.8, 0.6]), np.array([-0.6, 0.8]) + t = np.array([[-8000.0], [8000.0]]) + lines = [ + (np.array([x0, y0]) + normal * (k * 50.0) + t * along, 500.0 + 5.0 * k) + for k in range(-40, 160) + ] + hx, hy = hill + lines.append((_circle(hx, hy, 60.0), 1000.0)) + lines.append((_circle(hx, hy, 30.0), 1005.0)) + return lines + + +def test_burn_ignores_direction_and_split(): + spec = GridSpec(x_min=0.0, y_max=200.0, cell_m=2.0, n_rows=100, n_cols=100) + rng = np.random.default_rng(7) + pts = np.cumsum(rng.normal(0, 3.3, (300, 2)), axis=0) + 100.0 + a, b, z = surf.segments_from_lines([(pts, 10.0), (_circle(80, 90, 41.3), 15.0)]) + one = surf.burn_segments(spec, a, b, z) + flipped = surf.burn_segments(spec, b[::-1], a[::-1], z[::-1]) + assert np.array_equal(np.isfinite(one), np.isfinite(flipped)) + assert np.array_equal(one[np.isfinite(one)], flipped[np.isfinite(flipped)]) + # 선을 두 조각으로 나눠도 같음 + a2, b2, z2 = surf.segments_from_lines([(pts[:150], 10.0), (pts[149:], 10.0)]) + c = surf.segments_from_lines([(_circle(80, 90, 41.3), 15.0)]) + split = surf.burn_segments( + spec, np.vstack([a2, c[0]]), np.vstack([b2, c[1]]), np.concatenate([z2, c[2]]) + ) + assert np.array_equal(np.nan_to_num(one, nan=-1), np.nan_to_num(split, nan=-1)) + # 굵기 1 칸 — 가로 선은 한 줄 + flat = surf.burn_segments(spec, np.array([[10.0, 51.0]]), np.array([[150.0, 51.0]]), np.ones(1)) + assert int(np.isfinite(flat).sum()) == 71 and np.isfinite(flat).any(axis=1).sum() == 1 + + +def test_fill_direction_follows_own_ring(): + """같은 높이 고리 둘 — 하나는 봉우리(바깥 낮음) · 하나는 웅덩이(바깥 높음).""" + spec = GridSpec(x_min=0.0, y_max=400.0, cell_m=2.0, n_rows=200, n_cols=300) + lines = [ + (_circle(100, 200, 80), 95.0), + (_circle(100, 200, 40), 100.0), # 봉우리: 95 안에 100 + (_circle(300, 200, 80), 105.0), + (_circle(300, 200, 40), 100.0), # 웅덩이: 105 안에 100 + ] + a, b, z = surf.segments_from_lines(lines) + out = surf.build_contour_surface(spec, a, b, z) + row = 100 + hill, pit = out[row, 50], out[row, 150] + assert hill > 100.5 and pit < 99.5 + # 옛 B 판은 같은 높이 선 전체 둘레로 재서 뒤집힘(이 꼴은 둘 다) + from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import build_laplace + + burned = surf.burn_segments(spec, a, b, z) + present = sorted(np.unique(burned[np.isfinite(burned)]).tolist()) + old = build_laplace(spec, burned, None, 2.0) + _resolve_enclosed_interiors(burned, present, old, 2.0, 5.0) + assert not (old[row, 50] > 100.5 and old[row, 150] < 99.5) + + +def _tm(lon, lat): + return Transformer.from_crs("EPSG:4326", f"EPSG:{EPSG}", always_xy=True).transform(lon, lat) + + +@pytest.fixture +def fake(tmp_path, monkeypatch): + lon0, lat0, lon1, lat1 = sheet5k_to_bounds(WEST) + x0, y0 = _tm(lon0 - 0.01, lat0 - 0.01) + edge_x, mid_y = _tm(lon1, (lat0 + lat1) / 2) + lines = _field_lines(x0, y0, (edge_x + 10.0, mid_y)) + allseg = surf.segments_from_lines(lines) + gpkg = tmp_path / "national.gpkg" + gpkg.write_bytes(b"v1") + calls = {"build": 0, "fetch": 0} + real_build = store.build_contour_surface + + def segments(epsg, box): + a, b, z = allseg + hit = (np.maximum(a[:, 0], b[:, 0]) >= box[0]) & (np.minimum(a[:, 0], b[:, 0]) <= box[2]) + hit &= (np.maximum(a[:, 1], b[:, 1]) >= box[1]) & (np.minimum(a[:, 1], b[:, 1]) <= box[3]) + return a[hit], b[hit], z[hit] + + def counted(*args): + calls["build"] += 1 + return real_build(*args) + + monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path / "terrain") + monkeypatch.setattr(store, "NATIONAL_GPKG", gpkg) + monkeypatch.setattr(store, "_gpkg_segments", segments) + monkeypatch.setattr(store, "build_contour_surface", counted) + return {"gpkg": gpkg, "calls": calls, "allseg": allseg} + + +def test_store_tiles_equal_whole_and_compute_once(fake): + cell = 10.0 + seen = [] + got = store.load_terrain( + EPSG, sheets=[WEST, EAST], cell_m=cell, on_progress=lambda p, s: seen.append(p) + ) + assert got["computed"] == [WEST, EAST] and not got["failed"] + assert seen[-1] == 100 and seen == sorted(seen) + # 통째 셈(같은 계산 · 넉넉한 상자)과 칸마다 같음 + ix0, iy0 = int(round(got["x"][0] / cell - 0.5)), int(round(got["y"][0] / cell - 0.5)) + ny, nx = got["z"].shape + pad = 40 + spec = GridSpec( + x_min=(ix0 - pad) * cell, + y_max=(iy0 + ny + pad) * cell, + cell_m=cell, + n_rows=ny + 2 * pad, + n_cols=nx + 2 * pad, + ) + whole = surf.surface_grid(spec, surf.build_contour_surface(spec, *fake["allseg"]))["z"] + whole = whole[pad:-pad, pad:-pad] + both = got["valid_mask"] + assert both.mean() > 0.9 + assert float(np.max(np.abs(got["z"][both] - whole[both]))) < 1e-3 + # 두 번째 읽기 = 셈 없음 + builds = fake["calls"]["build"] + again = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=cell) + assert again["computed"] == [] and fake["calls"]["build"] == builds + assert np.array_equal(np.nan_to_num(again["z"], nan=-1), np.nan_to_num(got["z"], nan=-1)) + # 상자로 읽기 = 그 상자 칸만 + box = (got["x"][10] - 5, got["y"][20] - 5, got["x"][60] + 5, got["y"][90] + 5) + part = store.load_terrain(EPSG, box=box, cell_m=cell) + assert part["z"].shape == (71, 51) + assert np.array_equal( + np.nan_to_num(part["z"], nan=-1), np.nan_to_num(got["z"][20:91, 10:61], nan=-1) + ) + # 원천(전국 파일) 바뀜 = 다시 셈 + fake["gpkg"].write_bytes(b"v2-longer") + store.load_terrain(EPSG, sheets=[WEST], cell_m=cell) + assert fake["calls"]["build"] == builds + 1 + + +def test_same_sheet_twice_at_once_computes_once(fake): + results = [] + threads = [ + threading.Thread( + target=lambda: results.append(store.ensure_sheet_terrain(WEST, EPSG, 10.0)) + ) + for _ in range(3) + ] + for t in threads: + t.start() + for t in threads: + t.join() + assert fake["calls"]["build"] == 1 and len(set(results)) == 1 and results[0] is not None + assert not list((store.TERRAIN_ROOT / WEST).glob(".*.lock")) + + +def test_missing_sheet_uses_zip_and_failure_not_stored(fake, monkeypatch): + a, b, z = fake["allseg"] + zlines = [(np.vstack([a[i], b[i]]), z[i]) for i in range(len(z))] + monkeypatch.setattr(store, "_gpkg_segments", lambda epsg, box: store._NO_SEGMENTS) + zip_path = store.TERRAIN_ROOT / WEST / f"{WEST}.zip" + + def fetch_fail(sheet_no, report): + return {}, "PermissionError: 로그인" + + monkeypatch.setattr(store, "_fetch_zips", fetch_fail) + got = store.load_terrain(EPSG, sheets=[WEST], cell_m=20.0) + assert ( + WEST in got["failed"] and not (store.TERRAIN_ROOT / WEST / "dem_5187_20m.key.json").exists() + ) + + stages = [] + + def fetch_ok(sheet_no, report): + report(10, f"도엽 받기 1/9 · {sheet_no}") + zip_path.parent.mkdir(parents=True, exist_ok=True) + zip_path.write_bytes(b"zip") + return {WEST: store._sha256(zip_path)}, None + + monkeypatch.setattr(store, "_fetch_zips", fetch_ok) + monkeypatch.setattr(store, "_zip_lines", lambda nos, epsg: zlines) + got = store.load_terrain( + EPSG, sheets=[WEST], cell_m=20.0, on_progress=lambda p, s: stages.append(s) + ) + assert got["computed"] == [WEST] and got["valid_mask"].mean() > 0.9 + assert any("도엽 받기" in s for s in stages) + builds = fake["calls"]["build"] + assert store.load_terrain(EPSG, sheets=[WEST], cell_m=20.0)["computed"] == [] + zip_path.write_bytes(b"zip-new") # 도엽 zip 바뀜 = 다시 셈 + store.load_terrain(EPSG, sheets=[WEST], cell_m=20.0) + assert fake["calls"]["build"] == builds + 1 + + +def test_empty_sheet_marked_and_skipped(fake, monkeypatch): + monkeypatch.setattr(store, "_gpkg_segments", lambda epsg, box: store._NO_SEGMENTS) + monkeypatch.setattr(store, "_fetch_zips", lambda sheet_no, report: ({}, None)) + got = store.load_terrain(EPSG, sheets=[WEST], cell_m=20.0) + assert got["z"].size == 0 and not got["failed"] + assert store.ensure_sheet_terrain(WEST, EPSG, 20.0) is None + assert fake["calls"]["build"] == 0