Aislo/common_util/common_util_drainage_detail.py

"""배수유역 세부 설계 공용 엔진 (B04 관리자 화면 · B05 사용자 화면 공용).

**격자 해석은 하지 않는다.** B04 배수유역 분석이 미리 돌려 저장한 결과를 읽어, 사용자가
실제로 손대는 두 가지만 처리한다(2026-07-31 사용자 지시).

  ⑨ 관 간격이 최대치를 넘는 구간에 **최소 개수**로 관을 보충
  ⑩ 측구 흐름으로 도로 셀 → 담당 관을 정하고, 셀이 도달한 도로 셀의 담당 관을 그대로
     그 셀의 유역 번호로 삼아 세부유역을 나눈다
  ⑪ 사용자가 관을 옮기거나 추가하면 ⑩만 다시 돈다 — 격자 해석은 재사용한다

읽어 오는 것(`{배수유역 폴더}/`):
  · `03_road_routing.geojson` — 계획도로선 · 기본 배관 · 2차 전체 배수유역
  · `03_road_routing.npz` — 셀 → 도로 셀 귀속, 유하장, 강도, 도로 셀 제원, 셀 표고
화살표(방향 코드)나 밴드 표고 같은 관리자 확인용 배열은 읽지 않는다 — 여기서는 필요 없고
파일만 무거워진다.

**왜 common_util인가**: B04(관리자 트러블슈팅)와 B05(일반 사용자)가 같은 이름의 버튼을
누르면 같은 결과가 나와야 한다(2026-08-01 사용자 지시). 두 벌로 두면 언젠가 갈라진다.
읽는 폴더만 다르므로 폴더를 인자로 받고, B04/B05 각자의 어댑터가 경로를 정한다.
격자 산출물의 규격(`GridSpec`·`STAGES`·`polygonize_labels`)은 B04가 만든 것이므로
정의처인 B04 엔진을 그대로 참조한다(역방향 참조 없음).
"""

from __future__ import annotations

import json
import logging
from dataclasses import dataclass, field
from pathlib import Path
from typing import Any

import numpy as np

from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Export import STAGES
from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Flow import largest_ring, polygonize_labels
from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Grid import GridSpec
from common_util.common_util_route_geometry import (
    RouteVertex,
    StructureCandidate,
    interpolate_vertex,
    is_uphill_at,
)
from config.config_system import (
    DRAINAGE_DITCH_SAMPLE_M,
    DRAINAGE_PIPE_MAX_SPACING_M,
    DRAINAGE_PIPE_MIN_SPACING_M,
)

logger = logging.getLogger(__name__)

# 관 위치 선정 점수 배분 — 흐름 강도가 주, 종단 저점이 보조.
_SCORE_WEIGHT_STRENGTH = 0.7
_SCORE_WEIGHT_SAG = 0.3
# 성토부(내리막)는 물이 노선 밖으로 빠지므로 관 위치로 덜 선호한다.
_SCORE_FILL_PENALTY = 0.5


@dataclass
class DrainageDetail:
    """세부 설계 산출물 — 화면에 그릴 기하와 세부유역."""

    route_lonlat: list[list[float]] = field(default_factory=list)
    basin_lonlat: list[list[float]] = field(default_factory=list)
    pipes: list[StructureCandidate] = field(default_factory=list)
    basins: list[WatershedBasin] = field(default_factory=list)
    grid_cell_m: float = 1.0
    # B04가 계산해 둔 평균 흐름 화살표를 그대로 넘긴다 — 여기서 다시 계산하지 않는다.
    flow_arrows: list[list[Any]] = field(default_factory=list)
    arrow_spacing_m: float = 0.0
    # 유역 안쪽 상류 세류망(WGS84 lon/lat 조각들). 화면 강조 표시용 — 계산에는 쓰지 않는다.
    upstream_lonlat: list[list[list[float]]] = field(default_factory=list)


