"""map.wamis.go.kr 확률강우량 등우선도 수급·내삽 (표준 라이브러리만 사용).

「전국 하천수계 홍수량 및 확률강우량 정보」가 공개하는 등우선 GeoJSON에서
지점의 지속시간×재현기간 확률강우량표를 만든다. 인증 불필요.

    https://map.wamis.go.kr/data/gis/contour/{재현기간3자리}yr_{지속시간}.json

파일 안 등우선의 ``ELEV`` 값이 그 선의 확률강우량(mm), 좌표계는 EPSG:3857.
지점값은 값이 다른 인접 등우선 두 개까지의 거리비례 내부삽입으로 얻는다.

- 원본 96개(8빈도×12지속시간)는 전국 공통이라 전역 폴더에 1회 캐시한다.
- 1시간 미만 자료는 사이트에 없다(2026-08 실측). 5분 강우강도는 100년빈도
  12점을 General형 강우강도식으로 적합해 얻는다(설계 근거:
  docs/raw/guidelines/2026-08-05_임도_배수설계_규정_조사.md).

검증: 2026-08-05 실서버 실측 — 문경 100년 1시간 142.5mm / 24시간 387.5mm.
"""

from __future__ import annotations

import json
import logging
import math
import urllib.error
import urllib.request
from pathlib import Path
from typing import Any, Callable, Iterable, Sequence

logger = logging.getLogger(__name__)

BASE_URL = "https://map.wamis.go.kr"
CONTOUR_PATH = "/data/gis/contour"

#: 사이트가 제공하는 재현기간(년). 2026-08 실측 8종.
RETURN_PERIODS: tuple[int, ...] = (2, 10, 30, 50, 80, 100, 200, 500)

#: 사이트가 제공하는 지속시간 (파일명 표기, 분). 2026-08 실측 12종 — 1시간 미만 없음.
DURATIONS: tuple[tuple[str, int], ...] = (
    ("01hr", 60),
    ("02hr", 120),
    ("03hr", 180),
    ("04hr", 240),
    ("05hr", 300),
    ("06hr", 360),
    ("08hr", 480),
    ("09hr", 540),
    ("10hr", 600),
    ("12hr", 720),
    ("18hr", 1080),
    ("24hr", 1440),
)

REQUEST_TIMEOUT = 120.0
EARTH_RADIUS_M = 6378137.0

_HEADERS = {
    "User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64)",
    "Accept": "application/json, text/plain, */*",
}


# --------------------------------------------------------------------------
# 좌표·주소
# --------------------------------------------------------------------------
def to_web_mercator(lat: float, lon: float) -> tuple[float, float]:
    """위경도(WGS84) → 웹지도 좌표(EPSG:3857) 미터."""
    if not -90.0 < lat < 90.0:
        raise ValueError(f"위도 범위를 벗어났습니다: {lat}")
    if not -180.0 <= lon <= 180.0:
        raise ValueError(f"경도 범위를 벗어났습니다: {lon}")
    x = EARTH_RADIUS_M * math.radians(lon)
    y = EARTH_RADIUS_M * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0))
    return x, y


def contour_url(return_period: int, duration_label: str, base_url: str = BASE_URL) -> str:
    """``100yr_01hr.json`` 형태의 주소. 재현기간은 3자리 0채움."""
    return f"{base_url.rstrip('/')}{CONTOUR_PATH}/{return_period:03d}yr_{duration_label}.json"


def cache_filename(return_period: int, duration_label: str) -> str:
    return f"{return_period:03d}yr_{duration_label}.json"


# --------------------------------------------------------------------------
# 등우선 읽기·내삽
# --------------------------------------------------------------------------
def parse_contours(payload: dict) -> list[tuple[float, list[tuple[float, float]]]]:
    """GeoJSON에서 ``(강우량mm, 좌표열)`` 목록을 뽑는다."""
    if not isinstance(payload, dict) or "features" not in payload:
        raise RuntimeError("등우선 자료 형식이 아닙니다 (features 없음)")
    lines: list[tuple[float, list[tuple[float, float]]]] = []
    for feature in payload.get("features") or []:
        if not isinstance(feature, dict):
            continue
        props = feature.get("properties") or {}
        elev = props.get("ELEV", props.get("elev"))
        try:
            level = float(elev)
        except (TypeError, ValueError):
            continue
        geometry = feature.get("geometry") or {}
        for coords in _iter_linestrings(geometry):
            points = [
                (float(p[0]), float(p[1]))
                for p in coords
                if isinstance(p, (list, tuple)) and len(p) >= 2
            ]
            if len(points) >= 2:
                lines.append((level, points))
    if not lines:
        raise RuntimeError("등우선을 하나도 읽지 못했습니다 (ELEV 값 확인 필요)")
    return lines


