Aislo/common_util/common_util_wamis_station.py

"""map.wamis.go.kr 관측소별 확률강우량 수급·지점 선정 (본토 공사지용).

등우선도는 서버 원본이 제주도만 커버한다(2026-08-13 실측). 본토는 전국 621개
강수량 관측소의 확률강우량 xlsx(2017 지역빈도해석)를 받아 쓴다. 인증 불필요.

    관측소 목록: /data/gis/gis_obsvt.json  (EPSG:3857 Point)
    확률강우량 : /file/down/download.do?path=rainfall/probability&file={관측소코드}.xlsx

지점 선정 — Aislo 자체 로직 (사용자 확정 2026-08-13, 지식DB 설계유량.md §3-2):
  계획노선 대표점의 최근접 관측소 + 반경 10km 내 관측소를 모두 받아
  설계빈도(100년)/24시간 확률강우량이 가장 큰 지점을 선정한다(안전측).
  후보 전 지점 표를 함께 기록한다 — 성과품(수리계산서)의 실무 유사 표 재료.

내려받은 xlsx는 프로젝트 영구저장소(drainage/rainfall_stations/)에 보존한다.
원천이 갱신되는 자료라 프로젝트 시점 스냅샷을 남긴다(전역 공통 캐시 아님).

설계빈도 강우강도(IDF)는 2방식을 모두 기록하고 기본값은 Mononobe형이다:
  - mononobe(기본): rt = (R24/24)·(24/T)^0.557 — 실무 수리계산서 방식(울진·소광).
    t(분) 기준 General형 a/(t^n+b)에서 b=0, n=0.557, a=(R24/24)·1440^0.557과
    정확히 같아 하류 소비자(idf_intensity)를 그대로 쓴다.
  - general: 지속시간별 강우강도점의 General형 적합 — 등우선 방식과 동일.

검증: 울진군(전곡리 20014210) 100년/24hr = 320.0mm (2026-08-13 실측).
"""

from __future__ import annotations

import json
import logging
import math
from pathlib import Path
from typing import Any, Sequence

from common_util.common_util_wamis_rainfall import (
    BASE_URL,
    _download,
    fit_general_idf,
)

logger = logging.getLogger(__name__)

STATION_LIST_PATH = "/data/gis/gis_obsvt.json"
PROBABILITY_DOWNLOAD_PATH = "/file/down/download.do"

#: Mononobe형 지수 (실무 수리계산서 관측값 — 종합비교 02).
MONONOBE_EXPONENT = 0.557

#: 등우선 방식과 같은 지속시간 축(분)으로 표를 만든다 — rainfall_table.json 스키마 호환.
TABLE_DURATIONS_HR: tuple[int, ...] = (1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 18, 24)

#: 제주 판정 상자(WGS84) — 등우선 서버 데이터의 실측 커버리지와 일치.
JEJU_LAT_RANGE = (32.8, 33.7)
JEJU_LON_RANGE = (126.0, 127.1)

_EARTH_RADIUS_M = 6371000.0


def is_jeju(lat: float, lon: float) -> bool:
    """계획노선 좌표의 제주 여부 — 제주면 기존 등우선 로직을 쓴다."""
    return (
        JEJU_LAT_RANGE[0] <= lat <= JEJU_LAT_RANGE[1]
        and JEJU_LON_RANGE[0] <= lon <= JEJU_LON_RANGE[1]
    )


# --------------------------------------------------------------------------
# 관측소 목록·후보 탐색
# --------------------------------------------------------------------------
def _mercator_to_lonlat(x: float, y: float) -> tuple[float, float]:
    lon = math.degrees(x / 6378137.0)
    lat = math.degrees(2.0 * math.atan(math.exp(y / 6378137.0)) - math.pi / 2.0)
    return lon, lat


def _haversine_m(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp, dl = math.radians(lat2 - lat1), math.radians(lon2 - lon1)
    a = math.sin(dp / 2.0) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2.0) ** 2
    return 2.0 * _EARTH_RADIUS_M * math.asin(math.sqrt(a))


