Aislo/B04_PreProcess/B04_PreProcess_Engine_Watershed_Analyze.py

"""배수유역 분석 오케스트레이터 (B04 — 관리자 확인용 전처리).

계획 노선(B03 업로드 파일)과 도엽 등고선·세류선만으로 배수유역을 끝까지 분석해
영구저장소에 남긴다. 30초 안팎이 걸리는 무거운 작업이라 여기서 한 번만 돌리고,
일반 사용자가 쓰는 B05는 그 결과를 읽어 쓰기만 한다(2026-07-31 사용자 지시).

  ① 도로 교차 세류망 상류측 추출 → 반경 버퍼 = 1차 배수유역
  ② 도로 시작점 기준 격자 생성 (1차 영역에 걸치는 셀만)
  ③ 등고선 하강 방향 — 높은 등고 라인에서 낮은 등고 라인으로. 보간면을 쓰지 않으므로
     가짜 웅덩이·평탄면이 원리적으로 생기지 않는다
  ④ 세류망 흐름 새김 → 최외곽부터 사슬 추적 → 도로 도달 여부(적/청) 판정
  ⑤ 최외곽 적색 셀 주변 확장 — 새로 추가한 셀에 적색이 없을 때까지
  ⑥ 도로 셀별 흐름 강도
  ⑦ 2차 전체 배수유역 외곽선
  ⑧ 기본 관 매설 위치 (도로 × 세류선 교차점)

⑨ 이후(관 최소 개수 보충, 세부유역 분할)는 사용자가 관을 옮길 수 있어야 하므로
B05에 남긴다.
"""

from __future__ import annotations

import logging
import math
import time
from dataclasses import dataclass, field
from typing import Any

import numpy as np
from shapely.geometry import LineString

from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import ContourDescent
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Expand import expand_by_red_boundary
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import (
    FlowClassification,
    RoadRaster,
    largest_ring,
    outer_boundary,
    trace_flow,
)
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
    AZIMUTH_STEPS,
    GridSpec,
    TerrainGrid,
    build_contour_cloud,
    route_elevation_floor,
)
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Stream import (
    PrimaryRegion,
    build_primary_region,
)
from common_util.common_util_route_geometry import (
    RouteVertex,
    StructureCandidate,
    find_stream_crossings,
)
from config.config_system import (
    DRAINAGE_ARROW_BLOCK_M,
    DRAINAGE_ARROW_MIN_AGREEMENT,
    DRAINAGE_ARROW_MIN_COVERAGE,
    DRAINAGE_ARROW_SPACING_M,
    DRAINAGE_GRID_SIZE_M,
    DRAINAGE_INITIAL_RADIUS_M,
    DRAINAGE_PIPE_MAX_SPACING_M,
    DRAINAGE_PIPE_MIN_SPACING_M,
)

logger = logging.getLogger(__name__)

# 강도 곡선 응답 간격(m). 도로 위 흐름 강도 표기는 이 간격으로 내보낸다.
# 1m = 계산 원본 그대로다. 5m로 줄이면 도로 색이 뭉개져 어느 자리에 물이 모이는지 못 읽는다
# (2026-08-01 사용자 지시: 단순화·평균 금지). 350m 노선이면 약 351점, 수 KB 수준이라 부담 없다.
_STRENGTH_OUTPUT_STEP_M = 1.0

# 구역마다 뽑을 유입 집중점 개수에 더하는 여유 몫.
# 개수 = floor(구역 길이 / 관 최대간격) + 이 값. 최대간격을 지키는 데 필요한 자리 수만 뽑으면
# 견줄 대상이 없어 그 자리가 최선인지 판단할 수 없다(2026-08-01 사용자 지시).
_HOTSPOT_EXTRA_PER_ZONE = 3


# ── ①~② 1차 배수유역 ────────────────────────────────────────────────────────