def _iter_linestrings(geometry: dict) -> Iterable[list]:
    kind = geometry.get("type")
    coords = geometry.get("coordinates") or []
    if kind == "LineString":
        yield coords
    elif kind == "MultiLineString":
        yield from coords
    elif kind == "Polygon":
        yield from coords
    elif kind == "MultiPolygon":
        for polygon in coords:
            yield from polygon


#: 등우선 커버리지 밖 판정 여유(m). 서버 등우선은 제주도만 커버(2026-08-13 실측) —
#: 본토 좌표를 넣으면 수백 km 밖 보간으로 "조용히 틀린 값"이 나오므로 차단한다.
COVERAGE_MARGIN_M = 10_000.0


def value_at(
    lines: Sequence[tuple[float, Sequence[tuple[float, float]]]], x: float, y: float
) -> float:
    """등우선에서 지점 ``(x, y)``의 확률강우량(mm) — 거리비례 내부삽입.

    지점이 등우선 범위(+여유 30km) 밖이면 RuntimeError — 커버리지 밖 외삽 금지.
    """
    if not lines:
        raise RuntimeError("등우선이 비어 있습니다")
    xs = [p[0] for _elev, points in lines for p in points]
    ys = [p[1] for _elev, points in lines for p in points]
    if not (
        min(xs) - COVERAGE_MARGIN_M <= x <= max(xs) + COVERAGE_MARGIN_M
        and min(ys) - COVERAGE_MARGIN_M <= y <= max(ys) + COVERAGE_MARGIN_M
    ):
        raise RuntimeError("지점이 등우선 커버리지(제주도) 밖입니다 — 관측소 방식을 사용하세요")
    distances = sorted((_distance_to_polyline(x, y, points), elev) for elev, points in lines)
    d1, e1 = distances[0]
    if d1 == 0.0:
        return e1
    for d2, e2 in distances[1:]:
        if e2 != e1:
            total = d1 + d2
            return e1 if total == 0 else e1 + (e2 - e1) * (d1 / total)
    return e1


def _distance_to_polyline(x: float, y: float, points: Sequence[tuple[float, float]]) -> float:
    return min(_distance_to_segment(x, y, p0, p1) for p0, p1 in zip(points, points[1:]))


def _distance_to_segment(x, y, p0, p1) -> float:
    x0, y0 = p0
    x1, y1 = p1
    dx, dy = x1 - x0, y1 - y0
    if dx == 0.0 and dy == 0.0:
        return math.hypot(x - x0, y - y0)
    t = ((x - x0) * dx + (y - y0) * dy) / (dx * dx + dy * dy)
    t = max(0.0, min(1.0, t))
    return math.hypot(x - (x0 + t * dx), y - (y0 + t * dy))


# --------------------------------------------------------------------------
# 내려받기·전역 캐시
# --------------------------------------------------------------------------
def _download(url: str, timeout: float = REQUEST_TIMEOUT) -> bytes:
    request = urllib.request.Request(url, headers=_HEADERS)
    with urllib.request.urlopen(request, timeout=timeout) as response:
        return response.read()


def ensure_contour_cache(
    cache_dir: Path,
    return_periods: Sequence[int] = RETURN_PERIODS,
    durations: Sequence[tuple[str, int]] = DURATIONS,
    base_url: str = BASE_URL,
    progress: Callable[[str], None] | None = None,
) -> tuple[int, list[str]]:
    """전역 캐시 폴더에 없는 등우선 원본만 내려받는다.

    Returns (캐시 확보 파일 수, 실패 메시지 목록). 실패는 비치명 — 있는 것만으로 진행.
    """
    cache_dir = Path(cache_dir)
    cache_dir.mkdir(parents=True, exist_ok=True)
    ok, failures = 0, []
    for period in return_periods:
        for label, _minutes in durations:
            path = cache_dir / cache_filename(period, label)
            if path.exists() and path.stat().st_size > 0:
                ok += 1
                continue
            url = contour_url(period, label, base_url)
            if progress:
                progress(f"등우선 내려받기 {period}년/{label}")
            try:
                raw = _download(url)
                json.loads(raw)  # JSON 검증 후에만 저장
                tmp = path.with_suffix(".part")
                tmp.write_bytes(raw)
                tmp.replace(path)
                ok += 1
            except (urllib.error.URLError, OSError, ValueError) as exc:
                failures.append(f"{period}년/{label}: {exc}")
    return ok, failures


def _load_cached_lines(
    cache_dir: Path, period: int, label: str
) -> list[tuple[float, list[tuple[float, float]]]] | None:
    path = Path(cache_dir) / cache_filename(period, label)
    if not path.exists():
        return None
    try:
        return parse_contours(json.loads(path.read_text(encoding="utf-8")))
    except (OSError, ValueError, RuntimeError) as exc:
        logger.warning("등우선 캐시 파일을 읽지 못했습니다 (%s): %s", path, exc)
        return None


