From 8b853e8deca3eb30268993bac7216969ba9eecd4 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Wed, 29 Jul 2026 18:07:39 +0900 Subject: [PATCH] auto: 2026-07-29 18:07 (EOMSANGDON-HOME) --- .../B05_wf2_Route_Engine_Drainage.py | 43 +- ...B05_wf2_Route_Engine_Drainage_Watershed.py | 653 +++++++----------- .../B05_wf2_Route_Engine_Watershed_Trace.py | 271 ++++++++ .../B05_wf2_Route_Router_Drainage.py | 7 +- 4 files changed, 571 insertions(+), 403 deletions(-) create mode 100644 B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py index 99be8d85..87de9cce 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py @@ -228,12 +228,18 @@ def _fill_spacing( start_m: float, end_m: float, ) -> list[StructureCandidate]: - """[start, end] 구간이 300m를 넘으면 절토부 지점에 보충 측점을 만든다.""" + """[start, end] 구간이 300m를 넘으면 보충 측점을 만든다. + + 종단도상 상대적으로 물이 모일 것으로 예상되는 지점(절토부 내 종단 저점)을 + 우선 배치한다(2026-07-29 사용자 지시). 저점이 없으면 목표 인근 절토부로 대체한다. + """ added: list[StructureCandidate] = [] cursor = start_m while end_m - cursor > MAX_STRUCTURE_SPACING_M: target = cursor + MAX_STRUCTURE_SPACING_M - placed = _nearest_uphill(vertices, target, end_m) + placed = _gather_low_point(vertices, cursor, target, end_m) + if placed is None: + placed = _nearest_uphill(vertices, target, end_m) if placed is None: break x, y, _ = _interpolate_vertex(vertices, placed) @@ -242,6 +248,39 @@ def _fill_spacing( return added +def _gather_low_point( + vertices: list[RouteVertex], + cursor_m: float, + target_m: float, + limit_m: float, + step_m: float = 10.0, +) -> float | None: + """탐색창 [cursor+150, target] 안 절토부의 종단 국소 저점(사그) 중 가장 낮은 지점. + + 창 하한을 간격의 절반으로 두어 보충 측점이 과밀하게 몰리지 않게 하고, + 국소 저점만 인정해 일정 오르막에서는 None(300m 규칙 폴백)을 돌려준다. + """ + window_start = cursor_m + MAX_STRUCTURE_SPACING_M / 2.0 + probes: list[float] = [] + probe = window_start - step_m + while probe <= target_m + step_m: + probes.append(probe) + probe += step_m + heights = [_interpolate_vertex(vertices, position)[2] for position in probes] + best: tuple[float, float] | None = None # (계획고 z, 누가거리) + for i in range(1, len(probes) - 1): + position = probes[i] + if position >= limit_m or position > target_m or position < window_start: + continue + # 국소 저점(양쪽이 같거나 높음) = 물이 모여 더 못 흐르는 지점. 앞쪽이 오르막인 + # 조건을 내포하므로 별도의 절토부(is_uphill_at) 판정은 두지 않는다. + if heights[i] > heights[i - 1] or heights[i] > heights[i + 1]: + continue + if best is None or heights[i] < best[0]: + best = (heights[i], position) + return best[1] if best else None + + def _nearest_uphill( vertices: list[RouteVertex], target_m: float, 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 ebf41f7b..5b676eaf 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py @@ -1,56 +1,50 @@ -"""배수유역 도로(노선) 기준 산정 엔진 — 측구 흐름 모델. +"""배수유역 산정 엔진 — 세류 기반 등고선 기하 직접 분석 (2026-07-29 합의). -목적: 도로로 오는 물의 양 산정. 참고 도면(2026-07-28 사용자 제공)과 같이 **도로 산측 -사면 전체를 관(구조물 측점) 개수만큼 빈틈없이 분할**한다: +DEM 보간·D8 전역 흐름분석을 쓰지 않는다. 계산 순서(사용자 정의 7단계): +① 도로(노선)가 유역의 하측 경계 ② 도로를 가로지르는 세류 교차점에서 출발 +③ 세류 상류망을 추적하고 연관 등고선만 분석해 메인 유역 선정(하류 무의미) +④ 세류 교차점 = 관매설 지점 ⑤ 300m 초과 구간은 종단 저점에 보충(제안 엔진 담당) +⑥ 관 사이 물갈림 고개에서 오르는 분할선(능선 근사)으로 유역을 세분화하고 + 번호·면적·표고차·유하장을 산출 ⑦ 관 추가·경로 변경 시 재호출로 재분석. -- 도엽 등고선·표고점을 격자 DEM으로 보간 (3D 라이다는 산 전체를 계측하지 않아 미사용). -- 노선을 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로 도로 셀마다 - 담당 관을 배정 — "도로에 닿은 물은 측구를 타고 내리막의 첫 관으로 들어간다". -- 사면 각 셀은 D8 흐름으로 내려가 처음 닿는 도로 셀의 관을 물려받는다. 도로를 만나지 - 못하는 셀(도로 하측 사면, 능선 너머)은 자동 제외. -- 함몰 보정은 Whitebox `fill_depressions`(Engine_Skeleton과 동일 패턴), 실패 시 원본 진행. -- 유역 경계의 산측이 분수령(능선)·지능선이고 하측이 도로선이다. 프론트가 능선 파선 표시. +유역 폴리곤 = 도로 구간(하측) + 좌우 분할선 + 최상위 공통 등고선 아크(상측)로 폐합. +등고선은 STRtree에서 필요한 것만 꺼내므로 분석량이 유역 크기에 비례한다(도엽 매수 무관). """ from __future__ import annotations import logging import math -import tempfile from dataclasses import dataclass, field -from pathlib import Path from typing import Any -import numpy as np -from scipy.interpolate import griddata -from shapely.geometry import shape +from shapely.geometry import LineString, Point, Polygon +from shapely.ops import substring, unary_union from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( StructureCandidate, + _interpolate_vertex, estimate_pipe_diameter_mm, ) +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Trace import ( + LOCAL_MAX_STEPS, + LOCAL_SEARCH_RADIUS_M, + STREAM_JOIN_TOL_M, + ContourIndex, + DividerStep, + _explode_lines, + trace_divider, + trace_upstream_network, +) logger = logging.getLogger(__name__) -# 격자 해상도(m)와 최대 격자 크기. 도엽 9매 범위라도 이 상한 안에서 해상도를 낮춰 계산한다. -GRID_RES_M = 10.0 -MAX_GRID_CELLS = 1_400_000 -# 유역 계산 범위: 노선 bbox + 여유폭(m). 주변 8도엽까지 확보되어 있어 넉넉히 잡는다. -BBOX_MARGIN_M = 1500.0 -# DEM 보간 표본 상한(속도 확보용 간축). 초과 시 균등 간격으로 추린다. -MAX_SAMPLE_POINTS = 250_000 - -# D8 이웃: (행 오프셋, 열 오프셋, 거리 계수) -_D8 = ( - (-1, -1, math.sqrt(2.0)), - (-1, 0, 1.0), - (-1, 1, math.sqrt(2.0)), - (0, -1, 1.0), - (0, 1, 1.0), - (1, -1, math.sqrt(2.0)), - (1, 0, 1.0), - (1, 1, math.sqrt(2.0)), -) +# 유효 유역 최소 면적(m²)과 상류망 커버 보정 버퍼(m). +MIN_BASIN_AREA_M2 = 100.0 +STREAM_COVER_BUFFER_M = 20.0 +# 산측 판정: 도로 접선의 좌우 수직 오프셋(m)과 표고 조회 반경(m). +UPHILL_PROBE_OFFSET_M = 40.0 +UPHILL_PROBE_RADIUS_M = 60.0 @dataclass @@ -61,7 +55,7 @@ class WatershedBasin: chainage_m: float outlet_x: float outlet_y: float - # 유역 경계(사업지 좌표계 m). 외곽선이 곧 분수령(능선). + # 유역 경계(사업지 좌표계 m). 외곽선의 산측이 곧 분수령(능선), 하측이 도로선. boundary_xy: list[list[float]] = field(default_factory=list) area_m2: float = 0.0 relief_m: float = 0.0 @@ -69,164 +63,13 @@ class WatershedBasin: pipe_diameter_mm: float | None = None -def _collect_samples( - contour_features: list[dict[str, Any]], - spot_features: list[dict[str, Any]], - elevation_keys: tuple[str, ...], -) -> np.ndarray: - """등고선 정점·표고점을 (x, y, z) 표본 배열로 모은다.""" - xs: list[float] = [] - ys: list[float] = [] - zs: list[float] = [] +def _divide_chainages(vertices: list[Any], ordered: list[StructureCandidate]) -> list[float]: + """유역 분할 누가거리 목록(양끝 포함, 관 개수+1개). - def _walk(coordinates: Any, elevation: float) -> None: - if not isinstance(coordinates, list) or not coordinates: - return - if isinstance(coordinates[0], (int, float)): - xs.append(float(coordinates[0])) - ys.append(float(coordinates[1])) - zs.append(elevation) - return - for item in coordinates: - _walk(item, elevation) - - for feature in [*contour_features, *spot_features]: - properties = feature.get("properties") or {} - elevation: float | None = None - for key in elevation_keys: - value = properties.get(key) - if value is None: - continue - try: - elevation = float(value) - break - except (TypeError, ValueError): - continue - if elevation is None: - continue - geometry = feature.get("geometry") or {} - _walk(geometry.get("coordinates"), elevation) - - if not xs: - return np.empty((0, 3)) - samples = np.column_stack([xs, ys, zs]) - if len(samples) > MAX_SAMPLE_POINTS: - step = len(samples) // MAX_SAMPLE_POINTS + 1 - samples = samples[::step] - return samples - - -def _build_dem( - samples: np.ndarray, - anchors_x: list[float], - anchors_y: list[float], -) -> tuple[np.ndarray, np.ndarray, np.ndarray, float] | None: - """표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도). - - 범위는 앵커(노선 전체 정점) bbox + 여유폭 — 유역이 노선 전 연장을 덮어야 하므로 - 측점 bbox가 아니라 노선 bbox를 쓴다. + 인접한 두 관 사이 **종단 계획선의 최고점(물갈림 고개)**이 분할점이다 — + "도로에 닿은 물은 측구를 타고 내리막의 첫 관으로 들어간다". """ - if len(samples) < 10 or not anchors_x: - return None - 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())) - min_y = max(min_y, float(samples[:, 1].min())) - max_y = min(max_y, float(samples[:, 1].max())) - if max_x - min_x < GRID_RES_M * 4 or max_y - min_y < GRID_RES_M * 4: - return None - - resolution = GRID_RES_M - while ((max_x - min_x) / resolution) * ((max_y - min_y) / resolution) > MAX_GRID_CELLS: - resolution *= 1.5 - x_coords = np.arange(min_x, max_x + resolution, resolution) - y_coords = np.arange(min_y, max_y + resolution, resolution) - grid_x, grid_y = np.meshgrid(x_coords, y_coords) - - points = samples[:, :2] - values = samples[:, 2] - dem = griddata(points, values, (grid_x, grid_y), method="linear") - # linear 보간 밖(볼록 껍질 바깥)은 nearest로 메워 유역 추적이 끊기지 않게 한다. - holes = ~np.isfinite(dem) - if holes.any(): - dem[holes] = griddata(points, values, (grid_x[holes], grid_y[holes]), method="nearest") - return dem.astype(np.float64), x_coords, y_coords, resolution - - -def _fill_depressions(dem: np.ndarray, resolution: float) -> np.ndarray: - """Whitebox로 함몰을 메운다. 실패하면 원본 그대로 진행한다.""" - try: - import rasterio - from rasterio.transform import from_origin - from whitebox import WhiteboxTools - except Exception: # noqa: BLE001 - return dem - rows, cols = dem.shape - try: - with tempfile.TemporaryDirectory(prefix="wbt_drain_") as tmp: - tmp_path = Path(tmp) - transform = from_origin(0.0, rows * resolution, resolution, resolution) - with rasterio.open( - tmp_path / "dem.tif", - "w", - driver="GTiff", - height=rows, - width=cols, - count=1, - dtype="float32", - nodata=-9999.0, - crs="EPSG:3857", - transform=transform, - ) as dst: - dst.write(dem.astype(np.float32)[::-1, :], 1) - wbt = WhiteboxTools() - wbt.set_verbose_mode(False) - wbt.set_working_dir(str(tmp_path)) - if wbt.fill_depressions("dem.tif", "filled.tif") != 0: - raise RuntimeError("fill_depressions 실패") - with rasterio.open(tmp_path / "filled.tif") as src: - filled = src.read(1).astype(np.float64)[::-1, :] - return np.where(np.isfinite(filled), filled, dem) - except Exception: # noqa: BLE001 - logger.warning("Whitebox 함몰 보정 실패 — 원본 DEM으로 진행") - return dem - - -def _d8_pointer(dem: np.ndarray) -> np.ndarray: - """각 셀의 최급강하 이웃 인덱스(0~7, 배수구 없으면 -1).""" - rows, cols = dem.shape - pointer = np.full((rows, cols), -1, dtype=np.int8) - best_drop = np.zeros((rows, cols), dtype=np.float64) - for direction, (dr, dc, distance) in enumerate(_D8): - shifted = np.full_like(dem, np.inf) - r_src = slice(max(0, -dr), rows - max(0, dr)) - c_src = slice(max(0, -dc), cols - max(0, dc)) - r_dst = slice(max(0, dr), rows - max(0, -dr)) - c_dst = slice(max(0, dc), cols - max(0, -dc)) - shifted[r_src, c_src] = dem[r_dst, c_dst] - drop = (dem - shifted) / distance - better = drop > best_drop - pointer[better] = direction - best_drop[better] = drop[better] - return pointer - - -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] = [] + divides = [vertices[0].chainage_m] for left, right in zip(ordered, ordered[1:]): window = [ vertex for vertex in vertices if left.chainage_m < vertex.chainage_m < right.chainage_m @@ -235,179 +78,137 @@ def _ditch_intervals( 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 + divides.append(vertices[-1].chainage_m) + return divides -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( - pointer: np.ndarray, - outlets: dict[tuple[int, int], int], -) -> np.ndarray: - """각 셀이 흐름을 따라 처음 만나는 pour point의 라벨을 붙인다(경로 메모이제이션).""" - rows, cols = pointer.shape - labels = np.zeros((rows, cols), dtype=np.int32) # 0 = 미소속 - for (r, c), label in outlets.items(): - labels[r, c] = label - flat_pointer = pointer.ravel() - flat_labels = labels.ravel() - for start in range(flat_labels.size): - if flat_labels[start] != 0: - continue - path: list[int] = [] - current = start - label = 0 - while True: - if flat_labels[current] != 0: - label = flat_labels[current] - break - direction = flat_pointer[current] - if direction < 0: - label = -1 # 배수구 없음(격자 밖 유출) — 어떤 유역에도 속하지 않음 - break - path.append(current) - dr, dc, _ = _D8[direction] - r, c = divmod(current, cols) - nr, nc = r + dr, c + dc - if not (0 <= nr < rows and 0 <= nc < cols): - label = -1 - break - current = nr * cols + nc - for cell in path: - flat_labels[cell] = label - return labels - - -def _flow_lengths(pointer: np.ndarray, labels: np.ndarray, resolution: float) -> dict[int, float]: - """라벨별 최장 흐름 경로(셀→해당 pour point) 길이.""" - rows, cols = pointer.shape - distance = np.full((rows, cols), -1.0, dtype=np.float64) - # pour point 셀은 자기 라벨의 시작점이므로 거리 0. - longest: dict[int, float] = {} - flat_pointer = pointer.ravel() - flat_labels = labels.ravel() - flat_distance = distance.ravel() - - def _resolve(start: int) -> float: - chain: list[int] = [] - current = start - total = 0.0 - while True: - if flat_distance[current] >= 0: - total = flat_distance[current] - break - direction = flat_pointer[current] - if direction < 0: - total = 0.0 - break - r, c = divmod(current, cols) - dr, dc, factor = _D8[direction] - nr, nc = r + dr, c + dc - if not (0 <= nr < rows and 0 <= nc < cols): - total = 0.0 - break - next_cell = nr * cols + nc - # 다음 셀이 다른 라벨이면(=pour point 통과) 여기서 경로가 끝난 것으로 본다. - chain.append(current) - if flat_labels[next_cell] != flat_labels[current]: - total = 0.0 - break - current = next_cell - # 뒤에서부터 거리를 되채운다. - for cell in reversed(chain): - direction = flat_pointer[cell] - factor = _D8[direction][2] if direction >= 0 else 0.0 - total += factor * resolution - flat_distance[cell] = total - return total - - for start in range(flat_labels.size): - label = int(flat_labels[start]) - if label <= 0: - continue - length = _resolve(start) - if length > longest.get(label, 0.0): - longest[label] = length - return longest - - -def _vectorize_basin( - labels: np.ndarray, - label: int, - x_coords: np.ndarray, - y_coords: np.ndarray, - resolution: float, -) -> list[list[float]]: - """유역 셀 집합을 폴리곤 외곽 링(사업지 좌표계)으로 벡터화한다.""" - try: - from rasterio import features as rio_features - from rasterio.transform import from_origin - except Exception: # noqa: BLE001 - return [] - mask = (labels == label).astype(np.uint8) - if mask.sum() == 0: - return [] - transform = from_origin( - float(x_coords[0]) - resolution / 2.0, - float(y_coords[-1]) + resolution / 2.0, - resolution, - resolution, +def _uphill_sign_at(vertices: list[Any], chainage_m: float, contour_index: ContourIndex) -> int: + """해당 측점의 산측이 도로 진행방향 기준 좌(+1)인지 우(-1)인지. 불명이면 0.""" + x, y, _ = _interpolate_vertex(vertices, chainage_m) + back = _interpolate_vertex(vertices, max(0.0, chainage_m - 10.0)) + forward = _interpolate_vertex(vertices, chainage_m + 10.0) + dx, dy = forward[0] - back[0], forward[1] - back[1] + norm = math.hypot(dx, dy) + if norm < 1e-6: + return 0 + dx, dy = dx / norm, dy / norm + # 좌측 법선 (-dy, dx) 방향 오프셋이 side_sign +1에 대응한다. + left_z = contour_index.nearest_elevation( + Point(x - dy * UPHILL_PROBE_OFFSET_M, y + dx * UPHILL_PROBE_OFFSET_M), + UPHILL_PROBE_RADIUS_M, ) - shapes = rio_features.shapes(mask[::-1, :], mask=mask[::-1, :] > 0, transform=transform) - polygons = [shape(geometry) for geometry, value in shapes if value == 1] - if not polygons: + right_z = contour_index.nearest_elevation( + Point(x + dy * UPHILL_PROBE_OFFSET_M, y - dx * UPHILL_PROBE_OFFSET_M), + UPHILL_PROBE_RADIUS_M, + ) + if left_z is None or right_z is None or left_z == right_z: + return 0 + return 1 if left_z > right_z else -1 + + +def _road_segment_coords( + vertices: list[Any], start_m: float, end_m: float +) -> list[tuple[float, float]]: + """분할점 사이 도로 구간의 평면 좌표열(유역 폴리곤의 하측 경계).""" + sx, sy, _ = _interpolate_vertex(vertices, start_m) + ex, ey, _ = _interpolate_vertex(vertices, end_m) + coords = [(sx, sy)] + coords.extend( + (vertex.x, vertex.y) for vertex in vertices if start_m < vertex.chainage_m < end_m + ) + coords.append((ex, ey)) + return coords + + +def _contour_arc( + line: Any, p_from: Point, p_to: Point, road_line: LineString +) -> list[tuple[float, float]]: + """등고선에서 두 분할선 접점 사이 아크(상측 경계)를 뽑는다. + + 폐합 등고선은 두 방향 아크가 생기므로 도로와 교차하지 않는(=산측) 쪽을 고른다. + """ + t1, t2 = sorted((line.project(p_from), line.project(p_to))) + arcs = [] + inner = substring(line, t1, t2) + if inner.geom_type == "LineString" and len(inner.coords) >= 2: + arcs.append(inner) + if getattr(line, "is_closed", False): + head = substring(line, t2, line.length) + tail = substring(line, 0.0, t1) + coords = list(head.coords) + list(tail.coords)[1:] + if len(coords) >= 2: + arcs.append(LineString(coords)) + if not arcs: return [] - merged = max(polygons, key=lambda polygon: polygon.area) - simplified = merged.simplify(resolution, preserve_topology=True) - if simplified.is_empty or simplified.geom_type != "Polygon": - simplified = merged - return [[float(x), float(y)] for x, y in simplified.exterior.coords] + scored = [] + for arc in arcs: + crosses = arc.crosses(road_line) + midpoint = arc.interpolate(0.5, normalized=True) + scored.append((crosses, -road_line.distance(midpoint), arc)) + scored.sort(key=lambda item: (item[0], item[1])) + arc = scored[0][2] + coords = list(arc.coords) + if Point(coords[0]).distance(p_from) > Point(coords[-1]).distance(p_from): + coords.reverse() + return [(float(x), float(y)) for x, y in coords] + + +def _junction(left: list[DividerStep], right: list[DividerStep]) -> tuple[int, int, int] | None: + """두 분할선이 같은 등고선 지오메트리를 밟은 최고 표고 지점(좌 idx, 우 idx, geom idx).""" + left_keys = { + (step.z, step.geom_index): position + for position, step in enumerate(left) + if step.geom_index >= 0 + } + best: tuple[float, int, int, int] | None = None + for position, step in enumerate(right): + if step.geom_index < 0: + continue + left_position = left_keys.get((step.z, step.geom_index)) + if left_position is None: + continue + if best is None or step.z > best[0]: + best = (step.z, left_position, position, step.geom_index) + if best is None: + return None + return best[1], best[2], best[3] + + +def _assemble_polygon( + vertices: list[Any], + start_m: float, + end_m: float, + left: list[DividerStep], + right: list[DividerStep], + contour_index: ContourIndex, + road_line: LineString, +) -> Polygon | None: + """도로 구간 + 우측 분할선 + 상측 등고선 아크 + 좌측 분할선으로 폴리곤을 폐합한다.""" + ring = _road_segment_coords(vertices, start_m, end_m) + junction = _junction(left, right) + if junction is not None: + left_position, right_position, geom_index = junction + left_used = left[: left_position + 1] + right_used = right[: right_position + 1] + arc = _contour_arc( + contour_index.geoms[geom_index], + right_used[-1].point, + left_used[-1].point, + road_line, + ) + else: + left_used, right_used, arc = left, right, [] + ring.extend((step.point.x, step.point.y) for step in right_used[1:]) + ring.extend(arc) + ring.extend((step.point.x, step.point.y) for step in reversed(left_used[1:])) + if len(ring) < 4: + return None + polygon = Polygon(ring).buffer(0) + if polygon.geom_type == "MultiPolygon": + polygon = max(polygon.geoms, key=lambda part: part.area) + if polygon.is_empty or polygon.geom_type != "Polygon": + return None + return polygon def build_watershed_basins( @@ -416,64 +217,116 @@ def build_watershed_basins( contour_features: list[dict[str, Any]], spot_features: list[dict[str, Any]], elevation_keys: tuple[str, ...], + stream_features: list[dict[str, Any]] | None = None, ) -> list[WatershedBasin]: - """도로(노선) 기준 배수유역을 산정한다 — 측구 흐름 모델. + """관 지점 배치를 기준으로 메인 배수유역을 세분화해 산정한다. - "도로에 닿은 물은 측구를 타고 종단 내리막 방향으로 흘러 첫 번째 관으로 들어간다"를 - 전제로, 도로 산측 사면 전체를 관 개수만큼 빈틈없이 분할한다(참고 도면과 동일 개념): - ① 노선 전체를 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로 - 각 도로 셀에 담당 관 라벨을 붙인다. - ② 사면 각 셀은 D8 흐름을 따라 내려가 처음 닿는 도로 셀의 관 라벨을 물려받는다. - 도로를 만나지 못하고 격자 밖으로 빠지는 셀(도로 하측 성토부 등)은 제외된다. - 반환된 boundary_xy 외곽선의 산측이 곧 분수령(능선)이며, 하측은 도로선을 따른다. + 세류 교차 관("stream")은 상류망을 추적해 넓은 한계로, 세류 없는 관은 도로 상측 + 첫 능선까지 소범위 한계로 분할선을 올린다(작은 유역, 영역 선정 주의 — 사용자 지시). 번호는 노선 시점에 가까운 순. """ if not candidates or len(vertices) < 2: return [] - samples = _collect_samples(contour_features, spot_features, elevation_keys) - 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)) + contour_index = ContourIndex(contour_features, elevation_keys) + if contour_index.tree is None: + logger.warning("표고 속성이 있는 등고선이 없어 유역을 산정하지 못했습니다.") return [] - dem, x_coords, y_coords, resolution = built - dem = _fill_depressions(dem, resolution) - pointer = _d8_pointer(dem) + spot_index = ContourIndex(spot_features, elevation_keys) + road_line = LineString([(vertex.x, vertex.y) for vertex in vertices]) ordered = sorted(candidates, key=lambda item: item.chainage_m) - 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) - lengths = _flow_lengths(pointer, labels, resolution) + divides = _divide_chainages(vertices, ordered) + signs = [_uphill_sign_at(vertices, item.chainage_m, contour_index) for item in ordered] + majority = 1 if sum(signs) >= 0 else -1 + signs = [sign or majority for sign in signs] + # 사용자 확정("confirmed") 측점도 세류에 닿아 있으면 세류 유역으로 취급한다. + stream_lines = _explode_lines(stream_features) if stream_features else [] + is_stream = [ + item.reason == "stream" + or any(line.distance(Point(item.x, item.y)) <= STREAM_JOIN_TOL_M for line in stream_lines) + for item in ordered + ] + + # 분할선은 인접 유역과 공유하므로 분할점마다 1회만 추적한다. + dividers: list[list[DividerStep]] = [] + for position, chainage in enumerate(divides): + neighbor_streams = [] + if position > 0: + neighbor_streams.append(is_stream[position - 1]) + if position < len(ordered): + neighbor_streams.append(is_stream[position]) + wide = any(neighbor_streams) + sign = signs[position - 1] if position > 0 else signs[0] + x, y, _ = _interpolate_vertex(vertices, chainage) + if wide: + steps = trace_divider(Point(x, y), contour_index, road_line, sign) + else: + steps = trace_divider( + Point(x, y), + contour_index, + road_line, + sign, + radius_m=LOCAL_SEARCH_RADIUS_M, + max_steps=LOCAL_MAX_STEPS, + ) + dividers.append(steps) basins: list[WatershedBasin] = [] - for label, candidate in enumerate(ordered, start=1): - mask = labels == label - cell_count = int(mask.sum()) - if cell_count < 4: + for position, candidate in enumerate(ordered): + outlet = Point(candidate.x, candidate.y) + network: list[Any] = [] + flow_length = 0.0 + if is_stream[position] and stream_features: + network, flow_length = trace_upstream_network( + outlet, stream_features, road_line, signs[position] + ) + polygon = _assemble_polygon( + vertices, + divides[position], + divides[position + 1], + dividers[position], + dividers[position + 1], + contour_index, + road_line, + ) + if polygon is None: continue - boundary = _vectorize_basin(labels, label, x_coords, y_coords, resolution) - if len(boundary) < 4: + if network: + # 상류망이 폴리곤 밖으로 뻗은 경우까지 유역이 덮도록 보정한다. + covered = polygon.union(unary_union(network).buffer(STREAM_COVER_BUFFER_M)).buffer(0) + if covered.geom_type == "MultiPolygon": + covered = max(covered.geoms, key=lambda part: part.area) + if covered.geom_type == "Polygon" and not covered.is_empty: + polygon = covered + if polygon.area < MIN_BASIN_AREA_M2: continue - # 측점(관) 위치의 표고를 기준으로 낙차를 계산한다. - 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()) + + outlet_z = contour_index.nearest_elevation(outlet, UPHILL_PROBE_RADIUS_M) + if outlet_z is None: + outlet_z = _interpolate_vertex(vertices, candidate.chainage_m)[2] + top_z = max( + contour_index.max_elevation_within(polygon) or outlet_z, + spot_index.max_elevation_within(polygon) or outlet_z, + ) + boundary_line = polygon.simplify(5.0, preserve_topology=True) + if boundary_line.is_empty or boundary_line.geom_type != "Polygon": + boundary_line = polygon + boundary = [[float(x), float(y)] for x, y in boundary_line.exterior.coords] + if flow_length <= 0.0: + flow_length = max( + (math.dist((candidate.x, candidate.y), point) for point in boundary), + default=0.0, + ) + basin = WatershedBasin( - index=label, + index=len(basins) + 1, chainage_m=candidate.chainage_m, outlet_x=candidate.x, outlet_y=candidate.y, boundary_xy=boundary, - area_m2=cell_count * resolution * resolution, - relief_m=max(0.0, float(dem[mask].max()) - outlet_z), - flow_length_m=lengths.get(label, 0.0), + area_m2=float(polygon.area), + relief_m=max(0.0, float(top_z) - float(outlet_z)), + flow_length_m=float(flow_length), ) basin.pipe_diameter_mm = estimate_pipe_diameter_mm( basin.area_m2, basin.relief_m, basin.flow_length_m diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py new file mode 100644 index 00000000..842a477a --- /dev/null +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py @@ -0,0 +1,271 @@ +"""배수유역 추적 유틸 — 등고선 공간 인덱스·세류 상류망·분수령(능선) 추적. + +등고선 기하 직접 분석(2026-07-29 합의)의 하위 도구 모음. DEM 보간 없이: +- `ContourIndex`: 등고선(선)·표고점(점)을 STRtree에 1회 적재하고 필요한 것만 꺼낸다. +- `trace_upstream_network`: 세류 교차점에서 도로 산측 상류망만 추적한다(하류 무시). +- `trace_divider`: 물갈림 지점에서 상향 등고선을 한 겹씩 따라 오르는 유역 분할선(능선 근사). + +전체 등고선을 순회하는 연산을 두지 않아 분석량이 유역 크기에 비례한다(도엽 매수와 무관). +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any + +from shapely.geometry import LineString, Point, shape +from shapely.ops import nearest_points, substring +from shapely.strtree import STRtree + +# 분할선(능선 근사) 추적: 다음 상위 등고선을 찾는 탐색 반경(m)과 최대 단계 수. +DIVIDER_SEARCH_RADIUS_M = 120.0 +DIVIDER_MAX_STEPS = 60 +# 세류 없는 소규모 유역: 도로 상측 첫 능선까지만 오르도록 좁힌 한계(영역 선정 주의). +LOCAL_SEARCH_RADIUS_M = 80.0 +LOCAL_MAX_STEPS = 12 +# 세류 연결 판정 이격(m)과 상류망 총연장 상한(m). +STREAM_JOIN_TOL_M = 15.0 +MAX_UPSTREAM_TOTAL_M = 5000.0 + + +@dataclass +class DividerStep: + """분할선의 한 단계 — 어느 등고선(geom_index)의 어느 지점을 밟았는지.""" + + point: Point + z: float + geom_index: int + + +class ContourIndex: + """표고 속성이 있는 등고선·표고점 피처의 STRtree 래퍼.""" + + def __init__(self, features: list[dict[str, Any]], elevation_keys: tuple[str, ...]) -> None: + self.geoms: list[Any] = [] + self.zs: list[float] = [] + for feature in features: + elevation = _feature_elevation(feature, elevation_keys) + if elevation is None: + continue + geometry = feature.get("geometry") or {} + try: + geom = shape(geometry) + except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다 + continue + if geom.is_empty: + continue + parts = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] + for part in parts: + self.geoms.append(part) + self.zs.append(elevation) + self.tree = STRtree(self.geoms) if self.geoms else None + + def query(self, geometry: Any) -> list[int]: + """geometry 근방(bbox 교차) 피처의 인덱스만 돌려준다.""" + if self.tree is None: + return [] + return [int(i) for i in self.tree.query(geometry)] + + def nearest_elevation(self, point: Point, radius_m: float) -> float | None: + """point에서 radius 안 가장 가까운 피처의 표고. 없으면 None.""" + best_z: float | None = None + best_distance = radius_m + for index in self.query(point.buffer(radius_m)): + distance = self.geoms[index].distance(point) + if distance <= best_distance: + best_distance = distance + best_z = self.zs[index] + return best_z + + def max_elevation_within(self, polygon: Any) -> float | None: + """polygon과 실제로 교차하는 피처들의 최고 표고.""" + best: float | None = None + for index in self.query(polygon): + if not polygon.intersects(self.geoms[index]): + continue + if best is None or self.zs[index] > best: + best = self.zs[index] + return best + + +def _feature_elevation(feature: dict[str, Any], elevation_keys: tuple[str, ...]) -> float | None: + properties = feature.get("properties") or {} + 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 side_sign(road_line: LineString, point: Point) -> int: + """도로선 기준 point가 어느 쪽인지(+1/-1, 선상이면 0). 국소 접선과의 외적 부호.""" + t = road_line.project(point) + a = road_line.interpolate(max(0.0, t - 5.0)) + b = road_line.interpolate(min(road_line.length, t + 5.0)) + cross = (b.x - a.x) * (point.y - a.y) - (b.y - a.y) * (point.x - a.x) + if cross > 0: + return 1 + if cross < 0: + return -1 + return 0 + + +def trace_divider( + start: Point, + contour_index: ContourIndex, + road_line: LineString, + uphill_sign: int, + radius_m: float = DIVIDER_SEARCH_RADIUS_M, + max_steps: int = DIVIDER_MAX_STEPS, +) -> list[DividerStep]: + """물갈림 지점에서 상향 등고선을 한 겹씩 밟아 오르는 분할선을 만든다. + + 각 단계에서 반경 안의 "현재보다 높은 등고선 중 가장 낮은 것"의 최근접점으로 이동한다. + 도로 산측(uphill_sign)을 벗어나거나 도로 쪽으로 되돌아가는 이동은 막는다. + 더 높은 등고선이 반경 안에 없으면 능선(분수령)에 닿은 것으로 보고 멈춘다. + """ + z = contour_index.nearest_elevation(start, radius_m) + if z is None: + return [] + steps = [DividerStep(point=start, z=z, geom_index=-1)] + current = start + road_distance = road_line.distance(start) + for _ in range(max_steps): + best: tuple[float, float, Point, int] | None = None + for index in contour_index.query(current.buffer(radius_m)): + candidate_z = contour_index.zs[index] + if candidate_z <= z + 0.01: + continue + if best is not None and candidate_z > best[0]: + continue + point = nearest_points(contour_index.geoms[index], current)[0] + distance = current.distance(point) + if distance > radius_m: + continue + # 도로 반대편·도로 방향 후퇴 금지 — 분할선은 산측으로만 오른다. + if side_sign(road_line, point) == -uphill_sign: + continue + if road_line.distance(point) + 1.0 < road_distance: + continue + if ( + best is None + or candidate_z < best[0] + or (candidate_z == best[0] and distance < best[1]) + ): + best = (candidate_z, distance, point, index) + if best is None: + break + z, _, current, geom_index = best + road_distance = max(road_distance, road_line.distance(current)) + steps.append(DividerStep(point=current, z=z, geom_index=geom_index)) + return steps + + +def _explode_lines(stream_features: list[dict[str, Any]]) -> list[LineString]: + lines: list[LineString] = [] + for feature in stream_features: + geometry = feature.get("geometry") or {} + try: + geom = shape(geometry) + except Exception: # noqa: BLE001 + continue + if geom.is_empty: + continue + parts = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] + lines.extend(part for part in parts if part.geom_type == "LineString") + return lines + + +def _oriented_from(line: LineString, origin: Point) -> LineString: + """origin에 가까운 끝이 시작점이 되도록 방향을 맞춘다.""" + if Point(line.coords[0]).distance(origin) <= Point(line.coords[-1]).distance(origin): + return line + return LineString(list(line.coords)[::-1]) + + +def _clip_uphill(line: LineString, road_line: LineString, origin: Point) -> LineString | None: + """도로를 다시 가로지르면 교차 지점에서 잘라 origin 쪽 조각만 남긴다.""" + if not line.crosses(road_line): + return line + t = line.project(nearest_points(line.intersection(road_line), origin)[0]) + piece = substring(line, 0.0, t) if line.project(origin) < t else substring(line, t, line.length) + if piece.geom_type != "LineString" or piece.length < 1.0: + return None + return piece + + +def trace_upstream_network( + crossing: Point, + stream_features: list[dict[str, Any]], + road_line: LineString, + uphill_sign: int, +) -> tuple[list[LineString], float]: + """세류 교차점에서 도로 산측으로 뻗는 상류망을 추적한다. + + ① 교차한 세류를 교차점에서 잘라 산측 조각을 뿌리로 삼는다. + ② 끝점이 기존 망에 근접(STREAM_JOIN_TOL_M)한 세류를 반복 편입한다(분기 포함). + 도로를 다시 가로지르는 조각은 절단하고, 총연장 상한을 두어 폭주를 막는다. + 반환: (상류망 폴리라인 목록, 최장 유하 경로 길이 m). + """ + lines = _explode_lines(stream_features) + network: list[tuple[LineString, float]] = [] # (폴리라인, 뿌리에서 시작점까지 누적거리) + used: set[int] = set() + total = 0.0 + + # ① 뿌리: 교차점을 지나는 세류의 산측 조각. + for index, line in enumerate(lines): + if line.distance(crossing) > 1.0: + continue + t = line.project(crossing) + for piece in (substring(line, 0.0, t), substring(line, t, line.length)): + if piece.geom_type != "LineString" or piece.length < STREAM_JOIN_TOL_M: + continue + oriented = _oriented_from(piece, crossing) + probe = oriented.interpolate(min(10.0, oriented.length)) + if side_sign(road_line, probe) != uphill_sign: + continue + clipped = _clip_uphill(oriented, road_line, crossing) + if clipped is None: + continue + network.append((clipped, 0.0)) + total += clipped.length + used.add(index) + if not network: + return [], 0.0 + + # ② 편입 반복: 끝점이 망에 닿는 세류를 상류로 붙인다. + grew = True + while grew and total < MAX_UPSTREAM_TOTAL_M: + grew = False + for index, line in enumerate(lines): + if index in used: + continue + attach: tuple[float, Point, LineString, float] | None = None + for endpoint in (Point(line.coords[0]), Point(line.coords[-1])): + for parent, parent_cum in network: + distance = parent.distance(endpoint) + if distance > STREAM_JOIN_TOL_M: + continue + cum = parent_cum + parent.project(endpoint) + if attach is None or distance < attach[0]: + attach = (distance, endpoint, parent, cum) + if attach is None: + continue + used.add(index) + _, endpoint, _, cum = attach + if side_sign(road_line, line.interpolate(0.5, normalized=True)) == -uphill_sign: + continue + oriented = _oriented_from(line, endpoint) + clipped = _clip_uphill(oriented, road_line, endpoint) + if clipped is None: + continue + network.append((clipped, cum)) + total += clipped.length + grew = True + + flow_length = max((cum + line.length for line, cum in network), default=0.0) + return [line for line, _ in network], flow_length diff --git a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py index ecb3a54d..71979888 100644 --- a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py @@ -194,7 +194,12 @@ async def post_drainage_basins( candidates = propose_structure_stations(vertices, prepared["streams"]) basins = build_watershed_basins( - vertices, candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS + vertices, + candidates, + prepared["contours"], + prepared["spots"], + _ELEVATION_KEYS, + stream_features=prepared["streams"], ) to_lonlat = prepared["to_lonlat"] return {