Aislo/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Subdivide.py

"""메인 배수유역 내부 세분화 — 분할선으로 폴리곤을 쪼갠다 (2026-07-30).

불변 조건(사용자 지시): 전체(메인) 배수유역 경계는 절대 변경하지 않는다.
세분화는 확정된 메인 유역 폴리곤을 관 사이 분할선(물갈림 고개에서 오르는
능선 근사선)으로 **내부에서만** 쪼개는 방식이다 — 외곽 재추적 금지.
따라서 세부유역은 서로 배타적이고 합집합은 항상 메인 유역과 동일하다.
"""

from __future__ import annotations

import logging
import math

from shapely.geometry import LineString, Point, Polygon
from shapely.ops import split as shapely_split
from shapely.ops import substring, unary_union

from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Trace import DividerStep

logger = logging.getLogger(__name__)

# 분할 절단선 연장: 도로 하류측(m)과 능선 너머(폴리곤 대각선 배수).
CUT_ROAD_TAIL_M = 40.0
CUT_RIDGE_EXTEND_RATIO = 1.5


def _cut_line(
    road_point: Point,
    steps: list[DividerStep],
    road_line: LineString,
    main_polygon: Polygon,
) -> LineString | None:
    """분할선 스텝을 절단선으로 확장한다 — 도로 하류측과 능선 너머까지 관통.

    split()은 절단선이 폴리곤 경계를 완전히 넘어야 동작하므로 양끝을 연장한다.
    스텝이 없으면(등고선 공백) 도로 법선 직선으로 폴백한다.
    """
    bounds = main_polygon.bounds
    reach = CUT_RIDGE_EXTEND_RATIO * math.hypot(bounds[2] - bounds[0], bounds[3] - bounds[1])
    points = [(road_point.x, road_point.y)]
    points.extend((step.point.x, step.point.y) for step in steps if step.geom_index >= 0)
    if len(points) < 2:
        # 폴백: 도로 접선의 법선 방향으로 폴리곤을 관통하는 직선.
        t = road_line.project(road_point)
        a = road_line.interpolate(max(0.0, t - 5.0))
        b = road_line.interpolate(min(road_line.length, t + 5.0))
        dx, dy = b.x - a.x, b.y - a.y
        norm = math.hypot(dx, dy)
        if norm < 1e-6:
            return None
        nx, ny = -dy / norm, dx / norm
        head = (road_point.x + nx * reach, road_point.y + ny * reach)
        tail = (road_point.x - nx * reach, road_point.y - ny * reach)
        return LineString([tail, (road_point.x, road_point.y), head])
    # 능선 너머 연장: 마지막 진행 방향 유지.
    (px, py), (qx, qy) = points[-2], points[-1]
    dx, dy = qx - px, qy - py
    norm = math.hypot(dx, dy)
    if norm >= 1e-6:
        points.append((qx + dx / norm * reach, qy + dy / norm * reach))
    # 도로 하류측 연장: 첫 스텝 → 도로점 방향을 그대로 지나쳐 내려간다.
    (fx, fy) = points[1]
    dx, dy = road_point.x - fx, road_point.y - fy
    norm = math.hypot(dx, dy)
    if norm >= 1e-6:
        points.insert(
            0,
            (
                road_point.x + dx / norm * CUT_ROAD_TAIL_M,
                road_point.y + dy / norm * CUT_ROAD_TAIL_M,
            ),
        )
    return LineString(points)


