"""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, }