From 59ae735eda3cc3d77726913016d74c78578ac874 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Fri, 30 Oct 2026 23:09:57 +0900 Subject: [PATCH] =?UTF-8?q?feat(B04=20=EB=B0=B0=EC=88=98):=20=ED=95=98?= =?UTF-8?q?=EA=B0=95=20=EB=B0=A9=ED=96=A5=EC=9E=A5=20=EB=8F=84=EC=97=BD?= =?UTF-8?q?=EB=B3=84=20=EA=B3=B5=EC=9A=A9=20=EC=B0=BD=EA=B3=A0=20=E2=80=94?= =?UTF-8?q?=20=EB=8F=84=EC=97=BD=20+=20=EC=97=AC=EC=9C=A0=20100=20m=20?= =?UTF-8?q?=EB=A1=9C=20=ED=95=9C=20=EB=B2=88=20=EC=A7=93=EA=B3=A0=20?= =?UTF-8?q?=EC=9D=B4=EC=96=B4=20=EB=B6=99=EC=97=AC=20=EC=9D=BD=EA=B8=B0?= =?UTF-8?q?=EB=A7=8C=20(PLAN=20609-1=20=C2=B7=20=EC=9A=94=EA=B5=AC=20486?= =?UTF-8?q?=20=C2=B7=20493)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 새 파일 B04_PreProcess_Engine_Watershed_DescentShared.py — storage/_shared/terrain/<도엽>/descent_<좌표계>_1m.npz(칸 코드 int8 · 밴드 표고) + 열쇠(판 · 여유 · 규칙 · 전국 파일 지문 · zip 해시) · 짓는 법은 build_contour_descent 그대로 · 원천 = 전국 등고선 · 빠진 곳 = 607 도엽 zip - 배수 확장 source_for — 저장소 셋째 값(좌표계)이 있으면 공용 길 · 못 지으면 옛 프로젝트 길 · 2 값 저장소(시험)는 그대로 - 라우터 두 곳 — _prepare 가 좌표계를 함께 넘김 - 공용 지형 창고 — 좌표 변환 6 곳이 번호 없는 좌표계 글(PRJ WKT)도 받음 · 정수 입력 그대로 - 시험 test_609_1 — 도엽 넷 이어 붙임 = 통째 짓기 · 두 번째 짓기 0 · 원천 바뀌면 다시 · WKT · 확장 공용 길 = 프로젝트 길 · 실패 시 프로젝트 길 Co-Authored-By: Claude Opus 5.5 --- ...eProcess_Engine_Watershed_DescentShared.py | 321 ++++++++++++++++++ .../B04_PreProcess_Engine_Watershed_Expand.py | 18 +- .../B04_PreProcess_Router_Watershed.py | 10 +- .../B04_PreProcess_Router_Watershed_Live.py | 2 +- common_util/common_util_terrain_store.py | 27 +- resources/tester/test_609_1_descent_shared.py | 143 ++++++++ 6 files changed, 501 insertions(+), 20 deletions(-) create mode 100644 B04_PreProcess/B04_PreProcess_Engine_Watershed_DescentShared.py create mode 100644 resources/tester/test_609_1_descent_shared.py diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_DescentShared.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_DescentShared.py new file mode 100644 index 000000000..1057b796d --- /dev/null +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_DescentShared.py @@ -0,0 +1,321 @@ +"""도엽별 공용 하강장(609-1) — 배수 하강 방향장을 도엽마다 한 번 짓고 공용 창고에 영구 저장. + +자리 = `storage/_shared/terrain/<도엽>/descent_<좌표계>_<셀>m.npz` + `.key.json`(607 지형 창고 옆). +원천 = 전국 등고선 · 그 도엽에 없으면 607 길로 받은 도엽 zip. 도엽 + 여유 100 m 로 짓고(짓는 법은 +`build_contour_descent` 그대로 — 1 m · 같은 굽기 · 같은 EDT) 도엽 안쪽 칸(칸 중심 위경도)만 남김 — +칸마다 방향 코드(0~7 = `_NEIGHBOURS` · 8 = 제자리 · -1 = 무효) + 밴드 표고. 프로젝트는 격자를 덮는 +도엽을 이어 붙여 읽기만 — 코드라 도엽 경계를 넘는 물길도 그대로 이어짐. 짓는 규칙이 바뀌면 판을 +올림. +""" + +from __future__ import annotations + +import hashlib +import json +import logging +import math +import time +from pathlib import Path +from typing import Any + +import numpy as np +from pyproj import CRS, Transformer +from pyproj.exceptions import CRSError +from shapely.geometry import LineString, mapping, shape + +from B04_PreProcess.B04_PreProcess_Engine_MapSheet import sheet5k_to_bounds +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import ( + _NEIGHBOURS, + DESCENT_TILE_CELLS, + DESCENT_WINDOW_MARGIN_CELLS, + ContourDescent, + build_contour_descent, +) +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( + AZIMUTH_INVALID, + AZIMUTH_SINK, + AZIMUTH_STEPS, + GridSpec, +) +from common_util.common_util_atomic import atomic_write_npz +from common_util.common_util_crs import tm_epsg_for_lon +from common_util.common_util_json import atomic_write_json +from common_util.common_util_terrain_store import ( + NATIONAL_GPKG, + TERRAIN_ROOT, + _crs, + _gpkg_stamp, + _Lock, + _names, + _owned, + _sha256, + _zip_lines, + ensure_sheet_terrain, + sheets_for_box, +) +from config.config_system import DRAINAGE_CONTOUR_MIN_LENGTH_M + +logger = logging.getLogger(__name__) + +SHARED_DESCENT_VERSION = 1 +SHARED_DESCENT_MARGIN_M = 100.0 +_CODE_SINK = 8 +# 3×3 자리 (dr+1)*3 + (dc+1) → 코드 · 코드 → (dr, dc) · 방위 +_CODE_OF_SHIFT = np.array([0, 1, 2, 3, _CODE_SINK, 4, 5, 6, 7], dtype=np.int8) +_SHIFT_ROW = np.array([r for r, _ in _NEIGHBOURS] + [0]) +_SHIFT_COL = np.array([c for _, c in _NEIGHBOURS] + [0]) +_CODE_AZIMUTH = np.append( + np.rint(np.arctan2(_SHIFT_ROW[:8], _SHIFT_COL[:8]) / (2.0 * math.pi / AZIMUTH_STEPS)).astype( + np.int16 + ) + % AZIMUTH_STEPS, + AZIMUTH_SINK, +).astype(np.int16) + + +Crs = int | str # EPSG 번호 또는 번호 없는 좌표계 글(PRJ WKT — 용화 등 Korean 1985 + TOWGS84) + + +def descent_crs(crs: str) -> Crs | None: + """프로젝트 좌표계 → EPSG 번호(있으면) 또는 그 글 · 못 읽으면 None(프로젝트 저장본 길).""" + try: + return CRS.from_user_input(crs).to_epsg() or crs + except CRSError: + return None + + +def _crs_tag(crs: Crs) -> str: + """파일 이름 · 열쇠의 좌표계 — 번호 또는 WKT 지문(같은 PRJ 를 쓰는 프로젝트끼리 같음).""" + if isinstance(crs, int) or str(crs).isdigit(): + return str(int(crs)) + wkt = CRS.from_user_input(crs).to_wkt() + return "wkt" + hashlib.sha1(wkt.encode(), usedforsecurity=False).hexdigest()[:10] + + +def _shared_names(sheet_no: str, crs: Crs, cell_m: float) -> tuple[Path, Path, Path]: + stem = f"descent_{_crs_tag(crs)}_{cell_m:g}m" + folder = TERRAIN_ROOT / str(sheet_no) + return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock" + + +def _shared_key(crs: Crs, cell_m: float, zips: dict[str, str | None]) -> dict[str, Any]: + return { + "version": SHARED_DESCENT_VERSION, + "margin_m": SHARED_DESCENT_MARGIN_M, + "rule": [DRAINAGE_CONTOUR_MIN_LENGTH_M, DESCENT_WINDOW_MARGIN_CELLS, DESCENT_TILE_CELLS], + "crs": _crs_tag(crs), + "cell_m": float(cell_m), + "gpkg": _gpkg_stamp(), + "zips": zips, + } + + +def _shared_ok(sheet_no: str, crs: Crs, cell_m: float) -> dict[str, Any] | None: + """저장 열쇠가 지금 원천 · 규칙과 같으면 그 열쇠 — 아니면 None(다시 지음).""" + npz, key_path, _ = _shared_names(sheet_no, crs, 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 = _shared_key(crs, cell_m, zips) + if now["gpkg"] is None: # 전국 파일 없는 PC — 있는 저장본을 믿음(607 과 같음) + now["gpkg"] = stored.get("gpkg") + return stored if {k: stored.get(k) for k in now} == now else None + + +def _feature(line: Any, level: float) -> dict[str, Any]: + return {"geometry": mapping(line), "properties": {"elevation": float(level)}} + + +def _gpkg_features(crs: Crs, region: tuple[float, float, float, float]) -> list[dict[str, Any]]: + """전국 등고선에서 TM 상자에 닿는 선 — 통째 선(길이 거름이 프로젝트 길과 같게) · 사업 좌표.""" + if not NATIONAL_GPKG.is_file(): + return [] + import pyogrio + import shapely + + x0, y0, x1, y1 = region + t = np.linspace(0.0, 1.0, 21) + bx, by = Transformer.from_crs(_crs(crs), "EPSG:5179", always_xy=True).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]), + ) + frame = pyogrio.read_dataframe( + NATIONAL_GPKG, + bbox=(bx.min() - 20.0, by.min() - 20.0, bx.max() + 20.0, by.max() + 20.0), + columns=["CTRLN_HG"], + ) + if frame.empty: + return [] + to_tm = Transformer.from_crs("EPSG:5179", _crs(crs), always_xy=True) + lines = shapely.transform( + frame.geometry.values, lambda xy: np.column_stack(to_tm.transform(xy[:, 0], xy[:, 1])) + ) + levels = frame["CTRLN_HG"].to_numpy(dtype=np.float64) + return [ + _feature(line, level) + for line, level in zip(lines, levels) + if line is not None and np.isfinite(level) + ] + + +def _touches(features: list[dict[str, Any]], box: tuple[float, float, float, float]) -> bool: + x0, y0, x1, y1 = box + for feature in features: + bx0, by0, bx1, by1 = shape(feature["geometry"]).bounds + if bx1 >= x0 and bx0 <= x1 and by1 >= y0 and by0 <= y1: + return True + return False + + +def _codes(descent: ContourDescent) -> np.ndarray: + """방향장 → 칸 코드 (R, C) int8.""" + spec = descent.spec + index = np.arange(spec.size, dtype=np.int64) + receiver = descent.receiver.astype(np.int64) + shift = (receiver // spec.n_cols - index // spec.n_cols + 1) * 3 + ( + receiver % spec.n_cols - index % spec.n_cols + 1 + ) + code = _CODE_OF_SHIFT[shift] + code[~descent.valid.reshape(-1)] = -1 + return code.reshape(spec.n_rows, spec.n_cols) + + +def ensure_sheet_descent(sheet_no: str, crs: Crs, cell_m: float = 1.0) -> Path | None: + """도엽 하나의 공용 하강장 — 있으면 그대로 · 없거나 원천이 바뀌었으면 지음. 빈 도엽이면 None. + + 빠진 곳 받기 실패(로그인 · 망)는 RuntimeError — 저장 안 함(다음에 다시). + """ + sheet_no = str(sheet_no) + npz, key_path, lock_path = _shared_names(sheet_no, crs, cell_m) + if (stored := _shared_ok(sheet_no, crs, cell_m)) is not None: + return None if stored.get("empty") else npz + with _Lock(lock_path): + if (stored := _shared_ok(sheet_no, crs, cell_m)) is not None: + return None if stored.get("empty") else npz + started = time.monotonic() + ix0, iy0, own = _owned(sheet_no, crs, cell_m) + pad = math.ceil(SHARED_DESCENT_MARGIN_M / cell_m) + rows, cols = own.shape + spec = GridSpec( + (ix0 - pad) * cell_m, + (iy0 + rows + pad) * cell_m, + cell_m, + rows + 2 * pad, + cols + 2 * pad, + ) + x_max, y_min = spec.x_min + spec.n_cols * cell_m, spec.y_max - spec.n_rows * cell_m + features = _gpkg_features(crs, (spec.x_min, y_min, x_max, spec.y_max)) + zips: dict[str, str | None] = {} + own_box = (ix0 * cell_m, iy0 * cell_m, (ix0 + cols) * cell_m, (iy0 + rows) * cell_m) + if not _touches(features, own_box): + # 빠진 곳(섬 등) — 607 길로 도엽 zip 받기 · 그 열쇠의 zip 을 원천으로 + lon0, _, lon1, _ = sheet5k_to_bounds(sheet_no) + tm = crs if isinstance(crs, int) else tm_epsg_for_lon((lon0 + lon1) / 2) + ensure_sheet_terrain(sheet_no, tm) + terrain_key = json.loads(_names(sheet_no, tm, 2.0)[1].read_text(encoding="utf-8")) + zips = terrain_key.get("zips") or {} + features += [_feature(LineString(xy), z) for xy, z in _zip_lines(list(zips), crs)] + descent = build_contour_descent(spec, features) if features else None + key = _shared_key(crs, cell_m, zips) + key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게 + if descent is None or not descent.valid.any(): + npz.unlink(missing_ok=True) + atomic_write_json(key_path, {**key, "empty": True}) + return None + inner = (slice(pad, pad + rows), slice(pad, pad + cols)) + code = np.where(own[::-1, :], _codes(descent)[inner], -1).astype(np.int8) + band = np.where(code >= 0, descent.band_elevation[inner], np.nan).astype(np.float32) + atomic_write_npz( + npz, + code=code, + band=band, + origin=np.array([ix0 * cell_m, (iy0 + rows) * cell_m, cell_m]), + levels=np.array(descent.levels, dtype=np.float64), + ) + seconds = round(time.monotonic() - started, 1) + atomic_write_json( + key_path, {**key, "source": "zip" if zips else "gpkg", "seconds": seconds} + ) + logger.info("공용 하강장 %s (%s · %gm) %.1fs", sheet_no, _crs_tag(crs), cell_m, seconds) + return npz + + +def shared_descent(crs: Crs, spec: GridSpec) -> ContourDescent | None: + """`spec` 격자를 덮는 도엽 하강장을 이어 붙임(없는 도엽은 지음). + + 못 지은 도엽 · 칸 어긋남 = None(프로젝트 길로). + """ + + started = time.perf_counter() + rows, cols, cell_m = spec.n_rows, spec.n_cols, spec.cell_m + code = np.full((rows, cols), -1, dtype=np.int8) + band = np.full((rows, cols), np.nan, dtype=np.float32) + levels: set[float] = set() + box = (spec.x_min, spec.y_max - rows * cell_m, spec.x_min + cols * cell_m, spec.y_max) + for sheet_no in sheets_for_box(box, crs): + try: + path = ensure_sheet_descent(sheet_no, crs, cell_m) + except Exception: # noqa: BLE001 — 한 장이라도 못 지으면 프로젝트 길로 + logger.exception( + "배수유역: 공용 하강장 %s 을 짓지 못해 프로젝트 길로 갑니다.", sheet_no + ) + return None + if path is None: + continue + with np.load(path) as data: + x_min, y_max, _ = data["origin"].tolist() + row_shift = (spec.y_max - y_max) / cell_m + col_shift = (x_min - spec.x_min) / cell_m + row0, col0 = int(round(row_shift)), int(round(col_shift)) + if abs(row_shift - row0) > 1e-6 or abs(col_shift - col0) > 1e-6: + return None + piece, piece_band = data["code"], data["band"] + levels.update(data["levels"].tolist()) + r0, c0 = max(row0, 0), max(col0, 0) + r1, c1 = min(row0 + piece.shape[0], rows), min(col0 + piece.shape[1], cols) + if r0 >= r1 or c0 >= c1: + continue + src = (slice(r0 - row0, r1 - row0), slice(c0 - col0, c1 - col0)) + mine = piece[src] >= 0 + np.copyto(code[r0:r1, c0:c1], piece[src], where=mine) + np.copyto(band[r0:r1, c0:c1], piece_band[src], where=mine) + descent = _from_codes(spec, code, band, sorted(levels, reverse=True)) + logger.info( + "[계측] 공용 하강장 이어 붙임 %.2fs · %d×%d", time.perf_counter() - started, rows, cols + ) + return descent + + +def _from_codes( + spec: GridSpec, code: np.ndarray, band: np.ndarray, levels: list[float] +) -> ContourDescent: + """칸 코드 → 방향장. 격자 밖을 가리키는 칸은 제자리(`crop_descent` 와 같은 판정).""" + rows, cols = spec.n_rows, spec.n_cols + valid = code >= 0 + safe = np.where(valid, code, _CODE_SINK) + target_row = np.arange(rows)[:, None] + _SHIFT_ROW[safe] + target_col = np.arange(cols)[None, :] + _SHIFT_COL[safe] + moves = ( + valid + & (safe != _CODE_SINK) + & (target_row >= 0) + & (target_row < rows) + & (target_col >= 0) + & (target_col < cols) + ) + flat = np.arange(spec.size, dtype=np.int64).reshape(rows, cols) + receiver = np.where(moves, target_row * cols + target_col, flat) + length = np.hypot(_SHIFT_ROW[safe], _SHIFT_COL[safe]) * spec.cell_m + return ContourDescent( + spec=spec, + band_elevation=band, + valid=valid, + receiver=receiver.reshape(-1).astype(np.int32), + step_length=np.where(moves, length, 0.0).astype(np.float32).reshape(-1), + azimuth=np.where(valid, _CODE_AZIMUTH[safe], AZIMUTH_INVALID).astype(np.int16).reshape(-1), + levels=levels, + ) diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py index 8865fb76f..c9ea61078 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py @@ -15,6 +15,8 @@ (`DESCENT_SPARE_M` · 342-3). 격자가 커질 때마다 다시 짓던 것(8회 · 30 s)을 한 번으로 줄인다. 저장소(`descent_store` = 폴더 · 도엽 등고선 지문)를 주면 지은 것을 남기고, 덮는 저장본이 있으면 읽기만 한다(347-1). 방향장은 노선 최저점 하한 없이 짓는다 — 노선이 바뀌어도 같은 저장본(520-6). +저장소에 좌표계(셋째 값)가 있으면 도엽별 공용 하강장을 이어 붙여 읽기만 한다(609-1) — 못 지으면 +위 프로젝트 길로. 확장은 **큰 방향장 위 한 번 셈**으로 낸다(520-6) — 큰 격자 전체에서 도로 도달을 한 번 판정하고 「도로로 오는 셀 + 확장 폭 띠 ∪ 1차 영역」으로 잘라 돌려준다. 도로로 오는 셀이 큰 격자 @@ -41,6 +43,10 @@ from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import ( iter_saved_descents, save_descent, ) +from B04_PreProcess.B04_PreProcess_Engine_Watershed_DescentShared import ( + descent_crs, + shared_descent, +) from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import ( FlowClassification, RoadRaster, @@ -150,7 +156,7 @@ def expand_by_red_boundary( elevation_floor_m: float | None = None, band_m: float = DRAINAGE_RED_EXPAND_BAND_M, max_rounds: int = DRAINAGE_RED_EXPAND_MAX_ROUNDS, - descent_store: tuple[Path, str] | None = None, + descent_store: tuple[Path, str] | tuple[Path, str, str] | None = None, should_stop: Callable[[], bool] | None = None, ) -> RedExpansion | None: """큰 방향장 위 한 번 셈 — 가장자리에 닿으면 최외곽 적색 셀 주변으로 넓히며 반복한다. @@ -165,15 +171,21 @@ def expand_by_red_boundary( started_cells = int(domain.sum()) spare_cells = max(0, int(round(DESCENT_SPARE_M / spec.cell_m))) store = None + shared_crs = None if descent_store is not None: - folder, contour_key = descent_store + folder, contour_key, *crs = descent_store store = ( folder, descent_input_key(contour_key, contour_features, spec.cell_m), ) + shared_crs = descent_crs(crs[0]) if crs else None def source_for(grid: GridSpec, mask: np.ndarray) -> ContourDescent: - """`grid` 를 덮는 큰 방향장 — 저장본이 덮으면 읽기만 · 아니면 둘레 여유째 짓고 남김.""" + """`grid` 를 덮는 큰 방향장 — 공용 도엽 하강장(609-1) · 프로젝트 저장본 · 짓고 남김 차례.""" + if shared_crs is not None: + shared = shared_descent(shared_crs, _pad_spec(grid, mask, spare_cells)[0]) + if shared is not None: + return shared if store is not None: for saved in iter_saved_descents(*store): if crop_descent(saved, grid) is not None: diff --git a/B04_PreProcess/B04_PreProcess_Router_Watershed.py b/B04_PreProcess/B04_PreProcess_Router_Watershed.py index 992b055ec..8fde1465c 100644 --- a/B04_PreProcess/B04_PreProcess_Router_Watershed.py +++ b/B04_PreProcess/B04_PreProcess_Router_Watershed.py @@ -215,12 +215,13 @@ async def _prepare(project_id: UUID) -> dict[str, Any] | JSONResponse: ) # 설계 계통과 같은 프로젝트 좌표계로 맞춘다. 도엽 재투영도 이 좌표계로 간다. - def sheets() -> tuple[Transformer, list[dict[str, Any]], list[dict[str, Any]], str]: + def sheets() -> tuple[str, Transformer, list[dict[str, Any]], list[dict[str, Any]], str]: source_crs = resolve_project_crs( project_root, route_crs_input=planned.crs_input, db_epsg=epsg ) to_metric_transformer = Transformer.from_crs("EPSG:4326", source_crs, always_xy=True) return ( + source_crs, Transformer.from_crs(source_crs, "EPSG:4326", always_xy=True), *_metric_sheets( _sheet_dir(stored_path), @@ -231,7 +232,9 @@ async def _prepare(project_id: UUID) -> dict[str, Any] | JSONResponse: ) # 재투영(저장본 없으면 3.8 s) · 피클 읽기는 스레드로 — 이벤트 루프를 안 막음(505-1) - to_lonlat_transformer, streams, contour_features, contour_key = await asyncio.to_thread(sheets) + crs, to_lonlat_transformer, streams, contour_features, contour_key = await asyncio.to_thread( + sheets + ) return { "route_source": route_file.name, "vertices": planned.vertices, @@ -239,6 +242,7 @@ async def _prepare(project_id: UUID) -> dict[str, Any] | JSONResponse: "streams": streams, "contours": contour_features, "contour_key": contour_key, + "crs": crs, # 공용 하강장(609-1) "stored_path": stored_path, "to_lonlat": lambda x, y: to_lonlat_transformer.transform(x, y), } @@ -413,7 +417,7 @@ async def get_primary_region( prepared["contours"], prepared["streams"], # 흐름 방향 저장본(347-1) — 노선이 저장 범위 안이면 읽기만 - (drainage_dir(prepared["stored_path"]), prepared["contour_key"]), + (drainage_dir(prepared["stored_path"]), prepared["contour_key"], prepared["crs"]), ) if preview is None: return JSONResponse( diff --git a/B04_PreProcess/B04_PreProcess_Router_Watershed_Live.py b/B04_PreProcess/B04_PreProcess_Router_Watershed_Live.py index 336f6fb46..10df2a35b 100644 --- a/B04_PreProcess/B04_PreProcess_Router_Watershed_Live.py +++ b/B04_PreProcess/B04_PreProcess_Router_Watershed_Live.py @@ -134,7 +134,7 @@ async def post_drainage_live( vertices, prepared["contours"], prepared["streams"], - (drainage_dir(context.stored_path), prepared["contour_key"]), + (drainage_dir(context.stored_path), prepared["contour_key"], prepared["crs"]), superseded, ) except AnalysisStopped: diff --git a/common_util/common_util_terrain_store.py b/common_util/common_util_terrain_store.py index dad1b3263..62850abf8 100644 --- a/common_util/common_util_terrain_store.py +++ b/common_util/common_util_terrain_store.py @@ -54,6 +54,11 @@ Segments = tuple[np.ndarray, np.ndarray, np.ndarray] _NO_SEGMENTS: Segments = (np.zeros((0, 2)), np.zeros((0, 2)), np.zeros(0)) +def _crs(epsg: int | str) -> str: + """EPSG 번호 또는 좌표계 글(번호 없는 WKT · 배수 하강장 609-1) → pyproj 입력.""" + return f"EPSG:{int(epsg)}" if isinstance(epsg, int) or str(epsg).isdigit() else str(epsg) + + 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" @@ -130,15 +135,13 @@ class _Lock: self.path.unlink(missing_ok=True) -def sheets_for_box(box: Box, epsg: int) -> list[str]: +def sheets_for_box(box: Box, epsg: int | str) -> 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 - ) + lon, lat = Transformer.from_crs(_crs(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): @@ -149,10 +152,10 @@ def sheets_for_box(box: Box, epsg: int) -> list[str]: return sheets -def _owned(sheet_no: str, epsg: int, cell_m: float) -> tuple[int, int, np.ndarray]: +def _owned(sheet_no: str, epsg: int | str, 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) + to_tm = Transformer.from_crs("EPSG:4326", _crs(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)]) @@ -164,9 +167,7 @@ def _owned(sheet_no: str, epsg: int, cell_m: float) -> tuple[int, int, np.ndarra 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 - ) + clon, clat = Transformer.from_crs(_crs(epsg), "EPSG:4326", always_xy=True).transform(gx, gy) own = (clon >= lon0) & (clon < lon1) & (clat >= lat0) & (clat < lat1) return ix0, iy0, own @@ -179,7 +180,7 @@ def _gpkg_segments(epsg: int, box: Box) -> Segments: import shapely x0, y0, x1, y1 = box - to_5179 = Transformer.from_crs(f"EPSG:{int(epsg)}", "EPSG:5179", always_xy=True) + to_5179 = Transformer.from_crs(_crs(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)]), @@ -208,13 +209,13 @@ def _gpkg_segments(epsg: int, box: Box) -> Segments: 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) + to_tm = Transformer.from_crs("EPSG:5179", _crs(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]]: +def _zip_lines(sheet_nos: list[str], epsg: int | str) -> 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 @@ -230,7 +231,7 @@ def _zip_lines(sheet_nos: list[str], epsg: int) -> list[tuple[np.ndarray, float] 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)}") + frame = frame.to_crs(_crs(epsg)) column = next((k for k in ELEVATION_KEYS if k in frame.columns), None) if column is None: continue diff --git a/resources/tester/test_609_1_descent_shared.py b/resources/tester/test_609_1_descent_shared.py new file mode 100644 index 000000000..76181c4a6 --- /dev/null +++ b/resources/tester/test_609_1_descent_shared.py @@ -0,0 +1,143 @@ +"""609-1 도엽별 공용 하강장 — 도엽 넷 이어 붙임 = 통째 짓기 · 두 번째는 짓지 않음 · +원천 바뀌면 다시 · 번호 없는 좌표계 · 배수 확장이 공용 길 = 프로젝트 길 · 못 지으면 프로젝트 길. + +./venv/Scripts/python.exe -m pytest resources/tester/test_609_1_descent_shared.py -q +""" + +import numpy as np +import pytest +from pyproj import Transformer +from shapely.geometry import LineString, box, mapping + +from B04_PreProcess import B04_PreProcess_Engine_Watershed_DescentShared as shared +from B04_PreProcess import B04_PreProcess_Engine_Watershed_Expand as expand +from B04_PreProcess.B04_PreProcess_Engine_MapSheet import sheet5k_to_bounds +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import build_contour_descent +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec +from common_util import common_util_terrain_store as store + +EPSG = 5187 +CELL = 10.0 # 시험만 10 m(실제 1 m) — 도엽 넷을 빨리 지음 + + +def _corner(): + """도엽 36901074 남동 모서리(도엽 넷이 만나는 점) · 사업 좌표 10 m 눈금.""" + _, lat0, lon1, _ = sheet5k_to_bounds("36901074") + x, y = Transformer.from_crs("EPSG:4326", f"EPSG:{EPSG}", always_xy=True).transform(lon1, lat0) + return round(x / CELL) * CELL, round(y / CELL) * CELL + + +def _features(): + """모서리 곁 봉우리 — 30 m 마다 5 m 내려가는 고리(중심을 비켜 도엽 경계를 비스듬히 지남).""" + cx, cy = _corner() + out = [] + for i, radius in enumerate(range(30, 4500, 30)): + t = np.linspace(0.0, 2 * np.pi, max(64, radius // 5) + 1) + line = LineString( + np.column_stack([cx + 137 + radius * np.cos(t), cy - 91 + radius * np.sin(t)]) + ) + out.append({"geometry": mapping(line), "properties": {"elevation": 1000.0 - 5.0 * i}}) + return out + + +@pytest.fixture +def world(tmp_path, monkeypatch): + monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path) + monkeypatch.setattr(shared, "TERRAIN_ROOT", tmp_path) + monkeypatch.setattr(shared, "_gpkg_stamp", lambda: "gpkg-1") + features = _features() + + def from_gpkg(epsg, region): + area = box(*region) + return [f for f in features if LineString(f["geometry"]["coordinates"]).intersects(area)] + + monkeypatch.setattr(shared, "_gpkg_features", from_gpkg) + built = [] + real = shared.build_contour_descent + monkeypatch.setattr(shared, "build_contour_descent", lambda *a: built.append(1) or real(*a)) + return features, built + + +def _spec(half=1500.0): + cx, cy = _corner() + n = int(2 * half / CELL) + return GridSpec(cx - half, cy + half, CELL, n, n) + + +def test_stitched_equals_whole(world): + features, built = world + spec = _spec() + stitched = shared.shared_descent(EPSG, spec) + assert len(built) == 4 # 도엽 넷 + pad = 20 + big = GridSpec( + spec.x_min - pad * CELL, + spec.y_max + pad * CELL, + CELL, + spec.n_rows + 2 * pad, + spec.n_cols + 2 * pad, + ) + whole = expand.crop_descent(build_contour_descent(big, features), spec) + assert stitched.valid.sum() > 0.99 * spec.size + assert np.array_equal(stitched.valid, whole.valid) + assert np.array_equal(stitched.receiver, whole.receiver) + assert np.array_equal(stitched.azimuth, whole.azimuth) + assert np.array_equal(stitched.step_length, whole.step_length) + assert np.array_equal(stitched.band_elevation, whole.band_elevation, equal_nan=True) + + +def test_second_read_builds_nothing_and_source_change_rebuilds(world, monkeypatch): + _, built = world + shared.shared_descent(EPSG, _spec(500.0)) + first = len(built) + shared.shared_descent(EPSG, _spec(800.0)) # 다른 프로젝트 · 같은 도엽 + assert first == 4 and len(built) == first + monkeypatch.setattr(shared, "_gpkg_stamp", lambda: "gpkg-2") # 전국 파일 바뀜 + shared.shared_descent(EPSG, _spec(500.0)) + assert len(built) == first + 4 + + +def _expand(store_tuple): + cx, cy = _corner() + spec = GridSpec(cx - 400, cy + 400, CELL, 80, 80) + domain = np.zeros((80, 80), bool) + domain[20:60, 20:60] = True + route = LineString([(cx - 300, cy - 250), (cx + 300, cy - 250)]) + return expand.expand_by_red_boundary( + spec, domain, _features(), route, [], None, descent_store=store_tuple + ) + + +def test_expand_shared_equals_project_path(world, tmp_path): + _, built = world + project = _expand((tmp_path / "p1", "도엽-1")) + via_shared = _expand((tmp_path / "p2", "도엽-1", f"EPSG:{EPSG}")) + assert built # 공용 길로 지음 + assert not list((tmp_path / "p2").glob("descent_*.npz")) # 프로젝트엔 안 남김 + assert via_shared.analysis.spec == project.analysis.spec + assert np.array_equal(via_shared.analysis.flow.reaches_road, project.analysis.flow.reaches_road) + assert np.array_equal(via_shared.analysis.descent.receiver, project.analysis.descent.receiver) + + +def test_shared_failure_falls_back_to_project(world, tmp_path, monkeypatch): + def fail(*_): + raise RuntimeError("받기 실패") + + monkeypatch.setattr(shared, "ensure_sheet_descent", fail) + result = _expand((tmp_path / "p3", "도엽-1", f"EPSG:{EPSG}")) + assert result is not None and list((tmp_path / "p3").glob("descent_*.npz")) + + +def test_crs_without_epsg_number(world, tmp_path): + """PRJ WKT(번호 없음 · 용화 등)도 공용 길 — 좌표계 지문으로 따로 저장.""" + text = ( + "+proj=tmerc +lat_0=38 +lon_0=129 +k=1 +x_0=200000 +y_0=600000 +ellps=GRS80" + " +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs" + ) # 5187 과 같은 격자 · 번호 없음 + crs = shared.descent_crs(text) + assert crs == text and shared.descent_crs("EPSG:5187") == EPSG + spec = _spec(500.0) + by_text = shared.shared_descent(crs, spec) + by_number = shared.shared_descent(EPSG, spec) + assert np.array_equal(by_text.receiver, by_number.receiver) + assert {p.name.split("_")[1][:3] for p in tmp_path.glob("*/descent_*.npz")} == {"wkt", "518"}