def _interval_index(piece: Polygon, road_line: LineString, divide_ts: list[float]) -> int:
    """조각이 어느 관 구간(k)에 속하는지 — 구간 도로와 맞닿는 길이가 최대인 곳.

    대표점 투영은 대형 계곡 조각에서 오판한다(상류로 길게 뻗은 조각의 대표점이
    엉뚱한 구간에 떨어짐). 도로 접촉이 전혀 없는 조각만 대표점 투영으로 폴백.
    """
    strip = piece.buffer(1.0)
    best_k, best_length = -1, 0.0
    for k in range(len(divide_ts) - 1):
        segment = substring(road_line, divide_ts[k], divide_ts[k + 1])
        if segment.is_empty:
            continue
        length = segment.intersection(strip).length
        if length > best_length:
            best_k, best_length = k, length
    if best_k >= 0:
        return best_k
    t = road_line.project(piece.representative_point())
    for k in range(len(divide_ts) - 1):
        if divide_ts[k] <= t <= divide_ts[k + 1]:
            return k
    return 0 if t < divide_ts[0] else len(divide_ts) - 2


def subdivide_main_polygon(
    main_polygon: Polygon,
    divide_points: list[Point],
    dividers: list[list[DividerStep]],
    road_line: LineString,
) -> list[Polygon | None]:
    """메인 유역 폴리곤을 내부 분할선으로 쪼개 관 구간별 조각을 돌려준다.

    반환 길이 = 관 개수(구간 수). 조각이 없는 구간은 None.
    split() 기반이므로 조각들은 배타적이고 합집합 == 메인 폴리곤이 보장된다.
    구간에 여러 조각이 잡히면(절단선 재진입) 모두 합쳐 가장 큰 폴리곤을 쓴다.
    """
    interval_count = len(divide_points) - 1
    if interval_count <= 1:
        return [main_polygon]
    pieces: list[Polygon] = [main_polygon]
    for position in range(1, interval_count):
        cut = _cut_line(divide_points[position], dividers[position], road_line, main_polygon)
        if cut is None:
            logger.warning("분할 절단선 생성 실패(구간 %d) — 해당 분할을 건너뜁니다.", position)
            continue
        next_pieces: list[Polygon] = []
        for piece in pieces:
            try:
                parts = shapely_split(piece, cut)
            except Exception:  # noqa: BLE001 - 절단 실패 시 조각 유지
                next_pieces.append(piece)
                continue
            split_parts = [
                part
                for part in getattr(parts, "geoms", [parts])
                if part.geom_type == "Polygon" and not part.is_empty
            ]
            next_pieces.extend(split_parts if split_parts else [piece])
        pieces = next_pieces
    divide_ts = [road_line.project(point) for point in divide_points]
    assigned: list[list[Polygon]] = [[] for _ in range(interval_count)]
    for piece in pieces:
        assigned[_interval_index(piece, road_line, divide_ts)].append(piece)
    # 구간별 대표 조각 = 최대 조각과 그에 붙는 조각들. 비연결 잔여 조각은 버리지
    # 않고(합집합 불변 조건) 맞닿는 인접 구간으로 재배정한다.
    result: list[Polygon | None] = []
    leftovers: list[Polygon] = []
    for group in assigned:
        merged = _merge_touching(group)
        result.append(merged[0] if merged else None)
        leftovers.extend(merged[1:])
    for extra in leftovers:
        for position in sorted(
            range(interval_count),
            key=lambda k: extra.distance(result[k]) if result[k] is not None else math.inf,
        ):
            base = result[position]
            if base is None or not extra.touches(base):
                continue
            candidate = base.union(extra).buffer(0)
            if candidate.geom_type == "Polygon":
                result[position] = candidate
                break
        else:
            logger.warning("세분화 잔여 조각(%.0f m²)을 재배정하지 못해 제외합니다.", extra.area)
    return result


def _merge_touching(group: list[Polygon]) -> list[Polygon]:
    """조각 묶음을 서로 맞닿는 것끼리 합쳐 면적 내림차순으로 돌려준다."""
    if not group:
        return []
    merged = unary_union(group).buffer(0)
    parts = list(merged.geoms) if merged.geom_type == "MultiPolygon" else [merged]
    parts = [part for part in parts if part.geom_type == "Polygon" and not part.is_empty]
    return sorted(parts, key=lambda part: part.area, reverse=True)