340 lines
14 KiB
Python
340 lines
14 KiB
Python
"""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
|