eomsangdon
9fc4462d74
사용자 지시(2026-08-30): IDW·거리비례·TIN 곡면 제거, 5m 도엽등고선 기준 TIN 추가. TIN(도엽선)은 격자에 구운 라인 셀이 아니라 벡터 등고선 정점의 원좌표를 그대로 Delaunay로 잇는다. 격자 밖은 100m 여유만 물고 자르며, 길이 필터는 두지 않는다 (짧은 봉우리 폐합 링을 버리면 마루가 통째로 평평해진다). 거리비례는 버튼에서만 빼고 함수는 남긴다 — 라플라스·TPS·ANUDEM·다중해상도가 초기추정으로 계속 쓴다. 기본 방식은 distance가 빠져 multires로 옮겼다. 검증: 실데이터 767×840·1m 6종 93.5s. 서피스에서 다시 뽑은 5m 선의 평면 이탈이 TIN(도엽선) 중앙 0.12m·2m초과 8.1%로 6종 중 최소(직전 최선 다중해상도 0.34m/18.3%). pytest 8 passed(격자와 어긋난 링 위 표고 재현 테스트 추가).
514 lines
25 KiB
Python
514 lines
25 KiB
Python
"""도엽등고선 → 표고 격자 보간 방식 모음 (비교용).
|
|
|
|
문헌(Hutchinson 1988/89 ANUDEM; Chaplot 2006; Arun 2013 등)은 지형 복잡도·자료 밀도에
|
|
따라 우열이 갈리며 단일 최적해가 없다고 본다. 그래서 방식을 하나로 고르지 않고 여기
|
|
모아 두고 B04 화면에서 바꿔 가며 보게 한다(2026-08-30 사용자 지시).
|
|
|
|
각 builder는 `(spec, burned, features, cell_m) -> (R, C) float32` 격자를 돌려준다.
|
|
`burned`는 등고 라인이 구워진 격자(라인 셀 = 표고, 그 외 NaN)다. 폐합 링 안쪽 처리와
|
|
프리뷰·저장은 호출측(`_SheetSurface`)이 방식과 무관하게 똑같이 해 준다.
|
|
"""
|
|
|
|
import logging
|
|
import warnings
|
|
from typing import Any, Callable
|
|
|
|
import numpy as np
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
# TIN(도엽선)이 물어 오는 격자 밖 여유 — 테두리가 삼각망 밖으로 나가지 않을 만큼만.
|
|
SHEET_TIN_CLIP_MARGIN_M = 100.0
|
|
|
|
# 화면 버튼에 쓰는 이름 — 키는 surface_models.generation_params.source_filter 접미사다.
|
|
SHEET_METHOD_LABELS: dict[str, str] = {
|
|
"tin_sheet": "TIN(도엽선)",
|
|
"tin": "TIN 격자",
|
|
"biharmonic": "TPS(박판)",
|
|
"anudem": "ANUDEM형",
|
|
"multires": "다중해상도",
|
|
"laplace": "라플라스",
|
|
}
|
|
|
|
|
|
def _laplacian(values: np.ndarray) -> np.ndarray:
|
|
padded = np.pad(values, 1, mode="edge")
|
|
return (
|
|
padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:] - 4.0 * values
|
|
)
|
|
|
|
|
|
def _contour_vertices(burned: np.ndarray, spec: Any) -> tuple[np.ndarray, np.ndarray]:
|
|
"""등고 라인 셀을 (N,2) 세계좌표와 표고로 바꾼다."""
|
|
rows, cols = np.nonzero(np.isfinite(burned))
|
|
xs = spec.cell_centers_x()[cols]
|
|
ys = spec.cell_centers_y()[rows]
|
|
return np.column_stack([xs, ys]), burned[rows, cols].astype(np.float64)
|
|
|
|
|
|
def _grid_points(spec: Any) -> tuple[np.ndarray, np.ndarray]:
|
|
grid_x, grid_y = np.meshgrid(spec.cell_centers_x(), spec.cell_centers_y())
|
|
return grid_x, grid_y
|
|
|
|
|
|
# ── ① 거리 비례 (버튼에서는 뺐지만 다른 방식의 초기추정으로 계속 쓴다) ───────
|
|
def build_distance(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""가장 가까운 서로 다른 표고 두 라인 사이를 거리 비례로 나눈다.
|
|
|
|
z = (L1·d2 + L2·d1) / (d1 + d2)
|
|
|
|
지도 제작의 고전적 손보간을 그대로 옮긴 것이다. 원뿔·능선(z=r) 형상을 정확히
|
|
재현하고 계단이 생기지 않는다. 표고별 거리장을 돌며 가장 작은 두 값을 추적하므로
|
|
L1≠L2가 보장된다.
|
|
"""
|
|
from scipy.ndimage import distance_transform_edt
|
|
|
|
levels = np.unique(burned[np.isfinite(burned)])
|
|
if len(levels) < 2:
|
|
return np.full(burned.shape, levels[0] if len(levels) else np.nan, dtype=np.float32)
|
|
|
|
infinity = np.float32(np.inf)
|
|
first_d = np.full(burned.shape, infinity, dtype=np.float32)
|
|
first_z = np.zeros(burned.shape, dtype=np.float32)
|
|
second_d = np.full(burned.shape, infinity, dtype=np.float32)
|
|
second_z = np.zeros(burned.shape, dtype=np.float32)
|
|
for level in levels:
|
|
distance = distance_transform_edt(burned != level, sampling=cell_m).astype(np.float32)
|
|
beats_first = distance < first_d
|
|
second_d = np.where(beats_first, first_d, second_d)
|
|
second_z = np.where(beats_first, first_z, second_z)
|
|
first_d = np.where(beats_first, distance, first_d)
|
|
first_z = np.where(beats_first, np.float32(level), first_z)
|
|
beats_second = ~beats_first & (distance < second_d)
|
|
second_d = np.where(beats_second, distance, second_d)
|
|
second_z = np.where(beats_second, np.float32(level), second_z)
|
|
|
|
total = first_d + second_d
|
|
usable = np.isfinite(second_d) & (total > 1e-9)
|
|
surface = first_z.astype(np.float32)
|
|
surface[usable] = (
|
|
(
|
|
first_z[usable].astype(np.float64) * second_d[usable].astype(np.float64)
|
|
+ second_z[usable].astype(np.float64) * first_d[usable].astype(np.float64)
|
|
)
|
|
/ total[usable].astype(np.float64)
|
|
).astype(np.float32)
|
|
return surface
|
|
|
|
|
|
def relax_laplace(surface: np.ndarray, fixed: np.ndarray, iterations: int) -> None:
|
|
"""등고 라인을 고정한 채 이웃 평균으로 다듬는다 (in-place, red-black 순서)."""
|
|
if iterations <= 0:
|
|
return
|
|
free = ~fixed & np.isfinite(surface)
|
|
if not free.any():
|
|
return
|
|
rows, cols = np.indices(surface.shape)
|
|
red = free & (((rows + cols) & 1) == 0)
|
|
black = free & ~red
|
|
padded = np.zeros((surface.shape[0] + 2, surface.shape[1] + 2), dtype=np.float32)
|
|
for _ in range(iterations):
|
|
for colour in (red, black):
|
|
padded[1:-1, 1:-1] = surface
|
|
padded[0, 1:-1] = surface[0]
|
|
padded[-1, 1:-1] = surface[-1]
|
|
padded[1:-1, 0] = surface[:, 0]
|
|
padded[1:-1, -1] = surface[:, -1]
|
|
neighbours = (
|
|
padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:]
|
|
) * np.float32(0.25)
|
|
surface[colour] = neighbours[colour]
|
|
|
|
|
|
# ── ② 라플라스(조화) ─────────────────────────────────────────────────────────
|
|
def build_laplace(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""등고선을 경계값으로 두고 Δz=0을 푼다.
|
|
|
|
면이 매끈해지지만 z=r(원뿔·능선)은 조화함수가 아니라 마루가 눌린다. 비교 기준으로
|
|
남겨 둔다 — ANUDEM이 라플라스 대신 박판 스플라인을 쓰는 이유를 눈으로 보기 위함.
|
|
"""
|
|
surface = build_distance(spec, burned, features, cell_m)
|
|
relax_laplace(surface, np.isfinite(burned), 400)
|
|
return surface
|
|
|
|
|
|
# ── ③ 박판 스플라인(중조화) ──────────────────────────────────────────────────
|
|
def solve_min_curvature(constrained: np.ndarray, guess: np.ndarray) -> np.ndarray:
|
|
"""제약 셀을 고정하고 Δ²z=0(박판 스플라인)을 최소곡률 최소제곱으로 푼다.
|
|
|
|
ANUDEM/Topo to Raster가 쓰는 박판 스플라인과 같은 연산자다. 라플라스와 달리 z=r을
|
|
그대로 통과시켜 능선·마루가 눌리지 않고, 경사가 등고선 너머로 자연스럽게 이어진다.
|
|
|
|
`constrained`: 값이 고정된 셀(등고 라인, 필요하면 구조선 앵커) — 그 외는 NaN.
|
|
`guess`: 시작값이자 감쇠 기준(보통 거리 보간 결과).
|
|
"""
|
|
from scipy.sparse.linalg import LinearOperator, lsmr
|
|
|
|
guess = guess.astype(np.float64)
|
|
burned = constrained
|
|
# 제약 셀 + **격자 테두리**를 고정한다. 최외곽 등고선 바깥이 통째로 자유면
|
|
# 1차함수가 Δ²의 영공간에 남아 해가 하나로 정해지지 않고 켤레기울기가 발산한다
|
|
# (2026-08-30 실측: |Δz| 2092m). 테두리는 거리 보간값으로 묶는다.
|
|
fixed = np.isfinite(burned)
|
|
fixed[0, :] = fixed[-1, :] = True
|
|
fixed[:, 0] = fixed[:, -1] = True
|
|
free = ~fixed & np.isfinite(guess)
|
|
if not free.any():
|
|
return guess.astype(np.float32)
|
|
index = np.flatnonzero(free.ravel())
|
|
values = np.where(np.isfinite(burned), burned, guess)
|
|
base = np.where(fixed, np.nan_to_num(values), 0.0)
|
|
|
|
# Δ²z=0을 정규방정식(CG)으로 풀면 조건수가 격자변 4제곱이라 발산한다(실측).
|
|
# 대신 **최소곡률** 최소제곱으로 세운다 — 자유 셀에 대해 ‖Δz‖를 최소화하며,
|
|
# 그 정상해가 곧 Δ²z=0이다. 조건수가 제곱으로 줄어 LSMR이 안정적으로 푼다
|
|
# (Briggs 1974의 최소곡률 격자화와 같은 목적함수).
|
|
# 최소곡률만으로는 제약(등고선)에서 먼 영역이 정해지지 않아 해가 폭주한다.
|
|
# ANUDEM의 거칠기 벌점과 같은 취지로 감쇠항을 붙여 거리 보간값에 묶어 둔다:
|
|
# minimize ‖Δz‖² + λ‖z − 거리보간‖²
|
|
# λ가 작을수록 더 매끈하고 클수록 거리 보간에 가깝다.
|
|
total_cells = base.size
|
|
free_count = len(index)
|
|
damping = np.float64(np.sqrt(0.02))
|
|
anchor = guess.ravel()[index]
|
|
|
|
def forward(vector: np.ndarray) -> np.ndarray:
|
|
# 반드시 **선형**이어야 한다 — 고정 셀 기여(base)를 여기서 더하면 아핀이 되어
|
|
# LSMR의 전제가 깨지고 해가 폭주한다. base 몫은 우변으로만 넘긴다.
|
|
scattered = np.zeros_like(base)
|
|
scattered.ravel()[index] = vector
|
|
return np.concatenate([_laplacian(scattered).ravel(), damping * vector])
|
|
|
|
def adjoint(vector: np.ndarray) -> np.ndarray:
|
|
curvature = _laplacian(vector[:total_cells].reshape(base.shape)).ravel()[index]
|
|
return curvature + damping * vector[total_cells:]
|
|
|
|
rhs = np.concatenate([-_laplacian(base).ravel(), damping * anchor])
|
|
linear = LinearOperator(
|
|
(total_cells + free_count, free_count),
|
|
matvec=forward,
|
|
rmatvec=adjoint,
|
|
dtype=np.float64,
|
|
)
|
|
result = lsmr(linear, rhs, x0=anchor, maxiter=400, atol=1e-8, btol=1e-8)
|
|
solution, info = result[0], result[1]
|
|
surface = base.copy()
|
|
surface.ravel()[index] = solution
|
|
surface[fixed] = values[fixed]
|
|
# 안전장치 — 발산하면 조용히 틀린 지형을 넘기지 말고 거리 보간으로 되돌린다.
|
|
drift = float(np.nanmax(np.abs(surface - guess)))
|
|
span = float(np.nanmax(guess) - np.nanmin(guess))
|
|
if not np.isfinite(drift) or drift > max(span, 1.0):
|
|
logger.warning(
|
|
"도엽 서피스(TPS): 해가 발산해(최대 %.1fm) 거리 보간으로 되돌립니다 (info=%s).",
|
|
drift,
|
|
info,
|
|
)
|
|
return guess.astype(np.float32)
|
|
logger.info("도엽 서피스(TPS): 최소제곱 info=%s, 최대 변화 %.2fm", info, drift)
|
|
return surface.astype(np.float32)
|
|
|
|
|
|
def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""등고선만 제약으로 둔 박판 스플라인."""
|
|
return solve_min_curvature(burned, build_distance(spec, burned, features, cell_m))
|
|
|
|
|
|
# ── ④ TIN ───────────────────────────────────────────────────────────────────
|
|
def _triangulate(
|
|
spec: Any, shape_: tuple[int, int], points: np.ndarray, values: np.ndarray
|
|
) -> np.ndarray:
|
|
"""정점 구름을 Delaunay 삼각망 선형 보간해 격자로 편다. 삼각망 밖은 NaN."""
|
|
from scipy.interpolate import LinearNDInterpolator
|
|
|
|
if len(points) < 3:
|
|
return np.full(shape_, np.nan, dtype=np.float32)
|
|
interpolator = LinearNDInterpolator(points, values)
|
|
grid_x, grid_y = _grid_points(spec)
|
|
surface = np.empty(shape_, dtype=np.float32)
|
|
chunk = max(1, int(4_000_000 // max(spec.n_cols, 1)))
|
|
for start in range(0, spec.n_rows, chunk):
|
|
stop = min(start + chunk, spec.n_rows)
|
|
surface[start:stop] = interpolator(grid_x[start:stop], grid_y[start:stop]).astype(
|
|
np.float32
|
|
)
|
|
return surface
|
|
|
|
|
|
def build_tin(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""격자에 구운 등고 라인 셀을 Delaunay 삼각망 선형 보간한다.
|
|
|
|
정점이 셀 중심에 맞춰져 있어 1m 계단이 삼각망에 그대로 실린다. 같은 표고 정점
|
|
3개로 이루어진 평탄 삼각형이 굴곡부·마루에 계단을 만든다. 비교 기준으로 남긴다.
|
|
"""
|
|
points, values = _contour_vertices(burned, spec)
|
|
if len(points) > 120_000: # 삼각망 비용은 정점 수에 비례한다
|
|
step = int(np.ceil(len(points) / 120_000))
|
|
points, values = points[::step], values[::step]
|
|
return _triangulate(spec, burned.shape, points, values)
|
|
|
|
|
|
def build_tin_sheet(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""**원본 5m 도엽 등고선 정점**을 그대로 이은 고전 TIN (2026-08-30 사용자 지시).
|
|
|
|
`tin`은 격자에 구운 라인 셀(=1m 계단으로 뭉개진 정점)을 쓰지만, 이쪽은 벡터
|
|
등고선의 정점을 좌표 그대로 쓴다 — 도면 등고선을 삼각망으로 잇는 측량 관행 그대로다.
|
|
격자 밖 등고선은 물지 않는다(삼각망 비용만 커지고 결과는 같다). 다만 격자 테두리가
|
|
삼각망 밖으로 나가지 않도록 여유를 두고 자른다.
|
|
|
|
길이 필터는 두지 않는다 — 짧은 봉우리 폐합 등고선을 버리면 마루가 통째로 평평해진다.
|
|
"""
|
|
from shapely.geometry import shape as to_shape
|
|
|
|
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
|
|
ELEVATION_KEYS,
|
|
iter_linestrings,
|
|
)
|
|
|
|
xs, ys = spec.cell_centers_x(), spec.cell_centers_y()
|
|
margin = max(SHEET_TIN_CLIP_MARGIN_M, cell_m * 2.0)
|
|
x_lo, x_hi = xs[0] - margin, xs[-1] + margin
|
|
y_lo, y_hi = ys[-1] - margin, ys[0] + margin
|
|
|
|
coords: list[np.ndarray] = []
|
|
levels: list[np.ndarray] = []
|
|
for feature in features or []:
|
|
properties = feature.get("properties") or {}
|
|
elevation = next(
|
|
(float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None),
|
|
None,
|
|
)
|
|
geometry = feature.get("geometry")
|
|
if elevation is None or not geometry:
|
|
continue
|
|
try:
|
|
parsed = to_shape(geometry)
|
|
except Exception: # noqa: BLE001 — 손상된 피처는 건너뛴다
|
|
continue
|
|
for line in iter_linestrings(parsed):
|
|
point = np.asarray(line.coords, dtype=np.float64)[:, :2]
|
|
inside = (
|
|
(point[:, 0] >= x_lo)
|
|
& (point[:, 0] <= x_hi)
|
|
& (point[:, 1] >= y_lo)
|
|
& (point[:, 1] <= y_hi)
|
|
)
|
|
if not inside.any():
|
|
continue
|
|
coords.append(point[inside])
|
|
levels.append(np.full(int(inside.sum()), elevation, dtype=np.float64))
|
|
|
|
if not coords:
|
|
logger.warning("도엽 서피스: TIN(도엽선)에 쓸 등고선 정점이 없습니다.")
|
|
return np.full(burned.shape, np.nan, dtype=np.float32)
|
|
|
|
points = np.vstack(coords)
|
|
values = np.concatenate(levels)
|
|
# 도엽 이음매에서 같은 정점이 겹쳐 들어온다 — Qhull 비용만 늘어 미리 접는다.
|
|
_, unique = np.unique(np.round(points, 3), axis=0, return_index=True)
|
|
points, values = points[unique], values[unique]
|
|
logger.info("도엽 서피스: TIN(도엽선) 정점 %d개", len(points))
|
|
return _triangulate(spec, burned.shape, points, values)
|
|
|
|
|
|
# ── ⑦ ANUDEM형 (구조선 + 배수 강제) ─────────────────────────────────────────
|
|
def _contour_corner_anchors(
|
|
spec: Any, burned: np.ndarray, guess: np.ndarray, cell_m: float
|
|
) -> tuple[np.ndarray, np.ndarray] | None:
|
|
"""등고선의 국소 최대 곡률점(코너)에서 능선·계곡 구조선 앵커를 만든다.
|
|
|
|
ANUDEM은 등고선 자체의 곡률에서 능선·계곡망을 먼저 뽑아 흐름 구조를 세운다
|
|
(Hutchinson 1988/89). 여기서도 같은 순서를 따른다.
|
|
|
|
① 라인마다 정점 곡률을 재 국소 최대점(V자 꼭짓점)을 고른다
|
|
② 굽은 안쪽이 더 높으면 **계곡**(등고선 V가 상류를 가리킴), 낮으면 **능선**
|
|
③ 같은 종류의 코너를 이웃 표고끼리 이어 그 사이를 선형 보간해 앵커로 심는다
|
|
|
|
앵커는 박판 해의 제약으로 들어가 계곡 바닥이 이어져 내려가고 능선 마루가 선다.
|
|
"""
|
|
from scipy.ndimage import label
|
|
from scipy.spatial import cKDTree
|
|
|
|
levels = np.unique(burned[np.isfinite(burned)])
|
|
if len(levels) < 2:
|
|
return None
|
|
interval = float(np.diff(levels).min())
|
|
xs = spec.cell_centers_x()
|
|
ys = spec.cell_centers_y()
|
|
|
|
corners: list[tuple[float, float, float, int]] = [] # x, y, level, +1 계곡 / -1 능선
|
|
span = 6 # 곡률을 재는 정점 간격(px) — 짧으면 노이즈, 길면 꼭짓점을 놓친다
|
|
for level in levels:
|
|
labelled, count = label(burned == level)
|
|
for component_id in range(1, count + 1):
|
|
line_rows, line_cols = np.nonzero(labelled == component_id)
|
|
if len(line_rows) < 3 * span:
|
|
continue
|
|
# 라인 셀을 한 줄로 세운다 — 좌표 정렬로 근사한다(정밀 추적은 과하다).
|
|
order = np.argsort(line_cols + line_rows * 1e-3)
|
|
path = np.column_stack([line_cols[order], line_rows[order]]).astype(np.float64)
|
|
before = np.roll(path, span, axis=0)
|
|
after = np.roll(path, -span, axis=0)
|
|
first = path - before
|
|
second = after - path
|
|
first_len = np.hypot(first[:, 0], first[:, 1])
|
|
second_len = np.hypot(second[:, 0], second[:, 1])
|
|
valid = (first_len > 1e-6) & (second_len > 1e-6)
|
|
cosine = np.ones(len(path))
|
|
cosine[valid] = (first[valid] * second[valid]).sum(axis=1) / (
|
|
first_len[valid] * second_len[valid]
|
|
)
|
|
sharp = np.flatnonzero(valid & (cosine < 0.3)) # 70도 이상 꺾인 자리
|
|
for i in sharp[:: max(1, span)]:
|
|
# 굽은 안쪽 방향 = 두 변 단위벡터 합의 반대
|
|
inward = -(first[i] / first_len[i] + second[i] / second_len[i])
|
|
norm = float(np.hypot(inward[0], inward[1]))
|
|
if norm < 1e-6:
|
|
continue
|
|
probe = path[i] + inward / norm * 6.0
|
|
probe_col = int(round(probe[0]))
|
|
probe_row = int(round(probe[1]))
|
|
if not (0 <= probe_row < guess.shape[0] and 0 <= probe_col < guess.shape[1]):
|
|
continue
|
|
inside = float(guess[probe_row, probe_col])
|
|
if not np.isfinite(inside) or abs(inside - level) < interval * 0.15:
|
|
continue
|
|
corners.append(
|
|
(
|
|
float(xs[int(path[i, 0])]),
|
|
float(ys[int(path[i, 1])]),
|
|
float(level),
|
|
1 if inside > level else -1,
|
|
)
|
|
)
|
|
if len(corners) < 4:
|
|
logger.info("도엽 서피스(ANUDEM형): 등고선 코너가 부족해 구조선을 건너뜁니다.")
|
|
return None
|
|
|
|
array = np.asarray(corners, dtype=np.float64)
|
|
anchor_xy: list[np.ndarray] = []
|
|
anchor_z: list[np.ndarray] = []
|
|
reach = interval * 20.0 # 이보다 먼 코너는 같은 구조선으로 보지 않는다
|
|
for level in levels[:-1]:
|
|
upper = level + interval
|
|
lower_set = array[np.abs(array[:, 2] - level) < 1e-6]
|
|
upper_set = array[np.abs(array[:, 2] - upper) < 1e-6]
|
|
if not len(lower_set) or not len(upper_set):
|
|
continue
|
|
tree = cKDTree(upper_set[:, :2])
|
|
distance, index = tree.query(lower_set[:, :2], k=1)
|
|
for i in range(len(lower_set)):
|
|
j = int(index[i])
|
|
if distance[i] > reach or lower_set[i, 3] != upper_set[j, 3]:
|
|
continue
|
|
start = lower_set[i, :2]
|
|
end = upper_set[j, :2]
|
|
steps = max(2, int(distance[i] / max(cell_m, 1e-6) / 4))
|
|
fraction = np.linspace(0.0, 1.0, steps + 1)[1:-1]
|
|
if not len(fraction):
|
|
continue
|
|
anchor_xy.append(start + (end - start) * fraction[:, None])
|
|
anchor_z.append(level + interval * fraction)
|
|
if not anchor_xy:
|
|
return None
|
|
return np.vstack(anchor_xy), np.concatenate(anchor_z)
|
|
|
|
|
|
def _enforce_drainage(surface: np.ndarray, epsilon: float = 0.01) -> int:
|
|
"""가짜 웅덩이를 메운다 — ANUDEM의 배수 강제와 같은 목적.
|
|
|
|
등고선만으로 만든 면에는 흐름이 끊기는 웅덩이가 남는다. 형태학적 재구성(erosion)
|
|
으로 채우되 완전 평탄해지지 않게 아주 작은 값을 얹는다. 채운 셀 수를 반환한다.
|
|
"""
|
|
from skimage.morphology import reconstruction
|
|
|
|
if not np.isfinite(surface).all():
|
|
return 0
|
|
seed = np.full(surface.shape, float(surface.max()), dtype=np.float64)
|
|
seed[0, :] = surface[0, :]
|
|
seed[-1, :] = surface[-1, :]
|
|
seed[:, 0] = surface[:, 0]
|
|
seed[:, -1] = surface[:, -1]
|
|
filled = reconstruction(seed, surface.astype(np.float64), method="erosion")
|
|
raised = filled > surface + 1e-6
|
|
if not raised.any():
|
|
return 0
|
|
# 그냥 채우면 웅덩이가 통째로 평탄해져 흐름 방향이 없어진다. 채운 영역 안쪽으로
|
|
# 갈수록 아주 조금 높아지게 해서 물이 가장자리(넘침점)로 빠져나가게 둔다
|
|
# (Garbrecht·Martz의 평탄면 해소를 간단히 옮긴 것 — 표고 변화는 cm 단위다).
|
|
from scipy.ndimage import distance_transform_edt
|
|
|
|
inner = distance_transform_edt(raised)
|
|
surface[raised] = (filled[raised] + epsilon * inner[raised]).astype(surface.dtype)
|
|
return int(raised.sum())
|
|
|
|
|
|
def build_anudem(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""ANUDEM형 — 등고선 곡률에서 능선·계곡 구조선을 뽑아 제약에 더하고, 박판으로 풀고,
|
|
가짜 웅덩이를 메운다. Topo to Raster가 밟는 세 단계를 그대로 옮긴 것이다."""
|
|
guess = build_distance(spec, burned, features, cell_m).astype(np.float64)
|
|
constrained = burned.astype(np.float64).copy()
|
|
anchors = _contour_corner_anchors(spec, burned, guess, cell_m)
|
|
if anchors is not None:
|
|
xy, z = anchors
|
|
row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy())
|
|
inside = (row >= 0) & (col >= 0)
|
|
row, col, z = row[inside], col[inside], z[inside]
|
|
free = ~np.isfinite(constrained[row, col])
|
|
constrained[row[free], col[free]] = z[free]
|
|
logger.info("도엽 서피스(ANUDEM형): 구조선 앵커 %d셀", int(free.sum()))
|
|
surface = solve_min_curvature(constrained, guess).astype(np.float32)
|
|
logger.info("도엽 서피스(ANUDEM형): 가짜 웅덩이 %d셀 메움", _enforce_drainage(surface))
|
|
return surface
|
|
|
|
|
|
# ── ⑧ 다중해상도 (coarse → fine) ────────────────────────────────────────────
|
|
def build_multires(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
|
"""성긴 격자에서 풀고 점차 세밀화한다 — ANUDEM의 다중해상도 전략.
|
|
|
|
전체 형상은 성긴 격자에서 싸게 잡고, 세밀한 격자에서는 등고선 근처만 다듬는다.
|
|
한 해상도에서만 풀 때보다 넓은 밴드가 고르게 퍼지고 값싸게 수렴한다.
|
|
"""
|
|
from scipy.ndimage import zoom
|
|
|
|
surface: np.ndarray | None = None
|
|
for factor in (8, 4, 2, 1):
|
|
if factor == 1:
|
|
coarse = burned
|
|
else:
|
|
# 성긴 격자의 제약 — 블록 안 등고 라인의 평균 표고를 대표로 쓴다.
|
|
rows = burned.shape[0] // factor * factor
|
|
cols = burned.shape[1] // factor * factor
|
|
blocks = burned[:rows, :cols].reshape(rows // factor, factor, cols // factor, factor)
|
|
# 라인이 하나도 없는 블록은 NaN이 정상이라 경고를 삼킨다.
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
coarse = np.nanmean(blocks, axis=(1, 3)).astype(np.float32)
|
|
level = build_distance(spec, coarse, features, cell_m * factor)
|
|
if surface is not None:
|
|
# 앞 단계 해를 지금 해상도로 올려 절반씩 섞는다 — 성긴 단계의 넓은 추세를
|
|
# 이어받되 이번 해상도의 등고선 정보를 덮지 않는다.
|
|
scale = (level.shape[0] / surface.shape[0], level.shape[1] / surface.shape[1])
|
|
upscaled = zoom(surface, scale, order=1)
|
|
free = ~np.isfinite(coarse)
|
|
level[free] = (level[free] + upscaled[free]) * 0.5
|
|
relax_laplace(level, np.isfinite(coarse), 8)
|
|
surface = level
|
|
assert surface is not None
|
|
if surface.shape != burned.shape: # 블록 자르기로 남은 가장자리 보정
|
|
scale = (burned.shape[0] / surface.shape[0], burned.shape[1] / surface.shape[1])
|
|
surface = zoom(surface, scale, order=1)
|
|
line = np.isfinite(burned)
|
|
surface[line] = burned[line]
|
|
return surface.astype(np.float32)
|
|
|
|
|
|
SHEET_METHOD_BUILDERS: dict[str, Callable[[Any, np.ndarray, Any, float], np.ndarray]] = {
|
|
"tin_sheet": build_tin_sheet,
|
|
"tin": build_tin,
|
|
"biharmonic": build_biharmonic,
|
|
"anudem": build_anudem,
|
|
"multires": build_multires,
|
|
"laplace": build_laplace,
|
|
}
|