def resolve_primary_region(
    vertices: list[RouteVertex],
    route_line: LineString,
    contour_features: list[dict[str, Any]],
    stream_features: list[dict[str, Any]],
) -> PrimaryRegion | None:
    """도로 교차 세류선(상류측)과 노선을 반경 버퍼한 1차 배수유역과 격자 범위를 정한다.

    상·하류 판정에 쓸 등고선은 노선 주변만 있으면 된다(교차점이 전부 노선 위이므로).
    도엽 전체를 읽으면 이 단계에서만 수십 초가 날아간다.
    """
    floor = route_elevation_floor([vertex.z for vertex in vertices])
    near_bounds = route_line.buffer(DRAINAGE_INITIAL_RADIUS_M * 2.0).bounds
    cloud = build_contour_cloud(contour_features, floor, near_bounds)
    if cloud.is_empty:
        logger.warning("배수유역: 노선 주변에 등고선이 없어 1차 영역을 정할 수 없습니다.")
        return None
    return build_primary_region(
        route_line, stream_features, cloud, DRAINAGE_INITIAL_RADIUS_M, DRAINAGE_GRID_SIZE_M
    )


def preview_primary_region(
    vertices: list[RouteVertex],
    contour_features: list[dict[str, Any]],
    stream_features: list[dict[str, Any]],
) -> PrimaryRegion | None:
    """단계 검증용 — TIN·흐름 계산 없이 1차 배수유역 근거만 뽑는다."""
    if len(vertices) < 2:
        return None
    route_line = LineString([(vertex.x, vertex.y) for vertex in vertices])
    return resolve_primary_region(vertices, route_line, contour_features, stream_features)


@dataclass
class StagePreview:
    """단계 검증 산출물 묶음. 기능을 붙일 때마다 여기에 항목이 하나씩 늘어난다.

    확장을 거치면 격자와 해석 영역이 1차 영역보다 커진다. 화면·저장은 `region.spec`이
    아니라 여기 `spec`/`domain`을 봐야 한다.
    """

    region: PrimaryRegion
    spec: GridSpec | None = None
    domain: np.ndarray | None = None
    terrain: TerrainGrid | None = None
    road: RoadRaster | None = None
    flow: FlowClassification | None = None
    descent: ContourDescent | None = None
    expand_rounds: int = 0
    expand_closed: bool = False
    expand_added_cells: int = 0
    # ⑥ 도로 위 흐름 강도 — (누가거리 m, 그 구간으로 모이는 상류 면적 ㎡).
    strength_profile: list[tuple[float, float]] = field(default_factory=list)
    # ⑥-1 유입 집중점 — (누가거리 m, 유입면적 ㎡, 구역 번호, 구역 내 순위). 화면 마커용.
    inflow_hotspots: list[tuple[float, float, int, int]] = field(default_factory=list)
    # ⑦ 2차 전체 배수유역 외곽선(= 분수령). 적색 셀 전체를 폴리곤화한 것.
    basin_boundary_xy: list[tuple[float, float]] = field(default_factory=list)
    basin_area_m2: float = 0.0
    # 셀 → 도로 셀 귀속. B05가 세부유역을 나눌 때 이 배열이 있어야 한다.
    routing: Any = None
    # B05용 평균 흐름 화살표 — (x, y, 방위 라디안, 도로 도달, 셀 수).
    flow_arrows: list[tuple[float, float, float, bool, int]] = field(default_factory=list)
    # ⑧ 기본 관 매설 위치 — 도로 × 세류선 교차점.
    pipes: list[StructureCandidate] = field(default_factory=list)


