diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py index bdcfe002..ebf41f7b 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py @@ -1,13 +1,15 @@ -"""배수유역 능선(분수령) 기반 산정 엔진. +"""배수유역 도로(노선) 기준 산정 엔진 — 측구 흐름 모델. -목적: 도로(관 매설 지점)로 모이는 물의 양을 알기 위한 유역 산정. 유역 경계는 반드시 -분수령(능선)을 따라야 하므로, 도엽 등고선·표고점을 격자 DEM으로 보간한 뒤 D8 흐름 -방향으로 "각 셀의 물이 어느 측점으로 흘러가는가"를 직접 추적한다. +목적: 도로로 오는 물의 양 산정. 참고 도면(2026-07-28 사용자 제공)과 같이 **도로 산측 +사면 전체를 관(구조물 측점) 개수만큼 빈틈없이 분할**한다: -- 3D 라이다는 산 전체를 계측하지 않으므로 쓰지 않는다. 도엽 데이터만 사용(사용자 확정). -- 함몰 보정은 Whitebox `fill_depressions`를 쓰고(기존 Engine_Skeleton과 동일 패턴), - 실패하면 원본 DEM으로 진행한다(결과 저하 가능하나 계산은 지속). -- 유역 폴리곤 외곽선이 곧 분수령(능선)이며 프론트가 능선 스타일로 표시한다. +- 도엽 등고선·표고점을 격자 DEM으로 보간 (3D 라이다는 산 전체를 계측하지 않아 미사용). +- 노선을 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로 도로 셀마다 + 담당 관을 배정 — "도로에 닿은 물은 측구를 타고 내리막의 첫 관으로 들어간다". +- 사면 각 셀은 D8 흐름으로 내려가 처음 닿는 도로 셀의 관을 물려받는다. 도로를 만나지 + 못하는 셀(도로 하측 사면, 능선 너머)은 자동 제외. +- 함몰 보정은 Whitebox `fill_depressions`(Engine_Skeleton과 동일 패턴), 실패 시 원본 진행. +- 유역 경계의 산측이 분수령(능선)·지능선이고 하측이 도로선이다. 프론트가 능선 파선 표시. """ from __future__ import annotations @@ -33,10 +35,8 @@ logger = logging.getLogger(__name__) # 격자 해상도(m)와 최대 격자 크기. 도엽 9매 범위라도 이 상한 안에서 해상도를 낮춰 계산한다. GRID_RES_M = 10.0 MAX_GRID_CELLS = 1_400_000 -# 유역 계산 범위: 측점 bbox + 여유폭(m). 주변 8도엽까지 확보되어 있어 넉넉히 잡는다. +# 유역 계산 범위: 노선 bbox + 여유폭(m). 주변 8도엽까지 확보되어 있어 넉넉히 잡는다. BBOX_MARGIN_M = 1500.0 -# pour point 스냅 반경(m): 측점을 주변 흐름 누적 최대 셀로 옮겨 세류 격자 정합 오차를 흡수. -SNAP_RADIUS_M = 50.0 # DEM 보간 표본 상한(속도 확보용 간축). 초과 시 균등 간격으로 추린다. MAX_SAMPLE_POINTS = 250_000 @@ -118,15 +118,20 @@ def _collect_samples( def _build_dem( samples: np.ndarray, - candidates: list[StructureCandidate], + anchors_x: list[float], + anchors_y: list[float], ) -> tuple[np.ndarray, np.ndarray, np.ndarray, float] | None: - """표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도).""" - if len(samples) < 10 or not candidates: + """표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도). + + 범위는 앵커(노선 전체 정점) bbox + 여유폭 — 유역이 노선 전 연장을 덮어야 하므로 + 측점 bbox가 아니라 노선 bbox를 쓴다. + """ + if len(samples) < 10 or not anchors_x: return None - min_x = min(candidate.x for candidate in candidates) - BBOX_MARGIN_M - max_x = max(candidate.x for candidate in candidates) + BBOX_MARGIN_M - min_y = min(candidate.y for candidate in candidates) - BBOX_MARGIN_M - max_y = max(candidate.y for candidate in candidates) + BBOX_MARGIN_M + min_x = min(anchors_x) - BBOX_MARGIN_M + max_x = max(anchors_x) + BBOX_MARGIN_M + min_y = min(anchors_y) - BBOX_MARGIN_M + max_y = max(anchors_y) + BBOX_MARGIN_M # 표본 범위 밖으로는 나가지 않는다(외삽 방지). min_x = max(min_x, float(samples[:, 0].min())) max_x = min(max_x, float(samples[:, 0].max())) @@ -210,49 +215,75 @@ def _d8_pointer(dem: np.ndarray) -> np.ndarray: return pointer -def _flow_accumulation(pointer: np.ndarray) -> np.ndarray: - """D8 포인터 기반 흐름 누적(자기 자신 포함 셀 수). 위상 순서로 한 번에 계산.""" - rows, cols = pointer.shape - accumulation = np.ones((rows, cols), dtype=np.float64) - indegree = np.zeros((rows, cols), dtype=np.int32) - for direction, (dr, dc, _) in enumerate(_D8): - sources = np.argwhere(pointer == direction) - for r, c in sources: - nr, nc = r + dr, c + dc - if 0 <= nr < rows and 0 <= nc < cols: - indegree[nr, nc] += 1 - stack = [tuple(cell) for cell in np.argwhere(indegree == 0)] - while stack: - r, c = stack.pop() - direction = pointer[r, c] - if direction < 0: - continue - dr, dc, _ = _D8[direction] - nr, nc = r + dr, c + dc - if not (0 <= nr < rows and 0 <= nc < cols): - continue - accumulation[nr, nc] += accumulation[r, c] - indegree[nr, nc] -= 1 - if indegree[nr, nc] == 0: - stack.append((nr, nc)) - return accumulation +def _ditch_intervals( + vertices: list[Any], + ordered: list[StructureCandidate], +) -> list[tuple[float, float, int]]: + """측구(도로변 배수로) 흐름 모델의 담당 구간을 만든다. + + 도로에 닿은 물은 종단 내리막 방향으로 흘러 첫 번째 관으로 들어간다. 따라서 인접한 + 두 관 사이 **종단 계획선의 최고점(물갈림 고개)**이 유역 분할점이 된다. + 반환: (구간 시작 chainage, 구간 끝 chainage, 관 라벨) 목록. + """ + total_length = vertices[-1].chainage_m if vertices else 0.0 + divides: list[float] = [] + for left, right in zip(ordered, ordered[1:]): + window = [ + vertex for vertex in vertices if left.chainage_m < vertex.chainage_m < right.chainage_m + ] + if window: + divides.append(max(window, key=lambda vertex: vertex.z).chainage_m) + else: + divides.append((left.chainage_m + right.chainage_m) / 2.0) + intervals: list[tuple[float, float, int]] = [] + start = 0.0 + for label, divide in enumerate(divides, start=1): + intervals.append((start, divide, label)) + start = divide + intervals.append((start, total_length + 1.0, len(ordered))) + return intervals -def _snap_outlet( - accumulation: np.ndarray, - row: int, - col: int, - radius_cells: int, -) -> tuple[int, int]: - """측점 주변 반경 안에서 흐름 누적이 가장 큰 셀로 옮긴다(물길 위로 스냅).""" - rows, cols = accumulation.shape - r0 = max(0, row - radius_cells) - r1 = min(rows, row + radius_cells + 1) - c0 = max(0, col - radius_cells) - c1 = min(cols, col + radius_cells + 1) - window = accumulation[r0:r1, c0:c1] - local = np.unravel_index(int(np.argmax(window)), window.shape) - return r0 + int(local[0]), c0 + int(local[1]) +def _rasterize_road( + vertices: list[Any], + intervals: list[tuple[float, float, int]], + x_coords: np.ndarray, + y_coords: np.ndarray, + resolution: float, + shape_rc: tuple[int, int], +) -> dict[tuple[int, int], int]: + """노선 폴리라인을 격자에 새겨 도로 셀마다 담당 관 라벨을 붙인다. + + 대각 누수(도로가 셀 사이 대각으로 지나가 물이 새는 것)를 막으려고 도로 셀 주변 + 3×3을 함께 같은 라벨로 칠한다. + """ + rows, cols = shape_rc + + def _label_of(chainage: float) -> int: + for start, end, label in intervals: + if start <= chainage < end: + return label + return intervals[-1][2] if intervals else 0 + + road: dict[tuple[int, int], int] = {} + step = resolution / 2.0 + for previous, current in zip(vertices, vertices[1:]): + span = math.dist((previous.x, previous.y), (current.x, current.y)) + count = max(1, int(span / step)) + for i in range(count + 1): + ratio = i / count + x = previous.x + (current.x - previous.x) * ratio + y = previous.y + (current.y - previous.y) * ratio + chainage = previous.chainage_m + (current.chainage_m - previous.chainage_m) * ratio + col = int(round((x - float(x_coords[0])) / resolution)) + row = int(round((y - float(y_coords[0])) / resolution)) + label = _label_of(chainage) + for dr in (-1, 0, 1): + for dc in (-1, 0, 1): + r, c = row + dr, col + dc + if 0 <= r < rows and 0 <= c < cols and (r, c) not in road: + road[(r, c)] = label + return road def _label_basins( @@ -380,40 +411,41 @@ def _vectorize_basin( def build_watershed_basins( + vertices: list[Any], candidates: list[StructureCandidate], contour_features: list[dict[str, Any]], spot_features: list[dict[str, Any]], elevation_keys: tuple[str, ...], ) -> list[WatershedBasin]: - """능선(분수령) 기반 배수유역을 산정한다. + """도로(노선) 기준 배수유역을 산정한다 — 측구 흐름 모델. - 반환된 boundary_xy 외곽선이 곧 분수령(능선)이다. 번호는 노선 시점에 가까운 순. + "도로에 닿은 물은 측구를 타고 종단 내리막 방향으로 흘러 첫 번째 관으로 들어간다"를 + 전제로, 도로 산측 사면 전체를 관 개수만큼 빈틈없이 분할한다(참고 도면과 동일 개념): + ① 노선 전체를 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로 + 각 도로 셀에 담당 관 라벨을 붙인다. + ② 사면 각 셀은 D8 흐름을 따라 내려가 처음 닿는 도로 셀의 관 라벨을 물려받는다. + 도로를 만나지 못하고 격자 밖으로 빠지는 셀(도로 하측 성토부 등)은 제외된다. + 반환된 boundary_xy 외곽선의 산측이 곧 분수령(능선)이며, 하측은 도로선을 따른다. + 번호는 노선 시점에 가까운 순. """ - if not candidates: + if not candidates or len(vertices) < 2: return [] samples = _collect_samples(contour_features, spot_features, elevation_keys) - built = _build_dem(samples, candidates) + built = _build_dem( + samples, + [vertex.x for vertex in vertices], + [vertex.y for vertex in vertices], + ) if built is None: logger.warning("DEM 보간 실패 — 표본 %d개", len(samples)) return [] dem, x_coords, y_coords, resolution = built dem = _fill_depressions(dem, resolution) pointer = _d8_pointer(dem) - accumulation = _flow_accumulation(pointer) ordered = sorted(candidates, key=lambda item: item.chainage_m) - radius_cells = max(1, int(SNAP_RADIUS_M / resolution)) - outlets: dict[tuple[int, int], int] = {} - outlet_cells: dict[int, tuple[int, int]] = {} - for label, candidate in enumerate(ordered, start=1): - col = int(round((candidate.x - float(x_coords[0])) / resolution)) - row = int(round((candidate.y - float(y_coords[0])) / resolution)) - if not (0 <= row < dem.shape[0] and 0 <= col < dem.shape[1]): - continue - snapped = _snap_outlet(accumulation, row, col, radius_cells) - outlets[snapped] = label - outlet_cells[label] = snapped - + intervals = _ditch_intervals(vertices, ordered) + outlets = _rasterize_road(vertices, intervals, x_coords, y_coords, resolution, dem.shape) if not outlets: return [] labels = _label_basins(pointer, outlets) @@ -421,9 +453,6 @@ def build_watershed_basins( basins: list[WatershedBasin] = [] for label, candidate in enumerate(ordered, start=1): - cell = outlet_cells.get(label) - if cell is None: - continue mask = labels == label cell_count = int(mask.sum()) if cell_count < 4: @@ -431,7 +460,11 @@ def build_watershed_basins( boundary = _vectorize_basin(labels, label, x_coords, y_coords, resolution) if len(boundary) < 4: continue - outlet_z = float(dem[cell[0], cell[1]]) + # 측점(관) 위치의 표고를 기준으로 낙차를 계산한다. + col = int(round((candidate.x - float(x_coords[0])) / resolution)) + row = int(round((candidate.y - float(y_coords[0])) / resolution)) + in_grid = 0 <= row < dem.shape[0] and 0 <= col < dem.shape[1] + outlet_z = float(dem[row, col]) if in_grid else float(dem[mask].min()) basin = WatershedBasin( index=label, chainage_m=candidate.chainage_m, diff --git a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py index 3d137c3f..ecb3a54d 100644 --- a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py @@ -194,7 +194,7 @@ async def post_drainage_basins( candidates = propose_structure_stations(vertices, prepared["streams"]) basins = build_watershed_basins( - candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS + vertices, candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS ) to_lonlat = prepared["to_lonlat"] return {