def fetch_station_list(stations_dir: Path) -> list[dict[str, Any]]:
    """관측소 목록(코드·이름·위경도). 프로젝트 폴더 스냅샷 우선, 없으면 다운로드."""
    stations_dir = Path(stations_dir)
    stations_dir.mkdir(parents=True, exist_ok=True)
    snapshot = stations_dir / "gis_obsvt.json"
    if snapshot.exists() and snapshot.stat().st_size > 0:
        raw = snapshot.read_bytes()
    else:
        raw = _download(f"{BASE_URL}{STATION_LIST_PATH}")
        json.loads(raw)  # JSON 검증 후에만 저장
        snapshot.write_bytes(raw)
    payload = json.loads(raw)
    stations: list[dict[str, Any]] = []
    for feature in payload.get("features") or []:
        props = feature.get("properties") or {}
        geometry = feature.get("geometry") or {}
        code = str(props.get("OBSCD") or "").strip()
        coords = geometry.get("coordinates")
        if not code or geometry.get("type") != "Point" or not coords:
            continue
        lon, lat = _mercator_to_lonlat(float(coords[0]), float(coords[1]))
        stations.append(
            {"code": code, "name": str(props.get("OBSNM") or code), "lat": lat, "lon": lon}
        )
    if not stations:
        raise RuntimeError("관측소 목록을 읽지 못했습니다 (gis_obsvt.json)")
    return stations


def find_candidates(
    lat: float, lon: float, stations: Sequence[dict[str, Any]], radius_m: float
) -> list[dict[str, Any]]:
    """최근접 관측소 + 반경 내 관측소 (거리 오름차순, 최근접은 반경 밖이어도 포함)."""
    ranked = sorted(
        ({**s, "distance_m": _haversine_m(lat, lon, s["lat"], s["lon"])} for s in stations),
        key=lambda s: s["distance_m"],
    )
    candidates = [s for s in ranked if s["distance_m"] <= radius_m]
    if not candidates:
        candidates = ranked[:1]
    return candidates


# --------------------------------------------------------------------------
# 확률강우량 xlsx 수급·파싱
# --------------------------------------------------------------------------
def download_probability_xlsx(code: str, stations_dir: Path) -> Path:
    """관측소 확률강우량 xlsx — 프로젝트 폴더에 있으면 재사용, 없으면 다운로드."""
    stations_dir = Path(stations_dir)
    stations_dir.mkdir(parents=True, exist_ok=True)
    path = stations_dir / f"{code}.xlsx"
    if path.exists() and path.stat().st_size > 0:
        return path
    query = f"path=rainfall/probability&file={code}.xlsx&name={code}.xlsx"
    url = f"{BASE_URL}{PROBABILITY_DOWNLOAD_PATH}?{query}"
    raw = _download(url)
    if not raw.startswith(b"PK"):  # xlsx(zip) 시그니처 — 오류 페이지 저장 방지
        raise RuntimeError(f"관측소 {code} 확률강우량 xlsx 응답이 올바르지 않습니다")
    tmp = path.with_suffix(".part")
    tmp.write_bytes(raw)
    tmp.replace(path)
    return path


def parse_probability_xlsx(path: Path) -> dict[int, dict[int, float]]:
    """xlsx → {재현기간(년): {지속시간(시간): 강우량(mm)}}.

    시트 구성(2026-08-13 실측): '확률강우량_*' 시트에 헤더행(관측소코드/재현기간 +
    지속기간 1~72) 뒤로 재현기간별 행이 온다.
    """
    from openpyxl import load_workbook  # 무거운 의존 — 사용 시점에만 로드

    workbook = load_workbook(path, read_only=True, data_only=True)
    try:
        sheet = next((s for s in workbook.sheetnames if "확률강우량" in s), None)
        if sheet is None:
            raise RuntimeError(f"확률강우량 시트를 찾지 못했습니다 ({path.name})")
        rows = workbook[sheet].iter_rows(values_only=True)
        durations: list[tuple[int, int]] = []  # (열 번호, 지속시간 hr)
        table: dict[int, dict[int, float]] = {}
        for row in rows:
            if not durations:
                cells = list(row)
                found = [
                    (idx, int(v))
                    for idx, v in enumerate(cells)
                    if isinstance(v, (int, float)) and float(v).is_integer() and 1 <= v <= 72
                ]
                # 헤더행: 지속기간 1,2,3…이 연속으로 나열된 행
                if len(found) >= 12 and found[0][1] == 1:
                    durations = found
                continue
            # 값 행: [관측소코드, 재현기간, 값1hr, 값2hr, …] — 헤더와 같은 열 정렬
            if len(row) < 3 or not isinstance(row[0], (int, float)):
                continue
            period = row[1]
            if not isinstance(period, (int, float)):
                continue
            per_duration = {}
            for idx, hours in durations:
                v = row[idx] if idx < len(row) else None
                if isinstance(v, (int, float)):
                    per_duration[hours] = float(v)
            if per_duration:
                table[int(period)] = per_duration
        if not table:
            raise RuntimeError(f"확률강우량 값을 읽지 못했습니다 ({path.name})")
        return table
    finally:
        workbook.close()