@dataclass
class WatershedBasin:
    """관 하나가 받는 세부 배수유역."""

    index: int
    chainage_m: float
    outlet_x: float
    outlet_y: float
    boundary_xy: list[tuple[float, float]] = field(default_factory=list)
    area_m2: float = 0.0
    relief_m: float = 0.0
    flow_length_m: float = 0.0
    pipe_diameter_mm: float | None = None


@dataclass
class RoadRouting:
    """B04가 남긴 배수유역 분석 결과 — 세부유역을 나누는 데 필요한 최소 묶음."""

    spec: GridSpec
    # (R*C,) int32 — 셀이 물길을 따라 도달하는 도로 셀 슬롯(−1 = 미도달).
    road_slot: np.ndarray
    path_length: np.ndarray  # (R*C,) float32 — 그 도로 셀까지 물길 길이(m)
    elevation: np.ndarray  # (R*C,) float32 — 셀 표고(유역 낙차 계산용)
    road_cell_index: np.ndarray  # (K,) int32 — 도로 셀의 평탄 인덱스
    road_chainage: np.ndarray  # (K,) float64 — 도로 셀의 누가거리(m)
    strength: np.ndarray  # (K,) int64 — 도로 셀별 상류 셀 수
    # 화면에 그대로 그릴 기하(WGS84 lon/lat).
    route_lonlat: list[list[float]] = field(default_factory=list)
    basin_lonlat: list[list[float]] = field(default_factory=list)
    base_pipes: list[StructureCandidate] = field(default_factory=list)
    # 평균 흐름 화살표 — [x, y, 방위(도), 도로도달, 셀 수]. B04가 계산해 둔 그대로.
    flow_arrows: list[list[Any]] = field(default_factory=list)
    arrow_spacing_m: float = 0.0

    @property
    def strength_curve(self) -> np.ndarray:
        """누가거리 1m 구간별 유입 면적(㎡) 곡선 — 관 보충 위치 점수의 근거."""
        if self.road_chainage.size == 0:
            return np.zeros(1)
        bins = max(1, int(np.ceil(self.road_chainage.max())) + 1)
        index = np.clip(np.round(self.road_chainage).astype(np.int64), 0, bins - 1)
        weights = self.strength.astype(np.float64) * self.spec.cell_area_m2
        return np.bincount(index, weights=weights, minlength=bins)


def build_detail(
    directory: Path,
    vertices: list[RouteVertex],
    confirmed_chainages: list[float] | None = None,
) -> DrainageDetail | None:
    """B04 분석 결과를 읽어 관을 보충하고 세부유역을 나눈다.

    `confirmed_chainages`를 주면 그 위치를 관으로 확정하고(사용자 편집), 비우면 B04의
    기본 관에 최대 간격 규칙으로 최소 개수만 보충한다. 어느 쪽이든 격자 해석은 하지 않는다.
    """
    routing = read_road_routing(directory)
    if routing is None or len(vertices) < 2:
        return None

    if confirmed_chainages:
        pipes = pipes_from_chainages(vertices, confirmed_chainages)
    else:
        # 저장분의 기본 관은 누가거리만 신뢰한다 — 좌표는 현재 노선 위로 다시 찍는다.
        base = [
            StructureCandidate(
                chainage_m=pipe.chainage_m,
                x=interpolate_vertex(vertices, pipe.chainage_m)[0],
                y=interpolate_vertex(vertices, pipe.chainage_m)[1],
                reason=pipe.reason,
            )
            for pipe in routing.base_pipes
        ]
        pipes = place_pipes(vertices, base, routing.strength_curve)

    detail = DrainageDetail(
        route_lonlat=routing.route_lonlat,
        basin_lonlat=routing.basin_lonlat,
        pipes=pipes,
        grid_cell_m=routing.spec.cell_m,
        flow_arrows=routing.flow_arrows,
        arrow_spacing_m=routing.arrow_spacing_m,
        upstream_lonlat=read_upstream_lines(directory),
    )
    if not pipes:
        return detail
    pipe_of_slot = assign_road_cells_to_pipes(vertices, pipes, routing.road_chainage)
    detail.basins = assemble_basins(routing, pipes, pipe_of_slot)
    logger.info(
        "배수유역: 세부 설계 — 관 %d개(기본 %d + 보충 %d), 세부유역 %d개",
        len(pipes),
        sum(1 for pipe in pipes if pipe.reason != "spacing"),
        sum(1 for pipe in pipes if pipe.reason == "spacing"),
        len(detail.basins),
    )
    return detail