def preview_stages(
    vertices: list[RouteVertex],
    contour_features: list[dict[str, Any]],
    stream_features: list[dict[str, Any]],
) -> StagePreview | None:
    """지금까지 구현·검증된 단계를 순서대로 돌려 결과를 모은다.

    현재 포함: ① 1차 배수유역 ② 격자 생성 ③ **등고선 하강 방향** ④ 도로 도달 판정
    ⑤ **최외곽 적색 셀 주변 확장**.

    ③은 보간면(TIN)을 쓰지 않는다. 높은 등고 라인에서 낮은 등고 라인으로 내려가며 방향을
    세우므로 가짜 웅덩이·평탄면이 원리적으로 생기지 않는다(2026-07-31 사용자 지시로 방식 교체).

    ⑤는 최외곽에 적색이 남아 있으면 그 주변으로 넓혀 다시 분석하고, **새로 추가한 셀에
    적색이 없으면** 멈춘다.
    """
    if len(vertices) < 2:
        return None
    from B03_FileInput.B03_FileInput_Service_Chain import _log_steps

    marks = [("시작", time.perf_counter())]
    route_line = LineString([(vertex.x, vertex.y) for vertex in vertices])
    region = resolve_primary_region(vertices, route_line, contour_features, stream_features)
    marks.append(("1차 영역(resolve_primary_region)", time.perf_counter()))
    if region is None:
        return None

    started = time.perf_counter()
    floor = route_elevation_floor([vertex.z for vertex in vertices])
    expansion = expand_by_red_boundary(
        region.spec,
        region.cell_mask,
        contour_features,
        route_line,
        region.split.upstream,
        floor,
    )
    if expansion is None:
        logger.warning("배수유역: 등고선 하강 방향을 세우지 못해 흐름 판정을 건너뜁니다.")
        return StagePreview(region=region)

    marks.append(("격자 확장·흐름 판정(expand_by_red_boundary)", time.perf_counter()))
    analysis = expansion.analysis
    spec = analysis.spec
    red = analysis.flow.reaches_road & analysis.flow.analyzed

    # ⑥ 흐름 강도 — 셀마다 물이 실제로 들어가는 도로 셀을 구해 도로 셀별로 센다.
    # 색 판정은 세류 셀에서 멈추지만(거기서 도달이 확정되므로), 강도는 그 물이 세류를 타고
    # 최종적으로 어느 도로 셀로 들어가는지를 봐야 하므로 도로만 흡수점으로 두고 다시 따라간다.
    routing = trace_flow(analysis.terrain, analysis.road) if analysis.road.count else None
    strength_curve = _preview_strength(analysis, routing, red, route_line.length)
    marks.append(("흐름 강도(trace_flow)", time.perf_counter()))

    # ⑦ 2차 전체 배수유역 외곽선 = 적색 셀 전체의 외곽.
    boundary = outer_boundary(spec, red.reshape(spec.n_rows, spec.n_cols))
    basin_ring = largest_ring(boundary) if boundary is not None else []

    # ⑧ 기본 관 매설 위치 = 도로 × 세류선 교차점(너무 가까운 것은 하나로 합침).
    pipes = _base_pipes(vertices, stream_features)

    # ⑥-1 유입 집중점 — 기본 관으로 나눈 구역마다 물이 많이 모이는 자리를 뽑는다.
    hotspots = find_inflow_hotspots(
        strength_curve, [pipe.chainage_m for pipe in pipes], route_line.length
    )

    # B05에 얹을 평균 흐름 화살표 — 셀 화살표는 도면 배율에서 안 보인다.
    flow_arrows = build_flow_arrows(analysis, analysis.flow)
    marks.append(("외곽선·기본 관·화살표", time.perf_counter()))
    _log_steps("배수유역 preview_stages", marks)

    logger.info(
        "배수유역: 단계 분석 %.1fs (확장 %d회, 최종 셀 %d개) — "
        "2차 유역 %.0f㎡, 기본 관 %d개, 강도 곡선 %d점",
        time.perf_counter() - started,
        expansion.rounds,
        spec.size,
        int(red.sum()) * spec.cell_area_m2,
        len(pipes),
        int((strength_curve > 0).sum()),
    )
    return StagePreview(
        region=region,
        spec=spec,
        domain=analysis.domain,
        terrain=analysis.terrain,
        road=analysis.road,
        flow=analysis.flow,
        descent=analysis.descent,
        expand_rounds=expansion.rounds,
        expand_closed=expansion.closed,
        expand_added_cells=expansion.added_cells,
        strength_profile=_downsample_strength(strength_curve),
        inflow_hotspots=hotspots,
        basin_boundary_xy=basin_ring,
        basin_area_m2=int(red.sum()) * spec.cell_area_m2,
        routing=routing,
        pipes=pipes,
        flow_arrows=flow_arrows,
    )


