feat(공용 지형): 공용 지형 창고 · 흠 고친 지형 계산 새 모듈 둘 (PLAN 607-1~3 · 요구 486~491)
- common_util_terrain_surface — 방향 무관 굽기(선분 끝 정렬 · 큰 축 칸 중심) · 폐합 안쪽 채움을 그 고리 둘레로 · 보간은 B04 build_laplace 그대로 부름 - common_util_terrain_store — storage/_shared/terrain/<도엽>/dem_<epsg>_<격자>m.npz · 전국 등고선 + 여유 200 m · 도엽 안쪽 칸만 · 열쇠 = 방식 판 + 전국 파일 크기·수정 시각 + 도엽 zip 해시 · 잠금 파일 · 빠진 곳만 도엽 zip(B04 ensure_sheets · prune 안 부름) · load_terrain 진행 % · 단계 글 - 시험 test_607_terrain_store — 뒤집기·조각 굽기 같음 · 봉우리·웅덩이 방향 · 조각 = 통째 · 한 번만 셈 · 원천 바뀜 · zip 길 · 실패 저장 안 함 - B04 · B05 파일 고침 없음 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -0,0 +1,417 @@
|
||||
"""공용 지형 창고 — 도엽(1:5,000)마다 2 m 지형을 처음 쓸 때 한 번 셈 · 영구 저장
|
||||
(PLAN 607-2 · 607-3).
|
||||
|
||||
자리 = `storage/_shared/terrain/<도엽>/dem_<epsg>_<격자>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,
|
||||
}
|
||||
@@ -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),
|
||||
}
|
||||
Reference in New Issue
Block a user