def read_road_routing(directory: Path) -> RoadRouting | None:
    """`03_road_routing` 산출물을 읽는다. 없으면 None."""
    prefix = STAGES["road_routing"]
    array_path = directory / f"{prefix}_road_routing.npz"
    if not array_path.exists():
        logger.warning("배수유역: B04 분석 결과가 없습니다 (%s).", array_path)
        return None
    try:
        with np.load(array_path, allow_pickle=False) as data:
            spec = GridSpec(
                x_min=float(data["x_min"]),
                y_max=float(data["y_max"]),
                cell_m=float(data["cell_m"]),
                n_rows=int(data["n_rows"]),
                n_cols=int(data["n_cols"]),
            )
            routing = RoadRouting(
                spec=spec,
                road_slot=data["road_slot"].reshape(-1),
                path_length=data["path_length"].reshape(-1),
                elevation=data["elevation"].reshape(-1),
                road_cell_index=data["road_cell_index"],
                road_chainage=data["road_chainage"],
                strength=data["strength"],
            )
    except (OSError, KeyError, ValueError):
        logger.warning("배수유역: B04 분석 결과를 읽지 못했습니다 (%s).", array_path)
        return None

    _read_geometry(directory / f"{prefix}_road_routing.geojson", routing)
    logger.info(
        "배수유역: B04 결과 로드 — 격자 %d×%d, 도로 셀 %d, 기본 관 %d",
        spec.n_rows,
        spec.n_cols,
        routing.road_cell_index.size,
        len(routing.base_pipes),
    )
    return routing


def _read_geometry(path: Path, routing: RoadRouting) -> None:
    """계획도로선·2차 유역 외곽선·기본 관을 GeoJSON에서 읽어 채운다."""
    if not path.exists():
        logger.warning("배수유역: B04 기하 산출물이 없습니다 (%s).", path)
        return
    try:
        with path.open("r", encoding="utf-8") as file:
            document = json.load(file)
    except (OSError, json.JSONDecodeError):
        logger.warning("배수유역: B04 기하 산출물을 읽지 못했습니다 (%s).", path)
        return
    routing.arrow_spacing_m = float(
        (document.get("properties") or {}).get("arrow_spacing_m") or 0.0
    )
    for feature in document.get("features", []):
        properties = feature.get("properties") or {}
        geometry = feature.get("geometry") or {}
        coordinates = geometry.get("coordinates")
        kind = properties.get("kind")
        if kind == "route" and geometry.get("type") == "LineString":
            routing.route_lonlat = coordinates
        elif kind == "basin_boundary" and geometry.get("type") == "Polygon" and coordinates:
            routing.basin_lonlat = coordinates[0]
        elif kind == "flow_arrow" and geometry.get("type") == "Point":
            # 화면이 미터로 그리므로 속성의 x·y를 쓴다(기하는 저장 규약상 lon/lat).
            routing.flow_arrows.append(
                [
                    float(properties.get("x") or 0.0),
                    float(properties.get("y") or 0.0),
                    float(properties.get("azimuth_deg") or 0.0),
                    bool(properties.get("reaches_road")),
                    int(properties.get("cells") or 0),
                ]
            )
        elif kind == "pipe" and geometry.get("type") == "Point":
            routing.base_pipes.append(
                StructureCandidate(
                    chainage_m=float(properties.get("chainage_m") or 0.0),
                    x=0.0,
                    y=0.0,
                    reason=str(properties.get("reason") or "stream"),
                )
            )


