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 5593952a..c7697a8e 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py @@ -462,7 +462,7 @@ def build_watershed_basins( valley_top_z, ) valley = ( - valley_region_polygon(network_union, stream_features or [], outlet) + valley_region_polygon(network_union, stream_features or [], outlet, contour_index) if network_union is not None else None ) diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py index 25110185..7f880eab 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py @@ -13,7 +13,7 @@ from __future__ import annotations from dataclasses import dataclass from typing import Any -from shapely.geometry import LineString, Point, box, shape +from shapely.geometry import LineString, Point, Polygon, box, shape from shapely.ops import nearest_points, substring, unary_union from shapely.strtree import STRtree @@ -29,6 +29,9 @@ MAX_UPSTREAM_TOTAL_M = 5000.0 # 계곡 유역 근사 격자: 셀 크기(m)와 상류망에서의 최대 이격(m). VALLEY_CELL_M = 12.0 VALLEY_CAP_M = 350.0 +# 능선 스냅: 경계 정점에서 등고선 탐색 반경(m)과 등고선 위 능선 꼭짓점 탐색 폭(m). +RIDGE_SNAP_M = 40.0 +RIDGE_WALK_M = 80.0 @dataclass @@ -205,12 +208,17 @@ def valley_region_polygon( network_union: Any, stream_features: list[dict[str, Any]], crossing: Point, + contour_index: Any = None, ) -> Any | None: - """상류망 계곡의 유역 영역 — 분수령(능선)을 인접 세류와의 등거리선으로 근사한다. + """상류망 계곡의 유역 영역 — 경계는 등고선을 참고한 능선(분수령)으로 긋는다. - 격자 셀 중심이 ① 인접 계곡 세류보다 우리 상류망에 가깝고 ② 상한 거리 이내이면 - 이 유역에 속한 것으로 본다. 지류 사이 사면(지능선)은 자연히 포함되고, 능선 너머 - 인접 계곡은 자동 제외된다. 반환: 폴리곤(간략화됨) 또는 None. + ① 등거리 1차 근사: 격자 셀 중심이 인접 계곡 세류보다 우리 상류망에 가깝고 상한 + 거리 이내이면 유역 소속. 지류 사이 사면(지능선) 포함, 인접 계곡 자동 제외. + ② 등고선 스냅(2026-07-29 사용자 지시): 경계는 세류들 사이 중간 어딘가가 아니라 + **등고선을 참고해** 그어야 한다 — 각 경계 정점을 근처 등고선 위에서 두 세류 + 모두로부터 가장 먼 지점(능선 꼭짓점)으로 이동. 단 우리 세류를 침범하거나 + 인접 세류 너머로 나가지 않는다. + 반환: 폴리곤(간략화됨) 또는 None. """ lines = _explode_lines(stream_features) others = [line for line in lines if line.distance(network_union) > STREAM_JOIN_TOL_M] @@ -247,9 +255,70 @@ def valley_region_polygon( region = region.simplify(VALLEY_CELL_M, preserve_topology=True) if region.is_empty or region.geom_type != "Polygon": return None + if contour_index is not None: + region = _snap_boundary_to_ridge(region, contour_index, network_union, others, other_tree) return region +def _snap_boundary_to_ridge( + region: Any, + contour_index: Any, + network_union: Any, + others: list[LineString], + other_tree: STRtree | None, +) -> Any: + """등거리 경계 정점을 등고선 위 능선 꼭짓점으로 스냅한다. + + 능선 꼭짓점 = 근처 등고선을 따라 걸었을 때 두 세류(우리 상류망·인접 세류) + 모두로부터의 최소거리가 최대가 되는 지점. 우리 세류 침범(근접)과 인접 세류 + 이탈(인접이 더 가까워짐)은 금지한다. 실패 시 원본 경계를 유지한다. + """ + + def _score(point: Point) -> tuple[float, float]: + to_network = network_union.distance(point) + to_other = ( + others[int(other_tree.nearest(point))].distance(point) + if other_tree is not None + else float("inf") + ) + return to_network, to_other + + snapped: list[tuple[float, float]] = [] + for x, y in list(region.exterior.coords)[:-1]: + vertex = Point(x, y) + best = vertex + best_net, best_other = _score(vertex) + best_value = min(best_net, best_other) + for index in contour_index.query(vertex.buffer(RIDGE_SNAP_M)): + geom = contour_index.geoms[index] + if geom.distance(vertex) > RIDGE_SNAP_M: + continue + t0 = geom.project(vertex) + steps = int(RIDGE_WALK_M / 10.0) + for offset in [0.0] + [s * 10.0 for k in range(1, steps + 1) for s in (k, -k)]: + t = min(max(t0 + offset, 0.0), geom.length) + candidate = geom.interpolate(t) + if candidate.distance(vertex) > RIDGE_SNAP_M + RIDGE_WALK_M: + continue + to_network, to_other = _score(candidate) + if to_network < STREAM_JOIN_TOL_M: + continue # 우리 세류 침범 금지 + if to_other < to_network: + continue # 인접 세류 쪽으로 이탈 금지 + if min(to_network, to_other) > best_value: + best_value = min(to_network, to_other) + best = candidate + snapped.append((best.x, best.y)) + if len(snapped) < 4: + return region + polygon = Polygon(snapped).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 region + return polygon + + def trace_upstream_network( crossing: Point, stream_features: list[dict[str, Any]],