"""배수유역 흐름 해석 — 도로 굽기 · 상류 추적(포인터 더블링) · 유역 폴리곤화. 핵심은 하나다: **물길을 따라가 도로에 닿는 셀만 유역이다.** 셀마다 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)