# --------------------------------------------------------------------------
# General형 강우강도식 적합 — I(t) = a / (t^n + b)  [t: 분, I: mm/hr]
# --------------------------------------------------------------------------
def fit_general_idf(points: Sequence[tuple[float, float]]) -> dict[str, float] | None:
    """``(지속시간 분, 강우강도 mm/hr)`` 점들로 General형 계수를 적합한다.

    scipy가 있으면 curve_fit, 없거나 실패하면 b·n 격자탐색 + a 최소제곱 폴백.
    1시간 미만(도달시간 5분~) 외삽의 근거식이므로 적합 잔차(rmse)를 함께 기록한다.
    """
    pts = [(float(t), float(i)) for t, i in points if t > 0 and i > 0]
    if len(pts) < 4:
        return None
    t = [p[0] for p in pts]
    i = [p[1] for p in pts]

    def rmse(a: float, b: float, n: float) -> float:
        return math.sqrt(sum((a / (tk**n + b) - ik) ** 2 for tk, ik in zip(t, i)) / len(t))

    best: tuple[float, float, float, float] | None = None
    try:
        from scipy.optimize import curve_fit  # 프로젝트 GIS 스택에 이미 포함

        popt, _ = curve_fit(
            lambda tt, a, b, n: a / (tt**n + b),
            t,
            i,
            p0=(3000.0, 15.0, 0.7),
            bounds=((1.0, 0.0, 0.05), (1e6, 1e4, 1.5)),
            maxfev=20000,
        )
        a, b, n = (float(v) for v in popt)
        best = (rmse(a, b, n), a, b, n)
    except Exception:  # noqa: BLE001 — 폴백 격자탐색으로 계속
        pass
    if best is None:
        for n in (x / 100.0 for x in range(30, 121, 5)):
            for b in (0.0, 1.0, 2.0, 5.0, 10.0, 20.0, 40.0, 80.0):
                denom = [tk**n + b for tk in t]
                # 최소제곱: min Σ(a/d − i)² → a = Σ(i/d) / Σ(1/d²)
                a = sum(ik / dk for ik, dk in zip(i, denom)) / sum(1.0 / (dk * dk) for dk in denom)
                err = rmse(a, b, n)
                if best is None or err < best[0]:
                    best = (err, a, b, n)
    if best is None:
        return None
    err, a, b, n = best
    return {"a": round(a, 4), "b": round(b, 4), "n": round(n, 5), "rmse_mm_hr": round(err, 3)}


def idf_intensity(idf: dict[str, float], duration_min: float) -> float:
    """General형 계수로 지속시간(분)의 강우강도(mm/hr)를 계산한다."""
    return float(idf["a"]) / (duration_min ** float(idf["n"]) + float(idf["b"]))


# --------------------------------------------------------------------------
# 프로젝트 강우량표 생성
# --------------------------------------------------------------------------
def build_rainfall_table(
    lat: float,
    lon: float,
    cache_dir: Path,
    design_return_period: int = 100,
    return_periods: Sequence[int] = RETURN_PERIODS,
    durations: Sequence[tuple[str, int]] = DURATIONS,
) -> dict[str, Any]:
    """전역 캐시된 등우선에서 지점 강우량표 + 설계빈도 IDF 적합계수를 만든다."""
    x, y = to_web_mercator(lat, lon)
    values: list[dict[str, float]] = []
    failures: list[str] = []
    design_points: list[tuple[float, float]] = []
    for period in return_periods:
        for label, minutes in durations:
            lines = _load_cached_lines(cache_dir, period, label)
            if lines is None:
                failures.append(f"{period}년/{label}: 캐시 없음")
                continue
            try:
                depth = value_at(lines, x, y)
            except RuntimeError as exc:
                failures.append(f"{period}년/{label}: {exc}")
                continue
            values.append(
                {"return_period_yr": period, "duration_min": minutes, "depth_mm": round(depth, 2)}
            )
            if period == design_return_period:
                design_points.append((minutes, depth / (minutes / 60.0)))
    idf = fit_general_idf(design_points) if design_points else None
    return {
        "source": f"{BASE_URL} 확률강우량 등우선도 (거리비례 내부삽입)",
        "lat": lat,
        "lon": lon,
        "design_return_period_yr": design_return_period,
        "values": values,
        # 설계빈도 IDF: I(t)=a/(t^n+b), t분 → mm/hr. 도달시간(≥5분) 강우강도 산출용.
        "idf_design": idf,
        "failures": failures,
    }


def design_depth_mm(table: dict[str, Any], duration_min: int) -> float | None:
    """저장된 표에서 설계빈도·지속시간 강우량(mm)을 찾는다."""
    period = table.get("design_return_period_yr")
    for row in table.get("values") or []:
        if row.get("return_period_yr") == period and row.get("duration_min") == duration_min:
            return float(row["depth_mm"])
    return None