# --------------------------------------------------------------------------
# 강우강도(IDF) — Mononobe(기본)·General형
# --------------------------------------------------------------------------
def mononobe_idf(r24_mm: float) -> dict[str, float]:
    """R24 하나로 전 지속시간 강도 — rt=(R24/24)·(24/T)^0.557 (실무 수리계산서 방식).

    t(분) 기준 General형 a/(t^n+b) 표현: b=0, n=0.557, a=(R24/24)·1440^0.557.
    하류의 idf_intensity()가 그대로 계산한다.
    """
    a = (r24_mm / 24.0) * (1440.0**MONONOBE_EXPONENT)
    return {"a": round(a, 4), "b": 0.0, "n": MONONOBE_EXPONENT, "r24_mm": round(r24_mm, 2)}


# --------------------------------------------------------------------------
# 프로젝트 강우량표 생성 (등우선판 build_rainfall_table과 같은 스키마 + 확장 필드)
# --------------------------------------------------------------------------
def build_station_rainfall_table(
    lat: float,
    lon: float,
    stations_dir: Path,
    design_return_period: int = 100,
    radius_m: float = 10000.0,
    idf_method: str = "mononobe",
) -> dict[str, Any]:
    """관측소 방식 강우량표. 실패는 명시적 RuntimeError — 제주 등우선 폴백 금지."""
    stations = fetch_station_list(stations_dir)
    candidates = find_candidates(lat, lon, stations, radius_m)
    failures: list[str] = []
    evaluated: list[dict[str, Any]] = []
    for cand in candidates:
        try:
            xlsx = download_probability_xlsx(cand["code"], stations_dir)
            data = parse_probability_xlsx(xlsx)
            r24 = data.get(design_return_period, {}).get(24)
            if r24 is None:
                raise RuntimeError(f"{design_return_period}년/24hr 값 없음")
            evaluated.append({**cand, "data": data, "r24_design_mm": round(r24, 2)})
        except (RuntimeError, OSError, ValueError) as exc:
            failures.append(f"관측소 {cand['code']}({cand['name']}): {exc}")
    if not evaluated:
        raise RuntimeError(
            "후보 관측소의 확률강우량을 하나도 받지 못했습니다: " + "; ".join(failures)
        )
    # Aislo 선정 규칙: 설계빈도/24hr 최대 지점 (안전측)
    selected = max(evaluated, key=lambda s: s["r24_design_mm"])
    data = selected["data"]
    values = [
        {
            "return_period_yr": period,
            "duration_min": hours * 60,
            "depth_mm": round(per_duration[hours], 2),
        }
        for period, per_duration in sorted(data.items())
        for hours in TABLE_DURATIONS_HR
        if hours in per_duration
    ]
    design_points = [
        (hours * 60.0, data[design_return_period][hours] / hours)
        for hours in TABLE_DURATIONS_HR
        if hours in data.get(design_return_period, {})
    ]
    idf_general = fit_general_idf(design_points) if design_points else None
    idf_mononobe = mononobe_idf(selected["r24_design_mm"])
    idf_design = idf_mononobe if idf_method == "mononobe" else (idf_general or idf_mononobe)
    return {
        "source": f"{BASE_URL} 관측소 확률강우량 (2017 지역빈도해석, 지점 {selected['name']})",
        "lat": lat,
        "lon": lon,
        "design_return_period_yr": design_return_period,
        "values": values,
        "idf_design": idf_design,
        "idf_method": idf_method,
        "idf_mononobe": idf_mononobe,
        "idf_general": idf_general,
        "station": {
            "code": selected["code"],
            "name": selected["name"],
            "distance_m": round(selected["distance_m"], 1),
            "r24_design_mm": selected["r24_design_mm"],
        },
        # 성과품용 후보 표 — 실무 수리계산서의 "인근 지점 확률강우량 표" 재현 재료
        "candidates": [
            {
                "code": s["code"],
                "name": s["name"],
                "distance_m": round(s["distance_m"], 1),
                "r24_design_mm": s["r24_design_mm"],
                "selected": s["code"] == selected["code"],
            }
            for s in evaluated
        ],
        "failures": failures,
    }