Aislo/B04_PreProcess/B04_PreProcess_Router_Inflow.py

"""도로 유입 셀 조회 API (B04 — 흐름 강도 검토용).

화면에서 **유입 집중점 마커**를 고르면, 그 지점으로 실제 물이 들어오는 셀들의 외곽선을
돌려준다. 계산은 하지 않는다 — 배수유역 분석이 이미 만들어 둔 `03_road_routing.npz`의
`road_slot`(셀 → 도달 도로 셀)을 읽어 해당 도로 셀에 귀속된 셀만 골라낼 뿐이다.

그래서 원본 격자(1m, 방향 지정 원본) 그대로이며 새로 근사하거나 평균 내지 않는다.
"""

from __future__ import annotations

import asyncio
import logging
from pathlib import Path
from typing import Any
from uuid import UUID

import numpy as np
from fastapi import APIRouter
from fastapi.responses import JSONResponse
from pyproj import Transformer

from B03_FileInput.B03_FileInput_Repository import get_project_storage_relative_path
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Export import STAGES, drainage_dir
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import polygonize_labels
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec
from B05_Profile.B05_Profile_Repository import get_surface_crs_epsg
from common_util.common_util_route_geometry import load_design_route
from common_util.common_util_storage import resolve_stored_project_path
from config.config_db import get_db_pool

logger = logging.getLogger(__name__)
router = APIRouter(prefix="/api/projects", tags=["B04 Surface Inflow"])

# 읽어 둔 배열을 프로세스 메모리에 들고 있는다. 마커를 누를 때마다 수십 MB npz를 다시 푸는
# 것을 막는다. 파일이 새로 쓰이면(재분석) mtime이 달라져 자동으로 버려진다.
_CACHE: dict[str, tuple[float, dict[str, Any]]] = {}
# 유입 셀 외곽선은 조각이 많을 수 있다. 화면에서 읽을 수 있는 수준까지만 보낸다.
_MAX_RINGS = 40


def _routing_path(stored_path: str) -> Path:
    return drainage_dir(stored_path) / f"{STAGES['road_routing']}_road_routing.npz"


def _load_routing(stored_path: str) -> dict[str, Any] | None:
    """`03_road_routing.npz`를 읽어 필요한 배열만 남긴다(mtime 기준 캐시)."""
    path = _routing_path(stored_path)
    if not path.exists():
        return None
    stamp = path.stat().st_mtime
    cached = _CACHE.get(stored_path)
    if cached is not None and cached[0] == stamp:
        return cached[1]
    with np.load(path, allow_pickle=False) as data:
        loaded = {
            "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"]),
            ),
            "road_slot": data["road_slot"],
            "road_chainage": data["road_chainage"],
            "path_length": data["path_length"],
        }
    _CACHE[stored_path] = (stamp, loaded)
    logger.info("배수유역: 도로 귀속 배열을 읽었습니다 (%s).", path.name)
    return loaded


async def _resolve_epsg(project_id: UUID, stored_path: str) -> str:
    """분석에 쓰인 좌표계를 그대로 되찾는다.

    격자 산출물은 `load_design_route()`가 맞춘 **사업지(.prj) 좌표계**에 있다. 노선 CSV의
    `crs_epsg` 열은 표시용 라벨이라 그 값을 쓰면 좌표가 딴 곳으로 간다
    (2026-09-01 실측: 라벨 5179, 실제 5176 — 유입 폴리곤이 1,500km 밖에 찍혔다).
    """
    from B04_PreProcess.B04_PreProcess_Engine_Extent import project_epsg_from_prj

    project_root = Path(resolve_stored_project_path(stored_path))
    planned = load_design_route(project_root)
    if planned is not None and planned.crs_input:
        return planned.crs_input
    prj_crs = project_epsg_from_prj(project_root)
    if prj_crs:
        return prj_crs
    pool = get_db_pool()
    async with pool.acquire() as connection:
        epsg = await get_surface_crs_epsg(connection, project_id, 0)
    return f"EPSG:{epsg or 5186}"