def _preview_strength(
    analysis: Any, routing: Any, red: np.ndarray, route_length_m: float
) -> np.ndarray:
    """적색 셀이 실제로 들어가는 도로 셀을 세어 누가거리별 유입 면적 곡선을 만든다."""
    road = analysis.road
    if routing is None or road.count == 0:
        return np.zeros(1)
    slots = routing.road_slot
    counted = red & (slots >= 0)
    strength = np.bincount(slots[counted], minlength=road.count).astype(np.float64)
    return _strength_by_chainage(
        road.chainage, strength * analysis.spec.cell_area_m2, route_length_m
    )


def build_flow_arrows(analysis: Any, flow: Any) -> list[tuple[float, float, float, bool, int]]:
    """셀 흐름을 블록 단위로 평균해 B05에 얹을 화살표를 뽑는다.

    셀 화살표는 1m라 도면 배율에서 경향이 안 보인다. 겹치지 않는 블록으로 나눠 방향을
    평균하고, 화살표끼리 최소 간격을 두어 솎아낸다(2026-07-31 사용자 지시).

    **세류선 셀과 도로 셀은 뺀다.** 그 자리 흐름은 지형 경사가 아니라 확정된 물길·노면을
    따르는 값이라 사면 경향을 왜곡한다.

    방향 평균은 산술평균이 아니라 **원형 평균**으로 낸다(0°와 359°의 평균은 180°가 아니라
    0°다). 평균 벡터 길이가 일치도이므로, 블록 안 방향이 제각각이면 그 블록은 버린다.
    """
    spec = analysis.spec
    rows, cols = spec.n_rows, spec.n_cols
    direction = flow.direction.reshape(rows, cols)
    usable = (
        analysis.domain
        & flow.analyzed.reshape(rows, cols)
        & (direction < AZIMUTH_STEPS)  # 싱크·무효 제외
        & ~analysis.road.mask
    )
    if flow.burned is not None:
        usable &= ~flow.burned.reshape(rows, cols)
    if not usable.any():
        return []

    block = max(1, int(round(DRAINAGE_ARROW_BLOCK_M / spec.cell_m)))
    stride = max(1, int(round(DRAINAGE_ARROW_SPACING_M / (block * spec.cell_m))))
    angle = direction.astype(np.float64) * (2.0 * math.pi / AZIMUTH_STEPS)
    reaches = flow.reaches_road.reshape(rows, cols)

    arrows: list[tuple[float, float, float, bool, int]] = []
    for row0 in range(0, rows - block + 1, block * stride):
        for col0 in range(0, cols - block + 1, block * stride):
            window = usable[row0 : row0 + block, col0 : col0 + block]
            count = int(window.sum())
            if count < DRAINAGE_ARROW_MIN_COVERAGE * block * block:
                continue
            local = angle[row0 : row0 + block, col0 : col0 + block][window]
            mean_x = float(np.cos(local).mean())
            mean_y = float(np.sin(local).mean())
            agreement = math.hypot(mean_x, mean_y)
            if agreement < DRAINAGE_ARROW_MIN_AGREEMENT:
                continue  # 방향이 제각각인 블록 — 평균이 경향을 대표하지 못한다
            centre_row = row0 + block / 2.0
            centre_col = col0 + block / 2.0
            arrows.append(
                (
                    spec.x_min + centre_col * spec.cell_m,
                    spec.y_max - centre_row * spec.cell_m,
                    math.atan2(mean_y, mean_x),
                    bool(reaches[row0 : row0 + block, col0 : col0 + block][window].mean() >= 0.5),
                    count,
                )
            )
    logger.info(
        "배수유역: 평균 흐름 화살표 %d개 (블록 %.0fm, 간격 %.0fm, 세류·도로 셀 제외)",
        len(arrows),
        block * spec.cell_m,
        block * stride * spec.cell_m,
    )
    return arrows