def read_upstream_lines(directory: Path) -> list[list[list[float]]]:
    """`01_primary_region`에서 상류 세류망만 읽는다(화면 강조용).

    유역 판정의 기준선이라 B04 오버레이에서도 같은 선을 굵게 그린다 — B05는 그 선을
    그대로 받아 표시만 한다.
    """
    path = directory / f"{STAGES['primary_region']}_primary_region.geojson"
    if not path.exists():
        return []
    try:
        with path.open("r", encoding="utf-8") as file:
            document = json.load(file)
    except (OSError, json.JSONDecodeError):
        logger.warning("배수유역: 상류 세류망을 읽지 못했습니다 (%s).", path)
        return []
    lines: list[list[list[float]]] = []
    for feature in document.get("features", []):
        properties = feature.get("properties") or {}
        geometry = feature.get("geometry") or {}
        if properties.get("kind") != "upstream":
            continue
        coordinates = geometry.get("coordinates")
        if geometry.get("type") == "LineString" and coordinates:
            lines.append(coordinates)
        elif geometry.get("type") == "MultiLineString" and coordinates:
            lines.extend(part for part in coordinates if part)
    return lines


# ── ⑨ 관 최소 개수 보충 ─────────────────────────────────────────────────────


def place_pipes(
    vertices: list[RouteVertex],
    base_pipes: list[StructureCandidate],
    strength_curve: np.ndarray,
) -> list[StructureCandidate]:
    """B04가 정한 기본 관(세류 교차점)에, 최대 간격을 넘는 구간만 최소 개수로 보충한다.

    기본 관은 여기서 다시 찾지 않는다 — B04 산출물에 이미 들어 있다.
    """
    total_length = vertices[-1].chainage_m
    base: list[StructureCandidate] = []
    for candidate in sorted(base_pipes, key=lambda item: item.chainage_m):
        if base and candidate.chainage_m - base[-1].chainage_m < DRAINAGE_PIPE_MIN_SPACING_M:
            continue
        base.append(candidate)

    filled: list[StructureCandidate] = []
    previous = 0.0
    for candidate in [*base, None]:
        boundary = candidate.chainage_m if candidate else total_length
        filled.extend(_fill_gap(vertices, strength_curve, previous, boundary))
        if candidate:
            filled.append(candidate)
            previous = candidate.chainage_m
        else:
            previous = boundary
    filled.sort(key=lambda item: item.chainage_m)
    return filled


def _fill_gap(
    vertices: list[RouteVertex],
    strength_curve: np.ndarray,
    start_m: float,
    end_m: float,
) -> list[StructureCandidate]:
    """[start, end] 구간에 최대 간격을 지키는 **최소 개수**의 관을 배치한다.

    필요 개수 n은 구간 길이로 정해지고(ceil(L/max) − 1), 각 관은 등분 위치를 중심으로
    허용 여유(slack) 안에서만 움직인다. 그래서 개수는 늘지 않으면서도 흐름 강도가 크고
    종단이 낮은 지점으로 붙는다.
    """
    span = end_m - start_m
    if span <= DRAINAGE_PIPE_MAX_SPACING_M:
        return []
    count = int(np.ceil(span / DRAINAGE_PIPE_MAX_SPACING_M)) - 1
    if count <= 0:
        return []
    spacing = span / (count + 1)
    slack = max(0.0, (DRAINAGE_PIPE_MAX_SPACING_M - spacing) / 2.0)
    placed: list[StructureCandidate] = []
    for order in range(1, count + 1):
        nominal = start_m + spacing * order
        low = max(start_m + DRAINAGE_PIPE_MIN_SPACING_M, nominal - slack)
        high = min(end_m - DRAINAGE_PIPE_MIN_SPACING_M, nominal + slack)
        chosen = _best_position(vertices, strength_curve, low, high, nominal)
        x, y, _ = interpolate_vertex(vertices, chosen)
        placed.append(StructureCandidate(chainage_m=chosen, x=x, y=y, reason="spacing"))
    return placed


