"""배수유역 분석 오케스트레이터 (B04 — 관리자 확인용 전처리). 계획 노선(B03 업로드 파일)과 도엽 등고선·세류선만으로 배수유역을 끝까지 분석해 영구저장소에 남긴다. 30초 안팎이 걸리는 무거운 작업이라 여기서 한 번만 돌리고, 일반 사용자가 쓰는 B05는 그 결과를 읽어 쓰기만 한다(2026-07-31 사용자 지시). ① 도로 교차 세류망 상류측 추출 → 반경 버퍼 = 1차 배수유역 ② 도로 시작점 기준 격자 생성 (1차 영역에 걸치는 셀만) ③ 등고선 하강 방향 — 높은 등고 라인에서 낮은 등고 라인으로. 보간면을 쓰지 않으므로 가짜 웅덩이·평탄면이 원리적으로 생기지 않는다 ④ 세류망 흐름 새김 → 최외곽부터 사슬 추적 → 도로 도달 여부(적/청) 판정 ⑤ 최외곽 적색 셀 주변 확장 — 새로 추가한 셀에 적색이 없을 때까지 ⑥ 도로 셀별 흐름 강도 ⑦ 2차 전체 배수유역 외곽선 ⑧ 기본 관 매설 위치 (도로 × 세류선 교차점) ⑨ 이후(관 최소 개수 보충, 세부유역 분할)는 사용자가 관을 옮길 수 있어야 하므로 B05에 남긴다. """ from __future__ import annotations import logging import math import time from dataclasses import dataclass, field from typing import Any import numpy as np from shapely.geometry import LineString from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import ContourDescent from B04_PreProcess.B04_PreProcess_Engine_Watershed_Expand import expand_by_red_boundary from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import ( FlowClassification, RoadRaster, largest_ring, outer_boundary, trace_flow, ) from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( AZIMUTH_STEPS, GridSpec, TerrainGrid, build_contour_cloud, route_elevation_floor, ) from B04_PreProcess.B04_PreProcess_Engine_Watershed_Stream import ( PrimaryRegion, build_primary_region, ) from common_util.common_util_route_geometry import ( RouteVertex, StructureCandidate, find_stream_crossings, ) from config.config_system import ( DRAINAGE_ARROW_BLOCK_M, DRAINAGE_ARROW_MIN_AGREEMENT, DRAINAGE_ARROW_MIN_COVERAGE, DRAINAGE_ARROW_SPACING_M, DRAINAGE_GRID_SIZE_M, DRAINAGE_INITIAL_RADIUS_M, DRAINAGE_PIPE_MAX_SPACING_M, DRAINAGE_PIPE_MIN_SPACING_M, ) logger = logging.getLogger(__name__) # 강도 곡선 응답 간격(m). 도로 위 흐름 강도 표기는 이 간격으로 내보낸다. # 1m = 계산 원본 그대로다. 5m로 줄이면 도로 색이 뭉개져 어느 자리에 물이 모이는지 못 읽는다 # (2026-08-01 사용자 지시: 단순화·평균 금지). 350m 노선이면 약 351점, 수 KB 수준이라 부담 없다. _STRENGTH_OUTPUT_STEP_M = 1.0 # 구역마다 뽑을 유입 집중점 개수에 더하는 여유 몫. # 개수 = floor(구역 길이 / 관 최대간격) + 이 값. 최대간격을 지키는 데 필요한 자리 수만 뽑으면 # 견줄 대상이 없어 그 자리가 최선인지 판단할 수 없다(2026-08-01 사용자 지시). _HOTSPOT_EXTRA_PER_ZONE = 3 # ── ①~② 1차 배수유역 ──────────────────────────────────────────────────────── def resolve_primary_region( vertices: list[RouteVertex], route_line: LineString, contour_features: list[dict[str, Any]], stream_features: list[dict[str, Any]], ) -> PrimaryRegion | None: """도로 교차 세류선(상류측)과 노선을 반경 버퍼한 1차 배수유역과 격자 범위를 정한다. 상·하류 판정에 쓸 등고선은 노선 주변만 있으면 된다(교차점이 전부 노선 위이므로). 도엽 전체를 읽으면 이 단계에서만 수십 초가 날아간다. """ floor = route_elevation_floor([vertex.z for vertex in vertices]) near_bounds = route_line.buffer(DRAINAGE_INITIAL_RADIUS_M * 2.0).bounds cloud = build_contour_cloud(contour_features, floor, near_bounds) if cloud.is_empty: logger.warning("배수유역: 노선 주변에 등고선이 없어 1차 영역을 정할 수 없습니다.") return None return build_primary_region( route_line, stream_features, cloud, DRAINAGE_INITIAL_RADIUS_M, DRAINAGE_GRID_SIZE_M ) def preview_primary_region( vertices: list[RouteVertex], contour_features: list[dict[str, Any]], stream_features: list[dict[str, Any]], ) -> PrimaryRegion | None: """단계 검증용 — TIN·흐름 계산 없이 1차 배수유역 근거만 뽑는다.""" if len(vertices) < 2: return None route_line = LineString([(vertex.x, vertex.y) for vertex in vertices]) return resolve_primary_region(vertices, route_line, contour_features, stream_features) @dataclass class StagePreview: """단계 검증 산출물 묶음. 기능을 붙일 때마다 여기에 항목이 하나씩 늘어난다. 확장을 거치면 격자와 해석 영역이 1차 영역보다 커진다. 화면·저장은 `region.spec`이 아니라 여기 `spec`/`domain`을 봐야 한다. """ region: PrimaryRegion spec: GridSpec | None = None domain: np.ndarray | None = None terrain: TerrainGrid | None = None road: RoadRaster | None = None flow: FlowClassification | None = None descent: ContourDescent | None = None expand_rounds: int = 0 expand_closed: bool = False expand_added_cells: int = 0 # ⑥ 도로 위 흐름 강도 — (누가거리 m, 그 구간으로 모이는 상류 면적 ㎡). strength_profile: list[tuple[float, float]] = field(default_factory=list) # ⑥-1 유입 집중점 — (누가거리 m, 유입면적 ㎡, 구역 번호, 구역 내 순위). 화면 마커용. inflow_hotspots: list[tuple[float, float, int, int]] = field(default_factory=list) # ⑦ 2차 전체 배수유역 외곽선(= 분수령). 적색 셀 전체를 폴리곤화한 것. basin_boundary_xy: list[tuple[float, float]] = field(default_factory=list) basin_area_m2: float = 0.0 # 셀 → 도로 셀 귀속. B05가 세부유역을 나눌 때 이 배열이 있어야 한다. routing: Any = None # B05용 평균 흐름 화살표 — (x, y, 방위 라디안, 도로 도달, 셀 수). flow_arrows: list[tuple[float, float, float, bool, int]] = field(default_factory=list) # ⑧ 기본 관 매설 위치 — 도로 × 세류선 교차점. pipes: list[StructureCandidate] = field(default_factory=list) def preview_stages( vertices: list[RouteVertex], contour_features: list[dict[str, Any]], stream_features: list[dict[str, Any]], ) -> StagePreview | None: """지금까지 구현·검증된 단계를 순서대로 돌려 결과를 모은다. 현재 포함: ① 1차 배수유역 ② 격자 생성 ③ **등고선 하강 방향** ④ 도로 도달 판정 ⑤ **최외곽 적색 셀 주변 확장**. ③은 보간면(TIN)을 쓰지 않는다. 높은 등고 라인에서 낮은 등고 라인으로 내려가며 방향을 세우므로 가짜 웅덩이·평탄면이 원리적으로 생기지 않는다(2026-07-31 사용자 지시로 방식 교체). ⑤는 최외곽에 적색이 남아 있으면 그 주변으로 넓혀 다시 분석하고, **새로 추가한 셀에 적색이 없으면** 멈춘다. """ if len(vertices) < 2: return None route_line = LineString([(vertex.x, vertex.y) for vertex in vertices]) region = resolve_primary_region(vertices, route_line, contour_features, stream_features) if region is None: return None started = time.perf_counter() floor = route_elevation_floor([vertex.z for vertex in vertices]) expansion = expand_by_red_boundary( region.spec, region.cell_mask, contour_features, route_line, region.split.upstream, floor, ) if expansion is None: logger.warning("배수유역: 등고선 하강 방향을 세우지 못해 흐름 판정을 건너뜁니다.") return StagePreview(region=region) analysis = expansion.analysis spec = analysis.spec red = analysis.flow.reaches_road & analysis.flow.analyzed # ⑥ 흐름 강도 — 셀마다 물이 실제로 들어가는 도로 셀을 구해 도로 셀별로 센다. # 색 판정은 세류 셀에서 멈추지만(거기서 도달이 확정되므로), 강도는 그 물이 세류를 타고 # 최종적으로 어느 도로 셀로 들어가는지를 봐야 하므로 도로만 흡수점으로 두고 다시 따라간다. routing = trace_flow(analysis.terrain, analysis.road) if analysis.road.count else None strength_curve = _preview_strength(analysis, routing, red, route_line.length) # ⑦ 2차 전체 배수유역 외곽선 = 적색 셀 전체의 외곽. boundary = outer_boundary(spec, red.reshape(spec.n_rows, spec.n_cols)) basin_ring = largest_ring(boundary) if boundary is not None else [] # ⑧ 기본 관 매설 위치 = 도로 × 세류선 교차점(너무 가까운 것은 하나로 합침). pipes = _base_pipes(vertices, stream_features) # ⑥-1 유입 집중점 — 기본 관으로 나눈 구역마다 물이 많이 모이는 자리를 뽑는다. hotspots = find_inflow_hotspots( strength_curve, [pipe.chainage_m for pipe in pipes], route_line.length ) # B05에 얹을 평균 흐름 화살표 — 셀 화살표는 도면 배율에서 안 보인다. flow_arrows = build_flow_arrows(analysis, analysis.flow) logger.info( "배수유역: 단계 분석 %.1fs (확장 %d회, 최종 셀 %d개) — " "2차 유역 %.0f㎡, 기본 관 %d개, 강도 곡선 %d점", time.perf_counter() - started, expansion.rounds, spec.size, int(red.sum()) * spec.cell_area_m2, len(pipes), int((strength_curve > 0).sum()), ) return StagePreview( region=region, spec=spec, domain=analysis.domain, terrain=analysis.terrain, road=analysis.road, flow=analysis.flow, descent=analysis.descent, expand_rounds=expansion.rounds, expand_closed=expansion.closed, expand_added_cells=expansion.added_cells, strength_profile=_downsample_strength(strength_curve), inflow_hotspots=hotspots, basin_boundary_xy=basin_ring, basin_area_m2=int(red.sum()) * spec.cell_area_m2, routing=routing, pipes=pipes, flow_arrows=flow_arrows, ) def _preview_strength( analysis: Any, routing: Any, red: np.ndarray, route_length_m: float ) -> np.ndarray: """적색 셀이 실제로 들어가는 도로 셀을 세어 누가거리별 유입 면적 곡선을 만든다.""" road = analysis.road if routing is None or road.count == 0: return np.zeros(1) slots = routing.road_slot counted = red & (slots >= 0) strength = np.bincount(slots[counted], minlength=road.count).astype(np.float64) return _strength_by_chainage( road.chainage, strength * analysis.spec.cell_area_m2, route_length_m ) def build_flow_arrows(analysis: Any, flow: Any) -> list[tuple[float, float, float, bool, int]]: """셀 흐름을 블록 단위로 평균해 B05에 얹을 화살표를 뽑는다. 셀 화살표는 1m라 도면 배율에서 경향이 안 보인다. 겹치지 않는 블록으로 나눠 방향을 평균하고, 화살표끼리 최소 간격을 두어 솎아낸다(2026-07-31 사용자 지시). **세류선 셀과 도로 셀은 뺀다.** 그 자리 흐름은 지형 경사가 아니라 확정된 물길·노면을 따르는 값이라 사면 경향을 왜곡한다. 방향 평균은 산술평균이 아니라 **원형 평균**으로 낸다(0°와 359°의 평균은 180°가 아니라 0°다). 평균 벡터 길이가 일치도이므로, 블록 안 방향이 제각각이면 그 블록은 버린다. """ spec = analysis.spec rows, cols = spec.n_rows, spec.n_cols direction = flow.direction.reshape(rows, cols) usable = ( analysis.domain & flow.analyzed.reshape(rows, cols) & (direction < AZIMUTH_STEPS) # 싱크·무효 제외 & ~analysis.road.mask ) if flow.burned is not None: usable &= ~flow.burned.reshape(rows, cols) if not usable.any(): return [] block = max(1, int(round(DRAINAGE_ARROW_BLOCK_M / spec.cell_m))) stride = max(1, int(round(DRAINAGE_ARROW_SPACING_M / (block * spec.cell_m)))) angle = direction.astype(np.float64) * (2.0 * math.pi / AZIMUTH_STEPS) reaches = flow.reaches_road.reshape(rows, cols) arrows: list[tuple[float, float, float, bool, int]] = [] for row0 in range(0, rows - block + 1, block * stride): for col0 in range(0, cols - block + 1, block * stride): window = usable[row0 : row0 + block, col0 : col0 + block] count = int(window.sum()) if count < DRAINAGE_ARROW_MIN_COVERAGE * block * block: continue local = angle[row0 : row0 + block, col0 : col0 + block][window] mean_x = float(np.cos(local).mean()) mean_y = float(np.sin(local).mean()) agreement = math.hypot(mean_x, mean_y) if agreement < DRAINAGE_ARROW_MIN_AGREEMENT: continue # 방향이 제각각인 블록 — 평균이 경향을 대표하지 못한다 centre_row = row0 + block / 2.0 centre_col = col0 + block / 2.0 arrows.append( ( spec.x_min + centre_col * spec.cell_m, spec.y_max - centre_row * spec.cell_m, math.atan2(mean_y, mean_x), bool(reaches[row0 : row0 + block, col0 : col0 + block][window].mean() >= 0.5), count, ) ) logger.info( "배수유역: 평균 흐름 화살표 %d개 (블록 %.0fm, 간격 %.0fm, 세류·도로 셀 제외)", len(arrows), block * spec.cell_m, block * stride * spec.cell_m, ) return arrows def _base_pipes( vertices: list[RouteVertex], stream_features: list[dict[str, Any]] ) -> list[StructureCandidate]: """도로 × 세류선 교차점을 기본 관 위치로 삼는다. 300m 보충 배치는 다음 단계다.""" pipes: list[StructureCandidate] = [] for candidate in find_stream_crossings(vertices, stream_features): if pipes and candidate.chainage_m - pipes[-1].chainage_m < DRAINAGE_PIPE_MIN_SPACING_M: continue pipes.append(candidate) return pipes # ── 흐름 강도 곡선 ────────────────────────────────────────────────────────── def _strength_by_chainage( road_chainage: np.ndarray, strength_area: np.ndarray, total_length: float ) -> np.ndarray: """도로 셀 강도를 1m 누가거리 구간으로 합산한 곡선(㎡/m 구간 합). 구간 k는 **[k, k+1)** 이다. `round`를 쓰면 안 된다 — 도로 셀 누가거리는 노선을 0.5m 간격으로 샘플해 붙인 값이라 전부 0.5의 배수인데, numpy의 `round`는 정확히 .5인 값을 짝수 쪽으로 보낸다(0.5→0, 1.5→2). 그러면 짝수 칸에는 세 샘플({k−0.5, k, k+0.5}), 홀수 칸에는 한 샘플만 들어가 강도가 2m 주기로 흔들리고, 도로가 격자축과 나란한 구간은 홀수 칸이 전부 0이 되어 화면에서 1m 간격 점선으로 보인다(2026-08-01 검증에서 발견). """ bins = max(1, int(np.ceil(total_length)) + 1) if road_chainage.size == 0: return np.zeros(bins) index = np.clip(np.floor(road_chainage).astype(np.int64), 0, bins - 1) return np.bincount(index, weights=strength_area, minlength=bins) def find_inflow_hotspots( curve: np.ndarray, pipe_chainages: list[float], route_length_m: float, ) -> list[tuple[float, float, int, int]]: """도로 위 **유입 집중점**을 뽑는다 — (누가거리 m, 유입면적 ㎡, 구역 번호, 구역 내 순위). 구역 번호 `-1`은 **기본 관(세류 교차점) 자리**를 뜻한다. 나머지는 그 구역에서 보충 후보로 뽑은 자리다. 왜 필요한가(2026-08-01 사용자 지시): 기본 관은 국가 수치지형도의 세류선이 도로를 가로지르는 자리라 근거가 확실하다. 그런데 관 **최대 간격**을 지키려면 세류가 없는 긴 구간에도 관을 넣어야 하는데, 이때 등간격으로 기계적으로 꽂는 대신 **물이 실제로 많이 모이는 자리**를 골라야 한다. 여기서 뽑은 지점이 그 후보이며, B05가 관을 옮겨도 세부유역을 다시 나누는 기준으로 쓸 수 있다. 규칙 · 구역 = 시점 → 기본 관들 → 종점 으로 자른 구간 · 구역별 개수 = `floor(구역 길이 / 관 최대간격) + 3` (관 최대간격을 지키는 데 꼭 필요한 자리 수에 검토 여유 몫을 더한 것이다. 한 자리만 보여 주면 그 자리가 실제로 최선인지 견줄 대상이 없다 — 2026-08-01 사용자 지시) · 고르는 법 = 강도가 가장 큰 1m 칸을 고르고 **그 지점 양 편측으로 관 최소간격만큼**을 후보에서 뺀 뒤 다음으로 큰 칸을 고른다. 빼지 않으면 같은 계곡의 이웃 칸들이 연달아 뽑혀 한쪽에 쏠린다. · **기본 관 쪽 경계에서도 관 최소간격만큼 물러나 본다.** 강도 최대점은 원리상 세류가 도로를 가로지르는 자리, 곧 기본 관 자신이다. 관 옆 1m만 비워 두면 그 관이 받는 물 전부(전체 유역의 대부분)가 관에서 두어 걸음 떨어진 칸에서 다시 최대값으로 잡혀, 정작 관이 없어 보충이 필요한 자리는 뽑히지 못한다. 노선 시·종점은 관이 아니므로 물러나지 않는다. · **기본 관 자리 자체도 집중점으로 넣는다(구역 번호 −1).** 세류가 도로를 가로지르는 자리가 대개 그 노선에서 가장 큰 유역을 받는데, 그 자리는 구역 경계라 위 규칙으로는 영영 뽑히지 않는다(2026-08-01 사용자 지시). 위치는 관 누가거리 그대로가 아니라 관 ±최소간격 안에서 강도가 가장 큰 칸으로 잡는다 — 도로를 폭으로 굽는 과정에서 실제 유입 봉우리가 관에서 두어 걸음 옆에 놓이기 때문이다. · 강도 0인 칸은 뽑지 않는다(개수를 못 채워도 그대로 둔다). """ if curve.size == 0 or route_length_m <= 0: return [] # 구역 경계 = 0, 기본 관들, 노선 끝. 순서·중복을 정리해 둔다. edges = sorted({0.0, float(route_length_m), *(float(p) for p in pipe_chainages)}) exclusion = max(1, int(round(DRAINAGE_PIPE_MIN_SPACING_M))) hotspots: list[tuple[float, float, int, int]] = [] # ① 기본 관 자리 — 관마다 그 근처 봉우리 하나. for pipe_m in sorted({float(p) for p in pipe_chainages}): center = int(round(pipe_m)) low = max(0, center - exclusion) high = min(curve.size - 1, center + exclusion) if high < low: continue window = curve[low : high + 1] best = int(np.argmax(window)) if window[best] <= 0: continue hotspots.append((float(low + best), float(window[best]), -1, 0)) # 관 자리와 겹치는 칸은 구역 쪽에서 다시 뽑지 않는다 — 같은 자리에 마커가 두 개 겹친다. taken = {int(chainage) for chainage, _, _, _ in hotspots} # ② 구역별 보충 후보 — 관과 관 사이에서 물이 많이 모이는 자리. for zone_index, (start_m, end_m) in enumerate(zip(edges, edges[1:])): span = end_m - start_m if span <= 0: continue wanted = int(span // DRAINAGE_PIPE_MAX_SPACING_M) + _HOTSPOT_EXTRA_PER_ZONE # 이 구역에서 볼 곡선 조각. 관 쪽 경계는 최소간격만큼, 노선 끝 쪽은 한 칸만 물린다. start_is_pipe = zone_index > 0 end_is_pipe = zone_index < len(edges) - 2 low = int(math.floor(start_m)) + (exclusion if start_is_pipe else 1) high = int(math.ceil(end_m)) - (exclusion if end_is_pipe else 1) high = min(high, curve.size - 1) if high < low: continue window = curve[low : high + 1].copy() for index in taken: if low <= index <= high: window[index - low] = 0.0 for rank in range(wanted): best = int(np.argmax(window)) if window[best] <= 0: break # 남은 칸이 전부 0 — 뽑을 것이 없다 hotspots.append((float(low + best), float(window[best]), zone_index, rank)) # 고른 지점 양 편측으로 최소 간격만큼 지운다(한쪽 쏠림 방지). window[max(0, best - exclusion) : best + exclusion + 1] = 0.0 hotspots.sort(key=lambda item: item[0]) return hotspots def _downsample_strength(curve: np.ndarray) -> list[tuple[float, float]]: """응답용으로 강도 곡선을 일정 간격으로 줄인다(구간 합 유지). 끝자락을 잘라내면 종점 부근 유입 면적이 통째로 사라지므로 0으로 채워 맞춘다. """ step = max(1, int(_STRENGTH_OUTPUT_STEP_M)) if curve.size == 0: return [] padding = (-curve.size) % step padded = np.append(curve, np.zeros(padding)) if padding else curve summed = padded.reshape(-1, step).sum(axis=1) return [ (float(position * step), float(value)) for position, value in enumerate(summed) if value > 0 ]