diff --git a/common_util/common_util_terrain_store.py b/common_util/common_util_terrain_store.py index 62850abf8..b4302f407 100644 --- a/common_util/common_util_terrain_store.py +++ b/common_util/common_util_terrain_store.py @@ -8,6 +8,8 @@ 도엽 + 여유 200 m 등고선으로 셈 → 도엽 안쪽 칸(칸 중심 위경도가 그 도엽)만 저장 → 이어 붙이면 통째 셈과 같음. 격자 원점 = 좌표 격자 배수(칸 모서리 = k·cell_m). 열쇠 = 계산 방식 판 + 여유 + 전국 파일 지문(크기 · 수정 시각) + 쓴 도엽 zip 해시. +띠 밖 좌표계(5176 · 띠 아닌 5186 등)는 그 도엽 띠(5185~5187) 지형을 쌍선형으로 옮겨 담음(614) — +열쇠에 `via` = 띠 · 이미 직접 센 저장본은 그대로 씀. """ from __future__ import annotations @@ -32,6 +34,7 @@ from B04_PreProcess.B04_PreProcess_Engine_MapSheet import ( ) 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_crs import tm_epsg_for_lon from common_util.common_util_json import atomic_write_json from common_util.common_util_terrain_surface import ( TERRAIN_ALGO, @@ -47,6 +50,9 @@ NATIONAL_GPKG = Path(PROJECT_ROOT) / "resources" / "data_global_contours" / "nat TERRAIN_MARGIN_M = 200.0 # ponytail: 잠금 파일이 이보다 오래되면 죽은 셈으로 보고 깸 · 도엽 하나 셈 1~5 분 LOCK_STALE_S = 15 * 60 +# 띠 밖 좌표계 = 띠 지형 옮겨 담기(614) · 직접 셈 대비 문턱 못 넘으면 False(직접 셈) +# 문턱 = 참값(등고선 위 점) 대비 오차 — 옮김이 더하는 95 % ≤ 0.10 m · 최대 ≤ 직접 셈 최대 +REPROJECT_FROM_BAND = True Progress = Callable[[int, str], None] | None Box = tuple[float, float, float, float] @@ -65,6 +71,15 @@ def _names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock" +def _via(sheet_no: str, epsg: int | str) -> int | None: + """옮겨 담을 띠 좌표계 — 이미 띠 · 옮기기 꺼짐 · 번호 없는 좌표계면 None(직접 셈).""" + if not REPROJECT_FROM_BAND or not str(epsg).isdigit(): + return None + lon0, _, lon1, _ = sheet5k_to_bounds(sheet_no) + band = tm_epsg_for_lon((lon0 + lon1) / 2) + return None if band == int(epsg) else band + + def _gpkg_stamp() -> str | None: try: stat = NATIONAL_GPKG.stat() @@ -100,6 +115,8 @@ def _stored_ok(sheet_no: str, epsg: int, cell_m: float) -> dict[str, Any] | None return None if not stored.get("empty") and not npz.is_file(): return None + if stored.get("via") is not None and stored["via"] != _via(sheet_no, epsg): + 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 — 있는 저장본을 믿음(시놀로지 덮어쓰기 핑퐁 막음) @@ -282,6 +299,8 @@ def ensure_sheet_terrain( stored = _stored_ok(sheet_no, epsg, cell_m) # 기다리는 동안 남이 셌을 수 있음 if stored is not None: return None if stored.get("empty") else npz + if (band := _via(sheet_no, epsg)) is not None: + return _reproject(sheet_no, epsg, band, cell_m, report) started = time.monotonic() report(2, f"전국 등고선 읽기 · {sheet_no}") ix0, iy0, own = _owned(sheet_no, epsg, cell_m) @@ -348,6 +367,75 @@ def ensure_sheet_terrain( return npz +def _extend_edges(z: np.ndarray) -> np.ndarray: + """빈 칸 중 값 칸 바로 옆 한 칸 = 안쪽 두 칸으로 직선 연장(가로 뒤 세로 · 모서리 포함) — + 이웃 띠 지형이 없는 도엽 가장자리에서도 쌍선형 네 칸을 채움.""" + for axis in (1, 0): + src = np.moveaxis(z, axis, -1) + out = src.copy() + inner, outer = src[..., 1:-1], (src[..., 2:], src[..., :-2]) + for side, far in ((slice(None, -2), outer[0]), (slice(2, None), outer[1])): + gap = np.isnan(out[..., side]) & np.isfinite(inner) & np.isfinite(far) + out[..., side][gap] = (2 * inner - far)[gap] + z = np.moveaxis(out, -1, axis) + return z + + +def _bilinear(grid: dict[str, Any], px: np.ndarray, py: np.ndarray) -> np.ndarray: + """격자(셀 중심 · y 오름차순) 위 점들의 쌍선형 높이 — 네 칸 중 빈 칸이 있으면 NaN.""" + z, cell = _extend_edges(grid["z"]), grid["cell_m"] + fx, fy = (px - grid["x"][0]) / cell, (py - grid["y"][0]) / cell + i = np.clip(np.floor(fx).astype(np.int64), 0, z.shape[1] - 2) + j = np.clip(np.floor(fy).astype(np.int64), 0, z.shape[0] - 2) + tx, ty = fx - i, fy - j + return (z[j, i] * (1 - tx) + z[j, i + 1] * tx) * (1 - ty) + ( + z[j + 1, i] * (1 - tx) + z[j + 1, i + 1] * tx + ) * ty + + +def _reproject( + sheet_no: str, epsg: int, band: int, cell_m: float, report: Callable[[int, str], None] +) -> Path | None: + """띠 지형(없으면 띠로 먼저 셈) → 이 좌표계 격자로 쌍선형 옮겨 담아 저장 — 잠금 안에서 부름. + + 가장자리 칸은 창고에 있는 이웃 띠 지형만 씀(없으면 직선 연장 · 새로 안 셈). + """ + started = time.monotonic() + npz, key_path, _ = _names(sheet_no, epsg, cell_m) + src = ensure_sheet_terrain(sheet_no, band, cell_m, lambda p, s: report(p * 80 // 100, s)) + band_key = json.loads(_names(sheet_no, band, cell_m)[1].read_text(encoding="utf-8")) + key = {**_key(epsg, cell_m, band_key.get("zips") or {}), "via": band} + key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게 + if src is None: + npz.unlink(missing_ok=True) + atomic_write_json(key_path, {**key, "empty": True}) + return None + report(85, f"좌표계 옮겨 담기 · {sheet_no}") + ix0, iy0, own = _owned(sheet_no, epsg, cell_m) + rows, cols = np.nonzero(own) + to_band = Transformer.from_crs(_crs(epsg), _crs(band), always_xy=True) + px, py = to_band.transform((ix0 + cols + 0.5) * cell_m, (iy0 + rows + 0.5) * cell_m) + pad = 2 * cell_m + box = (px.min() - pad, py.min() - pad, px.max() + pad, py.max() + pad) + near = [n for n in sheets_for_box(box, band) if n == sheet_no or _stored_ok(n, band, cell_m)] + core = np.full(own.shape, np.nan, dtype=np.float32) + core[rows, cols] = _bilinear(load_terrain(band, sheets=near, box=box, cell_m=cell_m), px, py) + 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": "reproject", "seconds": round(time.monotonic() - started, 1)}, + ) + report(100, f"완료 · {sheet_no}") + logger.info( + "공용 지형 %s 옮겨 담음 EPSG:%s → %s %.1fs", + sheet_no, + band, + epsg, + time.monotonic() - started, + ) + return npz + + def load_terrain( epsg: int, *, diff --git a/resources/tester/test_614_terrain_reproject.py b/resources/tester/test_614_terrain_reproject.py new file mode 100644 index 000000000..653cc12d0 --- /dev/null +++ b/resources/tester/test_614_terrain_reproject.py @@ -0,0 +1,125 @@ +"""띠 밖 좌표계 = 띠 지형 옮겨 담기(PLAN 614-1) — 띠로 먼저 셈 · 옮긴 것도 그 좌표계 열쇠로 저장 · +둘째 읽기 셈 없음 · 참값(등고선 위 점) 대비 옮김 오차 ≤ 직접 + 0.10 m · +이미 직접 센 저장본 그대로 · 끄면 직접 셈 · +빈 띠 = 빈 도엽 · 가장자리 직선 연장. + +./venv/Scripts/python.exe -m pytest resources/tester/test_614_terrain_reproject.py -q +""" + +import json + +import numpy as np +import pytest +from pyproj import Transformer + +from B04_PreProcess.B04_PreProcess_Engine_MapSheet import sheet5k_to_bounds +from common_util import common_util_terrain_store as store +from common_util import common_util_terrain_surface as surf + +BAND, OTHER = 5187, 5186 +SHEET = "36901074" +CELL = 10.0 + + +@pytest.fixture +def fake(tmp_path, monkeypatch): + """띠(5187) 좌표의 굽은 비탈(남서 먼 곳 중심 동심 호 · 간격 점점 넓음) — 요청 좌표계로 옮겨 + 상자 안 선분만 돌려줌.""" + lon0, lat0, _, _ = sheet5k_to_bounds(SHEET) + to_band = Transformer.from_crs("EPSG:4326", f"EPSG:{BAND}", always_xy=True) + cx, cy = to_band.transform(lon0 - 0.04, lat0 - 0.04) + arc = np.linspace(0.0, np.pi / 2, 401)[:, None] + lines = [] + for k in range(250): + r = 3000.0 + 30.0 * k + 0.1 * k * k + lines.append((np.hstack([cx + r * np.cos(arc), cy + r * np.sin(arc)]), 500.0 + 5.0 * k)) + a, b, z = surf.segments_from_lines(lines) + gpkg = tmp_path / "national.gpkg" + gpkg.write_bytes(b"v1") + calls = {"build": 0} + real_build = store.build_contour_surface + + def segments(epsg, box): + tr = Transformer.from_crs(f"EPSG:{BAND}", store._crs(epsg), always_xy=True) + pa, pb = np.column_stack(tr.transform(*a.T)), np.column_stack(tr.transform(*b.T)) + hit = (np.maximum(pa[:, 0], pb[:, 0]) >= box[0]) & ( + np.minimum(pa[:, 0], pb[:, 0]) <= box[2] + ) + hit &= (np.maximum(pa[:, 1], pb[:, 1]) >= box[1]) & ( + np.minimum(pa[:, 1], pb[:, 1]) <= box[3] + ) + return pa[hit], pb[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) + monkeypatch.setattr(store, "REPROJECT_FROM_BAND", True) + return calls + + +def _key(epsg): + path = store._names(SHEET, epsg, CELL)[1] + return json.loads(path.read_text(encoding="utf-8")) + + +def test_other_crs_reprojects_band_and_stores(fake): + assert store._via(SHEET, OTHER) == BAND and store._via(SHEET, BAND) is None + got = store.load_terrain(OTHER, sheets=[SHEET], cell_m=CELL) + assert got["computed"] == [SHEET] and not got["failed"] + assert fake["build"] == 1 # 띠 한 번만 셈 · 옮기기는 셈 아님 + assert _key(BAND).get("via") is None and _key(OTHER)["via"] == BAND + assert _key(OTHER)["source"] == "reproject" + again = store.load_terrain(OTHER, sheets=[SHEET], cell_m=CELL) + assert again["computed"] == [] and fake["build"] == 1 + assert np.array_equal(np.nan_to_num(again["z"], nan=-1), np.nan_to_num(got["z"], nan=-1)) + # 직접 셈(끔)과 견줌 — 같은 칸 거의 다 · 높이 차 작음 + store.REPROJECT_FROM_BAND = False + assert store._stored_ok(SHEET, OTHER, CELL) is None # 끄면 옮긴 저장본 버림 + direct = store.load_terrain(OTHER, sheets=[SHEET], cell_m=CELL) + assert direct["computed"] == [SHEET] and fake["build"] == 2 + assert _key(OTHER).get("via") is None + assert direct["z"].shape == got["z"].shape + both = direct["valid_mask"] & got["valid_mask"] + assert both.sum() > 0.98 * direct["valid_mask"].sum() + diff = np.abs(direct["z"][both] - got["z"][both]) + assert float(np.percentile(diff, 95)) < 0.5 and float(diff.max()) < 2.0 # 10 m 칸 견줌 + # 문턱(브레인) = 참값(등고선 위 점) 대비 오차 — 옮김이 더하는 95 % ≤ 0.10 m · 최대 ≤ 직접 최대 + gx, gy = direct["x"], direct["y"] + a, b, z = store._gpkg_segments(OTHER, (gx[2], gy[2], gx[-3], gy[-3])) + px, py = (a[:, 0] + b[:, 0]) / 2, (a[:, 1] + b[:, 1]) / 2 + inside = (px > gx[2]) & (px < gx[-3]) & (py > gy[2]) & (py < gy[-3]) + err_direct = np.abs(store._bilinear(direct, px[inside], py[inside]) - z[inside]) + err_moved = np.abs(store._bilinear(got, px[inside], py[inside]) - z[inside]) + ok = np.isfinite(err_direct) & np.isfinite(err_moved) + assert ok.sum() > 1000 + added = np.percentile(err_moved[ok], 95) - np.percentile(err_direct[ok], 95) + assert added <= 0.10 and err_moved[ok].max() <= err_direct[ok].max() + # 다시 켜도 직접 센 저장본은 그대로 씀 + store.REPROJECT_FROM_BAND = True + assert store._stored_ok(SHEET, OTHER, CELL) is not None + store.load_terrain(OTHER, sheets=[SHEET], cell_m=CELL) + assert fake["build"] == 2 + + +def test_empty_band_is_empty(fake, monkeypatch): + monkeypatch.setattr(store, "_gpkg_segments", lambda epsg, box: store._NO_SEGMENTS) + monkeypatch.setattr(store, "_fetch_zips", lambda sheet_no, report: ({}, None)) + assert store.ensure_sheet_terrain(SHEET, OTHER, CELL) is None + assert _key(BAND)["empty"] and _key(OTHER)["empty"] and _key(OTHER)["via"] == BAND + assert store.ensure_sheet_terrain(SHEET, OTHER, CELL) is None + + +def test_extend_edges_and_bilinear_exact_on_plane(): + y, x = np.mgrid[0:12, 0:15].astype(float) + plane = 3.0 * x - 2.0 * y + 7.0 + z = np.full_like(plane, np.nan) + z[1:-1, 1:-1] = plane[1:-1, 1:-1] + assert np.allclose(store._extend_edges(z), plane) + grid = {"z": z, "x": np.arange(15) + 0.5, "y": np.arange(12) + 0.5, "cell_m": 1.0} + px, py = np.array([0.5, 7.3, 14.4]), np.array([0.6, 5.9, 11.5]) + assert np.allclose(store._bilinear(grid, px, py), 3.0 * (px - 0.5) - 2.0 * (py - 0.5) + 7.0)