From 067de3a25f221d76de3d3c30c7fc9fc826b665fd Mon Sep 17 00:00:00 2001 From: umsangdon Date: Wed, 29 Jul 2026 18:40:45 +0900 Subject: [PATCH] auto: 2026-07-29 18:40 (EOMSANGDON-HOME) --- ...B05_wf2_Route_Engine_Drainage_Watershed.py | 152 ++++++++++++++---- .../B05_wf2_Route_Engine_Watershed_Trace.py | 18 ++- 2 files changed, 137 insertions(+), 33 deletions(-) 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 70ed3d31..33327c66 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py @@ -19,7 +19,7 @@ from dataclasses import dataclass, field from typing import Any from shapely.geometry import LineString, Point, Polygon -from shapely.ops import split, substring, unary_union +from shapely.ops import nearest_points, split, substring, unary_union from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( StructureCandidate, @@ -45,6 +45,8 @@ MIN_BASIN_AREA_M2 = 100.0 STREAM_COVER_BUFFER_M = 20.0 # 도로 절단용 노선 양끝 연장 길이(m). 노선 끝 너머로 새는 하류측 조각까지 잘라낸다. ROAD_EXTEND_M = 500.0 +# 폐합 등고선 기하 탐색: 분할선 경로와 등고선의 근접 허용 거리(m). +CLOSING_SEARCH_M = 30.0 # 산측 판정: 도로 접선의 좌우 수직 오프셋(m)과 표고 조회 반경(m). UPHILL_PROBE_OFFSET_M = 40.0 UPHILL_PROBE_RADIUS_M = 60.0 @@ -124,11 +126,16 @@ def _road_segment_coords( def _contour_arc( - line: Any, p_from: Point, p_to: Point, road_line: LineString + line: Any, + p_from: Point, + p_to: Point, + road_line: LineString, + network_union: Any = None, ) -> list[tuple[float, float]]: """등고선에서 두 분할선 접점 사이 아크(상측 경계)를 뽑는다. - 폐합 등고선은 두 방향 아크가 생기므로 도로와 교차하지 않는(=산측) 쪽을 고른다. + 폐합 등고선은 두 방향 아크가 생기므로, 세류 상류망을 가로지르지 않고(계곡을 + 자르지 않고) 도로와도 교차하지 않는(=산측) 쪽을 고른다. """ t1, t2 = sorted((line.project(p_from), line.project(p_to))) arcs = [] @@ -145,36 +152,93 @@ def _contour_arc( return [] scored = [] for arc in arcs: - crosses = arc.crosses(road_line) + crosses_stream = bool(network_union is not None and arc.crosses(network_union)) + crosses_road = 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] + scored.append((crosses_stream, crosses_road, -road_line.distance(midpoint), arc)) + scored.sort(key=lambda item: (item[0], item[1], item[2])) + arc = scored[0][3] 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).""" +def _junction( + left: list[DividerStep], + right: list[DividerStep], + min_z: float | None = None, +) -> tuple[int, int, int] | None: + """두 분할선이 같은 등고선 지오메트리를 밟은 폐합 지점(좌 idx, 우 idx, geom idx). + + min_z(세류 상류망 최고 표고)가 있으면 그 **이상인 가장 낮은** 공통 등고선을 고른다 — + "세류로 영역을 지정한 뒤 가까운 등고선으로 바로 올려치면 안 된다"(2026-07-29 사용자 + 지시). 계곡 발원부를 넘긴 첫 등고선이 유역 상측 경계가 된다. 없으면 최고 공통 등고선. + """ 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 + matches: list[tuple[float, int, int, int]] = [] 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) + matches.append((step.z, left_position, position, step.geom_index)) + if not matches: + return None + if min_z is not None: + above = [match for match in matches if match[0] >= min_z] + if above: + best = min(above) + return best[1], best[2], best[3] + best = max(matches) + return best[1], best[2], best[3] + + +def _closing_contour( + left: list[DividerStep], + right: list[DividerStep], + contour_index: ContourIndex, + min_z: float, +) -> tuple[int, int, int, Point, Point] | None: + """두 분할선 경로에 모두 근접한 등고선 중 min_z 이상 최저를 찾는다. + + 분할선이 같은 스텝에서 같은 지오메트리를 밟지 못해도(도엽 분할 등) 계곡 발원부 + 위를 지나는 폐합 등고선을 기하적으로 찾아낸다. + 반환: (좌 절단 idx, 우 절단 idx, 등고선 geom idx, 좌 접점, 우 접점). + """ + if len(left) < 2 or len(right) < 2: + return None + left_line = LineString([step.point for step in left]) + right_line = LineString([step.point for step in right]) + shared = set(contour_index.query(left_line.buffer(CLOSING_SEARCH_M))) & set( + contour_index.query(right_line.buffer(CLOSING_SEARCH_M)) + ) + best: tuple[float, int] | None = None + for index in shared: + z = contour_index.zs[index] + if z < min_z: + continue + geom = contour_index.geoms[index] + if ( + geom.distance(left_line) > CLOSING_SEARCH_M + or geom.distance(right_line) > CLOSING_SEARCH_M + ): + continue + if best is None or z < best[0]: + best = (z, index) if best is None: return None - return best[1], best[2], best[3] + geom = contour_index.geoms[best[1]] + left_touch = nearest_points(geom, left_line)[0] + right_touch = nearest_points(geom, right_line)[0] + left_position = min(range(len(left)), key=lambda i: left[i].point.distance(left_touch)) + right_position = min(range(len(right)), key=lambda i: right[i].point.distance(right_touch)) + return left_position, right_position, best[1], left_touch, right_touch def _assemble_polygon( @@ -185,22 +249,41 @@ def _assemble_polygon( right: list[DividerStep], contour_index: ContourIndex, road_line: LineString, + network_union: Any = None, + valley_top_z: float | None = None, ) -> Polygon | None: - """도로 구간 + 우측 분할선 + 상측 등고선 아크 + 좌측 분할선으로 폴리곤을 폐합한다.""" + """도로 구간 + 우측 분할선 + 상측 등고선 아크 + 좌측 분할선으로 폴리곤을 폐합한다. + + 세류 유역(valley_top_z 지정)은 발원부 위를 지나는 폐합 등고선을 기하 탐색으로 + 먼저 찾고, 실패 시 같은 스텝 매칭(_junction)으로 폐합한다. + """ 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, right_used, arc = left, right, [] + closure = ( + _closing_contour(left, right, contour_index, valley_top_z) + if valley_top_z is not None + else None + ) + if closure is not None: + left_position, right_position, geom_index, left_touch, right_touch = closure 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, + contour_index.geoms[geom_index], right_touch, left_touch, road_line, network_union ) else: - left_used, right_used, arc = left, right, [] + junction = _junction(left, right, min_z=valley_top_z) + 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, + network_union, + ) 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:])) @@ -230,13 +313,15 @@ def _extended_road_line(vertices: list[Any]) -> LineString: def _clip_to_uphill( polygon: Polygon, - road_line: LineString, extended_road: LineString, uphill_sign: int, + keep_geom: Any = None, ) -> Polygon | None: """폴리곤을 도로선으로 절단해 산측 조각만 남긴다 — 도로가 유역의 한쪽 경계. 도로 하류측(성토부 아래)은 유역에 포함하지 않는다(2026-07-29 사용자 지시). + keep_geom(세류 상류망)이 걸린 조각은 부호와 무관하게 유지한다 — 상류망은 정의상 + 산측 배수인데, 노선 끝을 감아 도는 계곡은 연장 도로선 기준 부호가 뒤집힐 수 있다. """ try: pieces = split(polygon, extended_road) @@ -246,7 +331,11 @@ def _clip_to_uphill( piece for piece in getattr(pieces, "geoms", [pieces]) if piece.geom_type == "Polygon" - and side_sign(road_line, piece.representative_point()) == uphill_sign + and ( + side_sign(extended_road, piece.representative_point()) == uphill_sign + # 교차점(도로 위 시작점) 접촉만으로는 부족 — 상류망이 실제로 지나가야 한다. + or (keep_geom is not None and piece.intersection(keep_geom).length > 5.0) + ) ] if not kept: return None @@ -328,6 +417,13 @@ def build_watershed_basins( network, flow_length = trace_upstream_network( outlet, stream_features, road_line, signs[position] ) + # 계곡 발원부(상류망 최고 표고) — 유역 상측 폐합 등고선은 이보다 높아야 한다. + network_union = unary_union(network) if network else None + valley_top_z = ( + contour_index.max_elevation_within(network_union.buffer(10.0)) + if network_union is not None + else None + ) polygon = _assemble_polygon( vertices, divides[position], @@ -336,18 +432,20 @@ def build_watershed_basins( dividers[position + 1], contour_index, road_line, + network_union, + valley_top_z, ) if polygon is None: continue - if network: + if network_union is not None: # 상류망이 폴리곤 밖으로 뻗은 경우까지 유역이 덮도록 보정한다. - covered = polygon.union(unary_union(network).buffer(STREAM_COVER_BUFFER_M)).buffer(0) + covered = polygon.union(network_union.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 # 마지막에 도로선으로 절단해 하류측을 제거한다 — 도로가 유역의 한쪽 경계. - polygon = _clip_to_uphill(polygon, road_line, extended_road, signs[position]) + polygon = _clip_to_uphill(polygon, extended_road, signs[position], network_union) if polygon is None or polygon.area < MIN_BASIN_AREA_M2: continue 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 842a477a..2eb965ee 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py @@ -216,13 +216,20 @@ def trace_upstream_network( used: set[int] = set() total = 0.0 - # ① 뿌리: 교차점을 지나는 세류의 산측 조각. + # ① 뿌리: 교차점을 지나는 세류의 산측 조각. 도엽 세류는 교차점 부근에서 별도 + # 피처로 조각나 있는 경우가 많아, 근접(STREAM_JOIN_TOL_M) 조각도 뿌리로 받는다. for index, line in enumerate(lines): - if line.distance(crossing) > 1.0: + distance = line.distance(crossing) + if distance > STREAM_JOIN_TOL_M: 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: + used.add(index) + if distance <= 1.0: + t = line.project(crossing) + pieces = [substring(line, 0.0, t), substring(line, t, line.length)] + else: + pieces = [line] + for piece in pieces: + if piece.geom_type != "LineString" or piece.length < 1.0: continue oriented = _oriented_from(piece, crossing) probe = oriented.interpolate(min(10.0, oriented.length)) @@ -233,7 +240,6 @@ def trace_upstream_network( continue network.append((clipped, 0.0)) total += clipped.length - used.add(index) if not network: return [], 0.0