"""도로 유입 셀 조회 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_wf2_Route.B05_wf2_Route_Repository import get_surface_crs_epsg from common_util.common_util_route_geometry import ( find_planned_route_file, read_planned_route_csv, ) 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: """분석에 쓰인 좌표계를 그대로 되찾는다(노선 파일이 명시하면 그 값 우선).""" pool = get_db_pool() async with pool.acquire() as connection: epsg = await get_surface_crs_epsg(connection, project_id, 0) project_root = Path(resolve_stored_project_path(stored_path)) route_file = find_planned_route_file(project_root / "B03_FileInput" / "input") planned = read_planned_route_csv(route_file) if route_file else None source = (planned.epsg if planned else None) or epsg or 5186 return f"EPSG:{source}" 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), }