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>
583 lines
24 KiB
Python
583 lines
24 KiB
Python
"""배수유역 흐름 해석 — 도로 굽기 · 상류 추적(포인터 더블링) · 유역 폴리곤화.
|
|
|
|
핵심은 하나다: **물길을 따라가 도로에 닿는 셀만 유역이다.**
|
|
셀마다 D8 수신 셀을 따라가 종착점(root)을 구하고, 그 종착점이 도로 셀이면 활성이다.
|
|
|
|
이 방식은 능선을 따로 찾지 않는다. 능선 너머 셀의 물은 다른 계곡으로 빠져 도로에
|
|
닿지 못하므로 자동으로 비활성이 되고, 그 경계선이 곧 능선이다. 유역 안쪽 봉우리는
|
|
물이 결국 도로로 흘러 자동으로 포함된다.
|
|
|
|
종착점은 세부유역 라벨의 근거로도 그대로 쓴다 — 셀이 도달한 도로 셀이 정해지면
|
|
그 도로 셀을 담당하는 관이 곧 그 셀의 유역 번호다. 관 배치가 바뀌어도 격자 해석을
|
|
다시 돌릴 필요 없이 "도로 셀 → 관" 대응만 다시 계산하면 된다.
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import logging
|
|
from dataclasses import dataclass
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
from rasterio.features import rasterize, shapes
|
|
from scipy.spatial import cKDTree
|
|
from shapely.geometry import LineString, MultiPolygon, Polygon, shape
|
|
from shapely.ops import unary_union
|
|
|
|
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
|
|
AZIMUTH_STEPS,
|
|
ContourCloud,
|
|
GridSpec,
|
|
TerrainGrid,
|
|
build_cell_mask,
|
|
build_terrain_grid,
|
|
descent_azimuth,
|
|
expand_grid_spec,
|
|
grid_transform,
|
|
)
|
|
from config.config_system import (
|
|
DRAINAGE_EXPAND_STEP_M,
|
|
DRAINAGE_MAX_EXPAND_ROUNDS,
|
|
DRAINAGE_MIN_BASIN_AREA_M2,
|
|
DRAINAGE_POLYGON_SIMPLIFY_M,
|
|
DRAINAGE_ROAD_WIDTH_M,
|
|
)
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
# 포인터 더블링 반복 상한. 한 번에 경로 길이가 2배가 되므로 2^40 스텝이면 어떤 격자도 덮는다.
|
|
_MAX_DOUBLING_ROUNDS = 40
|
|
|
|
# 8이웃 (행 증분, 열 증분) — 최외곽 판정용. 거리는 쓰지 않는다.
|
|
_NEIGHBOR_SHIFTS = (
|
|
(-1, 0, 1.0),
|
|
(1, 0, 1.0),
|
|
(0, -1, 1.0),
|
|
(0, 1, 1.0),
|
|
(-1, -1, 1.0),
|
|
(-1, 1, 1.0),
|
|
(1, -1, 1.0),
|
|
(1, 1, 1.0),
|
|
)
|
|
|
|
|
|
@dataclass
|
|
class RoadRaster:
|
|
"""격자에 구운 도로. 도로 셀은 흐름을 흡수하는 종착점이 된다."""
|
|
|
|
mask: np.ndarray # (R, C) bool
|
|
cell_index: np.ndarray # (K,) int32 — 도로 셀의 평탄 인덱스
|
|
chainage: np.ndarray # (K,) float64 — 도로 셀의 누가거리(m)
|
|
slot_of_cell: np.ndarray # (R*C,) int32 — 도로 셀이면 K 내 위치, 아니면 −1
|
|
|
|
@property
|
|
def count(self) -> int:
|
|
return int(self.cell_index.size)
|
|
|
|
|
|
@dataclass
|
|
class FlowResult:
|
|
"""상류 추적 결과."""
|
|
|
|
root: np.ndarray # (R*C,) int32 — 흐름 종착 셀의 평탄 인덱스
|
|
road_slot: np.ndarray # (R*C,) int32 — 도달한 도로 셀 슬롯, 도달 못하면 −1
|
|
active: np.ndarray # (R, C) bool — 도로에 물이 닿는 셀
|
|
path_length: np.ndarray # (R*C,) float32 — 종착점까지 물길 길이(m)
|
|
strength: np.ndarray # (K,) int64 — 도로 셀별 상류 셀 수(흐름 강도)
|
|
|
|
|
|
# ── 도로 굽기 ───────────────────────────────────────────────────────────────
|
|
|
|
|
|
def rasterize_road(
|
|
spec: GridSpec,
|
|
route_line: LineString,
|
|
width_m: float = DRAINAGE_ROAD_WIDTH_M,
|
|
) -> RoadRaster:
|
|
"""노선을 노폭만큼 두껍게 격자에 굽고, 각 도로 셀에 누가거리를 붙인다.
|
|
|
|
폭을 주는 이유는 실제 노면이 물을 받기 때문이기도 하지만, 1셀 선으로 구우면 D8
|
|
대각 이동이 도로를 건너뛰어 상류 물이 도로를 지나쳐 버리기 때문이다. 3셀 이상 두께면
|
|
내리막 물길이 반드시 도로 셀을 한 번은 밟는다.
|
|
"""
|
|
half_width = max(width_m / 2.0, spec.cell_m)
|
|
burned = rasterize(
|
|
[(route_line.buffer(half_width), 1)],
|
|
out_shape=(spec.n_rows, spec.n_cols),
|
|
transform=grid_transform(spec),
|
|
fill=0,
|
|
dtype="uint8",
|
|
all_touched=True,
|
|
).astype(bool)
|
|
|
|
slot_of_cell = np.full(spec.size, -1, dtype=np.int32)
|
|
cell_index = np.flatnonzero(burned.reshape(-1)).astype(np.int32)
|
|
if cell_index.size == 0:
|
|
logger.warning("배수유역: 노선이 격자 범위 밖입니다 — 도로 셀 0개.")
|
|
return RoadRaster(burned, cell_index, np.zeros(0), slot_of_cell)
|
|
|
|
# 도로 셀 누가거리는 노선을 촘촘히 샘플링해 가장 가까운 샘플의 누가거리로 준다.
|
|
step = max(spec.cell_m / 2.0, 0.25)
|
|
positions = np.arange(0.0, route_line.length + step, step)
|
|
samples = np.array([list(route_line.interpolate(p).coords)[0] for p in positions])
|
|
rows = (cell_index // spec.n_cols).astype(np.float64)
|
|
cols = (cell_index % spec.n_cols).astype(np.float64)
|
|
centers = np.column_stack(
|
|
(
|
|
spec.x_min + (cols + 0.5) * spec.cell_m,
|
|
spec.y_max - (rows + 0.5) * spec.cell_m,
|
|
)
|
|
)
|
|
_, nearest = cKDTree(samples).query(centers)
|
|
chainage = np.minimum(positions[nearest], route_line.length)
|
|
slot_of_cell[cell_index] = np.arange(cell_index.size, dtype=np.int32)
|
|
return RoadRaster(burned, cell_index, chainage, slot_of_cell)
|
|
|
|
|
|
# ── 상류 추적 ───────────────────────────────────────────────────────────────
|
|
|
|
|
|
def trace_flow(terrain: TerrainGrid, road: RoadRaster) -> FlowResult:
|
|
"""모든 셀의 물길 종착점을 구하고 도로 도달 여부(=유역 포함 여부)를 판정한다.
|
|
|
|
포인터 더블링으로 한 번에 경로 길이를 2배씩 늘려 종착점을 찾는다. 채움·평탄해소를
|
|
거친 표고에서는 흐름을 따라 표고가 단조 감소하므로 순환이 없고, 반복은 항상 끝난다.
|
|
"""
|
|
spec = terrain.spec
|
|
receiver = terrain.receiver.astype(np.int32, copy=True)
|
|
step_length = terrain.step_length.astype(np.float32, copy=True)
|
|
|
|
# 도로 셀은 흐름을 흡수한다 — 물이 도로에 닿으면 거기서 끝난다.
|
|
receiver[road.cell_index] = road.cell_index
|
|
step_length[road.cell_index] = 0.0
|
|
|
|
jump = receiver
|
|
path_length = step_length
|
|
for _ in range(_MAX_DOUBLING_ROUNDS):
|
|
next_jump = jump[jump]
|
|
if np.array_equal(next_jump, jump):
|
|
break
|
|
path_length = path_length + path_length[jump]
|
|
jump = next_jump
|
|
|
|
road_slot = road.slot_of_cell[jump]
|
|
active_flat = road_slot >= 0
|
|
strength = (
|
|
np.bincount(road_slot[active_flat], minlength=max(road.count, 1)).astype(np.int64)
|
|
if road.count
|
|
else np.zeros(0, dtype=np.int64)
|
|
)
|
|
logger.info(
|
|
"배수유역: 활성 셀 %d / %d (도로 셀 %d)", int(active_flat.sum()), spec.size, road.count
|
|
)
|
|
return FlowResult(
|
|
root=jump,
|
|
road_slot=road_slot,
|
|
active=active_flat.reshape(spec.n_rows, spec.n_cols),
|
|
path_length=path_length,
|
|
strength=strength,
|
|
)
|
|
|
|
|
|
def burn_stream_flow(
|
|
terrain: TerrainGrid, road: RoadRaster, streams: list[LineString]
|
|
) -> tuple[TerrainGrid, np.ndarray]:
|
|
"""확정된 상류 세류망을 따라 격자 흐름 방향을 강제로 새긴다.
|
|
|
|
등고선 TIN 보간면은 실제 물골(thalweg)을 그대로 재현하지 못한다. 그래서 세류선 위
|
|
셀인데도 D8이 옆 사면으로 흘려보내 도로에 닿지 못하는 일이 생긴다. 세류선은 이미
|
|
"도로를 건너 하류로 빠지는 물길"로 확정된 자료이므로, 그 위 셀의 흐름 방향은 추정할
|
|
것이 아니라 **그대로 따라야 한다**(2026-07-31 사용자 지시).
|
|
|
|
세류선 위 셀은 물길 방향의 다음 셀을 수신 셀로 삼는다. 그러면 세류선 셀은 물론,
|
|
세류선으로 흘러드는 사면 셀까지 전부 도로에 도달한다. 도로 셀은 흡수점이므로 수신 셀은
|
|
건드리지 않되, 세류 셀 목록에는 포함한다.
|
|
|
|
**표고 유무를 따지지 않는다.** 세류선은 확정된 자료이므로 등고선 TIN 껍질 밖이라
|
|
표고가 없는 셀이라도 물이 지나간다는 사실은 변하지 않는다. 표고를 조건으로 걸면 그런
|
|
셀이 새김에서 빠져 파랑·회색으로 남는다.
|
|
|
|
돌려주는 값: (흐름이 새겨진 지형, 세류선이 지나는 셀 마스크).
|
|
"""
|
|
spec = terrain.spec
|
|
receiver = terrain.receiver.copy()
|
|
step_length = terrain.step_length.copy()
|
|
road_cells = road.mask.reshape(-1)
|
|
burned = np.zeros(spec.size, dtype=bool)
|
|
|
|
tails: list[int] = []
|
|
# 셀이 몇 갈래에 속하는지 — 하류 끝이 합류 지점인지 판단하는 근거.
|
|
visits = np.zeros(spec.size, dtype=np.int32)
|
|
for line in streams:
|
|
chain = _line_cell_chain(spec, line)
|
|
if not chain:
|
|
continue
|
|
# 사슬의 **모든** 셀을 세류 셀로 표시한다 — 마지막 셀도 물길 위다.
|
|
burned[chain] = True
|
|
np.add.at(visits, np.unique(chain), 1)
|
|
for current, following in zip(chain, chain[1:]):
|
|
if road_cells[current] or not _is_neighbour(spec, current, following):
|
|
continue
|
|
receiver[current] = following
|
|
step_length[current] = _cell_distance(spec, current, following)
|
|
tails.append(chain[-1])
|
|
|
|
# 하류 끝이 도로에도 닿지 않고 다른 갈래와도 겹치지 않으면 그 갈래는 떠 있는 것이다.
|
|
# 지류가 본류 중간에 합류하는 경우는 끝 셀을 두 갈래가 공유하므로 정상이다.
|
|
detached = sum(1 for tail in tails if not road_cells[tail] and visits[tail] < 2)
|
|
if detached:
|
|
logger.warning(
|
|
"배수유역: 세류망 %d갈래의 하류 끝이 도로·다른 세류 어디에도 닿지 않습니다.", detached
|
|
)
|
|
logger.info("배수유역: 세류망 셀 %d개 표시 (세류 %d갈래)", int(burned.sum()), len(streams))
|
|
return (
|
|
TerrainGrid(
|
|
spec=spec,
|
|
elevation=terrain.elevation,
|
|
valid=terrain.valid,
|
|
receiver=receiver,
|
|
step_length=step_length,
|
|
),
|
|
burned,
|
|
)
|
|
|
|
|
|
def _line_cell_chain(spec: GridSpec, line: LineString) -> list[int]:
|
|
"""선을 따라 지나가는 셀을 순서대로 뽑는다(연속 중복 제거)."""
|
|
step = max(spec.cell_m / 2.0, 0.1)
|
|
positions = np.arange(0.0, line.length + step, step)
|
|
chain: list[int] = []
|
|
for position in positions:
|
|
point = line.interpolate(float(position))
|
|
col = int((point.x - spec.x_min) // spec.cell_m)
|
|
row = int((spec.y_max - point.y) // spec.cell_m)
|
|
if not (0 <= row < spec.n_rows and 0 <= col < spec.n_cols):
|
|
continue
|
|
index = row * spec.n_cols + col
|
|
if not chain or chain[-1] != index:
|
|
chain.append(index)
|
|
return chain
|
|
|
|
|
|
def _is_neighbour(spec: GridSpec, first: int, second: int) -> bool:
|
|
row_delta = abs(first // spec.n_cols - second // spec.n_cols)
|
|
col_delta = abs(first % spec.n_cols - second % spec.n_cols)
|
|
return max(row_delta, col_delta) == 1
|
|
|
|
|
|
def _cell_distance(spec: GridSpec, first: int, second: int) -> float:
|
|
row_delta = abs(first // spec.n_cols - second // spec.n_cols)
|
|
col_delta = abs(first % spec.n_cols - second % spec.n_cols)
|
|
return spec.cell_m * float(np.hypot(row_delta, col_delta))
|
|
|
|
|
|
@dataclass
|
|
class FlowClassification:
|
|
"""셀별 흐름 방향과 도로 도달 여부 — 확장 없이 현재 격자만 본 결과."""
|
|
|
|
direction: np.ndarray # (R*C,) int16 — 32방위 코드(32=제자리, 33=무효)
|
|
reaches_road: np.ndarray # (R*C,) bool — 물길을 따라가면 도로에 닿는가
|
|
analyzed: np.ndarray # (R*C,) bool — 실제로 판정한 셀
|
|
outer_seeds: int # 최외곽에서 출발해 판정한 셀 수
|
|
interior_seeds: int # 최외곽 추적에 안 걸려 따로 출발시킨 내부 셀 수
|
|
burned: np.ndarray | None = None # (R*C,) bool — 세류망을 따라 흐름을 새긴 셀
|
|
|
|
|
|
def outermost_cells(domain: np.ndarray) -> np.ndarray:
|
|
"""해석 영역의 최외곽 셀 — 영역 밖(또는 격자 밖)에 8이웃이 하나라도 닿는 셀."""
|
|
padded = np.zeros((domain.shape[0] + 2, domain.shape[1] + 2), dtype=bool)
|
|
padded[1:-1, 1:-1] = domain
|
|
exposed = np.zeros_like(domain)
|
|
rows, cols = domain.shape
|
|
for row_shift, col_shift, _ in _NEIGHBOR_SHIFTS:
|
|
neighbour = padded[
|
|
1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols
|
|
]
|
|
exposed |= ~neighbour
|
|
return exposed & domain
|
|
|
|
|
|
def classify_flow(
|
|
terrain: TerrainGrid,
|
|
road: RoadRaster,
|
|
burned: np.ndarray | None = None,
|
|
azimuth: np.ndarray | None = None,
|
|
) -> FlowClassification:
|
|
"""최외곽 셀부터 물길을 따라가며 도로 도달 여부를 판정한다.
|
|
|
|
물이 흐르는 순서 그대로 따라간다 — 셀이 각자 도로를 바라보는 게 아니라, 한 셀에서
|
|
출발해 화살표가 가리키는 다음 셀, 그 셀이 가리키는 그다음 셀로 사슬처럼 이어 간다
|
|
(2026-07-31 사용자 지시).
|
|
|
|
판정 순서는 다음과 같다(2026-07-31 사용자 지시):
|
|
|
|
⓪ **세류선과 겹치는 셀을 먼저 적색으로 확정한다.** 세류망은 이미 "도로를 건너 하류로
|
|
빠지는 물길"로 확정된 자료다. 표고가 있든 없든 물이 지나간다는 사실은 변하지 않으므로
|
|
추적 결과를 기다릴 이유가 없다.
|
|
① 최외곽 셀에서 출발해 사슬을 따라간다.
|
|
② 사슬이 도로 셀이나 세류 셀에 닿으면 사슬 전체가 적색, 싱크에서 멈추거나 해석 영역
|
|
밖으로 나가면 전체가 미도달(파랑 채움 + 백색 화살표)이다.
|
|
③ 이미 색이 정해진 셀을 만나면 **그 셀의 색을 그대로 물려받고** 끝낸다. 판정된 셀은
|
|
다시 분석하지 않는다.
|
|
④ 최외곽 추적에 안 걸린 내부 셀을 그다음에 따로 출발시킨다.
|
|
|
|
`azimuth`를 주면 그 32방위 코드를 그대로 화살표로 쓴다(등고선 하강 방향). 주지 않으면
|
|
지표면 기울기에서 뽑는다. **화살표는 실제 수신 셀과 같은 방향이어야 한다** — 어긋나면
|
|
화살표로 사슬을 따라가는 눈 검증이 성립하지 않는다.
|
|
"""
|
|
spec = terrain.spec
|
|
valid = terrain.valid.reshape(-1)
|
|
receiver = terrain.receiver
|
|
# 도로 셀과 세류망 셀 둘 다 "여기 닿으면 적색"인 종결점이다.
|
|
stream_cells = np.zeros(spec.size, dtype=bool) if burned is None else burned
|
|
absorbing = (road.mask.reshape(-1) & valid) | stream_cells
|
|
|
|
if azimuth is None:
|
|
direction = descent_azimuth(spec, terrain.elevation, terrain.valid, receiver, burned)
|
|
else:
|
|
direction = _azimuth_with_burn(spec, azimuth, receiver, burned)
|
|
reaches = np.zeros(spec.size, dtype=bool)
|
|
# 0=미방문, 1=경로에 올라 있음, 2=판정 완료
|
|
state = np.zeros(spec.size, dtype=np.int8)
|
|
|
|
# ⓪ 세류선과 겹치는 셀을 먼저 적색으로 못박는다.
|
|
reaches[stream_cells] = True
|
|
state[stream_cells] = 2
|
|
stream_done = int(stream_cells.sum())
|
|
|
|
outer = np.flatnonzero(outermost_cells(terrain.valid))
|
|
remaining = np.flatnonzero(valid)
|
|
outer_done = _walk_from(outer, receiver, valid, absorbing, state, reaches)
|
|
interior_done = _walk_from(remaining, receiver, valid, absorbing, state, reaches)
|
|
|
|
analyzed = state == 2
|
|
logger.info(
|
|
"배수유역: 흐름 판정 %d셀 (세류 선확정 %d / 최외곽 출발 %d / 내부 보충 %d) — "
|
|
"도로 도달 %d, 미도달 %d",
|
|
int(analyzed.sum()),
|
|
stream_done,
|
|
outer_done,
|
|
interior_done,
|
|
int((reaches & analyzed).sum()),
|
|
int((~reaches & analyzed).sum()),
|
|
)
|
|
return FlowClassification(
|
|
direction=direction,
|
|
reaches_road=reaches,
|
|
analyzed=analyzed,
|
|
outer_seeds=outer_done,
|
|
interior_seeds=interior_done,
|
|
burned=burned,
|
|
)
|
|
|
|
|
|
def _azimuth_with_burn(
|
|
spec: GridSpec, azimuth: np.ndarray, receiver: np.ndarray, burned: np.ndarray | None
|
|
) -> np.ndarray:
|
|
"""세류망을 따라 흐름을 새긴 셀은 그 수신 셀 방향으로 화살표를 덮어쓴다."""
|
|
direction = azimuth.astype(np.int16, copy=True)
|
|
if burned is None or not burned.any():
|
|
return direction
|
|
index = np.arange(spec.size, dtype=np.int64)
|
|
moved = burned & (receiver != index)
|
|
if not moved.any():
|
|
return direction
|
|
row_delta = (receiver[moved] // spec.n_cols - index[moved] // spec.n_cols).astype(np.float64)
|
|
col_delta = (receiver[moved] % spec.n_cols - index[moved] % spec.n_cols).astype(np.float64)
|
|
angle = np.arctan2(row_delta, col_delta)
|
|
direction[moved] = (
|
|
np.rint(angle / (2.0 * np.pi / AZIMUTH_STEPS)).astype(np.int16) % AZIMUTH_STEPS
|
|
)
|
|
return direction
|
|
|
|
|
|
def _walk_from(
|
|
starts: np.ndarray,
|
|
receiver: np.ndarray,
|
|
valid: np.ndarray,
|
|
absorbing: np.ndarray,
|
|
state: np.ndarray,
|
|
reaches: np.ndarray,
|
|
) -> int:
|
|
"""출발 셀에서 물길 사슬을 따라가며 판정하고, 사슬 전체에 같은 색을 적는다.
|
|
|
|
`absorbing`은 도로 셀과 확정된 세류망 셀 — 여기 닿으면 사슬 전체가 도로 도달이다.
|
|
새로 판정한 셀 수를 돌려준다.
|
|
"""
|
|
resolved = 0
|
|
path: list[int] = []
|
|
for start in starts.tolist():
|
|
if state[start] == 2:
|
|
continue
|
|
path.clear()
|
|
node = start
|
|
while True:
|
|
if state[node] == 2:
|
|
verdict = bool(reaches[node]) # 이미 색이 정해진 셀 — 그 색을 물려받는다
|
|
break
|
|
if state[node] == 1: # 방어: 채움·평탄해소 후에는 순환이 없어야 한다
|
|
verdict = False
|
|
break
|
|
state[node] = 1
|
|
path.append(node)
|
|
if absorbing[node]:
|
|
verdict = True # 도로 또는 세류망에 합류 — 여기서 하류로 빠진다
|
|
break
|
|
following = int(receiver[node])
|
|
if following == node:
|
|
verdict = False # 싱크에 갇힘
|
|
break
|
|
if not valid[following] and state[following] != 2:
|
|
# 해석 영역 밖으로 빠짐. 단, 이미 색이 정해진 셀(세류 선확정 등)이면 따라간다.
|
|
verdict = False
|
|
break
|
|
node = following
|
|
for visited in path:
|
|
reaches[visited] = verdict
|
|
state[visited] = 2
|
|
resolved += len(path)
|
|
return resolved
|
|
|
|
|
|
# ── 격자 확장 (별도 단계) ───────────────────────────────────────────────────
|
|
|
|
|
|
@dataclass
|
|
class ExpansionResult:
|
|
"""확장 루프 결과. 확장을 쓰지 않는 경로에서는 이 모듈을 부르지 않는다."""
|
|
|
|
spec: GridSpec
|
|
terrain: TerrainGrid
|
|
road: RoadRaster
|
|
flow: FlowResult
|
|
rounds: int # 실제로 넓힌 횟수 (0 = 처음부터 닫혀 있었음)
|
|
closed: bool # 경계 링이 전부 비활성이 되어 스스로 멈췄는가
|
|
|
|
|
|
def expand_until_closed(
|
|
spec: GridSpec,
|
|
cloud: ContourCloud,
|
|
route_line: LineString,
|
|
region_area: Any = None,
|
|
road_width_m: float = DRAINAGE_ROAD_WIDTH_M,
|
|
max_rounds: int = DRAINAGE_MAX_EXPAND_ROUNDS,
|
|
step_m: float = DRAINAGE_EXPAND_STEP_M,
|
|
) -> ExpansionResult:
|
|
"""활성 셀이 격자 최외곽에 닿은 방향으로만 넓히며 유역이 닫힐 때까지 반복한다.
|
|
|
|
**본 계산 경로에서만 쓰는 별도 단계다**(2026-07-31 사용자 지시로 분리). 단계 검증
|
|
미리보기는 확장 없이 현재 격자만 본다 — 확장 로직 자체가 아직 검증 대상이기 때문이다.
|
|
|
|
종료 조건은 반경 상한이 아니라 **경계 링 전체가 비활성**이 되는 것이다. 비활성 셀이
|
|
최외곽에 띠로 완성되면 그 바깥은 볼 필요가 없다. `max_rounds`는 무한 반복 방지용이다.
|
|
|
|
`region_area`를 주면 확장된 격자에서도 그 영역에 걸치는 셀만 해석 대상으로 삼는다.
|
|
"""
|
|
terrain = road = flow = None
|
|
rounds = 0
|
|
closed = False
|
|
for attempt in range(max_rounds + 1):
|
|
domain = build_cell_mask(spec, region_area) if region_area is not None else None
|
|
terrain = build_terrain_grid(spec, cloud, domain)
|
|
road = rasterize_road(spec, route_line, road_width_m)
|
|
flow = trace_flow(terrain, road)
|
|
contact = border_contact(flow.active)
|
|
if not any(contact.values()):
|
|
closed = True
|
|
break
|
|
if attempt == max_rounds:
|
|
logger.warning(
|
|
"배수유역: 확장 상한(%d회) 도달 — 경계 %s가 아직 활성입니다.",
|
|
max_rounds,
|
|
[side for side, touched in contact.items() if touched],
|
|
)
|
|
break
|
|
widened = expand_grid_spec(spec, contact, step_m)
|
|
if widened == spec:
|
|
break
|
|
logger.info(
|
|
"배수유역: 경계 %s 활성 — %.0fm 확장 (%d회차)",
|
|
[side for side, touched in contact.items() if touched],
|
|
step_m,
|
|
attempt + 1,
|
|
)
|
|
spec = widened
|
|
rounds += 1
|
|
|
|
assert terrain is not None and road is not None and flow is not None
|
|
return ExpansionResult(
|
|
spec=spec, terrain=terrain, road=road, flow=flow, rounds=rounds, closed=closed
|
|
)
|
|
|
|
|
|
def border_contact(active: np.ndarray) -> dict[str, bool]:
|
|
"""활성 셀이 격자 최외곽에 닿은 방향. 전부 False면 유역이 능선 안에서 닫힌 것이다."""
|
|
return {
|
|
"north": bool(active[0, :].any()),
|
|
"south": bool(active[-1, :].any()),
|
|
"west": bool(active[:, 0].any()),
|
|
"east": bool(active[:, -1].any()),
|
|
}
|
|
|
|
|
|
# ── 폴리곤화 ────────────────────────────────────────────────────────────────
|
|
|
|
|
|
def polygonize_labels(
|
|
spec: GridSpec,
|
|
labels: np.ndarray,
|
|
min_area_m2: float = DRAINAGE_MIN_BASIN_AREA_M2,
|
|
simplify_m: float = DRAINAGE_POLYGON_SIMPLIFY_M,
|
|
) -> dict[int, Polygon | MultiPolygon]:
|
|
"""라벨 격자를 라벨별 폴리곤으로 바꾼다. 음수 라벨은 배경으로 무시한다.
|
|
|
|
`simplify_m`은 격자 계단 경계를 줄여 응답 크기를 낮추는 값이다. 셀 경계를 눈으로 대조해야
|
|
하는 용도(유입 셀 확인 등)에서는 0을 줘서 원래 계단 그대로 받는다.
|
|
"""
|
|
label_grid = np.ascontiguousarray(labels.reshape(spec.n_rows, spec.n_cols), dtype=np.int32)
|
|
valid_mask = label_grid >= 0
|
|
if not valid_mask.any():
|
|
return {}
|
|
collected: dict[int, list[Polygon]] = {}
|
|
for geometry, value in shapes(
|
|
label_grid, mask=valid_mask, transform=grid_transform(spec), connectivity=4
|
|
):
|
|
polygon = shape(geometry)
|
|
if polygon.is_empty or polygon.area < min_area_m2:
|
|
continue
|
|
collected.setdefault(int(value), []).append(polygon)
|
|
|
|
merged: dict[int, Polygon | MultiPolygon] = {}
|
|
for label, parts in collected.items():
|
|
union = unary_union(parts)
|
|
if union.is_empty:
|
|
continue
|
|
if simplify_m <= 0:
|
|
merged[label] = union
|
|
continue
|
|
simplified = union.simplify(simplify_m, preserve_topology=True)
|
|
merged[label] = simplified if not simplified.is_empty else union
|
|
return merged
|
|
|
|
|
|
def largest_ring(geometry: Polygon | MultiPolygon) -> list[tuple[float, float]]:
|
|
"""폴리곤(또는 멀티폴리곤)에서 가장 큰 조각의 외곽 링 좌표를 뽑는다."""
|
|
if geometry.is_empty:
|
|
return []
|
|
if geometry.geom_type == "MultiPolygon":
|
|
geometry = max(geometry.geoms, key=lambda part: part.area)
|
|
return [(float(x), float(y)) for x, y in geometry.exterior.coords]
|
|
|
|
|
|
def outer_boundary(
|
|
spec: GridSpec, active: np.ndarray, min_area_m2: float = DRAINAGE_MIN_BASIN_AREA_M2
|
|
) -> Polygon | MultiPolygon | None:
|
|
"""활성 셀 전체의 외곽 = 2차 전체 배수유역 경계.
|
|
|
|
비활성 셀이 격자 최외곽에 띠로 완성되면 활성 영역이 그 안에 닫힌다. 그 닫힌 영역의
|
|
바깥선이 곧 분수령이므로 능선을 따로 그릴 필요가 없다.
|
|
"""
|
|
labels = np.where(active.reshape(-1), 0, -1).astype(np.int32)
|
|
polygons = polygonize_labels(spec, labels, min_area_m2)
|
|
return polygons.get(0)
|