eomsangdonandClaude Fable 5
f7528a4aa4
- 폴더·내부 파일 51개 접두사 개명 (git mv, 이력 보존) - 저장소 전체 참조 치환 67파일: import 경로, 라우트 슬러그(b04-preprocess), 라우트 키(B04_PREPROCESS), storage 경로 상수, locale, SQL 주석 - 로직 변경 없음 (기계적 치환). typecheck·백엔드 import 검증 통과 Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
560 lines
23 KiB
Python
560 lines
23 KiB
Python
"""배수유역 해석 격자 생성 — 등고선 TIN 보간 · 웅덩이 채움 · D8 물 방향.
|
||
|
||
지형 근거는 **도엽 등고선**뿐이다. 라이다 DEM은 노선 주변만 커버해 유역 산정에 필요한
|
||
상류 범위를 담지 못하므로 쓰지 않는다(2026-07-31 사용자 지시). 표고점도 쓰지 않는다.
|
||
|
||
처리 순서
|
||
① 계획선 최저점 아래 등고선·짧은 파편 제거
|
||
② 세류선을 노선 교차점에서 잘라 상류측만 남김
|
||
③ 남은 세류선 + 노선을 반경 버퍼한 범위의 bbox로 격자 생성
|
||
④ 등고선 정점 Delaunay TIN 선형보간으로 셀 표고 산출
|
||
⑤ 웅덩이 채움(형태학적 재구성) + 평탄면 미세경사 부여
|
||
⑥ D8(8방향 최급강하) 수신 셀 인덱스 산출
|
||
|
||
⑤가 없으면 등고선 TIN 특유의 가짜 웅덩이·평탄 삼각형에서 흐름이 끊겨 상류 추적이
|
||
도중에 멈춘다. 능선 탐지는 하지 않는다 — 흐름이 도로에 닿는지 여부만으로 유역이 정해진다.
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import logging
|
||
import math
|
||
from dataclasses import dataclass
|
||
from typing import Any
|
||
|
||
import numpy as np
|
||
from affine import Affine
|
||
from rasterio.features import rasterize
|
||
from rasterio.transform import from_origin
|
||
from scipy.interpolate import LinearNDInterpolator
|
||
from scipy.ndimage import distance_transform_edt
|
||
from shapely import segmentize
|
||
from shapely.geometry import LineString, shape
|
||
from shapely.geometry.base import BaseGeometry
|
||
from skimage.morphology import reconstruction
|
||
|
||
from config.config_system import (
|
||
DRAINAGE_CONTOUR_CLIP_MARGIN_M,
|
||
DRAINAGE_CONTOUR_MARGIN_M,
|
||
DRAINAGE_CONTOUR_MIN_LENGTH_M,
|
||
DRAINAGE_CONTOUR_RESAMPLE_M,
|
||
DRAINAGE_FLAT_EPSILON_M,
|
||
DRAINAGE_MAX_GRID_CELLS,
|
||
)
|
||
|
||
logger = logging.getLogger(__name__)
|
||
|
||
# 표고 속성 키: 도엽 등고선(등고수치)·gpkg 등고선(CTRLN_HG) 통합.
|
||
ELEVATION_KEYS = ("등고수치", "CTRLN_HG", "수치", "표고", "높이", "elevation", "ELEV")
|
||
|
||
# D8 이웃 (행 증분, 열 증분, 거리계수). 행은 아래로 증가(북 → 남).
|
||
_NEIGHBORS = (
|
||
(-1, 0, 1.0),
|
||
(1, 0, 1.0),
|
||
(0, -1, 1.0),
|
||
(0, 1, 1.0),
|
||
(-1, -1, math.sqrt(2.0)),
|
||
(-1, 1, math.sqrt(2.0)),
|
||
(1, -1, math.sqrt(2.0)),
|
||
(1, 1, math.sqrt(2.0)),
|
||
)
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class GridSpec:
|
||
"""해석 격자 기하. (0,0) 셀 중심이 (x_min + cell/2, y_max − cell/2)에 놓인다."""
|
||
|
||
x_min: float
|
||
y_max: float
|
||
cell_m: float
|
||
n_rows: int
|
||
n_cols: int
|
||
|
||
@property
|
||
def size(self) -> int:
|
||
return self.n_rows * self.n_cols
|
||
|
||
@property
|
||
def cell_area_m2(self) -> float:
|
||
return self.cell_m * self.cell_m
|
||
|
||
def world_to_rc(self, x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
|
||
"""세계좌표(m)를 격자 행·열로 바꾼다. 범위를 벗어나면 −1을 돌려준다."""
|
||
col = np.floor((x - self.x_min) / self.cell_m).astype(np.int64)
|
||
row = np.floor((self.y_max - y) / self.cell_m).astype(np.int64)
|
||
outside = (col < 0) | (col >= self.n_cols) | (row < 0) | (row >= self.n_rows)
|
||
col[outside] = -1
|
||
row[outside] = -1
|
||
return row, col
|
||
|
||
def cell_centers_x(self) -> np.ndarray:
|
||
return self.x_min + (np.arange(self.n_cols, dtype=np.float64) + 0.5) * self.cell_m
|
||
|
||
def cell_centers_y(self) -> np.ndarray:
|
||
return self.y_max - (np.arange(self.n_rows, dtype=np.float64) + 0.5) * self.cell_m
|
||
|
||
|
||
@dataclass
|
||
class ContourCloud:
|
||
"""등고선에서 뽑은 정점 구름. TIN 보간과 세류 상·하류 판정에 함께 쓴다."""
|
||
|
||
xy: np.ndarray # (N, 2) float64
|
||
z: np.ndarray # (N,) float64
|
||
|
||
@property
|
||
def is_empty(self) -> bool:
|
||
return self.xy.shape[0] < 3
|
||
|
||
|
||
@dataclass
|
||
class TerrainGrid:
|
||
"""격자 지형 해석 결과."""
|
||
|
||
spec: GridSpec
|
||
elevation: np.ndarray # (R, C) float32 — 채움·평탄해소 후 표고, 무효 셀은 NaN
|
||
valid: np.ndarray # (R, C) bool — 등고선 TIN 내부 여부
|
||
receiver: np.ndarray # (R*C,) int32 — D8 수신 셀의 평탄 인덱스, 싱크는 자기 자신
|
||
step_length: np.ndarray # (R*C,) float32 — 수신 셀까지 거리(m), 싱크는 0
|
||
|
||
|
||
# ── ① 등고선 정리 ────────────────────────────────────────────────────────────
|
||
|
||
|
||
def _feature_elevation(properties: dict[str, Any]) -> float | None:
|
||
for key in ELEVATION_KEYS:
|
||
value = properties.get(key)
|
||
if value is None:
|
||
continue
|
||
try:
|
||
return float(value)
|
||
except (TypeError, ValueError):
|
||
continue
|
||
return None
|
||
|
||
|
||
def iter_linestrings(geometry: Any) -> list[LineString]:
|
||
if geometry.geom_type == "LineString":
|
||
return [geometry]
|
||
if geometry.geom_type in {"MultiLineString", "GeometryCollection"}:
|
||
lines: list[LineString] = []
|
||
for part in geometry.geoms:
|
||
lines.extend(iter_linestrings(part))
|
||
return lines
|
||
return []
|
||
|
||
|
||
def build_contour_cloud(
|
||
contour_features: list[dict[str, Any]],
|
||
elevation_floor_m: float | None = None,
|
||
clip_bounds: tuple[float, float, float, float] | None = None,
|
||
) -> ContourCloud:
|
||
"""등고선 피처를 표고가 붙은 정점 구름으로 바꾼다.
|
||
|
||
`elevation_floor_m` 아래 등고선은 계획선 최저점보다 낮아 상류 기여가 불가능하므로
|
||
버린다. 길이가 짧은 파편도 노이즈로 보고 버리되, 임계 이상인 봉우리 폐합 등고선은
|
||
남긴다(봉우리 표고가 사라지면 그 일대 흐름 방향이 통째로 틀어진다).
|
||
|
||
`clip_bounds`(x_min, y_min, x_max, y_max)를 주면 그 밖 등고선은 읽지 않는다. 도엽
|
||
전체 등고선을 다 물고 가면 TIN 삼각망 비용만 커지고 결과는 같다.
|
||
"""
|
||
xs: list[np.ndarray] = []
|
||
ys: list[np.ndarray] = []
|
||
zs: list[np.ndarray] = []
|
||
dropped_low = 0
|
||
dropped_short = 0
|
||
dropped_outside = 0
|
||
for feature in contour_features:
|
||
geometry = feature.get("geometry")
|
||
if not geometry:
|
||
continue
|
||
elevation = _feature_elevation(feature.get("properties") or {})
|
||
if elevation is None:
|
||
continue
|
||
if elevation_floor_m is not None and elevation < elevation_floor_m:
|
||
dropped_low += 1
|
||
continue
|
||
try:
|
||
parsed = shape(geometry)
|
||
except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다
|
||
continue
|
||
if clip_bounds is not None and _outside_bounds(parsed.bounds, clip_bounds):
|
||
dropped_outside += 1
|
||
continue
|
||
for line in iter_linestrings(parsed):
|
||
if line.length < DRAINAGE_CONTOUR_MIN_LENGTH_M:
|
||
dropped_short += 1
|
||
continue
|
||
coords = np.asarray(segmentize(line, DRAINAGE_CONTOUR_RESAMPLE_M).coords)
|
||
if coords.shape[0] < 2:
|
||
continue
|
||
xs.append(coords[:, 0])
|
||
ys.append(coords[:, 1])
|
||
zs.append(np.full(coords.shape[0], elevation, dtype=np.float64))
|
||
if not xs:
|
||
logger.warning(
|
||
"배수유역: 사용할 등고선이 없습니다(저지대 %d, 파편 %d, 범위밖 %d 제외).",
|
||
dropped_low,
|
||
dropped_short,
|
||
dropped_outside,
|
||
)
|
||
return ContourCloud(np.zeros((0, 2)), np.zeros(0))
|
||
xy = np.column_stack((np.concatenate(xs), np.concatenate(ys)))
|
||
z = np.concatenate(zs)
|
||
logger.info(
|
||
"배수유역: 등고선 정점 %d개 (저지대 %d, 파편 %d, 범위밖 %d 제외)",
|
||
xy.shape[0],
|
||
dropped_low,
|
||
dropped_short,
|
||
dropped_outside,
|
||
)
|
||
return ContourCloud(xy, z)
|
||
|
||
|
||
def _outside_bounds(
|
||
bounds: tuple[float, float, float, float], clip: tuple[float, float, float, float]
|
||
) -> bool:
|
||
return bounds[2] < clip[0] or bounds[0] > clip[2] or bounds[3] < clip[1] or bounds[1] > clip[3]
|
||
|
||
|
||
def grid_spec_from_bounds(
|
||
x_min: float,
|
||
y_min: float,
|
||
x_max: float,
|
||
y_max: float,
|
||
cell_m: float,
|
||
anchor_xy: tuple[float, float] | None = None,
|
||
) -> GridSpec:
|
||
"""범위를 덮는 격자를 만든다. `anchor_xy`를 주면 그 점에 셀 모서리를 맞춘다.
|
||
|
||
격자 원점은 **도로 시작점**에 고정한다(2026-07-31 사용자 지시). bbox 좌상단에 맞추면
|
||
1차 영역이 조금만 달라져도 격자가 통째로 밀려 이전 결과와 셀이 대응되지 않는다.
|
||
도로 시작점에 맞추면 반경·영역을 바꿔도 같은 자리의 셀은 같은 자리에 남는다.
|
||
|
||
격자 크기는 절대 자동으로 바꾸지 않는다 — config 값이 그대로 쓰인다(사용자 지시).
|
||
셀 수가 많으면 경고만 남기고 그대로 진행한다.
|
||
"""
|
||
if anchor_xy is not None:
|
||
anchor_x, anchor_y = anchor_xy
|
||
# 앵커에서 셀 정수배만큼 밖으로 나가 범위를 덮는다(넓어질 뿐 좁아지지 않는다).
|
||
x_min = anchor_x - math.ceil((anchor_x - x_min) / cell_m) * cell_m
|
||
y_max = anchor_y + math.ceil((y_max - anchor_y) / cell_m) * cell_m
|
||
n_cols = max(1, int(math.ceil((x_max - x_min) / cell_m)))
|
||
n_rows = max(1, int(math.ceil((y_max - y_min) / cell_m)))
|
||
if n_cols * n_rows > DRAINAGE_MAX_GRID_CELLS:
|
||
logger.warning(
|
||
"배수유역: 격자 %d×%d = %d셀 (%.0fm × %.0fm, 셀 %.2fm) — 권장 상한 %d셀 초과. "
|
||
"그대로 진행합니다. 느리면 DRAINAGE_GRID_SIZE_M을 올리세요.",
|
||
n_rows,
|
||
n_cols,
|
||
n_rows * n_cols,
|
||
n_cols * cell_m,
|
||
n_rows * cell_m,
|
||
cell_m,
|
||
DRAINAGE_MAX_GRID_CELLS,
|
||
)
|
||
return GridSpec(x_min=x_min, y_max=y_max, cell_m=cell_m, n_rows=n_rows, n_cols=n_cols)
|
||
|
||
|
||
def expand_grid_spec(spec: GridSpec, sides: dict[str, bool], step_m: float) -> GridSpec:
|
||
"""활성 셀이 닿은 방향으로만 격자를 넓힌다.
|
||
|
||
셀 정수배로만 넓혀 격자 격자점(도로 시작점 기준)이 그대로 유지되게 한다.
|
||
"""
|
||
steps = max(1, int(math.ceil(step_m / spec.cell_m)))
|
||
west = steps if sides.get("west") else 0
|
||
east = steps if sides.get("east") else 0
|
||
north = steps if sides.get("north") else 0
|
||
south = steps if sides.get("south") else 0
|
||
return GridSpec(
|
||
x_min=spec.x_min - west * spec.cell_m,
|
||
y_max=spec.y_max + north * spec.cell_m,
|
||
cell_m=spec.cell_m,
|
||
n_rows=spec.n_rows + north + south,
|
||
n_cols=spec.n_cols + west + east,
|
||
)
|
||
|
||
|
||
def grid_transform(spec: GridSpec) -> Affine:
|
||
"""rasterio 아핀 변환. 행 0이 북쪽(y_max)이다."""
|
||
return from_origin(spec.x_min, spec.y_max, spec.cell_m, spec.cell_m)
|
||
|
||
|
||
def build_cell_mask(spec: GridSpec, geometry: BaseGeometry) -> np.ndarray:
|
||
"""영역에 **조금이라도 걸치는** 셀만 True인 (rows, cols) 마스크.
|
||
|
||
`all_touched=True`라서 셀이 영역과 한 점만 스쳐도 생성 대상이 된다(사용자 지시).
|
||
"""
|
||
if geometry is None or geometry.is_empty:
|
||
return np.zeros((spec.n_rows, spec.n_cols), dtype=bool)
|
||
burned = rasterize(
|
||
[(geometry, 1)],
|
||
out_shape=(spec.n_rows, spec.n_cols),
|
||
transform=grid_transform(spec),
|
||
fill=0,
|
||
dtype="uint8",
|
||
all_touched=True,
|
||
)
|
||
return burned.astype(bool)
|
||
|
||
|
||
def mask_row_spans(mask: np.ndarray) -> list[tuple[int, int, int]]:
|
||
"""마스크를 행별 연속 구간 [행, 시작열, 끝열(포함)]으로 압축한다.
|
||
|
||
셀 수십만 개를 그대로 내보낼 수 없으니 구간으로 줄인다 — 프론트는 이 구간만 받아
|
||
실제 셀 사각형을 그린다.
|
||
"""
|
||
spans: list[tuple[int, int, int]] = []
|
||
for row in range(mask.shape[0]):
|
||
line = mask[row]
|
||
if not line.any():
|
||
continue
|
||
padded = np.concatenate(([False], line, [False]))
|
||
edges = np.flatnonzero(padded[1:] != padded[:-1])
|
||
for start, stop in zip(edges[0::2], edges[1::2]):
|
||
spans.append((row, int(start), int(stop) - 1))
|
||
return spans
|
||
|
||
|
||
# ── ④ TIN 보간 ──────────────────────────────────────────────────────────────
|
||
|
||
|
||
def interpolate_elevation(spec: GridSpec, cloud: ContourCloud) -> np.ndarray:
|
||
"""등고선 정점 Delaunay TIN으로 셀 표고를 선형보간한다. 외부는 NaN.
|
||
|
||
삼각망 비용은 정점 수에 비례한다. 격자 밖 정점으로 만든 삼각형은 어차피 쓰이지 않으므로
|
||
격자 범위 + 여유만큼만 남기고 잘라낸다 — 결과 표고는 그대로고 속도만 는다.
|
||
"""
|
||
surface = np.full((spec.n_rows, spec.n_cols), np.nan, dtype=np.float32)
|
||
if cloud.is_empty:
|
||
return surface
|
||
margin = DRAINAGE_CONTOUR_CLIP_MARGIN_M
|
||
inside = (
|
||
(cloud.xy[:, 0] >= spec.x_min - margin)
|
||
& (cloud.xy[:, 0] <= spec.x_min + spec.n_cols * spec.cell_m + margin)
|
||
& (cloud.xy[:, 1] >= spec.y_max - spec.n_rows * spec.cell_m - margin)
|
||
& (cloud.xy[:, 1] <= spec.y_max + margin)
|
||
)
|
||
if inside.sum() < 3:
|
||
logger.warning("배수유역: 격자 범위 안에 등고선 정점이 없습니다.")
|
||
return surface
|
||
logger.info("배수유역: TIN 정점 %d개 사용 (전체 %d개)", int(inside.sum()), cloud.xy.shape[0])
|
||
interpolator = LinearNDInterpolator(cloud.xy[inside], cloud.z[inside])
|
||
xs = spec.cell_centers_x()
|
||
ys = spec.cell_centers_y()
|
||
# 행 묶음 단위로 평가해 (행×열) 좌표 배열을 한 번에 들고 있지 않게 한다.
|
||
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)
|
||
grid_x, grid_y = np.meshgrid(xs, ys[start:stop])
|
||
surface[start:stop] = interpolator(grid_x, grid_y).astype(np.float32)
|
||
return surface
|
||
|
||
|
||
# ── ⑤ 웅덩이 채움 + 평탄면 해소 ─────────────────────────────────────────────
|
||
|
||
|
||
def condition_surface(
|
||
surface: np.ndarray, domain: np.ndarray | None = None
|
||
) -> tuple[np.ndarray, np.ndarray]:
|
||
"""가짜 웅덩이를 채우고 평탄면에 미세 경사를 준다.
|
||
|
||
등고선 TIN은 같은 표고 정점 3개로 이루어진 평탄 삼각형과 계단형 가짜 웅덩이를
|
||
필연적으로 만든다. 그대로 D8을 돌리면 흐름이 거기서 끊겨 상류 추적이 멈춘다.
|
||
|
||
채움은 형태학적 재구성(erosion)으로 한다. 배출구는 격자 최외곽과 유효 영역 경계(무효
|
||
셀에 맞닿은 유효 셀)로 둔다 — 그래야 유효 영역 전체가 하나의 평탄면으로 잠기지 않는다.
|
||
|
||
`domain`을 주면 그 안쪽만 해석 대상으로 삼는다(1차 영역에 걸쳐 실제 생성된 셀 마스크).
|
||
"""
|
||
valid = np.isfinite(surface)
|
||
if domain is not None:
|
||
valid &= domain
|
||
if not valid.any():
|
||
return surface, valid
|
||
ceiling = float(np.nanmax(surface)) + 1000.0
|
||
mask = np.where(valid, surface, ceiling).astype(np.float32)
|
||
|
||
open_boundary = np.zeros_like(valid)
|
||
open_boundary[0, :] = True
|
||
open_boundary[-1, :] = True
|
||
open_boundary[:, 0] = True
|
||
open_boundary[:, -1] = True
|
||
open_boundary |= _dilate(~valid) & valid
|
||
open_boundary &= valid
|
||
if not open_boundary.any():
|
||
open_boundary = valid & _dilate(~valid)
|
||
|
||
seed = np.full_like(mask, ceiling)
|
||
seed[open_boundary] = mask[open_boundary]
|
||
filled = reconstruction(seed, mask, method="erosion", footprint=np.ones((3, 3), dtype=bool))
|
||
filled = filled.astype(np.float32)
|
||
|
||
# 채움 뒤 더 낮은 이웃이 없는 셀 = 평탄면. 가장 가까운 비평탄 셀 쪽으로 미세 경사를 준다.
|
||
flat = valid & ~_has_lower_neighbour(filled, valid)
|
||
if flat.any():
|
||
distance = distance_transform_edt(flat).astype(np.float32)
|
||
filled = filled + distance * np.float32(DRAINAGE_FLAT_EPSILON_M)
|
||
filled[~valid] = np.nan
|
||
return filled, valid
|
||
|
||
|
||
def _dilate(mask: np.ndarray) -> np.ndarray:
|
||
"""8이웃 1스텝 팽창(외부는 False)."""
|
||
padded = np.zeros((mask.shape[0] + 2, mask.shape[1] + 2), dtype=bool)
|
||
padded[1:-1, 1:-1] = mask
|
||
result = np.zeros_like(mask)
|
||
for row_shift in (0, 1, 2):
|
||
for col_shift in (0, 1, 2):
|
||
result |= padded[
|
||
row_shift : row_shift + mask.shape[0], col_shift : col_shift + mask.shape[1]
|
||
]
|
||
return result
|
||
|
||
|
||
def _has_lower_neighbour(surface: np.ndarray, valid: np.ndarray) -> np.ndarray:
|
||
"""8이웃 중 자기보다 낮은 셀이 하나라도 있는지. 무효 셀은 +∞로 보아 제외한다."""
|
||
rows, cols = surface.shape
|
||
padded = np.full((rows + 2, cols + 2), np.inf, dtype=np.float32)
|
||
padded[1:-1, 1:-1] = np.where(valid, surface, np.inf)
|
||
result = np.zeros((rows, cols), dtype=bool)
|
||
center = padded[1:-1, 1:-1]
|
||
for row_shift, col_shift, _ in _NEIGHBORS:
|
||
neighbour = padded[
|
||
1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols
|
||
]
|
||
result |= neighbour < center
|
||
return result & valid
|
||
|
||
|
||
# ── ⑥ D8 물 방향 ────────────────────────────────────────────────────────────
|
||
|
||
|
||
def compute_receivers(
|
||
spec: GridSpec, surface: np.ndarray, valid: np.ndarray
|
||
) -> tuple[np.ndarray, np.ndarray]:
|
||
"""셀마다 8방향 최급강하 이웃(수신 셀)을 정한다.
|
||
|
||
돌려주는 `receiver`는 평탄 인덱스(row * n_cols + col)다. 더 낮은 이웃이 없는 셀(싱크)과
|
||
무효 셀은 자기 자신을 가리켜 흐름이 그 자리에서 멈춘다.
|
||
"""
|
||
rows, cols = spec.n_rows, spec.n_cols
|
||
padded = np.full((rows + 2, cols + 2), np.inf, dtype=np.float32)
|
||
padded[1:-1, 1:-1] = np.where(valid, surface, np.inf)
|
||
center = padded[1:-1, 1:-1]
|
||
|
||
flat_index = np.arange(rows * cols, dtype=np.int32).reshape(rows, cols)
|
||
padded_index = np.full((rows + 2, cols + 2), -1, dtype=np.int32)
|
||
padded_index[1:-1, 1:-1] = flat_index
|
||
|
||
best_slope = np.zeros((rows, cols), dtype=np.float32)
|
||
receiver = flat_index.copy()
|
||
step = np.zeros((rows, cols), dtype=np.float32)
|
||
|
||
for row_shift, col_shift, factor in _NEIGHBORS:
|
||
neighbour = padded[
|
||
1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols
|
||
]
|
||
distance = np.float32(factor * spec.cell_m)
|
||
# 무효 셀끼리는 ∞−∞ = NaN이 되지만 아래 isfinite에서 걸러진다.
|
||
with np.errstate(invalid="ignore"):
|
||
slope = (center - neighbour) / distance
|
||
better = np.isfinite(slope) & (slope > best_slope)
|
||
if not better.any():
|
||
continue
|
||
best_slope = np.where(better, slope, best_slope)
|
||
neighbour_index = padded_index[
|
||
1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols
|
||
]
|
||
receiver = np.where(better, neighbour_index, receiver)
|
||
step = np.where(better, distance, step)
|
||
return receiver.reshape(-1), step.reshape(-1)
|
||
|
||
|
||
# ── 오케스트레이션 ──────────────────────────────────────────────────────────
|
||
|
||
|
||
def build_terrain_grid(
|
||
spec: GridSpec, cloud: ContourCloud, domain: np.ndarray | None = None
|
||
) -> TerrainGrid:
|
||
"""격자 범위와 등고선 구름으로 지형 해석 격자를 만든다.
|
||
|
||
`domain`은 실제 해석할 셀 마스크(1차 영역에 걸친 셀). 주면 그 밖은 무효로 둔다.
|
||
"""
|
||
surface = interpolate_elevation(spec, cloud)
|
||
conditioned, valid = condition_surface(surface, domain)
|
||
receiver, step = compute_receivers(spec, conditioned, valid)
|
||
logger.info(
|
||
"배수유역: 격자 %d×%d (%.2fm), 해석 대상 셀 %d개",
|
||
spec.n_rows,
|
||
spec.n_cols,
|
||
spec.cell_m,
|
||
int(valid.sum()),
|
||
)
|
||
return TerrainGrid(
|
||
spec=spec, elevation=conditioned, valid=valid, receiver=receiver, step_length=step
|
||
)
|
||
|
||
|
||
# 화살표 방위 분해능. 0 = 화면상 오른쪽(+열), 시계방향으로 증가(행이 아래로 증가하므로).
|
||
AZIMUTH_STEPS = 32
|
||
# 방위 코드 특수값.
|
||
AZIMUTH_SINK = AZIMUTH_STEPS # 32 = 제자리(더 낮은 이웃 없음)
|
||
AZIMUTH_INVALID = AZIMUTH_STEPS + 1 # 33 = 표고 없음(해석 불가)
|
||
|
||
|
||
def descent_azimuth(
|
||
spec: GridSpec,
|
||
surface: np.ndarray,
|
||
valid: np.ndarray,
|
||
receiver: np.ndarray,
|
||
forced: np.ndarray | None = None,
|
||
) -> np.ndarray:
|
||
"""셀별 물 흐름 방위를 32방위 코드로 낸다.
|
||
|
||
D8은 연결(도로 도달 판정)에는 충분하지만 화면에 8방위밖에 못 그린다. 실제 지표수는
|
||
지형 최급강하 방향으로 흐르고 그 방향은 연속값이므로, **표시는 지표면 기울기에서 뽑은
|
||
연속 방위를 32단계로 양자화**해 보여 준다(2026-07-31 사용자 지시).
|
||
|
||
`forced`(세류망을 따라 흐름을 새긴 셀)는 기울기 대신 실제 수신 셀 방향을 쓴다 — 그
|
||
셀들은 지형 추정이 아니라 확정된 물길을 따르기 때문이다. 기울기가 0에 가까운 셀도
|
||
수신 셀 방향으로 대체한다.
|
||
"""
|
||
rows, cols = spec.n_rows, spec.n_cols
|
||
filled = np.where(valid, surface, np.nan)
|
||
# np.gradient는 NaN이 번지므로 무효 셀을 주변 유효값으로 임시 대체한 뒤 기울기를 잡는다.
|
||
working = np.where(np.isfinite(filled), filled, np.nanmean(filled) if valid.any() else 0.0)
|
||
grad_row, grad_col = np.gradient(working.astype(np.float64), spec.cell_m)
|
||
# 내리막 방향 = 기울기 반대. 행은 아래로 증가하므로 화면 좌표와 부호가 같다.
|
||
move_row = -grad_row
|
||
move_col = -grad_col
|
||
magnitude = np.hypot(move_row, move_col)
|
||
|
||
index = np.arange(receiver.size, dtype=np.int64)
|
||
receiver_row = (receiver // cols - index // cols).reshape(rows, cols).astype(np.float64)
|
||
receiver_col = (receiver % cols - index % cols).reshape(rows, cols).astype(np.float64)
|
||
use_receiver = magnitude < 1e-9
|
||
if forced is not None:
|
||
use_receiver |= forced.reshape(rows, cols)
|
||
move_row = np.where(use_receiver, receiver_row, move_row)
|
||
move_col = np.where(use_receiver, receiver_col, move_col)
|
||
|
||
angle = np.arctan2(move_row, move_col)
|
||
code = np.rint(angle / (2.0 * math.pi / AZIMUTH_STEPS)).astype(np.int16) % AZIMUTH_STEPS
|
||
# 수신 셀이 자기 자신이거나 이동량이 없는 셀은 방향이 없다.
|
||
is_sink = (receiver.reshape(rows, cols) == index.reshape(rows, cols)) | (
|
||
(np.abs(move_row) < 1e-12) & (np.abs(move_col) < 1e-12)
|
||
)
|
||
code = np.where(is_sink, AZIMUTH_SINK, code)
|
||
# 세류망을 따라 흐름을 새긴 셀은 표고가 없어도 방향이 확정돼 있다.
|
||
known = valid if forced is None else (valid | forced.reshape(rows, cols))
|
||
code = np.where(known, code, AZIMUTH_INVALID)
|
||
return code.reshape(-1).astype(np.int16)
|
||
|
||
|
||
def route_elevation_floor(route_z_values: list[float]) -> float | None:
|
||
"""계획선 최저점에서 여유를 뺀 등고선 하한. 값이 없으면 None(필터 미적용)."""
|
||
finite = [value for value in route_z_values if math.isfinite(value) and value != 0.0]
|
||
if not finite:
|
||
return None
|
||
return min(finite) - DRAINAGE_CONTOUR_MARGIN_M
|