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 33327c66..25d20dd8 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py +++ b/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)