def _base_pipes(
    vertices: list[RouteVertex], stream_features: list[dict[str, Any]]
) -> list[StructureCandidate]:
    """도로 × 세류선 교차점을 기본 관 위치로 삼는다. 300m 보충 배치는 다음 단계다."""
    pipes: list[StructureCandidate] = []
    for candidate in find_stream_crossings(vertices, stream_features):
        if pipes and candidate.chainage_m - pipes[-1].chainage_m < DRAINAGE_PIPE_MIN_SPACING_M:
            continue
        pipes.append(candidate)
    return pipes


# ── 흐름 강도 곡선 ──────────────────────────────────────────────────────────


def _strength_by_chainage(
    road_chainage: np.ndarray, strength_area: np.ndarray, total_length: float
) -> np.ndarray:
    """도로 셀 강도를 1m 누가거리 구간으로 합산한 곡선(㎡/m 구간 합).

    구간 k는 **[k, k+1)** 이다. `round`를 쓰면 안 된다 — 도로 셀 누가거리는 노선을 0.5m
    간격으로 샘플해 붙인 값이라 전부 0.5의 배수인데, numpy의 `round`는 정확히 .5인 값을
    짝수 쪽으로 보낸다(0.5→0, 1.5→2). 그러면 짝수 칸에는 세 샘플({k−0.5, k, k+0.5}),
    홀수 칸에는 한 샘플만 들어가 강도가 2m 주기로 흔들리고, 도로가 격자축과 나란한 구간은
    홀수 칸이 전부 0이 되어 화면에서 1m 간격 점선으로 보인다(2026-08-01 검증에서 발견).
    """
    bins = max(1, int(np.ceil(total_length)) + 1)
    if road_chainage.size == 0:
        return np.zeros(bins)
    index = np.clip(np.floor(road_chainage).astype(np.int64), 0, bins - 1)
    return np.bincount(index, weights=strength_area, minlength=bins)


