auto: 2026-07-29 19:13 (EOMSANGDON-HOME)

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 nearest_points, split, substring, unary_union
+from shapely.ops import nearest_points, substring, unary_union
 
 from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import (
     StructureCandidate,
@@ -36,6 +36,7 @@ from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Trace import (
     side_sign,
     trace_divider,
     trace_upstream_network,
+    valley_region_polygon,
 )
 
 logger = logging.getLogger(__name__)
@@ -43,12 +44,11 @@ logger = logging.getLogger(__name__)
 # 유효 유역 최소 면적(m²)과 상류망 커버 보정 버퍼(m).
 MIN_BASIN_AREA_M2 = 100.0
 STREAM_COVER_BUFFER_M = 20.0
-# 도로 절단용 노선 양끝 연장 길이(m). 노선 끝 너머로 새는 하류측 조각까지 잘라낸다.
-ROAD_EXTEND_M = 500.0
+# 도로 양끝에서 하류측으로 뻗는 절단 차단선 길이(m).
+DOWNHILL_BARRIER_M = 800.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
@@ -244,103 +244,6 @@ 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,
@@ -397,49 +300,72 @@ def _assemble_polygon(
     return polygon
 
 
-def _extended_road_line(vertices: list[Any]) -> LineString:
-    """양끝 접선 방향으로 연장한 도로선 — 폴리곤 절단이 끝에서 끊기지 않게 한다."""
-    coords = [(vertex.x, vertex.y) for vertex in vertices]
-    hx, hy = coords[0]
-    dx, dy = coords[1][0] - hx, coords[1][1] - hy
-    norm = math.hypot(dx, dy) or 1.0
-    head = (hx - dx / norm * ROAD_EXTEND_M, hy - dy / norm * ROAD_EXTEND_M)
-    tx, ty = coords[-1]
-    dx, dy = tx - coords[-2][0], ty - coords[-2][1]
-    norm = math.hypot(dx, dy) or 1.0
-    tail = (tx + dx / norm * ROAD_EXTEND_M, ty + dy / norm * ROAD_EXTEND_M)
-    return LineString([head, *coords, tail])
-
-
 def _clip_to_uphill(
     polygon: Polygon,
-    extended_road: LineString,
+    road_line: LineString,
     uphill_sign: int,
+    contour_index: ContourIndex,
     keep_geom: Any = None,
 ) -> Polygon | None:
-    """폴리곤을 도로선으로 절단해 산측 조각만 남긴다 — 도로가 유역의 한쪽 경계.
+    """도로 하류측 조각을 잘라낸다 — 도로가 유역의 한쪽 경계(2026-07-29 사용자 지시).
 
-    도로 하류측(성토부 아래)은 유역에 포함하지 않는다(2026-07-29 사용자 지시).
-    keep_geom(세류 상류망)이 걸린 조각은 부호와 무관하게 유지한다 — 상류망은 정의상
-    산측 배수인데, 노선 끝을 감아 도는 계곡은 연장 도로선 기준 부호가 뒤집힐 수 있다.
+    절단선 = 실제 도로선 + 양끝에서 **하류측으로 뻗는 수직 차단선** 2개. (후방 접선
+    연장은 곡선 노선에서 계곡 내부를 관통해 오절단을 일으켰다.) 노선 끝을 감아 도는
+    산측 사면은 남고, 도로 하류측 주머니만 분리된다.
+    조각 분류는 표고 기반: 대표점의 등고선 표고 > 최근접 도로 지점 표고 → 산측.
+    keep_geom(세류 상류망)이 실제로 지나가는 조각은 무조건 유지한다.
     """
+    length = road_line.length
+    cutters: list[LineString] = [road_line]
+    for t_end, t_inner in ((0.0, min(30.0, length)), (length, max(0.0, length - 30.0))):
+        end = road_line.interpolate(t_end)
+        inner = road_line.interpolate(t_inner)
+        dx, dy = end.x - inner.x, end.y - inner.y
+        norm = math.hypot(dx, dy) or 1.0
+        for nx, ny in ((-dy / norm, dx / norm), (dy / norm, -dx / norm)):
+            probe = Point(end.x + nx * 30.0, end.y + ny * 30.0)
+            if side_sign(road_line, probe) == -uphill_sign:
+                cutters.append(
+                    LineString(
+                        [
+                            (end.x, end.y),
+                            (end.x + nx * DOWNHILL_BARRIER_M, end.y + ny * DOWNHILL_BARRIER_M),
+                        ]
+                    )
+                )
+                break
+    # split()은 절단선이 폴리곤 경계와 겹치면(도로 = 유역 하측 경계) 동작하지 않는다.
+    # 얇은 스트립을 차감해 조각을 분리하고, 분류 후 buffer-교집합으로 원형을 복원한다.
     try:
-        pieces = split(polygon, extended_road)
+        strip = unary_union([cutter.buffer(0.5) for cutter in cutters])
+        separated = polygon.difference(strip)
     except Exception:  # noqa: BLE001 - 절단 실패 시 원본 유지
         return polygon
-    kept = [
-        piece
-        for piece in getattr(pieces, "geoms", [pieces])
-        if piece.geom_type == "Polygon"
-        and (
-            side_sign(extended_road, piece.representative_point()) == uphill_sign
-            # 교차점(도로 위 시작점) 접촉만으로는 부족 — 상류망이 실제로 지나가야 한다.
-            or (keep_geom is not None and piece.intersection(keep_geom).length > 5.0)
-        )
+    pieces = [
+        part
+        for part in (separated.geoms if separated.geom_type.startswith("Multi") else [separated])
+        if part.geom_type == "Polygon" and not part.is_empty
     ]
+    kept = []
+    for piece in pieces:
+        if keep_geom is not None and piece.intersection(keep_geom).length > 5.0:
+            kept.append(piece)
+            continue
+        representative = piece.representative_point()
+        piece_z = contour_index.nearest_elevation(representative, 2.0 * UPHILL_PROBE_RADIUS_M)
+        foot = road_line.interpolate(road_line.project(representative))
+        road_z = contour_index.nearest_elevation(foot, 2.0 * UPHILL_PROBE_RADIUS_M)
+        if piece_z is not None and road_z is not None and piece_z != road_z:
+            if piece_z > road_z:
+                kept.append(piece)
+            continue
+        # 표고 판정 불가(등고선 공백·동일 표고) 시에만 좌우 부호로 판정한다.
+        if side_sign(road_line, representative) == uphill_sign:
+            kept.append(piece)
     if not kept:
         return None
-    merged = unary_union(kept)
+    # 스트립 차감으로 깎인 0.5m를 되붙이되 원본 폴리곤 밖으로는 나가지 않는다.
+    merged = unary_union(kept).buffer(0.7).intersection(polygon).buffer(0)
     if merged.geom_type == "MultiPolygon":
         merged = max(merged.geoms, key=lambda part: part.area)
     if merged.is_empty or merged.geom_type != "Polygon":
@@ -469,7 +395,6 @@ def build_watershed_basins(
         return []
     spot_index = ContourIndex(spot_features, elevation_keys)
     road_line = LineString([(vertex.x, vertex.y) for vertex in vertices])
-    extended_road = _extended_road_line(vertices)
 
     ordered = sorted(candidates, key=lambda item: item.chainage_m)
     divides = _divide_chainages(vertices, ordered)
@@ -524,45 +449,43 @@ def build_watershed_basins(
             if network_union is not None
             else None
         )
-        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:
+        # 도로변 스트립(분할선 폐합) + 세류 계곡 영역(분수령=인접 세류 등거리)을 합친다.
+        base = _assemble_polygon(
+            vertices,
+            divides[position],
+            divides[position + 1],
+            dividers[position],
+            dividers[position + 1],
+            contour_index,
+            road_line,
+            network_union,
+            valley_top_z,
+        )
+        valley = (
+            valley_region_polygon(network_union, stream_features or [], outlet)
+            if network_union is not None
+            else None
+        )
+        if base is None and valley is None:
             continue
-        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)
-            if covered.geom_type == "Polygon" and not covered.is_empty:
-                polygon = covered
+        if base is not None and valley is not None:
+            merged = base.union(valley).buffer(0)
+            if merged.geom_type == "MultiPolygon":
+                merged = max(merged.geoms, key=lambda part: part.area)
+            polygon = merged if merged.geom_type == "Polygon" and not merged.is_empty else base
+        else:
+            polygon = base if base is not None else valley
+        if network_union is not None:
+            coverage = network_union.intersection(polygon).length / max(network_union.length, 1.0)
+            if 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)
+                if covered.geom_type == "Polygon" and not covered.is_empty:
+                    polygon = covered
         # 마지막에 도로선으로 절단해 하류측을 제거한다 — 도로가 유역의 한쪽 경계.
-        polygon = _clip_to_uphill(polygon, extended_road, signs[position], network_union)
+        polygon = _clip_to_uphill(polygon, road_line, signs[position], contour_index, network_union)
         if polygon is None or polygon.area < MIN_BASIN_AREA_M2:
             continue
 

B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py

@@ -13,8 +13,8 @@ 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.geometry import LineString, Point, box, shape
+from shapely.ops import nearest_points, substring, unary_union
 from shapely.strtree import STRtree
 
 # 분할선(능선 근사) 추적: 다음 상위 등고선을 찾는 탐색 반경(m)과 최대 단계 수.
@@ -26,6 +26,9 @@ LOCAL_MAX_STEPS = 12
 # 세류 연결 판정 이격(m)과 상류망 총연장 상한(m).
 STREAM_JOIN_TOL_M = 15.0
 MAX_UPSTREAM_TOTAL_M = 5000.0
+# 계곡 유역 근사 격자: 셀 크기(m)와 상류망에서의 최대 이격(m).
+VALLEY_CELL_M = 12.0
+VALLEY_CAP_M = 350.0
 
 
 @dataclass
@@ -198,6 +201,55 @@ def _clip_uphill(line: LineString, road_line: LineString, origin: Point) -> Line
     return piece
 
 
+def valley_region_polygon(
+    network_union: Any,
+    stream_features: list[dict[str, Any]],
+    crossing: Point,
+) -> Any | None:
+    """상류망 계곡의 유역 영역 — 분수령(능선)을 인접 세류와의 등거리선으로 근사한다.
+
+    격자 셀 중심이 ① 인접 계곡 세류보다 우리 상류망에 가깝고 ② 상한 거리 이내이면
+    이 유역에 속한 것으로 본다. 지류 사이 사면(지능선)은 자연히 포함되고, 능선 너머
+    인접 계곡은 자동 제외된다. 반환: 폴리곤(간략화됨) 또는 None.
+    """
+    lines = _explode_lines(stream_features)
+    others = [line for line in lines if line.distance(network_union) > STREAM_JOIN_TOL_M]
+    other_tree = STRtree(others) if others else None
+    min_x, min_y, max_x, max_y = network_union.buffer(VALLEY_CAP_M).bounds
+    cells = []
+    y = min_y
+    while y < max_y:
+        x = min_x
+        while x < max_x:
+            center = Point(x + VALLEY_CELL_M / 2.0, y + VALLEY_CELL_M / 2.0)
+            distance = network_union.distance(center)
+            if distance <= VALLEY_CAP_M:
+                if other_tree is not None:
+                    nearest = others[int(other_tree.nearest(center))]
+                    if nearest.distance(center) < distance:
+                        x += VALLEY_CELL_M
+                        continue
+                cells.append(box(x, y, x + VALLEY_CELL_M, y + VALLEY_CELL_M))
+            x += VALLEY_CELL_M
+        y += VALLEY_CELL_M
+    if not cells:
+        return None
+    region = unary_union(cells).buffer(0)
+    if region.geom_type == "MultiPolygon":
+        touching = [
+            part
+            for part in region.geoms
+            if part.intersects(network_union) or part.distance(crossing) < VALLEY_CELL_M * 2
+        ]
+        region = unary_union(touching) if touching else max(region.geoms, key=lambda p: p.area)
+        if region.geom_type == "MultiPolygon":
+            region = max(region.geoms, key=lambda p: p.area)
+    region = region.simplify(VALLEY_CELL_M, preserve_topology=True)
+    if region.is_empty or region.geom_type != "Polygon":
+        return None
+    return region
+
+
 def trace_upstream_network(
     crossing: Point,
     stream_features: list[dict[str, Any]],
@@ -263,7 +315,12 @@ def trace_upstream_network(
                 continue
             used.add(index)
             _, endpoint, _, cum = attach
-            if side_sign(road_line, line.interpolate(0.5, normalized=True)) == -uphill_sign:
+            # 좌우(산측) 판정은 도로 근처에서만 신뢰 — 노선에서 먼 상류는 투영 기준이
+            # 뒤틀려 부호가 뒤집히므로 연결성과 도로 재교차 절단만으로 판단한다.
+            near_road = road_line.distance(endpoint) < 2.0 * STREAM_JOIN_TOL_M
+            if near_road and 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)