def _collect_inflow(
    routing: dict[str, Any], chainage_m: float, span_m: float
) -> tuple[np.ndarray, Any]:
    """지정한 누가거리 구간의 도로 셀에 귀속된 셀 마스크와 외곽 폴리곤을 만든다."""
    spec: GridSpec = routing["spec"]
    road_slot: np.ndarray = routing["road_slot"]
    road_chainage: np.ndarray = routing["road_chainage"]

    half = max(span_m, spec.cell_m) / 2.0
    slots = np.flatnonzero(np.abs(road_chainage - chainage_m) <= half)
    if slots.size == 0:
        return np.zeros(0, dtype=bool), None

    # 슬롯 번호 조회표로 한 번에 거른다(도로 셀 수가 적어 표가 작다).
    selected = np.zeros(road_chainage.size, dtype=bool)
    selected[slots] = True
    mask = (road_slot >= 0) & selected[np.maximum(road_slot, 0)]
    if not mask.any():
        return mask, None

    labels = np.where(mask, 0, -1).astype(np.int32)
    # 최소 면적 걸러내기와 단순화를 **둘 다 끈다**. 이 외곽선은 "어느 셀이 이 지점으로
    # 들어오는가"를 눈으로 대조하는 용도라, 몇 셀이 빠지거나 경계가 2m 뭉개지면 표시된
    # 면적·셀 수와 그림이 어긋난다.
    polygons = polygonize_labels(spec, labels, min_area_m2=0.0, simplify_m=0.0)
    return mask, polygons.get(0)


def _rings_lonlat(geometry: Any, to_lonlat: Transformer) -> list[list[list[float]]]:
    """폴리곤(멀티 포함)의 바깥 링들을 WGS84 좌표 배열로 바꾼다. 큰 조각부터."""
    if geometry is None or geometry.is_empty:
        return []
    parts = list(geometry.geoms) if geometry.geom_type == "MultiPolygon" else [geometry]
    parts.sort(key=lambda part: part.area, reverse=True)
    rings: list[list[list[float]]] = []
    for part in parts[:_MAX_RINGS]:
        ring = [list(to_lonlat.transform(x, y)) for x, y in part.exterior.coords]
        if len(ring) >= 4:
            rings.append(ring)
    return rings


@router.get("/{project_id}/drainage/road-inflow", response_model=None)
async def get_road_inflow(
    project_id: UUID, chainage_m: float, span_m: float = 1.0
) -> dict[str, Any] | JSONResponse:
    """도로 위 한 지점(누가거리)으로 물이 들어오는 셀들의 외곽선을 돌려준다.

    `span_m`은 그 지점을 중심으로 한 도로 구간 길이다. 기본 1m로, 화면에서 색을 칠하는
    단위와 같다(2026-08-01 사용자 지시).
    """
    pool = get_db_pool()
    async with pool.acquire() as connection:
        stored_path = await get_project_storage_relative_path(connection, project_id)

    routing = await asyncio.to_thread(_load_routing, stored_path)
    if routing is None:
        return JSONResponse(
            status_code=404,
            content={
                "status": "error",
                "message": "배수유역 분석 결과가 없습니다. [유역 분석]을 먼저 실행하세요.",
            },
        )

    mask, geometry = await asyncio.to_thread(_collect_inflow, routing, chainage_m, span_m)
    spec: GridSpec = routing["spec"]
    cell_count = int(mask.sum()) if mask.size else 0
    path_length = routing["path_length"]
    max_path = float(path_length[mask].max()) if cell_count else 0.0

    epsg = await _resolve_epsg(project_id, stored_path)
    to_lonlat = Transformer.from_crs(epsg, "EPSG:4326", always_xy=True)
    return {
        "status": "success",
        "chainage_m": round(chainage_m, 2),
        "span_m": span_m,
        "cell_count": cell_count,
        "area_m2": round(cell_count * spec.cell_area_m2, 1),
        "max_path_length_m": round(max_path, 1),
        "rings_lonlat": _rings_lonlat(geometry, to_lonlat),
    }