def find_inflow_hotspots(
    curve: np.ndarray,
    pipe_chainages: list[float],
    route_length_m: float,
) -> list[tuple[float, float, int, int]]:
    """도로 위 **유입 집중점**을 뽑는다 — (누가거리 m, 유입면적 ㎡, 구역 번호, 구역 내 순위).

    구역 번호 `-1`은 **기본 관(세류 교차점) 자리**를 뜻한다. 나머지는 그 구역에서 보충 후보로
    뽑은 자리다.

    왜 필요한가(2026-08-01 사용자 지시): 기본 관은 국가 수치지형도의 세류선이 도로를 가로지르는
    자리라 근거가 확실하다. 그런데 관 **최대 간격**을 지키려면 세류가 없는 긴 구간에도 관을
    넣어야 하는데, 이때 등간격으로 기계적으로 꽂는 대신 **물이 실제로 많이 모이는 자리**를
    골라야 한다. 여기서 뽑은 지점이 그 후보이며, B05가 관을 옮겨도 세부유역을 다시 나누는
    기준으로 쓸 수 있다.

    규칙
      · 구역 = 시점 → 기본 관들 → 종점 으로 자른 구간
      · 구역별 개수 = `floor(구역 길이 / 관 최대간격) + 3`
        (관 최대간격을 지키는 데 꼭 필요한 자리 수에 검토 여유 몫을 더한 것이다.
         한 자리만 보여 주면 그 자리가 실제로 최선인지 견줄 대상이 없다 — 2026-08-01 사용자 지시)
      · 고르는 법 = 강도가 가장 큰 1m 칸을 고르고 **그 지점 양 편측으로 관 최소간격만큼**을
        후보에서 뺀 뒤 다음으로 큰 칸을 고른다. 빼지 않으면 같은 계곡의 이웃 칸들이 연달아
        뽑혀 한쪽에 쏠린다.
      · **기본 관 쪽 경계에서도 관 최소간격만큼 물러나 본다.** 강도 최대점은 원리상 세류가
        도로를 가로지르는 자리, 곧 기본 관 자신이다. 관 옆 1m만 비워 두면 그 관이 받는 물
        전부(전체 유역의 대부분)가 관에서 두어 걸음 떨어진 칸에서 다시 최대값으로 잡혀,
        정작 관이 없어 보충이 필요한 자리는 뽑히지 못한다. 노선 시·종점은 관이 아니므로
        물러나지 않는다.
      · **기본 관 자리 자체도 집중점으로 넣는다(구역 번호 −1).** 세류가 도로를 가로지르는
        자리가 대개 그 노선에서 가장 큰 유역을 받는데, 그 자리는 구역 경계라 위 규칙으로는
        영영 뽑히지 않는다(2026-08-01 사용자 지시). 위치는 관 누가거리 그대로가 아니라
        관 ±최소간격 안에서 강도가 가장 큰 칸으로 잡는다 — 도로를 폭으로 굽는 과정에서
        실제 유입 봉우리가 관에서 두어 걸음 옆에 놓이기 때문이다.
      · 강도 0인 칸은 뽑지 않는다(개수를 못 채워도 그대로 둔다).
    """
    if curve.size == 0 or route_length_m <= 0:
        return []
    # 구역 경계 = 0, 기본 관들, 노선 끝. 순서·중복을 정리해 둔다.
    edges = sorted({0.0, float(route_length_m), *(float(p) for p in pipe_chainages)})
    exclusion = max(1, int(round(DRAINAGE_PIPE_MIN_SPACING_M)))
    hotspots: list[tuple[float, float, int, int]] = []

    # ① 기본 관 자리 — 관마다 그 근처 봉우리 하나.
    for pipe_m in sorted({float(p) for p in pipe_chainages}):
        center = int(round(pipe_m))
        low = max(0, center - exclusion)
        high = min(curve.size - 1, center + exclusion)
        if high < low:
            continue
        window = curve[low : high + 1]
        best = int(np.argmax(window))
        if window[best] <= 0:
            continue
        hotspots.append((float(low + best), float(window[best]), -1, 0))
    # 관 자리와 겹치는 칸은 구역 쪽에서 다시 뽑지 않는다 — 같은 자리에 마커가 두 개 겹친다.
    taken = {int(chainage) for chainage, _, _, _ in hotspots}
    # ② 구역별 보충 후보 — 관과 관 사이에서 물이 많이 모이는 자리.
    for zone_index, (start_m, end_m) in enumerate(zip(edges, edges[1:])):
        span = end_m - start_m
        if span <= 0:
            continue
        wanted = int(span // DRAINAGE_PIPE_MAX_SPACING_M) + _HOTSPOT_EXTRA_PER_ZONE
        # 이 구역에서 볼 곡선 조각. 관 쪽 경계는 최소간격만큼, 노선 끝 쪽은 한 칸만 물린다.
        start_is_pipe = zone_index > 0
        end_is_pipe = zone_index < len(edges) - 2
        low = int(math.floor(start_m)) + (exclusion if start_is_pipe else 1)
        high = int(math.ceil(end_m)) - (exclusion if end_is_pipe else 1)
        high = min(high, curve.size - 1)
        if high < low:
            continue
        window = curve[low : high + 1].copy()
        for index in taken:
            if low <= index <= high:
                window[index - low] = 0.0
        for rank in range(wanted):
            best = int(np.argmax(window))
            if window[best] <= 0:
                break  # 남은 칸이 전부 0 — 뽑을 것이 없다
            hotspots.append((float(low + best), float(window[best]), zone_index, rank))
            # 고른 지점 양 편측으로 최소 간격만큼 지운다(한쪽 쏠림 방지).
            window[max(0, best - exclusion) : best + exclusion + 1] = 0.0
    hotspots.sort(key=lambda item: item[0])
    return hotspots


def _downsample_strength(curve: np.ndarray) -> list[tuple[float, float]]:
    """응답용으로 강도 곡선을 일정 간격으로 줄인다(구간 합 유지).

    끝자락을 잘라내면 종점 부근 유입 면적이 통째로 사라지므로 0으로 채워 맞춘다.
    """
    step = max(1, int(_STRENGTH_OUTPUT_STEP_M))
    if curve.size == 0:
        return []
    padding = (-curve.size) % step
    padded = np.append(curve, np.zeros(padding)) if padding else curve
    summed = padded.reshape(-1, step).sum(axis=1)
    return [
        (float(position * step), float(value)) for position, value in enumerate(summed) if value > 0
    ]