def _best_position(
    vertices: list[RouteVertex],
    strength_curve: np.ndarray,
    low_m: float,
    high_m: float,
    fallback_m: float,
) -> float:
    """허용 구간 안에서 흐름 강도가 크고 종단이 낮은 위치를 고른다."""
    if high_m <= low_m:
        return fallback_m
    positions = np.arange(low_m, high_m + 1.0, 1.0)
    if positions.size == 0:
        return fallback_m
    index = np.clip(np.round(positions).astype(np.int64), 0, strength_curve.size - 1)
    strength = strength_curve[index]
    heights = np.array([interpolate_vertex(vertices, float(p))[2] for p in positions])

    strength_score = strength / strength.max() if strength.max() > 0 else np.zeros_like(strength)
    height_span = float(heights.max() - heights.min())
    sag_score = (
        (heights.max() - heights) / height_span if height_span > 1e-6 else np.zeros_like(heights)
    )
    score = _SCORE_WEIGHT_STRENGTH * strength_score + _SCORE_WEIGHT_SAG * sag_score
    for order, position in enumerate(positions):
        if not is_uphill_at(vertices, float(position)):
            score[order] *= _SCORE_FILL_PENALTY
    return float(positions[int(np.argmax(score))])


def pipes_from_chainages(
    vertices: list[RouteVertex], chainages: list[float]
) -> list[StructureCandidate]:
    """사용자가 확정·편집한 누가거리 목록을 관 후보로 되돌린다.

    노선 밖 값은 시·종점으로 당긴다. 그대로 두면 마커는 끝점에 찍히는데 라벨만 −50m처럼
    나와 좌표와 표기가 어긋난다.
    """
    total_length = vertices[-1].chainage_m
    clamped = {min(max(round(float(item), 2), 0.0), total_length) for item in chainages}
    pipes: list[StructureCandidate] = []
    for value in sorted(clamped):
        x, y, _ = interpolate_vertex(vertices, value)
        pipes.append(StructureCandidate(chainage_m=value, x=x, y=y, reason="confirmed"))
    return pipes


# ── ⑩ 측구 흐름으로 도로 셀 → 담당 관 ───────────────────────────────────────


def assign_road_cells_to_pipes(
    vertices: list[RouteVertex],
    pipes: list[StructureCandidate],
    road_chainage: np.ndarray,
) -> np.ndarray:
    """도로 셀마다 물이 실제로 흘러가는 담당 관 번호를 정한다.

    노면 물은 측구를 타고 종단 내리막으로 흐르므로, 종단 계획선을 1차원 지형으로 보고
    같은 방식(내리막 추적 + 관에서 흡수)으로 푼다. 관이 없는 사그(저점)에 갇힌 구간은
    가장 가까운 관이 받는 것으로 본다.
    """
    total_length = vertices[-1].chainage_m
    step = max(DRAINAGE_DITCH_SAMPLE_M, 0.5)
    stations = np.arange(0.0, total_length + step, step)
    heights = np.array([interpolate_vertex(vertices, float(s))[2] for s in stations])
    pipe_chainages = np.array([pipe.chainage_m for pipe in pipes])
    pipe_station = np.clip(np.round(pipe_chainages / step).astype(np.int64), 0, stations.size - 1)

    # 앞뒤 이웃 중 더 낮은 쪽으로 흘려보낸다(양쪽 다 높으면 사그 = 제자리).
    back_z = np.full(stations.size, np.inf)
    back_z[1:] = heights[:-1]
    forward_z = np.full(stations.size, np.inf)
    forward_z[:-1] = heights[1:]
    go_back = (back_z < heights) & (back_z <= forward_z)
    go_forward = (forward_z < heights) & ~go_back
    receiver = np.arange(stations.size, dtype=np.int64)
    receiver[go_back] -= 1
    receiver[go_forward] += 1
    receiver[pipe_station] = pipe_station  # 관은 물을 흡수한다

    owner = np.full(stations.size, -1, dtype=np.int64)
    owner[pipe_station] = np.arange(pipe_chainages.size)
    jump = receiver
    for _ in range(40):
        next_jump = jump[jump]
        if np.array_equal(next_jump, jump):
            break
        jump = next_jump
    resolved = owner[jump]
    # 관 없는 사그에 갇힌 구간은 가장 가까운 관에 붙인다.
    orphan = resolved < 0
    if orphan.any() and pipe_chainages.size:
        nearest = np.abs(stations[orphan, None] - pipe_chainages[None, :]).argmin(axis=1)
        resolved[orphan] = nearest

    slot_station = np.clip(np.round(road_chainage / step).astype(np.int64), 0, stations.size - 1)
    return resolved[slot_station].astype(np.int32)


