auto: 2026-07-29 18:44 (EOMSANGDON-HOME)

B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py

@@ -47,6 +47,9 @@ STREAM_COVER_BUFFER_M = 20.0
 ROAD_EXTEND_M = 500.0
 # 폐합 등고선 기하 탐색: 분할선 경로와 등고선의 근접 허용 거리(m).
 CLOSING_SEARCH_M = 30.0
+# 세류 포함 폐합 등고선 탐색: 상류망과 등고선의 근접 허용 거리(m)와 커버리지 목표.
+CLOSING_NEAR_M = 200.0
+CLOSING_COVERAGE_GOAL = 0.95
 # 산측 판정: 도로 접선의 좌우 수직 오프셋(m)과 표고 조회 반경(m).
 UPHILL_PROBE_OFFSET_M = 40.0
 UPHILL_PROBE_RADIUS_M = 60.0
@@ -241,6 +244,103 @@ def _closing_contour(
     return left_position, right_position, best[1], left_touch, right_touch
 
 
+def _both_arcs(line: Any, p_from: Point, p_to: Point) -> list[list[tuple[float, float]]]:
+    """등고선 위 두 접점 사이의 가능한 아크(양방향) 좌표열. p_from에서 시작하도록 정렬."""
+    t1, t2 = sorted((line.project(p_from), line.project(p_to)))
+    raw: list[list[Any]] = []
+    inner = substring(line, t1, t2)
+    if inner.geom_type == "LineString" and len(inner.coords) >= 2:
+        raw.append(list(inner.coords))
+    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:
+            raw.append(coords)
+    oriented: list[list[tuple[float, float]]] = []
+    for coords in raw:
+        if Point(coords[0]).distance(p_from) > Point(coords[-1]).distance(p_from):
+            coords = coords[::-1]
+        oriented.append([(float(x), float(y)) for x, y in coords])
+    return oriented
+
+
+def _steps_until_near(steps: list[DividerStep], point: Point) -> list[DividerStep]:
+    """분할선을 point에 가장 가까운 스텝까지 자른다."""
+    if not steps:
+        return []
+    position = min(range(len(steps)), key=lambda i: steps[i].point.distance(point))
+    return steps[: position + 1]
+
+
+def _assemble_stream_polygon(
+    vertices: list[Any],
+    start_m: float,
+    end_m: float,
+    left: list[DividerStep],
+    right: list[DividerStep],
+    contour_index: ContourIndex,
+    network_union: Any,
+    valley_top_z: float,
+) -> tuple[Polygon, float] | None:
+    """세류 상류망을 **포함하는 등고선**을 유역 상측 경계로 폐합한다.
+
+    "세류가 중심이 되면 안 되고, 세류를 포함하는 등고라인이 경계가 되어야 한다"
+    (2026-07-29 사용자 지시). 발원부(valley_top_z) 이상이고 세류망을 가로지르지 않으며
+    세류망에 근접한 등고선을 낮은 것부터 시도해, 상류망을 가장 잘 덮는 아크 폐합을 고른다.
+    반환: (폴리곤, 상류망 커버리지 0~1). 후보가 없으면 None.
+    """
+    candidates: list[tuple[float, int]] = []
+    for index in contour_index.query(network_union.buffer(CLOSING_NEAR_M)):
+        z = contour_index.zs[index]
+        if z < valley_top_z:
+            continue
+        geom = contour_index.geoms[index]
+        if geom.distance(network_union) > CLOSING_NEAR_M or geom.intersects(network_union):
+            continue
+        candidates.append((z, index))
+    if not candidates:
+        return None
+    candidates.sort()
+    road_coords = _road_segment_coords(vertices, start_m, end_m)
+    left_anchor: Any = (
+        LineString([step.point for step in left]) if len(left) > 1 else Point(road_coords[0])
+    )
+    right_anchor: Any = (
+        LineString([step.point for step in right]) if len(right) > 1 else Point(road_coords[-1])
+    )
+    best: tuple[float, Polygon] | None = None
+    for _, index in candidates[:12]:
+        geom = contour_index.geoms[index]
+        left_touch = nearest_points(geom, left_anchor)[0]
+        right_touch = nearest_points(geom, right_anchor)[0]
+        for arc in _both_arcs(geom, right_touch, left_touch):
+            ring: list[tuple[float, float]] = list(road_coords)
+            ring.extend(
+                (step.point.x, step.point.y) for step in _steps_until_near(right, right_touch)[1:]
+            )
+            ring.extend(arc)
+            ring.extend(
+                (step.point.x, step.point.y)
+                for step in reversed(_steps_until_near(left, left_touch)[1:])
+            )
+            if len(ring) < 4:
+                continue
+            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":
+                continue
+            coverage = network_union.intersection(polygon).length / max(network_union.length, 1.0)
+            if best is None or coverage > best[0]:
+                best = (coverage, polygon)
+        if best is not None and best[0] >= CLOSING_COVERAGE_GOAL:
+            break
+    if best is None:
+        return None
+    return best[1], best[0]
+
+
 def _assemble_polygon(
     vertices: list[Any],
     start_m: float,
@@ -424,21 +524,38 @@ def build_watershed_basins(
             if network_union is not None
             else None
         )
-        polygon = _assemble_polygon(
-            vertices,
-            divides[position],
-            divides[position + 1],
-            dividers[position],
-            dividers[position + 1],
-            contour_index,
-            road_line,
-            network_union,
-            valley_top_z,
-        )
+        polygon = None
+        coverage = 0.0
+        if network_union is not None and valley_top_z is not None:
+            # 세류 유역: 상류망을 포함하는 등고선 아크로 폐합(사용자 지시).
+            closed = _assemble_stream_polygon(
+                vertices,
+                divides[position],
+                divides[position + 1],
+                dividers[position],
+                dividers[position + 1],
+                contour_index,
+                network_union,
+                valley_top_z,
+            )
+            if closed is not None:
+                polygon, coverage = closed
+        if polygon is None:
+            polygon = _assemble_polygon(
+                vertices,
+                divides[position],
+                divides[position + 1],
+                dividers[position],
+                dividers[position + 1],
+                contour_index,
+                road_line,
+                network_union,
+                valley_top_z,
+            )
         if polygon is None:
             continue
-        if network_union is not None:
-            # 상류망이 폴리곤 밖으로 뻗은 경우까지 유역이 덮도록 보정한다.
+        if network_union is not None and coverage < CLOSING_COVERAGE_GOAL:
+            # 폐합 등고선이 못 덮은 상류망만 버퍼로 보정한다(최후 폴백).
             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)