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 25d20dd8..5593952a 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 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 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 2eb965ee..25110185 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py +++ b/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)