"""배수유역 세부 설계 공용 엔진 (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 from dataclasses import dataclass, field from pathlib import Path from typing import Any import numpy as np from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Analyze import find_inflow_hotspots from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Export import STAGES from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Flow import largest_ring, polygonize_labels from B04_wf1_Surface.B04_wf1_Surface_Engine_Watershed_Grid import GridSpec from common_util.common_util_route_geometry import ( RouteVertex, StructureCandidate, interpolate_vertex, is_uphill_at, ) from config.config_system import ( DRAINAGE_DITCH_SAMPLE_M, DRAINAGE_PIPE_MAX_SPACING_M, DRAINAGE_PIPE_MIN_SPACING_M, ) 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 boundary_xy: 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 @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) 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차원 지형으로 보고 같은 방식(내리막 추적 + 관에서 흡수)으로 푼다. 관이 없는 사그(저점)에 갇힌 구간은 가장 가까운 관이 받는 것으로 본다. """ total_length = vertices[-1].chainage_m step = max(DRAINAGE_DITCH_SAMPLE_M, 0.5) stations = np.arange(0.0, total_length + step, step) heights = np.array([interpolate_vertex(vertices, float(s))[2] for s in stations]) pipe_chainages = np.array([pipe.chainage_m for pipe in pipes]) pipe_station = np.clip(np.round(pipe_chainages / step).astype(np.int64), 0, stations.size - 1) # 앞뒤 이웃 중 더 낮은 쪽으로 흘려보낸다(양쪽 다 높으면 사그 = 제자리). back_z = np.full(stations.size, np.inf) back_z[1:] = heights[:-1] forward_z = np.full(stations.size, np.inf) forward_z[:-1] = heights[1:] go_back = (back_z < heights) & (back_z <= forward_z) go_forward = (forward_z < heights) & ~go_back receiver = np.arange(stations.size, dtype=np.int64) receiver[go_back] -= 1 receiver[go_forward] += 1 receiver[pipe_station] = pipe_station # 관은 물을 흡수한다 owner = np.full(stations.size, -1, dtype=np.int64) owner[pipe_station] = np.arange(pipe_chainages.size) jump = receiver for _ in range(40): next_jump = jump[jump] if np.array_equal(next_jump, jump): break jump = next_jump resolved = owner[jump] # 관 없는 사그에 갇힌 구간은 가장 가까운 관에 붙인다. orphan = resolved < 0 if orphan.any() and pipe_chainages.size: nearest = np.abs(stations[orphan, None] - pipe_chainages[None, :]).argmin(axis=1) resolved[orphan] = nearest slot_station = np.clip(np.round(road_chainage / step).astype(np.int64), 0, stations.size - 1) return resolved[slot_station].astype(np.int32) # ── ⑩ 세부유역 조립 ──────────────────────────────────────────────────────── def assemble_basins( routing: RoadRouting, pipes: list[StructureCandidate], pipe_of_slot: np.ndarray, ) -> 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]] polygons = polygonize_labels(spec, labels) 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()) basins.append( WatershedBasin( index=len(basins) + 1, chainage_m=pipe.chainage_m, outlet_x=pipe.x, outlet_y=pipe.y, boundary_xy=largest_ring(geometry) if geometry is not None else [], area_m2=area, relief_m=relief, flow_length_m=flow_length, pipe_diameter_mm=estimate_pipe_diameter_mm(area, relief, flow_length), ) ) 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 estimate_pipe_diameter_mm( area_m2: float, relief_m: float, flow_length_m: float, rainfall_mm_per_hour: 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