From 69df5333501278ba62b17b94da1b0ac5ad5599c6 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Thu, 13 Aug 2026 18:51:25 +0900 Subject: [PATCH] =?UTF-8?q?feat(drainage):=20=ED=99=95=EB=A5=A0=EA=B0=95?= =?UTF-8?q?=EC=9A=B0=EB=9F=89=20=EB=B3=B8=ED=86=A0=20=EC=98=A4=EB=A5=98=20?= =?UTF-8?q?=EA=B0=9C=EC=84=A0=20=E2=80=94=20WAMIS=20=EA=B4=80=EC=B8=A1?= =?UTF-8?q?=EC=86=8C=20=EB=B0=A9=EC=8B=9D=20=EB=8F=84=EC=9E=85,=20?= =?UTF-8?q?=EC=A0=9C=EC=A3=BC=EB=8A=94=20=EB=93=B1=EC=9A=B0=EC=84=A0=20?= =?UTF-8?q?=EC=9C=A0=EC=A7=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - common_util_wamis_station.py 신규: 관측소 목록/확률강우량 xlsx 수급(프로젝트별 영구저장), 최근접+반경 10km 후보 중 100년/24hr 최대 지점 선정(Aislo 자체 로직), 후보표 기록(성과품용) - 강우강도 이원화: Mononobe형(기본, 실무 방식 rt=(R24/24)(24/T)^0.557) + General형 적합 옵션 — Mononobe는 b=0 General형과 동일해 하류 idf_intensity/size_pipe 무수정 - _build_rainfall_sync: 계획노선 좌표로 제주/본토 분기, 본토 실패 시 예외 전파(제주 폴백 금지) - value_at: 등우선 커버리지(제주+30km) 밖 좌표 차단 — 조용히 틀린 값 재발 방지 - config: WAMIS_STATION_RADIUS_M/DRAINAGE_RAINFALL_STATION_DIRNAME/DRAINAGE_RAINFALL_IDF_METHOD - 검증: 제주 등우선 96값 정상, 본토 차단, 울진 관측소 R24=305.83, Mononobe 검산 실무 일치 Co-Authored-By: Claude Fable 5 --- .../B04_PreProcess_Router_Watershed.py | 44 ++- common_util/common_util_wamis_rainfall.py | 17 +- common_util/common_util_wamis_station.py | 286 ++++++++++++++++++ config/config_system.py | 17 +- 4 files changed, 349 insertions(+), 15 deletions(-) create mode 100644 common_util/common_util_wamis_station.py diff --git a/B04_PreProcess/B04_PreProcess_Router_Watershed.py b/B04_PreProcess/B04_PreProcess_Router_Watershed.py index 206046f6..64418f83 100644 --- a/B04_PreProcess/B04_PreProcess_Router_Watershed.py +++ b/B04_PreProcess/B04_PreProcess_Router_Watershed.py @@ -46,13 +46,20 @@ from common_util.common_util_wamis_rainfall import ( build_rainfall_table, ensure_contour_cache, ) +from common_util.common_util_wamis_station import ( + build_station_rainfall_table, + is_jeju, +) from config.config_db import get_db_pool from config.config_system import ( DRAINAGE_ARROW_SPACING_M, DRAINAGE_DESIGN_RETURN_PERIOD_YR, DRAINAGE_RAINFALL_FILENAME, + DRAINAGE_RAINFALL_IDF_METHOD, + DRAINAGE_RAINFALL_STATION_DIRNAME, DRAINAGE_RESPONSE_FILENAME, WAMIS_CONTOUR_CACHE_DIR, + WAMIS_STATION_RADIUS_M, ) logger = logging.getLogger(__name__) @@ -249,18 +256,35 @@ def _route_center_lonlat(stored_path: str, fallback_epsg: int | None) -> tuple[f return lat, lon -def _build_rainfall_sync(lat: float, lon: float) -> dict[str, Any]: - """등우선 전역 캐시 확보(없는 파일만 다운로드) 후 지점 강우량표를 만든다.""" - cached, failures = ensure_contour_cache(WAMIS_CONTOUR_CACHE_DIR) - table = build_rainfall_table( +def _build_rainfall_sync(lat: float, lon: float, stored_path: str) -> dict[str, Any]: + """계획노선 좌표로 지역을 분별해 강우량표를 만든다 (PLAN.md A, 2026-08-13). + + 제주: 등우선 내삽(서버 원본이 제주만 커버) / 본토: 관측소 방식(최근접+10km + 최대값). 본토에서 관측소 수급이 실패하면 그대로 예외를 올린다 — 제주 등우선으로 + 폴백하면 "조용히 틀린 값"이 재발하므로 금지. + """ + if is_jeju(lat, lon): + cached, failures = ensure_contour_cache(WAMIS_CONTOUR_CACHE_DIR) + table = build_rainfall_table( + lat, + lon, + WAMIS_CONTOUR_CACHE_DIR, + design_return_period=DRAINAGE_DESIGN_RETURN_PERIOD_YR, + ) + table["contour_cache_files"] = cached + if failures: + table["failures"] = (table.get("failures") or []) + failures + table["region_mode"] = "jeju_contour" + return table + table = build_station_rainfall_table( lat, lon, - WAMIS_CONTOUR_CACHE_DIR, + drainage_dir(stored_path) / DRAINAGE_RAINFALL_STATION_DIRNAME, design_return_period=DRAINAGE_DESIGN_RETURN_PERIOD_YR, + radius_m=WAMIS_STATION_RADIUS_M, + idf_method=DRAINAGE_RAINFALL_IDF_METHOD, ) - table["contour_cache_files"] = cached - if failures: - table["failures"] = (table.get("failures") or []) + failures + table["region_mode"] = "mainland_station" return table @@ -273,7 +297,7 @@ async def _ensure_rainfall_table(stored_path: str, fallback_epsg: int | None) -> if center is None: return try: - table = await asyncio.to_thread(_build_rainfall_sync, *center) + table = await asyncio.to_thread(_build_rainfall_sync, *center, stored_path) path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(table, ensure_ascii=False, indent=1), encoding="utf-8") logger.info("확률강우량표 저장: %s (값 %d개)", path, len(table.get("values") or [])) @@ -307,7 +331,7 @@ async def get_drainage_rainfall( content={"status": "error", "message": "계획 노선 파일이 업로드되지 않았습니다."}, ) try: - table = await asyncio.to_thread(_build_rainfall_sync, *center) + table = await asyncio.to_thread(_build_rainfall_sync, *center, stored_path) except Exception as exc: # noqa: BLE001 return JSONResponse( status_code=502, diff --git a/common_util/common_util_wamis_rainfall.py b/common_util/common_util_wamis_rainfall.py index 7831cdca..4ce8f31f 100644 --- a/common_util/common_util_wamis_rainfall.py +++ b/common_util/common_util_wamis_rainfall.py @@ -127,12 +127,27 @@ def _iter_linestrings(geometry: dict) -> Iterable[list]: yield from polygon +#: 등우선 커버리지 밖 판정 여유(m). 서버 등우선은 제주도만 커버(2026-08-13 실측) — +#: 본토 좌표를 넣으면 수백 km 밖 보간으로 "조용히 틀린 값"이 나오므로 차단한다. +COVERAGE_MARGIN_M = 30_000.0 + + def value_at( lines: Sequence[tuple[float, Sequence[tuple[float, float]]]], x: float, y: float ) -> float: - """등우선에서 지점 ``(x, y)``의 확률강우량(mm) — 거리비례 내부삽입.""" + """등우선에서 지점 ``(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: diff --git a/common_util/common_util_wamis_station.py b/common_util/common_util_wamis_station.py new file mode 100644 index 00000000..c7af45a7 --- /dev/null +++ b/common_util/common_util_wamis_station.py @@ -0,0 +1,286 @@ +"""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, + } diff --git a/config/config_system.py b/config/config_system.py index 8f1c54b4..0647bcea 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -347,13 +347,22 @@ DRAINAGE_ARROW_MIN_AGREEMENT = float(os.getenv("DRAINAGE_ARROW_MIN_AGREEMENT", " # 단계별 검증 산출물은 같은 폴더에 `{번호}_{단계}.geojson` + `manifest.json`으로 쌓인다. # 파일명 규칙은 B05_Profile_Engine_Watershed_Export.STAGES가 유일한 정의처다. -# ── 확률강우량 (map.wamis.go.kr 등우선도) ── -# 전국 공통 등우선 GeoJSON 원본을 이 폴더에 1회 내려받아 영구 캐시한다(약 20~30MB). -# 근거·검증: docs/raw/guidelines/2026-08-05_임도_배수설계_규정_조사.md 5절. +# ── 확률강우량 (map.wamis.go.kr 등우선도 — 제주 전용) ── +# 등우선 서버 원본은 제주도만 커버한다(2026-08-13 실측) — 제주 공사지에서만 쓴다. +# 근거·검증: docs/raw/guidelines/2026-08-05_임도_배수설계_규정_조사.md 5절, PLAN.md A. WAMIS_CONTOUR_CACHE_DIR = PROJECT_ROOT / "resources" / "data_rainfall_idf_cache" -# 프로젝트별 내삽 결과 파일명. drainage/ 아래에 놓인다. +# 프로젝트별 결과 파일명. drainage/ 아래에 놓인다. DRAINAGE_RAINFALL_FILENAME = "rainfall_table.json" +# ── 확률강우량 (관측소 방식 — 본토, 2026-08-13 사용자 확정 로직) ── +# 최근접 관측소 + 이 반경 안 관측소를 모두 받아 100년/24hr 최대 지점을 선정한다. +WAMIS_STATION_RADIUS_M = float(os.getenv("WAMIS_STATION_RADIUS_M", "10000")) +# 관측소 목록·확률강우량 xlsx의 프로젝트별 보존 폴더명. drainage/ 아래에 놓인다 +# (원천이 갱신되는 자료라 프로젝트 시점 스냅샷 — 전역 공통 캐시 아님). +DRAINAGE_RAINFALL_STATION_DIRNAME = "rainfall_stations" +# 설계빈도 강우강도 산출 방식: "mononobe"(기본, 실무 수리계산서 방식) | "general" +DRAINAGE_RAINFALL_IDF_METHOD = os.getenv("DRAINAGE_RAINFALL_IDF_METHOD", "mononobe") + # ── 배수 유효직경 계산 계수 (2026-08-05 사용자 확정, 위 규정 조사 문서 6절) ── # 합리식 유출계수 C. 산악지 준용 0.8 — 필요 시 사용자가 이 값을 고친다. DRAINAGE_RUNOFF_COEFFICIENT = float(os.getenv("DRAINAGE_RUNOFF_COEFFICIENT", "0.8"))