eomsangdon/Aislo
1
0
Fork 0
Code Issues Pull Requests 2 Actions Packages Projects Releases Wiki Activity

feat(B04/B05): 확률강우량 수급·배수 유효직경 합리식 계산 (계획서 Phase 1·2)

Browse Source 
- common_util_wamis_rainfall.py 신설: map.wamis.go.kr 등우선도 urllib 수급,
  전역 캐시(resources/wamis_contours), 거리비례 내삽, General형 IDF 적합
- /drainage/rainfall API + primary-region 백그라운드 강우량표 생성 훅
- estimate_pipe_diameter_mm TODO 구현: Kirpich 도달시간(하한5분) + 합리식(C=0.8,
  2.0배) + Manning(경사10도, n=0.024, V<=3.0) + 통수단면 70% -> 배수 유효직경
- basins 응답·B05 유역 제원에 유효면적/유효직경·산출근거(tc/I/Qd) 표기
- config_system: 계산 계수·구조물 옵션(관종/직경/대피로) 기본값 등록

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This commit is contained in:
eomsangdon
2026-08-05 18:19:48 +09:00
co-authored by Claude Fable 5
parent 5378b9d807
commit fa6953ce6d
8 changed files with 600 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common_util/common_util_drainage_detail.py
+95 -12
View File
@@ -25,6 +25,7 @@ from __future__ import annotations
import json
import logging
import math
from dataclasses import dataclass, field
from pathlib import Path
from typing import Any
@@ -41,10 +42,20 @@ from common_util.common_util_route_geometry import (
interpolate_vertex,
is_uphill_at,
)
from common_util.common_util_wamis_rainfall import idf_intensity
from config.config_system import (
DRAINAGE_DESIGN_FLOW_FACTOR,
DRAINAGE_DITCH_SAMPLE_M,
DRAINAGE_FLOW_AREA_RATIO,
DRAINAGE_MANNING_N,
DRAINAGE_PIPE_MAX_SPACING_M,
DRAINAGE_PIPE_MIN_SPACING_M,
DRAINAGE_PIPE_SLOPE_DEG,
DRAINAGE_RAINFALL_FILENAME,
DRAINAGE_RUNOFF_COEFFICIENT,
DRAINAGE_TC_MIN_MINUTES,
DRAINAGE_VELOCITY_MAX_MS,
DRAINAGE_VELOCITY_MIN_MS,
)
logger = logging.getLogger(__name__)
@@ -89,6 +100,10 @@ class WatershedBasin:
relief_m: float = 0.0
flow_length_m: float = 0.0
pipe_diameter_mm: float | None = None
# 유효직경 산출 근거(합리식) — 도달시간·설계강우강도·설계유량. 미산출이면 None.
tc_minutes: float | None = None
intensity_mm_hr: float | None = None
design_flow_m3s: float | None = None
@dataclass
@@ -172,7 +187,7 @@ def build_detail(
if not pipes:
return detail
pipe_of_slot = assign_road_cells_to_pipes(vertices, pipes, routing.road_chainage)
detail.basins = assemble_basins(routing, pipes, pipe_of_slot)
detail.basins = assemble_basins(routing, pipes, pipe_of_slot, load_rainfall_idf(directory))
logger.info(
"배수유역: 세부 설계 — 관 %d개(기본 %d + 보충 %d), 세부유역 %d개",
len(pipes),
@@ -476,6 +491,7 @@ def assemble_basins(
routing: RoadRouting,
pipes: list[StructureCandidate],
pipe_of_slot: np.ndarray,
idf: dict[str, float] | None = None,
) -> list[WatershedBasin]:
"""셀이 도달한 도로 셀의 담당 관을 그대로 유역 번호로 삼아 세부유역을 만든다."""
spec = routing.spec
@@ -498,6 +514,7 @@ def assemble_basins(
area = count * cell_area
relief = max(0.0, highest - outlet_z)
flow_length = float(routing.path_length[member].max())
sizing = size_pipe(area, relief, flow_length, idf)
basins.append(
WatershedBasin(
index=len(basins) + 1,
@@ -508,7 +525,10 @@ def assemble_basins(
area_m2=area,
relief_m=relief,
flow_length_m=flow_length,
pipe_diameter_mm=estimate_pipe_diameter_mm(area, relief, flow_length),
pipe_diameter_mm=sizing["diameter_mm"] if sizing else None,
tc_minutes=sizing["tc_minutes"] if sizing else None,
intensity_mm_hr=sizing["intensity_mm_hr"] if sizing else None,
design_flow_m3s=sizing["design_flow_m3s"] if sizing else None,
)
)
return basins
@@ -524,18 +544,81 @@ def _outlet_elevation(routing: RoadRouting, pipe_order: int, pipe_of_slot: np.nd
return float(finite.min()) if finite.size else 0.0
def load_rainfall_idf(directory: Path) -> dict[str, float] | None:
"""배수유역 폴더의 rainfall_table.json에서 설계빈도 IDF 적합계수를 읽는다.
파일은 B04 전처리(`/drainage/rainfall`)가 만든다. 없으면 None — 유효직경은
"미정"으로 남고, 강우량표가 생기는 순간 다음 세부유역 계산부터 채워진다.
"""
path = Path(directory) / DRAINAGE_RAINFALL_FILENAME
if not path.exists():
return None
try:
table = json.loads(path.read_text(encoding="utf-8"))
except (OSError, json.JSONDecodeError):
logger.warning("강우량표를 읽지 못했습니다 (%s) — 유효직경 미산출.", path)
return None
idf = table.get("idf_design")
if not isinstance(idf, dict) or not all(k in idf for k in ("a", "b", "n")):
return None
return idf
def size_pipe(
area_m2: float,
relief_m: float,
flow_length_m: float,
idf: dict[str, float] | None,
) -> dict[str, float] | None:
"""유역 제원 + 설계빈도 IDF로 배수 유효직경(mm)과 산출 근거를 계산한다.
시행규칙 별표2 (가) 1호 방법: 100년빈도 확률강우량 + 홍수도달시간 → 합리식.
전 과정·계수 근거: docs/raw/guidelines/2026-08-05_임도_배수설계_규정_조사.md 4·6절.
① 도달시간 tc = 0.0663·L^0.77·(H/L)^-0.385 [Kirpich, hr] — 하한 5분
② 강우강도 I = General형 적합식 I(tc) [mm/hr]
③ 합리식 Q = (1/3.6)·C·I·A[km²], 설계유량 = 2.0·Q (별표2)
④ 유속 V = Manning(경사 10도 고정, n=0.024) → 0.8~3.0m/s 클램프
⑤ 통수단면 70%만 유효 → D = √(4·(Q설계/V)/(0.7π))
"""
if idf is None or area_m2 <= 0:
return None
length_km = max(flow_length_m, 1.0) / 1000.0
slope = max(relief_m, 0.1) / max(flow_length_m, 1.0)
tc_hr = 0.0663 * (length_km**0.77) * (slope**-0.385)
tc_min = max(DRAINAGE_TC_MIN_MINUTES, tc_hr * 60.0)
intensity = idf_intensity(idf, tc_min)
flow = DRAINAGE_RUNOFF_COEFFICIENT * intensity * (area_m2 / 1e6) / 3.6
design_flow = DRAINAGE_DESIGN_FLOW_FACTOR * flow
if design_flow <= 0:
return None
pipe_slope = math.tan(math.radians(DRAINAGE_PIPE_SLOPE_DEG))
diameter = 0.5
for _ in range(20):
velocity = (1.0 / DRAINAGE_MANNING_N) * (diameter / 4.0) ** (2.0 / 3.0) * math.sqrt(
pipe_slope
)
velocity = min(DRAINAGE_VELOCITY_MAX_MS, max(DRAINAGE_VELOCITY_MIN_MS, velocity))
required_area = design_flow / velocity
updated = math.sqrt(4.0 * required_area / (DRAINAGE_FLOW_AREA_RATIO * math.pi))
if abs(updated - diameter) < 1e-4:
diameter = updated
break
diameter = updated
return {
"tc_minutes": round(tc_min, 1),
"intensity_mm_hr": round(intensity, 1),
"design_flow_m3s": round(design_flow, 4),
"diameter_mm": round(diameter * 1000.0, 1),
}
def estimate_pipe_diameter_mm(
area_m2: float,
relief_m: float,
flow_length_m: float,
rainfall_mm_per_hour: float | None = None,
idf: dict[str, float] | None = None,
) -> float | None:
"""유역 제원으로 배수 파이프 관경(mm)을 산정한다.
100년 강우빈도와 유역 경사면을 곱해 유출량을 구하고, 그 유량으로 관경을 정하는 것이
목적이다. **수식은 아직 확정되지 않았다** — 사용자가 로직을 제공하면 여기를 채운다.
그때까지는 None을 돌려 호출부가 "미정"으로 표기하게 한다.
"""
# TODO(사용자 로직 대기): 100년 강우강도 × 유역면적 × 유출계수 → 유량 Q → 관경 D 산정.
_ = (area_m2, relief_m, flow_length_m, rainfall_mm_per_hour)
return None
"""유역 제원으로 배수 유효직경(mm)만 돌려주는 축약형 — 상세는 size_pipe()."""
sizing = size_pipe(area_m2, relief_m, flow_length_m, idf)
return sizing["diameter_mm"] if sizing else None
common_util/common_util_wamis_rainfall.py
+324
View File
@@ -0,0 +1,324 @@
"""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
def value_at(
lines: Sequence[tuple[float, Sequence[tuple[float, float]]]], x: float, y: float
) -> float:
"""등우선에서 지점 ``(x, y)``의 확률강우량(mm) — 거리비례 내부삽입."""
if not lines:
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