From 5354e722580291c937fbd00ab5dd5fa767222126 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Wed, 29 Jul 2026 19:31:38 +0900 Subject: [PATCH] auto: 2026-07-29 19:31 (EOMSANGDON-HOME) --- .../B05_wf2_Route_Engine_Watershed_Trace.py | 81 ++++++++++++++----- 1 file changed, 59 insertions(+), 22 deletions(-) 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 7f880eab..c5877c23 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py @@ -10,6 +10,7 @@ from __future__ import annotations +import bisect from dataclasses import dataclass from typing import Any @@ -256,23 +257,42 @@ def valley_region_polygon( 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) + region = _contour_chain_boundary(region, contour_index, network_union, others, other_tree) return region -def _snap_boundary_to_ridge( +def _collect_intersection_points(geometry: Any) -> list[Point]: + """교차 결과에서 대표 점들을 뽑는다.""" + if geometry.is_empty: + return [] + if geometry.geom_type == "Point": + return [geometry] + if geometry.geom_type in {"MultiPoint", "GeometryCollection"}: + points: list[Point] = [] + for part in geometry.geoms: + points.extend(_collect_intersection_points(part)) + return points + if geometry.geom_type in {"LineString", "MultiLineString"}: + return [geometry.interpolate(0.5, normalized=True)] + return [] + + +def _contour_chain_boundary( region: Any, contour_index: Any, network_union: Any, others: list[LineString], other_tree: STRtree | None, ) -> Any: - """등거리 경계 정점을 등고선 위 능선 꼭짓점으로 스냅한다. + """유역 경계를 **등고선마다 능선 포인트 1개씩 찍어 연결**한 체인으로 재구성한다. - 능선 꼭짓점 = 근처 등고선을 따라 걸었을 때 두 세류(우리 상류망·인접 세류) - 모두로부터의 최소거리가 최대가 되는 지점. 우리 세류 침범(근접)과 인접 세류 - 이탈(인접이 더 가까워짐)은 금지한다. 실패 시 원본 경계를 유지한다. + (2026-07-29 사용자 지시: 정점 스냅은 점이 듬성듬성해 등고선을 건너뛴다.) + 등거리 1차 경계 링은 순서 뼈대로만 쓴다: 링을 가로지르는 모든 등고선 교차점마다 + 그 등고선 위에서 두 세류(우리 상류망·인접 세류) 모두로부터 가장 먼 지점(능선 + 꼭짓점)을 정제해 포인트를 얻고, 링 위 위치 순으로 연결한다. 등고선이 없는 구간은 + 원래 링 정점으로 메운다. 제약: 우리 세류 침범·인접 세류 이탈 금지. """ + ring = LineString(region.exterior.coords) def _score(point: Point) -> tuple[float, float]: to_network = network_union.distance(point) @@ -283,22 +303,26 @@ def _snap_boundary_to_ridge( ) 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) + # 링을 가로지르는 등고선 교차점마다 능선 꼭짓점 1개. + chained: list[tuple[float, float, float]] = [] # (링 위치 s, x, y) + for index in contour_index.query(ring.buffer(1.0)): + geom = contour_index.geoms[index] + try: + crossings = _collect_intersection_points(geom.intersection(ring)) + except Exception: # noqa: BLE001 + continue + for crossing in crossings: + s = ring.project(crossing) + t0 = geom.project(crossing) + best = None + best_value = -1.0 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)]: + for offset in [0.0] + [ + sign * k * 10.0 for k in range(1, steps + 1) for sign in (1, -1) + ]: t = min(max(t0 + offset, 0.0), geom.length) candidate = geom.interpolate(t) - if candidate.distance(vertex) > RIDGE_SNAP_M + RIDGE_WALK_M: + if candidate.distance(crossing) > RIDGE_WALK_M + RIDGE_SNAP_M: continue to_network, to_other = _score(candidate) if to_network < STREAM_JOIN_TOL_M: @@ -308,10 +332,23 @@ def _snap_boundary_to_ridge( 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: + if best is not None: + chained.append((s, best.x, best.y)) + if len(chained) < 4: return region - polygon = Polygon(snapped).buffer(0) + chained.sort() + # 등고선 공백 구간(교차점 사이가 먼 곳)은 원래 링 정점으로 메운다. + positions = [s for s, _, _ in chained] + filled: list[tuple[float, float, float]] = list(chained) + for x, y in list(region.exterior.coords)[:-1]: + s = ring.project(Point(x, y)) + slot = bisect.bisect_left(positions, s) + before = positions[slot - 1] if slot > 0 else positions[-1] - ring.length + after = positions[slot] if slot < len(positions) else positions[0] + ring.length + if min(s - before, after - s) > VALLEY_CELL_M * 2.5: + filled.append((s, x, y)) + filled.sort() + polygon = Polygon([(x, y) for _, x, y in filled]).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":