feat(공용 지형): 띠 밖 좌표계 = 띠 지형 쌍선형 옮겨 담기 (PLAN 614-1 · 요구 494)

- ensure_sheet_terrain: 도엽 띠(5185~5187) 아닌 좌표계면 띠 저장본(없으면 띠로 먼저 셈)을 그 좌표계 2 m 격자로 쌍선형 옮겨 담아 저장 · 열쇠 via = 띠
- 이미 직접 센 저장본 그대로 · REPROJECT_FROM_BAND 한 줄로 끔 · 가장자리 = 창고에 있는 이웃 띠만 · 없으면 직선 연장
- 문턱(브레인) = 참값(등고선 위 점) 대비 오차 · 옮김이 더하는 95 % ≤ 0.10 m · 최대 ≤ 직접 최대
- 시험 test_614 새로

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

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,
     *,

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)