# ── ⑩ 세부유역 조립 ────────────────────────────────────────────────────────


def assemble_basins(
    routing: RoadRouting,
    pipes: list[StructureCandidate],
    pipe_of_slot: np.ndarray,
) -> list[WatershedBasin]:
    """셀이 도달한 도로 셀의 담당 관을 그대로 유역 번호로 삼아 세부유역을 만든다."""
    spec = routing.spec
    labels = np.full(spec.size, -1, dtype=np.int32)
    reached = routing.road_slot >= 0
    labels[reached] = pipe_of_slot[routing.road_slot[reached]]

    polygons = polygonize_labels(spec, labels)
    cell_area = spec.cell_area_m2
    basins: list[WatershedBasin] = []
    for order, pipe in enumerate(pipes):
        member = labels == order
        count = int(member.sum())
        if count == 0:
            continue
        geometry = polygons.get(order)
        elevations = routing.elevation[member]
        highest = float(np.nanmax(elevations)) if np.isfinite(elevations).any() else 0.0
        outlet_z = _outlet_elevation(routing, order, pipe_of_slot)
        area = count * cell_area
        relief = max(0.0, highest - outlet_z)
        flow_length = float(routing.path_length[member].max())
        basins.append(
            WatershedBasin(
                index=len(basins) + 1,
                chainage_m=pipe.chainage_m,
                outlet_x=pipe.x,
                outlet_y=pipe.y,
                boundary_xy=largest_ring(geometry) if geometry is not None else [],
                area_m2=area,
                relief_m=relief,
                flow_length_m=flow_length,
                pipe_diameter_mm=estimate_pipe_diameter_mm(area, relief, flow_length),
            )
        )
    return basins


def _outlet_elevation(routing: RoadRouting, pipe_order: int, pipe_of_slot: np.ndarray) -> float:
    """관이 담당하는 도로 셀들의 최저 표고 = 유역 출구 표고."""
    slots = np.flatnonzero(pipe_of_slot == pipe_order)
    if slots.size == 0:
        return 0.0
    elevations = routing.elevation[routing.road_cell_index[slots]]
    finite = elevations[np.isfinite(elevations)]
    return float(finite.min()) if finite.size else 0.0


def estimate_pipe_diameter_mm(
    area_m2: float,
    relief_m: float,
    flow_length_m: float,
    rainfall_mm_per_hour: float | None = None,
) -> float | None:
    """유역 제원으로 배수 파이프 관경(mm)을 산정한다.

    100년 강우빈도와 유역 경사면을 곱해 유출량을 구하고, 그 유량으로 관경을 정하는 것이
    목적이다. **수식은 아직 확정되지 않았다** — 사용자가 로직을 제공하면 여기를 채운다.
    그때까지는 None을 돌려 호출부가 "미정"으로 표기하게 한다.
    """
    # TODO(사용자 로직 대기): 100년 강우강도 × 유역면적 × 유출계수 → 유량 Q → 관경 D 산정.
    _ = (area_m2, relief_m, flow_length_m, rainfall_mm_per_hour)
    return None