eomsangdonandClaude Opus 5
8393dfd619
100㎡ 미만 유역은 `polygonize_labels` 의 최소면적 필터에 걸려 폴리곤이 아예 만들어지지 않았다. 유역 목록에는 있는데 지도에는 없는 상태 — 실측에서 용화 22개 중 1개(64㎡), S자 48개 중 1개가 그랬다. 링 목록이 조각을 싣게 된 뒤로는 버릴 이유가 없다. 실측: 폴리곤 면적 오차(셀 면적 대비 평균) 용화 5.76%→1.86%, S자 3.56%→1.42%, 경계 없는 유역 1개→0개. 좌표점은 583→622 · 806→829 로 거의 그대로다. 1차 전체 배수유역 등 다른 호출부의 기본값(100㎡)은 그대로 둔다. S자 실사례 검증: 용화 지형 위에 스위치백 3구간 합성 노선(정점 439·연장 4,528m)으로 신규 프로젝트를 만들어 자동 체인까지 완주 — 관 48개·세부유역 48개, 경계 없는 유역 0개. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
693 lines
31 KiB
Python
693 lines
31 KiB
Python
"""배수유역 세부 설계 공용 엔진 (B04 관리자 화면 · B05 사용자 화면 공용).
|
||
|
||
**격자 해석은 하지 않는다.** B04 배수유역 분석이 미리 돌려 저장한 결과를 읽어, 사용자가
|
||
실제로 손대는 두 가지만 처리한다(2026-07-31 사용자 지시).
|
||
|
||
⑨ 관 간격이 최대치를 넘는 구간에 **최소 개수**로 관을 보충
|
||
⑩ 측구 흐름으로 도로 셀 → 담당 관을 정하고, 셀이 도달한 도로 셀의 담당 관을 그대로
|
||
그 셀의 유역 번호로 삼아 세부유역을 나눈다
|
||
⑪ 사용자가 관을 옮기거나 추가하면 ⑩만 다시 돈다 — 격자 해석은 재사용한다
|
||
|
||
읽어 오는 것(`{배수유역 폴더}/`):
|
||
· `03_road_routing.geojson` — 계획도로선 · 기본 배관 · 2차 전체 배수유역
|
||
· `03_road_routing.npz` — 셀 → 도로 셀 귀속, 유하장, 강도, 도로 셀 제원, 셀 표고
|
||
화살표(방향 코드)나 밴드 표고 같은 관리자 확인용 배열은 읽지 않는다 — 여기서는 필요 없고
|
||
파일만 무거워진다.
|
||
|
||
**왜 common_util인가**: B04(관리자 트러블슈팅)와 B05(일반 사용자)가 같은 이름의 버튼을
|
||
누르면 같은 결과가 나와야 한다(2026-08-01 사용자 지시). 두 벌로 두면 언젠가 갈라진다.
|
||
읽는 폴더만 다르므로 폴더를 인자로 받고, B04/B05 각자의 어댑터가 경로를 정한다.
|
||
격자 산출물의 규격(`GridSpec`·`STAGES`·`polygonize_labels`)은 B04가 만든 것이므로
|
||
정의처인 B04 엔진을 그대로 참조한다(역방향 참조 없음).
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import json
|
||
import logging
|
||
import math
|
||
from dataclasses import dataclass, field
|
||
from pathlib import Path
|
||
from typing import Any
|
||
|
||
import numpy as np
|
||
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Analyze import find_inflow_hotspots
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Export import STAGES
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import polygon_parts, polygonize_labels
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec
|
||
from common_util.common_util_drainage_pipes import (
|
||
PIPE_FACILITY_BOX,
|
||
PIPE_FACILITY_FORD_BRIDGE,
|
||
PIPE_FACILITY_PIPE,
|
||
)
|
||
from common_util.common_util_route_geometry import (
|
||
RouteVertex,
|
||
StructureCandidate,
|
||
interpolate_vertex,
|
||
is_uphill_at,
|
||
)
|
||
from common_util.common_util_wamis_rainfall import idf_intensity
|
||
from config.config_system import (
|
||
DRAINAGE_BOX_THRESHOLD_MM,
|
||
DRAINAGE_BRIDGE_THRESHOLD_MM,
|
||
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_RECOMMEND_DIAMETERS_MM,
|
||
DRAINAGE_RUNOFF_COEFFICIENT,
|
||
DRAINAGE_TC_MIN_MINUTES,
|
||
DRAINAGE_VELOCITY_MAX_MS,
|
||
DRAINAGE_VELOCITY_MIN_MS,
|
||
PIPE_DEFAULT_DIAMETER_MM,
|
||
)
|
||
|
||
logger = logging.getLogger(__name__)
|
||
|
||
# 관 위치 선정 점수 배분 — 흐름 강도가 주, 종단 저점이 보조.
|
||
_SCORE_WEIGHT_STRENGTH = 0.7
|
||
_SCORE_WEIGHT_SAG = 0.3
|
||
# 성토부(내리막)는 물이 노선 밖으로 빠지므로 관 위치로 덜 선호한다.
|
||
_SCORE_FILL_PENALTY = 0.5
|
||
|
||
|
||
@dataclass
|
||
class DrainageDetail:
|
||
"""세부 설계 산출물 — 화면에 그릴 기하와 세부유역."""
|
||
|
||
route_lonlat: list[list[float]] = field(default_factory=list)
|
||
basin_lonlat: list[list[float]] = field(default_factory=list)
|
||
pipes: list[StructureCandidate] = field(default_factory=list)
|
||
basins: list[WatershedBasin] = field(default_factory=list)
|
||
grid_cell_m: float = 1.0
|
||
# B04가 계산해 둔 평균 흐름 화살표를 그대로 넘긴다 — 여기서 다시 계산하지 않는다.
|
||
flow_arrows: list[list[Any]] = field(default_factory=list)
|
||
arrow_spacing_m: float = 0.0
|
||
# 유역 안쪽 상류 세류망(WGS84 lon/lat 조각들). 화면 강조 표시용 — 계산에는 쓰지 않는다.
|
||
upstream_lonlat: list[list[list[float]]] = field(default_factory=list)
|
||
# 도로 1m 구간별 유입 면적(㎡) — [누가거리, 면적]. 계획선을 색으로 칠하는 데 쓴다.
|
||
strength_profile: list[list[float]] = field(default_factory=list)
|
||
# 유입 집중점 — [누가거리, 유입면적, 구역번호, 구역 내 순위]. 관 자리를 판단하는 근거.
|
||
inflow_hotspots: list[list[float]] = field(default_factory=list)
|
||
|
||
|
||
@dataclass
|
||
class WatershedBasin:
|
||
"""관 하나가 받는 세부 배수유역."""
|
||
|
||
index: int
|
||
chainage_m: float
|
||
outlet_x: float
|
||
outlet_y: float
|
||
# 유역 경계 — 조각마다 [외곽 링, 구멍 링...]. 도넛(아래 유역이 위 유역을 감싼 경우)과
|
||
# 떨어진 조각을 그대로 싣는다. 단일 링만 쓰던 시절에는 이 둘이 소실돼 화면에서 중첩·
|
||
# 빈공간으로 보였다(2026-09-03).
|
||
boundary_parts: list[list[list[tuple[float, float]]]] = field(default_factory=list)
|
||
area_m2: float = 0.0
|
||
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
|
||
# 유효직경이 관 최대 규격을 넘는 계류 유역 — 관이 아니라 세월교·물넘이·교량 대상
|
||
# (임도설치규정 제12조). 화면은 "세월교 검토"로 표기하고 관경 자동 지정을 하지 않는다.
|
||
bridge_required: bool = False
|
||
# 필요 통수단면적(㎡) = 설계유량 / 유속. 물넘이·세월교 개략 단면의 출발값이다.
|
||
required_area_m2: float | None = None
|
||
# 유효직경으로 고른 추천 구조물·관경 (2026-08-17 사용자 확정: 유량 근거만).
|
||
# 배관이면 레지스트리 선택지로 스냅한 관경이 붙고, BOX암거·세월교면 None이다.
|
||
recommended_facility: str = "pipe"
|
||
recommended_diameter_mm: int | None = None
|
||
|
||
@property
|
||
def boundary_xy(self) -> list[tuple[float, float]]:
|
||
"""가장 넓은 조각의 외곽 링 — 링 하나만 받는 옛 소비처를 위한 자리."""
|
||
return self.boundary_parts[0][0] if self.boundary_parts else []
|
||
|
||
@property
|
||
def boundary_rings(self) -> list[list[tuple[float, float]]]:
|
||
"""조각 구분 없이 편 링 목록 — 캔버스는 even-odd로 한 번에 채운다."""
|
||
return [ring for part in self.boundary_parts for ring in part]
|
||
|
||
|
||
@dataclass
|
||
class RoadRouting:
|
||
"""B04가 남긴 배수유역 분석 결과 — 세부유역을 나누는 데 필요한 최소 묶음."""
|
||
|
||
spec: GridSpec
|
||
# (R*C,) int32 — 셀이 물길을 따라 도달하는 도로 셀 슬롯(−1 = 미도달).
|
||
road_slot: np.ndarray
|
||
path_length: np.ndarray # (R*C,) float32 — 그 도로 셀까지 물길 길이(m)
|
||
elevation: np.ndarray # (R*C,) float32 — 셀 표고(유역 낙차 계산용)
|
||
road_cell_index: np.ndarray # (K,) int32 — 도로 셀의 평탄 인덱스
|
||
road_chainage: np.ndarray # (K,) float64 — 도로 셀의 누가거리(m)
|
||
strength: np.ndarray # (K,) int64 — 도로 셀별 상류 셀 수
|
||
# 화면에 그대로 그릴 기하(WGS84 lon/lat).
|
||
route_lonlat: list[list[float]] = field(default_factory=list)
|
||
basin_lonlat: list[list[float]] = field(default_factory=list)
|
||
base_pipes: list[StructureCandidate] = field(default_factory=list)
|
||
# 평균 흐름 화살표 — [x, y, 방위(도), 도로도달, 셀 수]. B04가 계산해 둔 그대로.
|
||
flow_arrows: list[list[Any]] = field(default_factory=list)
|
||
arrow_spacing_m: float = 0.0
|
||
|
||
@property
|
||
def strength_curve(self) -> np.ndarray:
|
||
"""누가거리 1m 구간별 유입 면적(㎡) 곡선 — 관 보충 위치 점수의 근거."""
|
||
if self.road_chainage.size == 0:
|
||
return np.zeros(1)
|
||
bins = max(1, int(np.ceil(self.road_chainage.max())) + 1)
|
||
index = np.clip(np.round(self.road_chainage).astype(np.int64), 0, bins - 1)
|
||
weights = self.strength.astype(np.float64) * self.spec.cell_area_m2
|
||
return np.bincount(index, weights=weights, minlength=bins)
|
||
|
||
|
||
def build_detail(
|
||
directory: Path,
|
||
vertices: list[RouteVertex],
|
||
confirmed_chainages: list[float] | None = None,
|
||
) -> DrainageDetail | None:
|
||
"""B04 분석 결과를 읽어 관을 보충하고 세부유역을 나눈다.
|
||
|
||
`confirmed_chainages`를 주면 그 위치를 관으로 확정하고(사용자 편집), 비우면 B04의
|
||
기본 관에 최대 간격 규칙으로 최소 개수만 보충한다. 어느 쪽이든 격자 해석은 하지 않는다.
|
||
"""
|
||
routing = read_road_routing(directory)
|
||
if routing is None or len(vertices) < 2:
|
||
return None
|
||
|
||
if confirmed_chainages:
|
||
pipes = pipes_from_chainages(vertices, confirmed_chainages)
|
||
else:
|
||
# 저장분의 기본 관은 누가거리만 신뢰한다 — 좌표는 현재 노선 위로 다시 찍는다.
|
||
base = [
|
||
StructureCandidate(
|
||
chainage_m=pipe.chainage_m,
|
||
x=interpolate_vertex(vertices, pipe.chainage_m)[0],
|
||
y=interpolate_vertex(vertices, pipe.chainage_m)[1],
|
||
reason=pipe.reason,
|
||
)
|
||
for pipe in routing.base_pipes
|
||
]
|
||
pipes = place_pipes(vertices, base, routing.strength_curve)
|
||
|
||
detail = DrainageDetail(
|
||
route_lonlat=routing.route_lonlat,
|
||
basin_lonlat=routing.basin_lonlat,
|
||
pipes=pipes,
|
||
grid_cell_m=routing.spec.cell_m,
|
||
flow_arrows=routing.flow_arrows,
|
||
arrow_spacing_m=routing.arrow_spacing_m,
|
||
upstream_lonlat=read_upstream_lines(directory),
|
||
strength_profile=build_strength_profile(routing),
|
||
inflow_hotspots=[
|
||
[chainage, area, float(zone), float(rank)]
|
||
for chainage, area, zone, rank in find_inflow_hotspots(
|
||
routing.strength_curve,
|
||
[pipe.chainage_m for pipe in routing.base_pipes],
|
||
vertices[-1].chainage_m,
|
||
)
|
||
],
|
||
)
|
||
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, load_rainfall_idf(directory))
|
||
logger.info(
|
||
"배수유역: 세부 설계 — 관 %d개(기본 %d + 보충 %d), 세부유역 %d개",
|
||
len(pipes),
|
||
sum(1 for pipe in pipes if pipe.reason != "spacing"),
|
||
sum(1 for pipe in pipes if pipe.reason == "spacing"),
|
||
len(detail.basins),
|
||
)
|
||
return detail
|
||
|
||
|
||
def read_road_routing(directory: Path) -> RoadRouting | None:
|
||
"""`03_road_routing` 산출물을 읽는다. 없으면 None."""
|
||
prefix = STAGES["road_routing"]
|
||
array_path = directory / f"{prefix}_road_routing.npz"
|
||
if not array_path.exists():
|
||
logger.warning("배수유역: B04 분석 결과가 없습니다 (%s).", array_path)
|
||
return None
|
||
try:
|
||
with np.load(array_path, allow_pickle=False) as data:
|
||
spec = GridSpec(
|
||
x_min=float(data["x_min"]),
|
||
y_max=float(data["y_max"]),
|
||
cell_m=float(data["cell_m"]),
|
||
n_rows=int(data["n_rows"]),
|
||
n_cols=int(data["n_cols"]),
|
||
)
|
||
routing = RoadRouting(
|
||
spec=spec,
|
||
road_slot=data["road_slot"].reshape(-1),
|
||
path_length=data["path_length"].reshape(-1),
|
||
elevation=data["elevation"].reshape(-1),
|
||
road_cell_index=data["road_cell_index"],
|
||
road_chainage=data["road_chainage"],
|
||
strength=data["strength"],
|
||
)
|
||
except (OSError, KeyError, ValueError):
|
||
logger.warning("배수유역: B04 분석 결과를 읽지 못했습니다 (%s).", array_path)
|
||
return None
|
||
|
||
_read_geometry(directory / f"{prefix}_road_routing.geojson", routing)
|
||
logger.info(
|
||
"배수유역: B04 결과 로드 — 격자 %d×%d, 도로 셀 %d, 기본 관 %d",
|
||
spec.n_rows,
|
||
spec.n_cols,
|
||
routing.road_cell_index.size,
|
||
len(routing.base_pipes),
|
||
)
|
||
return routing
|
||
|
||
|
||
def _read_geometry(path: Path, routing: RoadRouting) -> None:
|
||
"""계획도로선·2차 유역 외곽선·기본 관을 GeoJSON에서 읽어 채운다."""
|
||
if not path.exists():
|
||
logger.warning("배수유역: B04 기하 산출물이 없습니다 (%s).", path)
|
||
return
|
||
try:
|
||
with path.open("r", encoding="utf-8") as file:
|
||
document = json.load(file)
|
||
except (OSError, json.JSONDecodeError):
|
||
logger.warning("배수유역: B04 기하 산출물을 읽지 못했습니다 (%s).", path)
|
||
return
|
||
routing.arrow_spacing_m = float(
|
||
(document.get("properties") or {}).get("arrow_spacing_m") or 0.0
|
||
)
|
||
for feature in document.get("features", []):
|
||
properties = feature.get("properties") or {}
|
||
geometry = feature.get("geometry") or {}
|
||
coordinates = geometry.get("coordinates")
|
||
kind = properties.get("kind")
|
||
if kind == "route" and geometry.get("type") == "LineString":
|
||
routing.route_lonlat = coordinates
|
||
elif kind == "basin_boundary" and geometry.get("type") == "Polygon" and coordinates:
|
||
routing.basin_lonlat = coordinates[0]
|
||
elif kind == "flow_arrow" and geometry.get("type") == "Point":
|
||
# 화면이 미터로 그리므로 속성의 x·y를 쓴다(기하는 저장 규약상 lon/lat).
|
||
routing.flow_arrows.append(
|
||
[
|
||
float(properties.get("x") or 0.0),
|
||
float(properties.get("y") or 0.0),
|
||
float(properties.get("azimuth_deg") or 0.0),
|
||
bool(properties.get("reaches_road")),
|
||
int(properties.get("cells") or 0),
|
||
]
|
||
)
|
||
elif kind == "pipe" and geometry.get("type") == "Point":
|
||
routing.base_pipes.append(
|
||
StructureCandidate(
|
||
chainage_m=float(properties.get("chainage_m") or 0.0),
|
||
x=0.0,
|
||
y=0.0,
|
||
reason=str(properties.get("reason") or "stream"),
|
||
)
|
||
)
|
||
|
||
|
||
def read_upstream_lines(directory: Path) -> list[list[list[float]]]:
|
||
"""`01_primary_region`에서 상류 세류망만 읽는다(화면 강조용).
|
||
|
||
유역 판정의 기준선이라 B04 오버레이에서도 같은 선을 굵게 그린다 — B05는 그 선을
|
||
그대로 받아 표시만 한다.
|
||
"""
|
||
path = directory / f"{STAGES['primary_region']}_primary_region.geojson"
|
||
if not path.exists():
|
||
return []
|
||
try:
|
||
with path.open("r", encoding="utf-8") as file:
|
||
document = json.load(file)
|
||
except (OSError, json.JSONDecodeError):
|
||
logger.warning("배수유역: 상류 세류망을 읽지 못했습니다 (%s).", path)
|
||
return []
|
||
lines: list[list[list[float]]] = []
|
||
for feature in document.get("features", []):
|
||
properties = feature.get("properties") or {}
|
||
geometry = feature.get("geometry") or {}
|
||
if properties.get("kind") != "upstream":
|
||
continue
|
||
coordinates = geometry.get("coordinates")
|
||
if geometry.get("type") == "LineString" and coordinates:
|
||
lines.append(coordinates)
|
||
elif geometry.get("type") == "MultiLineString" and coordinates:
|
||
lines.extend(part for part in coordinates if part)
|
||
return lines
|
||
|
||
|
||
def build_strength_profile(routing: RoadRouting) -> list[list[float]]:
|
||
"""도로 1m 구간별 유입 면적(㎡) 곡선을 응답용으로 정리한다.
|
||
|
||
값이 0인 구간은 빼고 보낸다 — 노선이 길면 대부분이 0이라 그대로 보내면 응답만 커진다.
|
||
화면은 받은 구간만 색칠하고 나머지는 계획선 원래 색을 남긴다.
|
||
"""
|
||
curve = routing.strength_curve
|
||
return [[float(index), float(value)] for index, value in enumerate(curve.tolist()) if value > 0]
|
||
|
||
|
||
# ── ⑨ 관 최소 개수 보충 ─────────────────────────────────────────────────────
|
||
|
||
|
||
def place_pipes(
|
||
vertices: list[RouteVertex],
|
||
base_pipes: list[StructureCandidate],
|
||
strength_curve: np.ndarray,
|
||
) -> list[StructureCandidate]:
|
||
"""B04가 정한 기본 관(세류 교차점)에, 최대 간격을 넘는 구간만 최소 개수로 보충한다.
|
||
|
||
기본 관은 여기서 다시 찾지 않는다 — B04 산출물에 이미 들어 있다.
|
||
"""
|
||
total_length = vertices[-1].chainage_m
|
||
base: list[StructureCandidate] = []
|
||
for candidate in sorted(base_pipes, key=lambda item: item.chainage_m):
|
||
if base and candidate.chainage_m - base[-1].chainage_m < DRAINAGE_PIPE_MIN_SPACING_M:
|
||
continue
|
||
base.append(candidate)
|
||
|
||
filled: list[StructureCandidate] = []
|
||
previous = 0.0
|
||
for candidate in [*base, None]:
|
||
boundary = candidate.chainage_m if candidate else total_length
|
||
filled.extend(_fill_gap(vertices, strength_curve, previous, boundary))
|
||
if candidate:
|
||
filled.append(candidate)
|
||
previous = candidate.chainage_m
|
||
else:
|
||
previous = boundary
|
||
filled.sort(key=lambda item: item.chainage_m)
|
||
return filled
|
||
|
||
|
||
def _fill_gap(
|
||
vertices: list[RouteVertex],
|
||
strength_curve: np.ndarray,
|
||
start_m: float,
|
||
end_m: float,
|
||
) -> list[StructureCandidate]:
|
||
"""[start, end] 구간에 최대 간격을 지키는 **최소 개수**의 관을 배치한다.
|
||
|
||
필요 개수 n은 구간 길이로 정해지고(ceil(L/max) − 1), 각 관은 등분 위치를 중심으로
|
||
허용 여유(slack) 안에서만 움직인다. 그래서 개수는 늘지 않으면서도 흐름 강도가 크고
|
||
종단이 낮은 지점으로 붙는다.
|
||
"""
|
||
span = end_m - start_m
|
||
if span <= DRAINAGE_PIPE_MAX_SPACING_M:
|
||
return []
|
||
count = int(np.ceil(span / DRAINAGE_PIPE_MAX_SPACING_M)) - 1
|
||
if count <= 0:
|
||
return []
|
||
spacing = span / (count + 1)
|
||
slack = max(0.0, (DRAINAGE_PIPE_MAX_SPACING_M - spacing) / 2.0)
|
||
placed: list[StructureCandidate] = []
|
||
for order in range(1, count + 1):
|
||
nominal = start_m + spacing * order
|
||
low = max(start_m + DRAINAGE_PIPE_MIN_SPACING_M, nominal - slack)
|
||
high = min(end_m - DRAINAGE_PIPE_MIN_SPACING_M, nominal + slack)
|
||
chosen = _best_position(vertices, strength_curve, low, high, nominal)
|
||
x, y, _ = interpolate_vertex(vertices, chosen)
|
||
placed.append(StructureCandidate(chainage_m=chosen, x=x, y=y, reason="spacing"))
|
||
return placed
|
||
|
||
|
||
def _best_position(
|
||
vertices: list[RouteVertex],
|
||
strength_curve: np.ndarray,
|
||
low_m: float,
|
||
high_m: float,
|
||
fallback_m: float,
|
||
) -> float:
|
||
"""허용 구간 안에서 흐름 강도가 크고 종단이 낮은 위치를 고른다."""
|
||
if high_m <= low_m:
|
||
return fallback_m
|
||
positions = np.arange(low_m, high_m + 1.0, 1.0)
|
||
if positions.size == 0:
|
||
return fallback_m
|
||
index = np.clip(np.round(positions).astype(np.int64), 0, strength_curve.size - 1)
|
||
strength = strength_curve[index]
|
||
heights = np.array([interpolate_vertex(vertices, float(p))[2] for p in positions])
|
||
|
||
strength_score = strength / strength.max() if strength.max() > 0 else np.zeros_like(strength)
|
||
height_span = float(heights.max() - heights.min())
|
||
sag_score = (
|
||
(heights.max() - heights) / height_span if height_span > 1e-6 else np.zeros_like(heights)
|
||
)
|
||
score = _SCORE_WEIGHT_STRENGTH * strength_score + _SCORE_WEIGHT_SAG * sag_score
|
||
for order, position in enumerate(positions):
|
||
if not is_uphill_at(vertices, float(position)):
|
||
score[order] *= _SCORE_FILL_PENALTY
|
||
return float(positions[int(np.argmax(score))])
|
||
|
||
|
||
def pipes_from_chainages(
|
||
vertices: list[RouteVertex], chainages: list[float]
|
||
) -> list[StructureCandidate]:
|
||
"""사용자가 확정·편집한 누가거리 목록을 관 후보로 되돌린다.
|
||
|
||
노선 밖 값은 시·종점으로 당긴다. 그대로 두면 마커는 끝점에 찍히는데 라벨만 −50m처럼
|
||
나와 좌표와 표기가 어긋난다.
|
||
"""
|
||
total_length = vertices[-1].chainage_m
|
||
clamped = {min(max(round(float(item), 2), 0.0), total_length) for item in chainages}
|
||
pipes: list[StructureCandidate] = []
|
||
for value in sorted(clamped):
|
||
x, y, _ = interpolate_vertex(vertices, value)
|
||
pipes.append(StructureCandidate(chainage_m=value, x=x, y=y, reason="confirmed"))
|
||
return pipes
|
||
|
||
|
||
# ── ⑩ 측구 흐름으로 도로 셀 → 담당 관 ───────────────────────────────────────
|
||
|
||
|
||
def assign_road_cells_to_pipes(
|
||
vertices: list[RouteVertex],
|
||
pipes: list[StructureCandidate],
|
||
road_chainage: np.ndarray,
|
||
) -> np.ndarray:
|
||
"""도로 셀마다 물이 실제로 흘러가는 담당 관 번호를 정한다.
|
||
|
||
노면 물은 측구를 타고 종단 내리막으로 흐르므로 종단 계획선을 1차원 지형으로 본다.
|
||
1차원에서는 물이 **마루(구간 최고점)를 넘지 못한다** — 이웃한 두 관 사이의 최고점이
|
||
곧 분수령이고, 그 왼쪽은 앞 관이, 오른쪽은 뒤 관이 받는다. 첫 관 앞과 마지막 관 뒤는
|
||
그 관이 받는다.
|
||
|
||
옛 방식(한 칸 이웃만 보는 국소 하강 + 관 없는 저점은 최근접 관)은 계획고의 미세
|
||
요철에 걸려 멈췄다. 용화 실측: 측점 2,138개 중 1,898개(88.8%)가 저점에 갇혀 흐름이
|
||
아니라 **누가거리 최근접**으로 배정됐고, 그 결과 도로 셀 43.2%가 자기보다 높은 관에
|
||
배정됐다(최대 6.82m 오르막). 마루 기준은 미세 요철을 타지 않으므로 오르막 배정이
|
||
구조적으로 생기지 않는다(2026-09-03 사용자 확정).
|
||
"""
|
||
total_length = vertices[-1].chainage_m
|
||
step = max(DRAINAGE_DITCH_SAMPLE_M, 0.5)
|
||
stations = np.arange(0.0, total_length + step, step)
|
||
pipe_chainages = np.array([pipe.chainage_m for pipe in pipes])
|
||
slot_station = np.clip(np.round(road_chainage / step).astype(np.int64), 0, stations.size - 1)
|
||
if pipe_chainages.size == 0:
|
||
return np.full(road_chainage.size, -1, dtype=np.int32)
|
||
|
||
heights = np.array([interpolate_vertex(vertices, float(s))[2] for s in stations])
|
||
# 관 순서는 호출자가 준 그대로 돌려줘야 한다 — 누가거리로 정렬해 풀고 끝에 되돌린다.
|
||
order = np.argsort(pipe_chainages, kind="stable")
|
||
pipe_station = np.clip(
|
||
np.round(pipe_chainages[order] / step).astype(np.int64), 0, stations.size - 1
|
||
)
|
||
|
||
owner = np.full(stations.size, pipe_station.size - 1, dtype=np.int64) # 마지막 관 뒤
|
||
owner[: pipe_station[0] + 1] = 0 # 첫 관 앞
|
||
for index in range(pipe_station.size - 1):
|
||
left = pipe_station[index]
|
||
right = pipe_station[index + 1]
|
||
if right <= left:
|
||
continue
|
||
ridge = left + int(np.argmax(heights[left : right + 1]))
|
||
owner[left : ridge + 1] = index
|
||
owner[ridge + 1 : right + 1] = index + 1
|
||
return order[owner][slot_station].astype(np.int32)
|
||
|
||
|
||
# ── ⑩ 세부유역 조립 ────────────────────────────────────────────────────────
|
||
|
||
|
||
def assemble_basins(
|
||
routing: RoadRouting,
|
||
pipes: list[StructureCandidate],
|
||
pipe_of_slot: np.ndarray,
|
||
idf: dict[str, float] | None = None,
|
||
) -> list[WatershedBasin]:
|
||
"""셀이 도달한 도로 셀의 담당 관을 그대로 유역 번호로 삼아 세부유역을 만든다."""
|
||
spec = routing.spec
|
||
labels = np.full(spec.size, -1, dtype=np.int32)
|
||
reached = routing.road_slot >= 0
|
||
labels[reached] = pipe_of_slot[routing.road_slot[reached]]
|
||
|
||
# 최소면적 필터를 끈다 — 그 필터가 곧 빈공간이었다. 떨어진 조각을 100㎡ 미만이라고
|
||
# 버리면 유역 면적과 그림이 어긋난다(실측: 용화 5.76%→1.86%, S자 3.56%→1.42%,
|
||
# 조각 유역 0→2·0→3 복원). 링 목록이 조각을 싣게 된 뒤로는 버릴 이유가 없고, 좌표점은
|
||
# 583→622·806→829로 거의 늘지 않는다(2026-09-03). 남은 오차는 simplify(2.0m) 몫.
|
||
polygons = polygonize_labels(spec, labels, min_area_m2=0.0)
|
||
cell_area = spec.cell_area_m2
|
||
basins: list[WatershedBasin] = []
|
||
for order, pipe in enumerate(pipes):
|
||
member = labels == order
|
||
count = int(member.sum())
|
||
if count == 0:
|
||
continue
|
||
geometry = polygons.get(order)
|
||
elevations = routing.elevation[member]
|
||
highest = float(np.nanmax(elevations)) if np.isfinite(elevations).any() else 0.0
|
||
outlet_z = _outlet_elevation(routing, order, pipe_of_slot)
|
||
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)
|
||
facility, recommended = recommend_structure(sizing["diameter_mm"] if sizing else None)
|
||
basins.append(
|
||
WatershedBasin(
|
||
index=len(basins) + 1,
|
||
chainage_m=pipe.chainage_m,
|
||
outlet_x=pipe.x,
|
||
outlet_y=pipe.y,
|
||
boundary_parts=polygon_parts(geometry) if geometry is not None else [],
|
||
area_m2=area,
|
||
relief_m=relief,
|
||
flow_length_m=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,
|
||
bridge_required=bool(
|
||
sizing and sizing["diameter_mm"] > DRAINAGE_BRIDGE_THRESHOLD_MM
|
||
),
|
||
required_area_m2=sizing["required_area_m2"] if sizing else None,
|
||
recommended_facility=facility,
|
||
recommended_diameter_mm=recommended,
|
||
)
|
||
)
|
||
return basins
|
||
|
||
|
||
def _outlet_elevation(routing: RoadRouting, pipe_order: int, pipe_of_slot: np.ndarray) -> float:
|
||
"""관이 담당하는 도로 셀들의 최저 표고 = 유역 출구 표고."""
|
||
slots = np.flatnonzero(pipe_of_slot == pipe_order)
|
||
if slots.size == 0:
|
||
return 0.0
|
||
elevations = routing.elevation[routing.road_cell_index[slots]]
|
||
finite = elevations[np.isfinite(elevations)]
|
||
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),
|
||
# 물넘이·세월교는 관이 아니라 개수로라 직경 대신 이 단면적에서 출발한다.
|
||
"required_area_m2": round(required_area, 4),
|
||
}
|
||
|
||
|
||
def recommend_structure(diameter_mm: float | None) -> tuple[str, int | None]:
|
||
"""유효직경으로 추천 구조물과 관경을 고른다 (2026-08-17 사용자 확정: 유량 근거만).
|
||
|
||
지식DB 선정 트리(개거_세월시설 §5·횡단배수관_암거 §5)에서 유량 조건만 취했다 —
|
||
계곡 횡단경사·하천 차수는 지형 계산이 필요해 화면이 "현장 확인"으로 안내한다.
|
||
|
||
D ≤ 1,500㎜ 배관 (레지스트리 선택지로 스냅, 하한 800㎜)
|
||
1,500 < D ≤ 2,000 BOX암거 후보
|
||
D > 2,000㎜ 세월교·물넘이 검토 (관 최대 규격 초과)
|
||
|
||
유효직경을 못 구한 유역(강우량표 없음)은 배관·관경 미정으로 둔다.
|
||
"""
|
||
if diameter_mm is None:
|
||
return PIPE_FACILITY_PIPE, None
|
||
if diameter_mm > DRAINAGE_BRIDGE_THRESHOLD_MM:
|
||
return PIPE_FACILITY_FORD_BRIDGE, None
|
||
if diameter_mm > DRAINAGE_BOX_THRESHOLD_MM:
|
||
return PIPE_FACILITY_BOX, None
|
||
# 별표2 (나) 예외 하한 800㎜ — 계산값이 더 작아도 그 아래로는 내리지 않는다.
|
||
need = max(float(PIPE_DEFAULT_DIAMETER_MM), diameter_mm)
|
||
size = next(
|
||
(value for value in DRAINAGE_RECOMMEND_DIAMETERS_MM if value >= need),
|
||
DRAINAGE_RECOMMEND_DIAMETERS_MM[-1],
|
||
)
|
||
return PIPE_FACILITY_PIPE, int(size)
|
||
|
||
|
||
def estimate_pipe_diameter_mm(
|
||
area_m2: float,
|
||
relief_m: float,
|
||
flow_length_m: float,
|
||
idf: dict[str, float] | None = None,
|
||
) -> float | None:
|
||
"""유역 제원으로 배수 유효직경(mm)만 돌려주는 축약형 — 상세는 size_pipe()."""
|
||
sizing = size_pipe(area_m2, relief_m, flow_length_m, idf)
|
||
return sizing["diameter_mm"] if sizing else None
|