From c1f139527d0f15898d111d37d267f937d9ff0f21 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Fri, 31 Jul 2026 16:51:47 +0900 Subject: [PATCH] =?UTF-8?q?feat(B05):=20=EB=B0=B0=EC=88=98=EC=9C=A0?= =?UTF-8?q?=EC=97=AD=EC=9D=84=20=EA=B2=A9=EC=9E=90=20=ED=9D=90=EB=A6=84=20?= =?UTF-8?q?=ED=95=B4=EC=84=9D=EC=9C=BC=EB=A1=9C=20=EC=A0=84=EB=A9=B4=20?= =?UTF-8?q?=EC=9E=AC=EC=84=A4=EA=B3=84?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 등고선 아크 추적 + 능선 행진 방식이 능선/계곡을 안정적으로 분리하지 못해 D8 물 방향 + 도로 기준 상류 추적 방식으로 교체한다. - 능선을 따로 탐지하지 않는다. 물길을 따라가 도로에 닿는 셀만 유역이고, 그 경계가 곧 능선이다. 유역 내부 봉우리는 자동으로 포함된다. - 등고선 TIN 보간 후 웅덩이 채움(형태학적 재구성) + 평탄면 미세경사로 가짜 웅덩이/평탄 삼각형에서 흐름이 끊기는 문제를 없앤다. - 상류 추적은 포인터 더블링으로 전 셀을 한 번에 푼다. 셀의 흐름 종착 도로 셀(root)이 유역 판정·흐름 강도·세부유역 라벨의 공통 근거가 되어, 관을 옮겨도 격자 해석 없이 측구 라우팅만 다시 돌면 된다(.npz 캐시). - 활성 셀이 격자 최외곽에 닿은 방향으로만 확장하고, 경계 링이 전부 비활성이 되면(띠 폐합) 멈춘다. 변경 사항 - 신규 엔진 3종: Engine_Watershed_Grid / _Flow / _Basin - 폐기 엔진 4종은 _legacy_watershed/ 로 원본 보관(ruff 제외) - config_system.py §5-3-1 에 DRAINAGE_* 파라미터 18개 (격자 1m, 반경 300m) - 표고점 데이터 사용 중단(유효 데이터 부족), 프론트 능선 토글 제거 (전체 유역 외곽선과 같은 선이므로 중복) - 응답에 main_polygon_lonlat / strength_profile 추가, 계획선 위 흐름 강도 표기 합성 지형 검증: 유역 179,919㎡ vs 이론 180,000㎡ (오차 0.04%), 능선 자동 검출, 확장 3회 후 자동 정지, 캐시 재사용 1.2s -> 0.1s Co-Authored-By: Claude Fable 5 --- B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts | 8 +- .../B05_wf2_Route_Engine_Drainage.py | 6 +- .../B05_wf2_Route_Engine_Watershed_Basin.py | 562 ++++++++++++++++++ .../B05_wf2_Route_Engine_Watershed_Flow.py | 224 +++++++ .../B05_wf2_Route_Engine_Watershed_Grid.py | 503 ++++++++++++++++ .../B05_wf2_Route_Router_Drainage.py | 79 +-- .../B05_wf2_Route_UI_Drainage_Panel.ts | 48 +- .../B05_wf2_Route_UI_Drainage_Pipes.ts | 40 ++ ...B05_wf2_Route_Engine_Drainage_Watershed.py | 0 ...B05_wf2_Route_Engine_Watershed_Assemble.py | 0 ...05_wf2_Route_Engine_Watershed_Subdivide.py | 0 .../B05_wf2_Route_Engine_Watershed_Trace.py | 0 B05_wf2_Route/_legacy_watershed/README.md | 16 + config/config_system.py | 43 ++ pyproject.toml | 2 + 15 files changed, 1460 insertions(+), 71 deletions(-) create mode 100644 B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py create mode 100644 B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py create mode 100644 B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py rename B05_wf2_Route/{ => _legacy_watershed}/B05_wf2_Route_Engine_Drainage_Watershed.py (100%) rename B05_wf2_Route/{ => _legacy_watershed}/B05_wf2_Route_Engine_Watershed_Assemble.py (100%) rename B05_wf2_Route/{ => _legacy_watershed}/B05_wf2_Route_Engine_Watershed_Subdivide.py (100%) rename B05_wf2_Route/{ => _legacy_watershed}/B05_wf2_Route_Engine_Watershed_Trace.py (100%) create mode 100644 B05_wf2_Route/_legacy_watershed/README.md diff --git a/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts b/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts index aff1d624..d28e705b 100644 --- a/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts +++ b/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts @@ -268,7 +268,7 @@ export interface DrainageCandidateResponse { candidates: DrainageCandidate[]; } -/** 배수유역 1개. 관경(pipe_diameter_mm)은 수식 미확정이라 당분간 항상 null이다. */ +/** 관 1개가 받는 세부 배수유역. 관경(pipe_diameter_mm)은 수식 미확정이라 당분간 항상 null이다. */ export interface DrainageBasin { index: number; chainage_m: number; @@ -287,6 +287,12 @@ export interface DrainageBasinResponse { route_id: number; /** 산정에 실제 사용된 배관 지점 — 유역이 없는 관도 포함(마커 동기화용). */ pipes: DrainageCandidate[]; + /** 2차 전체 배수유역 외곽선 = 분수령. 세부유역은 전부 이 안쪽이라 능선을 따로 그리지 않는다. */ + main_polygon_lonlat: Array<[number, number]>; + /** 도로 위 흐름 강도 — [누가거리 m, 그 구간으로 모이는 상류 면적 ㎡]. 관 추가 판단 근거. */ + strength_profile: Array<[number, number]>; + /** 해석에 실제 사용된 격자 한 변(m). 셀 수 상한에 걸리면 백엔드가 키워서 돌려준다. */ + grid_cell_m: number; basins: DrainageBasin[]; } diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py index 5121832c..c46b037a 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage.py @@ -1,8 +1,10 @@ """배수유역 산정 엔진. 관 매설 구조물 측점 후보를 제안하고, 각 측점이 받는 배수유역 경계를 산정한다. -지형 판단은 **도엽 등고선·세류선(하천중심선)·표고점**만 사용한다 — 3D 포인트클라우드나 -지형 메시는 쓰지 않는다(2026-07-28 사용자 지시). +지형 판단은 **도엽 등고선·세류선(하천중심선)**만 사용한다 — 3D 포인트클라우드나 지형 +메시는 쓰지 않고(2026-07-28 사용자 지시), 표고점도 유효 데이터가 적어 뺐다(2026-07-31). +유역 경계 산정 자체는 격자 흐름 해석(`..._Engine_Watershed_Basin`)이 맡고, 이 모듈은 +측점 후보 제안과 노선 정점·누가거리 보간만 담당한다. 유역을 나누는 최종 목적은 각 지점의 파이프 관경 결정이다. 유역 경사면에 100년 강우빈도를 적용해 모이는 물의 양을 산정하고 그 유량으로 관경을 정한다. 관경 수식은 아직 미확정이라 diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py new file mode 100644 index 00000000..b3a54450 --- /dev/null +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py @@ -0,0 +1,562 @@ +"""배수유역 산정 오케스트레이터 — 격자 해석 · 관 배치 · 세부유역 조립. + +전체 흐름 + ① 등고선 정리 → 상류 세류선 추출 → 1차 격자 범위(반경 버퍼 bbox) + ② 격자 지형 해석(TIN 보간 · 채움 · D8) → 도로에서 상류 추적 + ③ 활성 셀이 격자 최외곽에 닿으면 그 방향으로만 넓혀 다시 해석 (경계 링이 전부 + 비활성이 되면 정지 — 하드 반경 상한이 아니라 흐름 자체가 종료 조건이다) + ④ 도로 셀별 흐름 강도(상류 셀 수) 산출 → 2차 전체 배수유역 외곽선 확정 + ⑤ 관 배치: 세류 교차점이 기본, 간격이 최대치를 넘으면 흐름 강도·종단 저점을 보고 + **최소 개수**만 보충 + ⑥ 측구 흐름(종단 내리막)으로 도로 셀 → 담당 관을 정하고, 셀이 도달한 도로 셀의 + 담당 관을 그대로 그 셀의 유역 번호로 삼아 세부유역을 나눈다 + +②~④는 관 배치와 무관하므로 `.npz`로 캐시한다. 사용자가 관을 옮기거나 추가하면 +⑥만 다시 돌면 되고 격자 해석은 재사용한다. +""" + +from __future__ import annotations + +import hashlib +import logging +import time +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np +from shapely.geometry import LineString + +from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( + RouteVertex, + StructureCandidate, + _interpolate_vertex, + estimate_pipe_diameter_mm, + find_stream_crossings, + is_uphill_at, +) +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Flow import ( + border_contact, + largest_ring, + outer_boundary, + polygonize_labels, + rasterize_road, + trace_flow, +) +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import ( + GridSpec, + build_contour_cloud, + build_grid_spec, + build_terrain_grid, + expand_grid_spec, + route_elevation_floor, + select_upstream_streams, +) +from config.config_system import ( + DRAINAGE_DITCH_SAMPLE_M, + DRAINAGE_EXPAND_STEP_M, + DRAINAGE_GRID_SIZE_M, + DRAINAGE_INITIAL_RADIUS_M, + DRAINAGE_MAX_EXPAND_ROUNDS, + DRAINAGE_PIPE_MAX_SPACING_M, + DRAINAGE_PIPE_MIN_SPACING_M, + DRAINAGE_ROAD_WIDTH_M, +) + +logger = logging.getLogger(__name__) + +# 강도 곡선 응답 간격(m). 도로 위 흐름 강도 히트 표기는 이 간격으로 내보낸다. +_STRENGTH_OUTPUT_STEP_M = 5.0 +# 관 위치 선정 점수 배분 — 흐름 강도가 주, 종단 저점이 보조. +_SCORE_WEIGHT_STRENGTH = 0.7 +_SCORE_WEIGHT_SAG = 0.3 +# 성토부(내리막)는 물이 노선 밖으로 빠지므로 관 위치로 덜 선호한다. +_SCORE_FILL_PENALTY = 0.5 + + +@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 WatershedResult: + """배수유역 산정 결과 일체.""" + + basins: list[WatershedBasin] = field(default_factory=list) + pipes: list[StructureCandidate] = field(default_factory=list) + # 2차 전체 배수유역 외곽선(= 분수령). 세부유역 경계는 이 안쪽에서만 그어진다. + main_boundary_xy: list[tuple[float, float]] = field(default_factory=list) + # 도로 위 흐름 강도 곡선 — (누가거리 m, 그 지점으로 모이는 상류 면적 ㎡). + strength_profile: list[tuple[float, float]] = field(default_factory=list) + grid_cell_m: float = DRAINAGE_GRID_SIZE_M + + +@dataclass +class _GridSolution: + """관 배치와 무관한 격자 해석 결과 묶음(캐시 대상).""" + + spec: GridSpec + elevation: np.ndarray # (R*C,) float32 + road_cell_index: np.ndarray # (K,) int32 + road_chainage: np.ndarray # (K,) float64 + road_slot: np.ndarray # (R*C,) int32 — 셀이 도달한 도로 셀 슬롯(−1=미도달) + path_length: np.ndarray # (R*C,) float32 + strength: np.ndarray # (K,) int64 + active: np.ndarray # (R*C,) bool + signature: str + + +# ── 진입점 ────────────────────────────────────────────────────────────────── + + +def build_drainage_watershed( + vertices: list[RouteVertex], + contour_features: list[dict[str, Any]], + stream_features: list[dict[str, Any]], + confirmed_chainages: list[float] | None = None, + cache_path: Path | None = None, +) -> WatershedResult: + """배수유역과 관 배치를 산정한다. + + `confirmed_chainages`를 주면 그 위치를 관으로 확정하고(사용자 편집), 비우면 세류 + 교차 + 최소 보충으로 자동 배치한다. 두 경우 모두 격자 해석은 캐시를 재사용한다. + """ + if len(vertices) < 2: + return WatershedResult() + route_line = LineString([(vertex.x, vertex.y) for vertex in vertices]) + solution = _solve_grid(vertices, route_line, contour_features, stream_features, cache_path) + if solution is None or solution.road_cell_index.size == 0: + return WatershedResult() + + strength_area = solution.strength.astype(np.float64) * solution.spec.cell_area_m2 + strength_curve = _strength_by_chainage(solution.road_chainage, strength_area, route_line.length) + + if confirmed_chainages: + pipes = _pipes_from_chainages(vertices, confirmed_chainages) + else: + pipes = _place_pipes(vertices, stream_features, strength_curve) + if not pipes: + return WatershedResult( + main_boundary_xy=_main_boundary(solution), + strength_profile=_downsample_strength(strength_curve), + grid_cell_m=solution.spec.cell_m, + ) + + pipe_of_slot = _assign_road_cells_to_pipes(vertices, pipes, solution.road_chainage) + basins = _assemble_basins(solution, pipes, pipe_of_slot) + return WatershedResult( + basins=basins, + pipes=pipes, + main_boundary_xy=_main_boundary(solution), + strength_profile=_downsample_strength(strength_curve), + grid_cell_m=solution.spec.cell_m, + ) + + +# ── ①~④ 격자 해석 (캐시 대상) ─────────────────────────────────────────────── + + +def _solve_grid( + vertices: list[RouteVertex], + route_line: LineString, + contour_features: list[dict[str, Any]], + stream_features: list[dict[str, Any]], + cache_path: Path | None, +) -> _GridSolution | None: + signature = _signature(vertices, len(contour_features), len(stream_features)) + cached = _load_cache(cache_path, signature) + if cached is not None: + logger.info("배수유역: 격자 캐시 재사용 (%s)", cache_path) + return cached + + floor = route_elevation_floor([vertex.z for vertex in vertices]) + cloud = build_contour_cloud(contour_features, floor) + if cloud.is_empty: + logger.warning("배수유역: 등고선이 없어 격자 해석을 건너뜁니다.") + return None + streams = select_upstream_streams(route_line, stream_features, cloud) + spec = build_grid_spec(route_line, streams, DRAINAGE_INITIAL_RADIUS_M, DRAINAGE_GRID_SIZE_M) + + terrain = road = flow = None + for round_index in range(DRAINAGE_MAX_EXPAND_ROUNDS + 1): + started = time.perf_counter() + terrain = build_terrain_grid(spec, cloud) + road = rasterize_road(spec, route_line, DRAINAGE_ROAD_WIDTH_M) + flow = trace_flow(terrain, road) + # 격자 크기(config DRAINAGE_GRID_SIZE_M)를 조정할 근거가 되도록 회차별 소요를 남긴다. + logger.info( + "배수유역: %d회차 해석 %.1fs (셀 %d개)", + round_index + 1, + time.perf_counter() - started, + spec.size, + ) + contact = border_contact(flow.active) + if not any(contact.values()): + break + if round_index == DRAINAGE_MAX_EXPAND_ROUNDS: + logger.warning( + "배수유역: 확장 상한(%d회)에 도달 — 경계 %s가 아직 활성입니다.", + DRAINAGE_MAX_EXPAND_ROUNDS, + [side for side, touched in contact.items() if touched], + ) + break + widened = expand_grid_spec(spec, contact, DRAINAGE_EXPAND_STEP_M) + if widened == spec: + break + logger.info( + "배수유역: 경계 %s 활성 — %.0fm 확장", + [side for side, touched in contact.items() if touched], + DRAINAGE_EXPAND_STEP_M, + ) + spec = widened + + assert terrain is not None and road is not None and flow is not None + solution = _GridSolution( + spec=spec, + elevation=terrain.elevation.reshape(-1), + road_cell_index=road.cell_index, + road_chainage=road.chainage, + road_slot=flow.road_slot, + path_length=flow.path_length, + strength=flow.strength, + active=flow.active.reshape(-1), + signature=signature, + ) + _save_cache(cache_path, solution) + return solution + + +def _main_boundary(solution: _GridSolution) -> list[tuple[float, float]]: + boundary = outer_boundary( + solution.spec, solution.active.reshape(solution.spec.n_rows, solution.spec.n_cols) + ) + return largest_ring(boundary) if boundary is not None else [] + + +# ── 흐름 강도 곡선 ────────────────────────────────────────────────────────── + + +def _strength_by_chainage( + road_chainage: np.ndarray, strength_area: np.ndarray, total_length: float +) -> np.ndarray: + """도로 셀 강도를 1m 누가거리 구간으로 합산한 곡선(㎡/m 구간 합).""" + bins = max(1, int(np.ceil(total_length)) + 1) + if road_chainage.size == 0: + return np.zeros(bins) + index = np.clip(np.round(road_chainage).astype(np.int64), 0, bins - 1) + return np.bincount(index, weights=strength_area, minlength=bins) + + +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 + ] + + +# ── ⑤ 관 배치 ─────────────────────────────────────────────────────────────── + + +def _place_pipes( + vertices: list[RouteVertex], + stream_features: list[dict[str, Any]], + strength_curve: np.ndarray, +) -> list[StructureCandidate]: + """세류 교차점을 기본 관 위치로 두고, 최대 간격을 넘는 구간만 최소 개수로 보충한다. + + 교차점은 종단 절·성토를 가리지 않고 모두 관으로 둔다. 하류측 세류선은 이미 격자 + 해석 전에 제거되었으므로, 남은 교차점은 전부 상류에서 물이 실제로 들어오는 지점이다. + """ + total_length = vertices[-1].chainage_m + base: list[StructureCandidate] = [] + for candidate in find_stream_crossings(vertices, stream_features): + 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( + solution: _GridSolution, + pipes: list[StructureCandidate], + pipe_of_slot: np.ndarray, +) -> list[WatershedBasin]: + """셀이 도달한 도로 셀의 담당 관을 그대로 유역 번호로 삼아 세부유역을 만든다.""" + spec = solution.spec + labels = np.full(spec.size, -1, dtype=np.int32) + reached = solution.road_slot >= 0 + labels[reached] = pipe_of_slot[solution.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 = solution.elevation[member] + highest = float(np.nanmax(elevations)) if np.isfinite(elevations).any() else 0.0 + outlet_z = _outlet_elevation(solution, order, pipe_of_slot) + area = count * cell_area + relief = max(0.0, highest - outlet_z) + flow_length = float(solution.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(solution: _GridSolution, 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 = solution.elevation[solution.road_cell_index[slots]] + finite = elevations[np.isfinite(elevations)] + return float(finite.min()) if finite.size else 0.0 + + +# ── 캐시 ──────────────────────────────────────────────────────────────────── + + +def _signature(vertices: list[RouteVertex], contour_count: int, stream_count: int) -> str: + """노선 기하와 해석 파라미터가 바뀌면 캐시를 버리도록 하는 지문.""" + digest = hashlib.sha1() + for vertex in vertices: + digest.update(f"{vertex.x:.2f},{vertex.y:.2f},{vertex.z:.2f};".encode()) + digest.update( + f"|{contour_count}|{stream_count}|{DRAINAGE_GRID_SIZE_M}|{DRAINAGE_INITIAL_RADIUS_M}" + f"|{DRAINAGE_EXPAND_STEP_M}|{DRAINAGE_ROAD_WIDTH_M}".encode() + ) + return digest.hexdigest() + + +def _load_cache(cache_path: Path | None, signature: str) -> _GridSolution | None: + if cache_path is None or not cache_path.exists(): + return None + try: + with np.load(cache_path, allow_pickle=False) as data: + if str(data["signature"]) != signature: + return None + 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"]), + ) + return _GridSolution( + spec=spec, + elevation=data["elevation"], + road_cell_index=data["road_cell_index"], + road_chainage=data["road_chainage"], + road_slot=data["road_slot"], + path_length=data["path_length"], + strength=data["strength"], + active=data["active"], + signature=signature, + ) + except (OSError, KeyError, ValueError): + logger.warning("배수유역: 격자 캐시를 읽지 못해 다시 계산합니다 (%s).", cache_path) + return None + + +def _save_cache(cache_path: Path | None, solution: _GridSolution) -> None: + if cache_path is None: + return + try: + cache_path.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed( + cache_path, + signature=solution.signature, + x_min=solution.spec.x_min, + y_max=solution.spec.y_max, + cell_m=solution.spec.cell_m, + n_rows=solution.spec.n_rows, + n_cols=solution.spec.n_cols, + elevation=solution.elevation, + road_cell_index=solution.road_cell_index, + road_chainage=solution.road_chainage, + road_slot=solution.road_slot, + path_length=solution.path_length, + strength=solution.strength, + active=solution.active, + ) + except OSError: + logger.warning("배수유역: 격자 캐시를 저장하지 못했습니다 (%s).", cache_path) diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py new file mode 100644 index 00000000..bc567ca5 --- /dev/null +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py @@ -0,0 +1,224 @@ +"""배수유역 흐름 해석 — 도로 굽기 · 상류 추적(포인터 더블링) · 유역 폴리곤화. + +핵심은 하나다: **물길을 따라가 도로에 닿는 셀만 유역이다.** +셀마다 D8 수신 셀을 따라가 종착점(root)을 구하고, 그 종착점이 도로 셀이면 활성이다. + +이 방식은 능선을 따로 찾지 않는다. 능선 너머 셀의 물은 다른 계곡으로 빠져 도로에 +닿지 못하므로 자동으로 비활성이 되고, 그 경계선이 곧 능선이다. 유역 안쪽 봉우리는 +물이 결국 도로로 흘러 자동으로 포함된다. + +종착점은 세부유역 라벨의 근거로도 그대로 쓴다 — 셀이 도달한 도로 셀이 정해지면 +그 도로 셀을 담당하는 관이 곧 그 셀의 유역 번호다. 관 배치가 바뀌어도 격자 해석을 +다시 돌릴 필요 없이 "도로 셀 → 관" 대응만 다시 계산하면 된다. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass + +import numpy as np +from rasterio.features import rasterize, shapes +from rasterio.transform import from_origin +from scipy.spatial import cKDTree +from shapely.geometry import LineString, MultiPolygon, Polygon, shape +from shapely.ops import unary_union + +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import GridSpec, TerrainGrid +from config.config_system import ( + DRAINAGE_MIN_BASIN_AREA_M2, + DRAINAGE_POLYGON_SIMPLIFY_M, + DRAINAGE_ROAD_WIDTH_M, +) + +logger = logging.getLogger(__name__) + +# 포인터 더블링 반복 상한. 한 번에 경로 길이가 2배가 되므로 2^40 스텝이면 어떤 격자도 덮는다. +_MAX_DOUBLING_ROUNDS = 40 + + +@dataclass +class RoadRaster: + """격자에 구운 도로. 도로 셀은 흐름을 흡수하는 종착점이 된다.""" + + mask: np.ndarray # (R, C) bool + cell_index: np.ndarray # (K,) int32 — 도로 셀의 평탄 인덱스 + chainage: np.ndarray # (K,) float64 — 도로 셀의 누가거리(m) + slot_of_cell: np.ndarray # (R*C,) int32 — 도로 셀이면 K 내 위치, 아니면 −1 + + @property + def count(self) -> int: + return int(self.cell_index.size) + + +@dataclass +class FlowResult: + """상류 추적 결과.""" + + root: np.ndarray # (R*C,) int32 — 흐름 종착 셀의 평탄 인덱스 + road_slot: np.ndarray # (R*C,) int32 — 도달한 도로 셀 슬롯, 도달 못하면 −1 + active: np.ndarray # (R, C) bool — 도로에 물이 닿는 셀 + path_length: np.ndarray # (R*C,) float32 — 종착점까지 물길 길이(m) + strength: np.ndarray # (K,) int64 — 도로 셀별 상류 셀 수(흐름 강도) + + +def grid_transform(spec: GridSpec): + """rasterio 아핀 변환. 행 0이 북쪽(y_max)이다.""" + return from_origin(spec.x_min, spec.y_max, spec.cell_m, spec.cell_m) + + +# ── 도로 굽기 ─────────────────────────────────────────────────────────────── + + +def rasterize_road( + spec: GridSpec, + route_line: LineString, + width_m: float = DRAINAGE_ROAD_WIDTH_M, +) -> RoadRaster: + """노선을 노폭만큼 두껍게 격자에 굽고, 각 도로 셀에 누가거리를 붙인다. + + 폭을 주는 이유는 실제 노면이 물을 받기 때문이기도 하지만, 1셀 선으로 구우면 D8 + 대각 이동이 도로를 건너뛰어 상류 물이 도로를 지나쳐 버리기 때문이다. 3셀 이상 두께면 + 내리막 물길이 반드시 도로 셀을 한 번은 밟는다. + """ + half_width = max(width_m / 2.0, spec.cell_m) + burned = rasterize( + [(route_line.buffer(half_width), 1)], + out_shape=(spec.n_rows, spec.n_cols), + transform=grid_transform(spec), + fill=0, + dtype="uint8", + all_touched=True, + ).astype(bool) + + slot_of_cell = np.full(spec.size, -1, dtype=np.int32) + cell_index = np.flatnonzero(burned.reshape(-1)).astype(np.int32) + if cell_index.size == 0: + logger.warning("배수유역: 노선이 격자 범위 밖입니다 — 도로 셀 0개.") + return RoadRaster(burned, cell_index, np.zeros(0), slot_of_cell) + + # 도로 셀 누가거리는 노선을 촘촘히 샘플링해 가장 가까운 샘플의 누가거리로 준다. + step = max(spec.cell_m / 2.0, 0.25) + positions = np.arange(0.0, route_line.length + step, step) + samples = np.array([list(route_line.interpolate(p).coords)[0] for p in positions]) + rows = (cell_index // spec.n_cols).astype(np.float64) + cols = (cell_index % spec.n_cols).astype(np.float64) + centers = np.column_stack( + ( + spec.x_min + (cols + 0.5) * spec.cell_m, + spec.y_max - (rows + 0.5) * spec.cell_m, + ) + ) + _, nearest = cKDTree(samples).query(centers) + chainage = np.minimum(positions[nearest], route_line.length) + slot_of_cell[cell_index] = np.arange(cell_index.size, dtype=np.int32) + return RoadRaster(burned, cell_index, chainage, slot_of_cell) + + +# ── 상류 추적 ─────────────────────────────────────────────────────────────── + + +def trace_flow(terrain: TerrainGrid, road: RoadRaster) -> FlowResult: + """모든 셀의 물길 종착점을 구하고 도로 도달 여부(=유역 포함 여부)를 판정한다. + + 포인터 더블링으로 한 번에 경로 길이를 2배씩 늘려 종착점을 찾는다. 채움·평탄해소를 + 거친 표고에서는 흐름을 따라 표고가 단조 감소하므로 순환이 없고, 반복은 항상 끝난다. + """ + spec = terrain.spec + receiver = terrain.receiver.astype(np.int32, copy=True) + step_length = terrain.step_length.astype(np.float32, copy=True) + + # 도로 셀은 흐름을 흡수한다 — 물이 도로에 닿으면 거기서 끝난다. + receiver[road.cell_index] = road.cell_index + step_length[road.cell_index] = 0.0 + + jump = receiver + path_length = step_length + for _ in range(_MAX_DOUBLING_ROUNDS): + next_jump = jump[jump] + if np.array_equal(next_jump, jump): + break + path_length = path_length + path_length[jump] + jump = next_jump + + road_slot = road.slot_of_cell[jump] + active_flat = road_slot >= 0 + strength = ( + np.bincount(road_slot[active_flat], minlength=max(road.count, 1)).astype(np.int64) + if road.count + else np.zeros(0, dtype=np.int64) + ) + logger.info( + "배수유역: 활성 셀 %d / %d (도로 셀 %d)", int(active_flat.sum()), spec.size, road.count + ) + return FlowResult( + root=jump, + road_slot=road_slot, + active=active_flat.reshape(spec.n_rows, spec.n_cols), + path_length=path_length, + strength=strength, + ) + + +def border_contact(active: np.ndarray) -> dict[str, bool]: + """활성 셀이 격자 최외곽에 닿은 방향. 전부 False면 유역이 능선 안에서 닫힌 것이다.""" + return { + "north": bool(active[0, :].any()), + "south": bool(active[-1, :].any()), + "west": bool(active[:, 0].any()), + "east": bool(active[:, -1].any()), + } + + +# ── 폴리곤화 ──────────────────────────────────────────────────────────────── + + +def polygonize_labels( + spec: GridSpec, + labels: np.ndarray, + min_area_m2: float = DRAINAGE_MIN_BASIN_AREA_M2, +) -> dict[int, Polygon | MultiPolygon]: + """라벨 격자를 라벨별 폴리곤으로 바꾼다. 음수 라벨은 배경으로 무시한다.""" + label_grid = np.ascontiguousarray(labels.reshape(spec.n_rows, spec.n_cols), dtype=np.int32) + valid_mask = label_grid >= 0 + if not valid_mask.any(): + return {} + collected: dict[int, list[Polygon]] = {} + for geometry, value in shapes( + label_grid, mask=valid_mask, transform=grid_transform(spec), connectivity=4 + ): + polygon = shape(geometry) + if polygon.is_empty or polygon.area < min_area_m2: + continue + collected.setdefault(int(value), []).append(polygon) + + merged: dict[int, Polygon | MultiPolygon] = {} + for label, parts in collected.items(): + union = unary_union(parts) + if union.is_empty: + continue + simplified = union.simplify(DRAINAGE_POLYGON_SIMPLIFY_M, preserve_topology=True) + merged[label] = simplified if not simplified.is_empty else union + return merged + + +def largest_ring(geometry: Polygon | MultiPolygon) -> list[tuple[float, float]]: + """폴리곤(또는 멀티폴리곤)에서 가장 큰 조각의 외곽 링 좌표를 뽑는다.""" + if geometry.is_empty: + return [] + if geometry.geom_type == "MultiPolygon": + geometry = max(geometry.geoms, key=lambda part: part.area) + return [(float(x), float(y)) for x, y in geometry.exterior.coords] + + +def outer_boundary( + spec: GridSpec, active: np.ndarray, min_area_m2: float = DRAINAGE_MIN_BASIN_AREA_M2 +) -> Polygon | MultiPolygon | None: + """활성 셀 전체의 외곽 = 2차 전체 배수유역 경계. + + 비활성 셀이 격자 최외곽에 띠로 완성되면 활성 영역이 그 안에 닫힌다. 그 닫힌 영역의 + 바깥선이 곧 분수령이므로 능선을 따로 그릴 필요가 없다. + """ + labels = np.where(active.reshape(-1), 0, -1).astype(np.int32) + polygons = polygonize_labels(spec, labels, min_area_m2) + return polygons.get(0) diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py new file mode 100644 index 00000000..54f4a7ad --- /dev/null +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py @@ -0,0 +1,503 @@ +"""배수유역 해석 격자 생성 — 등고선 TIN 보간 · 웅덩이 채움 · D8 물 방향. + +지형 근거는 **도엽 등고선**뿐이다. 라이다 DEM은 노선 주변만 커버해 유역 산정에 필요한 +상류 범위를 담지 못하므로 쓰지 않는다(2026-07-31 사용자 지시). 표고점도 쓰지 않는다. + +처리 순서 + ① 계획선 최저점 아래 등고선·짧은 파편 제거 + ② 세류선을 노선 교차점에서 잘라 상류측만 남김 + ③ 남은 세류선 + 노선을 반경 버퍼한 범위의 bbox로 격자 생성 + ④ 등고선 정점 Delaunay TIN 선형보간으로 셀 표고 산출 + ⑤ 웅덩이 채움(형태학적 재구성) + 평탄면 미세경사 부여 + ⑥ D8(8방향 최급강하) 수신 셀 인덱스 산출 + +⑤가 없으면 등고선 TIN 특유의 가짜 웅덩이·평탄 삼각형에서 흐름이 끊겨 상류 추적이 +도중에 멈춘다. 능선 탐지는 하지 않는다 — 흐름이 도로에 닿는지 여부만으로 유역이 정해진다. +""" + +from __future__ import annotations + +import logging +import math +from dataclasses import dataclass +from typing import Any + +import numpy as np +from scipy.interpolate import LinearNDInterpolator +from scipy.ndimage import distance_transform_edt +from scipy.spatial import cKDTree +from shapely import segmentize +from shapely.geometry import LineString, shape +from shapely.ops import substring, unary_union +from skimage.morphology import reconstruction + +from config.config_system import ( + DRAINAGE_CONTOUR_MARGIN_M, + DRAINAGE_CONTOUR_MIN_LENGTH_M, + DRAINAGE_CONTOUR_RESAMPLE_M, + DRAINAGE_FLAT_EPSILON_M, + DRAINAGE_GRID_SIZE_M, + DRAINAGE_MAX_GRID_CELLS, +) + +logger = logging.getLogger(__name__) + +# 표고 속성 키: 도엽 등고선(등고수치)·gpkg 등고선(CTRLN_HG) 통합. +ELEVATION_KEYS = ("등고수치", "CTRLN_HG", "수치", "표고", "높이", "elevation", "ELEV") + +# D8 이웃 (행 증분, 열 증분, 거리계수). 행은 아래로 증가(북 → 남). +_NEIGHBORS = ( + (-1, 0, 1.0), + (1, 0, 1.0), + (0, -1, 1.0), + (0, 1, 1.0), + (-1, -1, math.sqrt(2.0)), + (-1, 1, math.sqrt(2.0)), + (1, -1, math.sqrt(2.0)), + (1, 1, math.sqrt(2.0)), +) + + +@dataclass(frozen=True) +class GridSpec: + """해석 격자 기하. (0,0) 셀 중심이 (x_min + cell/2, y_max − cell/2)에 놓인다.""" + + x_min: float + y_max: float + cell_m: float + n_rows: int + n_cols: int + + @property + def size(self) -> int: + return self.n_rows * self.n_cols + + @property + def cell_area_m2(self) -> float: + return self.cell_m * self.cell_m + + def world_to_rc(self, x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """세계좌표(m)를 격자 행·열로 바꾼다. 범위를 벗어나면 −1을 돌려준다.""" + col = np.floor((x - self.x_min) / self.cell_m).astype(np.int64) + row = np.floor((self.y_max - y) / self.cell_m).astype(np.int64) + outside = (col < 0) | (col >= self.n_cols) | (row < 0) | (row >= self.n_rows) + col[outside] = -1 + row[outside] = -1 + return row, col + + def cell_centers_x(self) -> np.ndarray: + return self.x_min + (np.arange(self.n_cols, dtype=np.float64) + 0.5) * self.cell_m + + def cell_centers_y(self) -> np.ndarray: + return self.y_max - (np.arange(self.n_rows, dtype=np.float64) + 0.5) * self.cell_m + + +@dataclass +class ContourCloud: + """등고선에서 뽑은 정점 구름. TIN 보간과 세류 상·하류 판정에 함께 쓴다.""" + + xy: np.ndarray # (N, 2) float64 + z: np.ndarray # (N,) float64 + + @property + def is_empty(self) -> bool: + return self.xy.shape[0] < 3 + + +@dataclass +class TerrainGrid: + """격자 지형 해석 결과.""" + + spec: GridSpec + elevation: np.ndarray # (R, C) float32 — 채움·평탄해소 후 표고, 무효 셀은 NaN + valid: np.ndarray # (R, C) bool — 등고선 TIN 내부 여부 + receiver: np.ndarray # (R*C,) int32 — D8 수신 셀의 평탄 인덱스, 싱크는 자기 자신 + step_length: np.ndarray # (R*C,) float32 — 수신 셀까지 거리(m), 싱크는 0 + + +# ── ① 등고선 정리 ──────────────────────────────────────────────────────────── + + +def _feature_elevation(properties: dict[str, Any]) -> float | None: + for key in ELEVATION_KEYS: + value = properties.get(key) + if value is None: + continue + try: + return float(value) + except (TypeError, ValueError): + continue + return None + + +def _iter_linestrings(geometry: Any) -> list[LineString]: + if geometry.geom_type == "LineString": + return [geometry] + if geometry.geom_type in {"MultiLineString", "GeometryCollection"}: + lines: list[LineString] = [] + for part in geometry.geoms: + lines.extend(_iter_linestrings(part)) + return lines + return [] + + +def build_contour_cloud( + contour_features: list[dict[str, Any]], + elevation_floor_m: float | None = None, +) -> ContourCloud: + """등고선 피처를 표고가 붙은 정점 구름으로 바꾼다. + + `elevation_floor_m` 아래 등고선은 계획선 최저점보다 낮아 상류 기여가 불가능하므로 + 버린다. 길이가 짧은 파편도 노이즈로 보고 버리되, 임계 이상인 봉우리 폐합 등고선은 + 남긴다(봉우리 표고가 사라지면 그 일대 흐름 방향이 통째로 틀어진다). + """ + xs: list[np.ndarray] = [] + ys: list[np.ndarray] = [] + zs: list[np.ndarray] = [] + dropped_low = 0 + dropped_short = 0 + for feature in contour_features: + geometry = feature.get("geometry") + if not geometry: + continue + elevation = _feature_elevation(feature.get("properties") or {}) + if elevation is None: + continue + if elevation_floor_m is not None and elevation < elevation_floor_m: + dropped_low += 1 + continue + try: + parsed = shape(geometry) + except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다 + continue + for line in _iter_linestrings(parsed): + if line.length < DRAINAGE_CONTOUR_MIN_LENGTH_M: + dropped_short += 1 + continue + coords = np.asarray(segmentize(line, DRAINAGE_CONTOUR_RESAMPLE_M).coords) + if coords.shape[0] < 2: + continue + xs.append(coords[:, 0]) + ys.append(coords[:, 1]) + zs.append(np.full(coords.shape[0], elevation, dtype=np.float64)) + if not xs: + logger.warning( + "배수유역: 사용할 등고선이 없습니다(저지대 %d, 파편 %d 제외).", + dropped_low, + dropped_short, + ) + return ContourCloud(np.zeros((0, 2)), np.zeros(0)) + xy = np.column_stack((np.concatenate(xs), np.concatenate(ys))) + z = np.concatenate(zs) + logger.info( + "배수유역: 등고선 정점 %d개 (저지대 %d, 파편 %d 제외)", + xy.shape[0], + dropped_low, + dropped_short, + ) + return ContourCloud(xy, z) + + +# ── ② 세류선 상류측만 남기기 ──────────────────────────────────────────────── + + +def select_upstream_streams( + route_line: LineString, + stream_features: list[dict[str, Any]], + cloud: ContourCloud, +) -> list[LineString]: + """노선과 교차하는 세류선을 교차점에서 잘라 상류(고지대)측만 돌려준다. + + 노선과 만나지 않는 세류선은 판단 근거가 없으므로 그대로 남긴다 — 어차피 격자 해석에서 + 도로에 물이 닿지 않으면 비활성 처리된다. 상·하류 판정은 가장 가까운 등고선 정점의 + 표고 평균으로 한다(TIN은 이 시점에 아직 없다). + """ + tree = cKDTree(cloud.xy) if not cloud.is_empty else None + kept: list[LineString] = [] + for feature in stream_features: + geometry = feature.get("geometry") + if not geometry: + continue + try: + parsed = shape(geometry) + except Exception: # noqa: BLE001 + continue + for line in _iter_linestrings(parsed): + if line.is_empty or line.length <= 0: + continue + if not line.intersects(route_line): + kept.append(line) + continue + kept.extend(_upstream_parts(line, route_line, tree, cloud)) + return kept + + +def _upstream_parts( + line: LineString, + route_line: LineString, + tree: cKDTree | None, + cloud: ContourCloud, +) -> list[LineString]: + """세류선을 노선 교차점에서 잘라 평균 표고가 높은 조각만 남긴다.""" + cuts = sorted( + { + line.project(point) + for point in _intersection_points(line.intersection(route_line)) + if 0.0 < line.project(point) < line.length + } + ) + if not cuts: + return [line] + bounds = [0.0, *cuts, line.length] + parts: list[tuple[float, LineString]] = [] + for start, end in zip(bounds, bounds[1:]): + if end - start < 1.0: + continue + piece = _substring(line, start, end) + if piece is None: + continue + parts.append((_mean_elevation(piece, tree, cloud), piece)) + if not parts: + return [] + highest = max(value for value, _ in parts) + # 최상류 조각과 표고가 비슷한(1m 이내) 조각까지 상류로 본다. 나머지는 하류이므로 버린다. + return [piece for value, piece in parts if highest - value <= 1.0] + + +def _intersection_points(geometry: Any) -> list[Any]: + if geometry.is_empty: + return [] + if geometry.geom_type == "Point": + return [geometry] + if geometry.geom_type in {"MultiPoint", "GeometryCollection", "MultiLineString"}: + points: list[Any] = [] + for part in geometry.geoms: + points.extend(_intersection_points(part)) + return points + if geometry.geom_type == "LineString": + return [geometry.interpolate(0.5, normalized=True)] + return [] + + +def _substring(line: LineString, start: float, end: float) -> LineString | None: + """선형 위 [start, end] 구간을 잘라낸다.""" + piece = substring(line, start, end) + if piece.is_empty or piece.geom_type != "LineString" or piece.length <= 0: + return None + return piece + + +def _mean_elevation(line: LineString, tree: cKDTree | None, cloud: ContourCloud) -> float: + if tree is None: + return 0.0 + samples = max(2, int(line.length // 10.0) + 1) + positions = np.linspace(0.0, line.length, samples) + points = np.array([list(line.interpolate(position).coords)[0] for position in positions]) + _, indices = tree.query(points) + return float(np.mean(cloud.z[indices])) + + +# ── ③ 격자 범위 ───────────────────────────────────────────────────────────── + + +def build_grid_spec( + route_line: LineString, + streams: list[LineString], + radius_m: float, + cell_m: float = DRAINAGE_GRID_SIZE_M, +) -> GridSpec: + """노선과 상류 세류선을 반경 버퍼한 범위의 bbox로 격자를 잡는다. + + 셀 수가 상한을 넘으면 셀 크기를 자동으로 키워 맞춘다(메모리 보호). 실제 유역 모양은 + 격자가 아니라 흐름 해석이 정한다 — 여기서는 넉넉한 사각 범위만 확보하면 된다. + """ + geometries = [route_line.buffer(radius_m)] + geometries.extend(line.buffer(radius_m) for line in streams) + x_min, y_min, x_max, y_max = unary_union(geometries).bounds + return _spec_from_bounds(x_min, y_min, x_max, y_max, cell_m) + + +def _spec_from_bounds( + x_min: float, y_min: float, x_max: float, y_max: float, cell_m: float +) -> GridSpec: + width = max(x_max - x_min, cell_m) + height = max(y_max - y_min, cell_m) + while (width / cell_m) * (height / cell_m) > DRAINAGE_MAX_GRID_CELLS: + cell_m *= 2.0 + logger.warning("배수유역: 셀 수 상한 초과 — 격자 크기를 %.1fm로 키웁니다.", cell_m) + n_cols = int(math.ceil(width / cell_m)) + n_rows = int(math.ceil(height / cell_m)) + return GridSpec( + x_min=x_min, y_max=y_min + n_rows * cell_m, cell_m=cell_m, n_rows=n_rows, n_cols=n_cols + ) + + +def expand_grid_spec(spec: GridSpec, sides: dict[str, bool], step_m: float) -> GridSpec: + """활성 셀이 닿은 방향으로만 격자를 넓힌다.""" + x_min = spec.x_min - (step_m if sides.get("west") else 0.0) + x_max = spec.x_min + spec.n_cols * spec.cell_m + (step_m if sides.get("east") else 0.0) + y_max = spec.y_max + (step_m if sides.get("north") else 0.0) + y_min = spec.y_max - spec.n_rows * spec.cell_m - (step_m if sides.get("south") else 0.0) + return _spec_from_bounds(x_min, y_min, x_max, y_max, spec.cell_m) + + +# ── ④ TIN 보간 ────────────────────────────────────────────────────────────── + + +def interpolate_elevation(spec: GridSpec, cloud: ContourCloud) -> np.ndarray: + """등고선 정점 Delaunay TIN으로 셀 표고를 선형보간한다. 외부는 NaN.""" + surface = np.full((spec.n_rows, spec.n_cols), np.nan, dtype=np.float32) + if cloud.is_empty: + return surface + interpolator = LinearNDInterpolator(cloud.xy, cloud.z) + xs = spec.cell_centers_x() + ys = spec.cell_centers_y() + # 행 묶음 단위로 평가해 (행×열) 좌표 배열을 한 번에 들고 있지 않게 한다. + chunk = max(1, int(4_000_000 // max(spec.n_cols, 1))) + for start in range(0, spec.n_rows, chunk): + stop = min(start + chunk, spec.n_rows) + grid_x, grid_y = np.meshgrid(xs, ys[start:stop]) + surface[start:stop] = interpolator(grid_x, grid_y).astype(np.float32) + return surface + + +# ── ⑤ 웅덩이 채움 + 평탄면 해소 ───────────────────────────────────────────── + + +def condition_surface(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """가짜 웅덩이를 채우고 평탄면에 미세 경사를 준다. + + 등고선 TIN은 같은 표고 정점 3개로 이루어진 평탄 삼각형과 계단형 가짜 웅덩이를 + 필연적으로 만든다. 그대로 D8을 돌리면 흐름이 거기서 끊겨 상류 추적이 멈춘다. + + 채움은 형태학적 재구성(erosion)으로 한다. 배출구는 격자 최외곽과 TIN 경계(무효 셀에 + 맞닿은 유효 셀)로 둔다 — 그래야 유효 영역 전체가 하나의 평탄면으로 잠기지 않는다. + """ + valid = np.isfinite(surface) + if not valid.any(): + return surface, valid + ceiling = float(np.nanmax(surface)) + 1000.0 + mask = np.where(valid, surface, ceiling).astype(np.float32) + + open_boundary = np.zeros_like(valid) + open_boundary[0, :] = True + open_boundary[-1, :] = True + open_boundary[:, 0] = True + open_boundary[:, -1] = True + open_boundary |= _dilate(~valid) & valid + open_boundary &= valid + if not open_boundary.any(): + open_boundary = valid & _dilate(~valid) + + seed = np.full_like(mask, ceiling) + seed[open_boundary] = mask[open_boundary] + filled = reconstruction(seed, mask, method="erosion", footprint=np.ones((3, 3), dtype=bool)) + filled = filled.astype(np.float32) + + # 채움 뒤 더 낮은 이웃이 없는 셀 = 평탄면. 가장 가까운 비평탄 셀 쪽으로 미세 경사를 준다. + flat = valid & ~_has_lower_neighbour(filled, valid) + if flat.any(): + distance = distance_transform_edt(flat).astype(np.float32) + filled = filled + distance * np.float32(DRAINAGE_FLAT_EPSILON_M) + filled[~valid] = np.nan + return filled, valid + + +def _dilate(mask: np.ndarray) -> np.ndarray: + """8이웃 1스텝 팽창(외부는 False).""" + padded = np.zeros((mask.shape[0] + 2, mask.shape[1] + 2), dtype=bool) + padded[1:-1, 1:-1] = mask + result = np.zeros_like(mask) + for row_shift in (0, 1, 2): + for col_shift in (0, 1, 2): + result |= padded[ + row_shift : row_shift + mask.shape[0], col_shift : col_shift + mask.shape[1] + ] + return result + + +def _has_lower_neighbour(surface: np.ndarray, valid: np.ndarray) -> np.ndarray: + """8이웃 중 자기보다 낮은 셀이 하나라도 있는지. 무효 셀은 +∞로 보아 제외한다.""" + rows, cols = surface.shape + padded = np.full((rows + 2, cols + 2), np.inf, dtype=np.float32) + padded[1:-1, 1:-1] = np.where(valid, surface, np.inf) + result = np.zeros((rows, cols), dtype=bool) + center = padded[1:-1, 1:-1] + for row_shift, col_shift, _ in _NEIGHBORS: + neighbour = padded[ + 1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols + ] + result |= neighbour < center + return result & valid + + +# ── ⑥ D8 물 방향 ──────────────────────────────────────────────────────────── + + +def compute_receivers( + spec: GridSpec, surface: np.ndarray, valid: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """셀마다 8방향 최급강하 이웃(수신 셀)을 정한다. + + 돌려주는 `receiver`는 평탄 인덱스(row * n_cols + col)다. 더 낮은 이웃이 없는 셀(싱크)과 + 무효 셀은 자기 자신을 가리켜 흐름이 그 자리에서 멈춘다. + """ + rows, cols = spec.n_rows, spec.n_cols + padded = np.full((rows + 2, cols + 2), np.inf, dtype=np.float32) + padded[1:-1, 1:-1] = np.where(valid, surface, np.inf) + center = padded[1:-1, 1:-1] + + flat_index = np.arange(rows * cols, dtype=np.int32).reshape(rows, cols) + padded_index = np.full((rows + 2, cols + 2), -1, dtype=np.int32) + padded_index[1:-1, 1:-1] = flat_index + + best_slope = np.zeros((rows, cols), dtype=np.float32) + receiver = flat_index.copy() + step = np.zeros((rows, cols), dtype=np.float32) + + for row_shift, col_shift, factor in _NEIGHBORS: + neighbour = padded[ + 1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols + ] + distance = np.float32(factor * spec.cell_m) + # 무효 셀끼리는 ∞−∞ = NaN이 되지만 아래 isfinite에서 걸러진다. + with np.errstate(invalid="ignore"): + slope = (center - neighbour) / distance + better = np.isfinite(slope) & (slope > best_slope) + if not better.any(): + continue + best_slope = np.where(better, slope, best_slope) + neighbour_index = padded_index[ + 1 + row_shift : 1 + row_shift + rows, 1 + col_shift : 1 + col_shift + cols + ] + receiver = np.where(better, neighbour_index, receiver) + step = np.where(better, distance, step) + return receiver.reshape(-1), step.reshape(-1) + + +# ── 오케스트레이션 ────────────────────────────────────────────────────────── + + +def build_terrain_grid(spec: GridSpec, cloud: ContourCloud) -> TerrainGrid: + """격자 범위와 등고선 구름으로 지형 해석 격자를 만든다.""" + surface = interpolate_elevation(spec, cloud) + conditioned, valid = condition_surface(surface) + receiver, step = compute_receivers(spec, conditioned, valid) + logger.info( + "배수유역: 격자 %d×%d (%.1fm), 유효 셀 %d개", + spec.n_rows, + spec.n_cols, + spec.cell_m, + int(valid.sum()), + ) + return TerrainGrid( + spec=spec, elevation=conditioned, valid=valid, receiver=receiver, step_length=step + ) + + +def route_elevation_floor(route_z_values: list[float]) -> float | None: + """계획선 최저점에서 여유를 뺀 등고선 하한. 값이 없으면 None(필터 미적용).""" + finite = [value for value in route_z_values if math.isfinite(value) and value != 0.0] + if not finite: + return None + return min(finite) - DRAINAGE_CONTOUR_MARGIN_M diff --git a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py index f489eed8..01b8b369 100644 --- a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py @@ -1,9 +1,11 @@ """배수유역도 API 라우터. -구조물 측점(관 매설) 후보 제안과 배수유역 산정을 제공한다. 지형 근거는 도엽 등고선·세류선· -표고점 GeoJSON뿐이며, 좌표는 사업지 CRS(m)에서 계산하고 응답만 WGS84로 바꿔 내보낸다. +구조물 측점(관 매설) 후보 제안과 배수유역 산정을 제공한다. 지형 근거는 **도엽 등고선과 +세류선 GeoJSON**뿐이며(표고점은 유효 데이터가 적어 2026-07-31 사용자 지시로 제외), +좌표는 사업지 CRS(m)에서 계산하고 응답만 WGS84로 바꿔 내보낸다. """ +import asyncio import json import logging from pathlib import Path @@ -20,7 +22,7 @@ from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( build_route_vertices, propose_structure_stations, ) -from B05_wf2_Route.B05_wf2_Route_Engine_Drainage_Watershed import build_watershed_basins +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Basin import build_drainage_watershed from B05_wf2_Route.B05_wf2_Route_Repository import ( get_latest_route, get_route_points, @@ -28,6 +30,7 @@ from B05_wf2_Route.B05_wf2_Route_Repository import ( ) from common_util.common_util_storage import resolve_stored_project_path from config.config_db import get_db_pool +from config.config_system import DRAINAGE_CACHE_DIRNAME, DRAINAGE_CACHE_FILENAME logger = logging.getLogger(__name__) router = APIRouter(prefix="/api/projects", tags=["B05 Route Drainage"]) @@ -35,15 +38,18 @@ router = APIRouter(prefix="/api/projects", tags=["B05 Route Drainage"]) # 도엽 레이어 파일명 (B04 전처리 산출물과 동일 위치) _CONTOUR_FILE = "도엽_등고선.geojson" _STREAM_FILE = "도엽_하천중심선.geojson" -_SPOT_FILE = "도엽_표고점.geojson" -# 표고 속성 키: 도엽 등고선(등고수치)·gpkg 등고선(CTRLN_HG)·표고점(수치/표고) 통합. -_ELEVATION_KEYS = ("등고수치", "CTRLN_HG", "수치", "표고", "높이", "elevation", "ELEV") def _sheet_dir(stored_path: str) -> Path: return Path(resolve_stored_project_path(stored_path)) / "B04_wf1_Surface" / "processed" +def _cache_path(stored_path: str) -> Path: + """격자 해석 캐시(.npz) 경로. 관을 옮겨도 격자를 다시 풀지 않게 여기에 남긴다.""" + root = Path(resolve_stored_project_path(stored_path)) / "B05_wf2_Route" + return root / DRAINAGE_CACHE_DIRNAME / DRAINAGE_CACHE_FILENAME + + def _load_features(directory: Path, filename: str) -> list[dict[str, Any]]: """도엽 GeoJSON을 읽어 피처 목록만 돌려준다. 없으면 빈 목록.""" path = directory / filename @@ -145,15 +151,12 @@ async def _prepare(project_id: UUID) -> dict[str, Any] | JSONResponse: contour_features = _reproject_features( _load_features(directory, _CONTOUR_FILE), to_metric_transformer ) - spot_features = _reproject_features( - _load_features(directory, _SPOT_FILE), to_metric_transformer - ) return { "route_id": int(route["id"]), "vertices": vertices, "streams": streams, "contours": contour_features, - "spots": spot_features, + "cache_path": _cache_path(stored_path), "to_lonlat": lambda x, y: to_lonlat_transformer.transform(x, y), } @@ -179,27 +182,26 @@ async def post_drainage_basins( project_id: UUID, payload: dict[str, Any] | None = None, ) -> dict[str, Any] | JSONResponse: - """확정된 구조물 측점별 배수유역을 산정한다. + """격자 흐름 해석으로 배수유역과 관 배치를 산정한다. - payload에 `chainages`(누가거리 목록)를 주면 그 위치로 확정하고, 없으면 자동 제안분을 쓴다. + payload에 `chainages`(누가거리 목록)를 주면 그 위치로 관을 확정하고, 없으면 세류 + 교차 + 최소 보충으로 자동 배치한다. 격자 해석은 `.npz` 캐시를 재사용하므로 관만 + 옮기는 재요청은 세부유역 분할만 다시 돈다. """ prepared = await _prepare(project_id) if isinstance(prepared, JSONResponse): return prepared - vertices = prepared["vertices"] - chainages = (payload or {}).get("chainages") - if isinstance(chainages, list) and chainages: - candidates = _candidates_from_chainages(vertices, chainages) - else: - candidates = propose_structure_stations(vertices, prepared["streams"]) + raw_chainages = (payload or {}).get("chainages") + confirmed = _parse_chainages(raw_chainages) if isinstance(raw_chainages, list) else [] - basins = build_watershed_basins( - vertices, - candidates, + # 격자 해석은 수백만 셀 numpy 연산이라 이벤트 루프를 막지 않도록 스레드로 뺀다. + result = await asyncio.to_thread( + build_drainage_watershed, + prepared["vertices"], prepared["contours"], - prepared["spots"], - _ELEVATION_KEYS, - stream_features=prepared["streams"], + prepared["streams"], + confirmed, + prepared["cache_path"], ) to_lonlat = prepared["to_lonlat"] return { @@ -207,17 +209,20 @@ async def post_drainage_basins( "project_id": str(project_id), "route_id": prepared["route_id"], # 계획선 위 배관(관 매설) 지점 — 유역이 없는 관도 마커로 표시해야 하므로 별도 목록. - "pipes": [ - _candidate_payload(candidate, to_lonlat) - for candidate in sorted(candidates, key=lambda item: item.chainage_m) + "pipes": [_candidate_payload(candidate, to_lonlat) for candidate in result.pipes], + # 2차 전체 배수유역 외곽선 = 분수령. 세부유역은 전부 이 안쪽에 들어간다. + "main_polygon_lonlat": [list(to_lonlat(x, y)) for x, y in result.main_boundary_xy], + # 도로 위 흐름 강도 — [누가거리 m, 그 지점으로 모이는 상류 면적 ㎡]. + "strength_profile": [ + [round(chainage, 1), round(area, 1)] for chainage, area in result.strength_profile ], + "grid_cell_m": result.grid_cell_m, "basins": [ { "index": basin.index, "chainage_m": round(basin.chainage_m, 2), # 관(배관) 매설 지점 좌표 — 계획선 위 마커 렌더용. "outlet_lonlat": list(to_lonlat(basin.outlet_x, basin.outlet_y)), - # 유역 경계 외곽선 = 분수령(능선). 프론트가 파스텔 채움 + 능선 파선으로 표시한다. "polygon_lonlat": [list(to_lonlat(x, y)) for x, y in basin.boundary_xy], "area_m2": round(basin.area_m2, 1), "relief_m": round(basin.relief_m, 2), @@ -225,21 +230,17 @@ async def post_drainage_basins( # 관경 수식 미확정 — 산정 함수가 None을 돌려주면 프론트가 "미정"으로 표기한다. "pipe_diameter_mm": basin.pipe_diameter_mm, } - for basin in basins + for basin in result.basins ], } -def _candidates_from_chainages(vertices: Any, chainages: list[Any]) -> list[StructureCandidate]: - """사용자가 확정한 누가거리 목록을 후보 구조로 되돌린다.""" - from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import _interpolate_vertex - - candidates: list[StructureCandidate] = [] - for value in chainages: +def _parse_chainages(values: list[Any]) -> list[float]: + """사용자가 확정·편집한 누가거리 목록을 숫자로 정리한다.""" + parsed: list[float] = [] + for value in values: try: - chainage = float(value) + parsed.append(float(value)) except (TypeError, ValueError): continue - x, y, _ = _interpolate_vertex(vertices, chainage) - candidates.append(StructureCandidate(chainage_m=chainage, x=x, y=y, reason="confirmed")) - return candidates + return parsed diff --git a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts index 8f6a6391..eba48bd3 100644 --- a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts +++ b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts @@ -31,20 +31,19 @@ import { createPipeEditor } from "./B05_wf2_Route_UI_Drainage_Pipes"; // 지도는 B04에서 분리한 렌더 엔진(B04_wf1_Surface_UI_MapRender)을 그대로 재사용해 // 사전 투영·LOD·뷰포트 컬링·커서 중심 줌 동작을 동일하게 얻는다. -/** 배수유역 산정의 근거가 되는 도엽 레이어. 3D는 쓰지 않는다(사용자 지시). */ -const DRAINAGE_LAYERS = ["도엽_등고선", "도엽_하천중심선", "도엽_표고점"] as const; +/** 배수유역 산정의 근거가 되는 도엽 레이어. 3D는 쓰지 않는다(사용자 지시). + * 표고점은 유효 데이터가 적어 산정에서 제외했으므로 배경에도 띄우지 않는다(2026-07-31). */ +const DRAINAGE_LAYERS = ["도엽_등고선", "도엽_하천중심선"] as const; type DrainageLayer = (typeof DRAINAGE_LAYERS)[number]; const LAYER_COLORS: Record = { 도엽_등고선: "#a5b4fc", 도엽_하천중심선: "#2563eb", - 도엽_표고점: "#f9a8d4", }; const LAYER_LABELS: Record = { 도엽_등고선: "등고선", 도엽_하천중심선: "세류", - 도엽_표고점: "표고점", }; const ROUTE_COLOR = "#f97316"; @@ -84,11 +83,13 @@ export function createDrainagePanel(): DrainagePanel { layerButtons.className = "b05-drainage__layers"; header.append(title, layerButtons); - // 유역 산정 실행 버튼 — 후보 제안·유역 산정을 한 번에 돌린다(자동 제안 + 사용자 확인 흐름). + // 배수유역 계산 실행 — 등고선 격자 흐름 해석으로 유역·관 위치를 한 번에 산정한다. + // (격자 해석 결과는 백엔드가 캐시하므로 관만 바꾼 재계산은 즉시 끝난다.) const analyzeButton = document.createElement("button"); analyzeButton.type = "button"; analyzeButton.className = "b05-drainage__analyze"; - analyzeButton.textContent = "유역 산정"; + analyzeButton.textContent = "배수유역 계산"; + analyzeButton.title = "등고선·세류선으로 배수유역과 관 매설 위치를 다시 계산합니다."; // 배관 편집 토글 — 켜면 계획선 클릭으로 배관 추가, 마커 드래그로 이동. const editButton = document.createElement("button"); editButton.type = "button"; @@ -141,8 +142,9 @@ export function createDrainagePanel(): DrainagePanel { syncPipeSelection(); scheduleDraw(); }); - // 유역 경계 외곽선 = 분수령(능선). 사용자 지시로 기본 표시. - let showRidge = true; + // 2차 전체 배수유역 외곽선 = 분수령. 별도 "능선" 레이어를 두지 않고 이 선 하나로 표시한다 + // (전체 유역 외곽선과 능선이 같은 선이므로 — 2026-07-31 사용자 지시). + let mainBoundary: Array<[number, number]> = []; let scale = 1; let offsetX = 0; let offsetY = 0; @@ -171,21 +173,6 @@ export function createDrainagePanel(): DrainagePanel { layerButtons.append(button); }); - // 능선(분수령) 표시 토글 — 유역 경계 파선. 기본 켜짐(사용자 지시). - const ridgeButton = document.createElement("button"); - ridgeButton.type = "button"; - ridgeButton.className = "b05-drainage__layer-button is-active"; - ridgeButton.textContent = "능선"; - ridgeButton.style.setProperty("--b05-layer-color", "#92400e"); - ridgeButton.setAttribute("aria-pressed", "true"); - ridgeButton.addEventListener("click", () => { - showRidge = !showRidge; - ridgeButton.classList.toggle("is-active", showRidge); - ridgeButton.setAttribute("aria-pressed", String(showRidge)); - scheduleDraw(); - }); - layerButtons.append(ridgeButton); - function updateImageTransform(): void { backgroundImage.style.transform = `translate(${offsetX}px, ${offsetY}px) scale(${scale})`; } @@ -210,7 +197,7 @@ export function createDrainagePanel(): DrainagePanel { context.clearRect(0, 0, width, height); const mapRect: MapRect = computeMapRect(meta, width, height); const view: ViewState = { width, height, scale, offsetX, offsetY, mapRect }; - // 유역 채움을 가장 아래에 깔아 등고선·세류 판독을 가리지 않게 한다. + // 세부유역 채움을 가장 아래에 깔아 등고선·세류 판독을 가리지 않게 한다. if (normalizer) { basins.forEach((basin) => { const color = BASIN_COLORS[(basin.index - 1) % BASIN_COLORS.length]; @@ -223,25 +210,26 @@ export function createDrainagePanel(): DrainagePanel { ? color : color.replace(/0\.45\)$/, "0.18)"), ); - // 유역 경계 = 분수령이므로 그 외곽선을 능선 파선으로 강조한다. - if (showRidge) drawRidgeRing(context, basin.polygon_lonlat, normalizer!, view); }); + // 전체 유역 외곽선 = 분수령(능선). 세부유역 경계와 구분되게 파선 한 겹만 얹는다. + if (mainBoundary.length > 2) drawRidgeRing(context, mainBoundary, normalizer, view); } - // 등고선을 얇게 깔고 세류·표고점을 그 위에, 노선을 맨 위에 둔다. + // 등고선을 얇게 깔고 세류를 그 위에, 노선을 맨 위에 둔다. DRAINAGE_LAYERS.forEach((layer) => { if (!activeLayers.has(layer)) return; const prepared = preparedLayers.get(layer); if (!prepared) return; context.lineWidth = layer === "도엽_등고선" ? 0.7 : 1.5; context.strokeStyle = LAYER_COLORS[layer]; - drawPreparedLayer(context, prepared, view, layer === "도엽_표고점" ? "x" : "dot"); + drawPreparedLayer(context, prepared, view, "dot"); }); if (routeLayer) { context.lineWidth = 2.4; context.strokeStyle = ROUTE_COLOR; drawPreparedLayer(context, routeLayer, view, "dot"); } - // 배관(관 매설) 마커 — 계획선 위 최상단. + // 계획선 위 흐름 강도 띠 → 그 위에 배관 마커. + pipeEditor.drawStrength(context, view); pipeEditor.draw(context, view, pipeColor); updateImageTransform(); } @@ -327,6 +315,7 @@ export function createDrainagePanel(): DrainagePanel { const chainages = !auto && pipeEditor.pipes().length > 0 ? pipeEditor.chainages() : undefined; const response = await fetchDrainageBasins(projectId, chainages); basins = response.basins; + mainBoundary = response.main_polygon_lonlat ?? []; // 실제 사용된 배관 목록으로 마커를 동기화한다(유역 없는 관 포함). pipeEditor.setPipes( (response.pipes ?? []).map((pipe) => ({ @@ -334,6 +323,7 @@ export function createDrainagePanel(): DrainagePanel { reason: pipe.reason, })), ); + pipeEditor.setStrength(response.strength_profile ?? []); selectedBasin = null; renderBasinList(); syncPipeSelection(); diff --git a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Pipes.ts b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Pipes.ts index a7ab7101..52c14e29 100644 --- a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Pipes.ts +++ b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Pipes.ts @@ -24,6 +24,8 @@ const ADD_SNAP_PX = 14; export interface PipeEditor { setContext(meta: VWorldMeta | null, points: ReadonlyArray): void; setPipes(pipes: ReadonlyArray): void; + /** 도로 위 흐름 강도 [누가거리 m, 상류 면적 ㎡]. 관 추가 판단 근거로 계획선에 덧그린다. */ + setStrength(profile: ReadonlyArray): void; pipes(): ReadonlyArray; chainages(): number[]; selected(): number | null; @@ -33,6 +35,8 @@ export interface PipeEditor { handleDown(view: ViewState, screenX: number, screenY: number, editMode: boolean): boolean; handleMove(view: ViewState, screenX: number, screenY: number): boolean; handleUp(): boolean; + /** 계획선 위 흐름 강도 띠. 마커보다 아래에 깔아야 하므로 draw()와 따로 호출한다. */ + drawStrength(context: CanvasRenderingContext2D, view: ViewState): void; draw( context: CanvasRenderingContext2D, view: ViewState, @@ -48,6 +52,9 @@ export function createPipeEditor(onChange: () => void): PipeEditor { let selectedIndex: number | null = null; let draggingIndex: number | null = null; let dragMoved = false; + let strength: Array = []; + let strengthPeak = 0; + let strengthSpan = 5; /** 화면 → 사업지 좌표계 m (MapRender affine의 역변환). */ function screenToMetric( @@ -146,6 +153,17 @@ export function createPipeEditor(onChange: () => void): PipeEditor { selectedIndex = null; draggingIndex = null; }, + setStrength(profile) { + strength = profile.map((entry) => [entry[0], entry[1]] as const); + strengthPeak = strength.reduce((peak, entry) => Math.max(peak, entry[1]), 0); + // 표본 간격은 백엔드 출력 간격을 그대로 따른다(값이 0인 구간은 빠져 있으므로 최소 간격 사용). + let span = Infinity; + for (let i = 1; i < strength.length; i += 1) { + const gap = strength[i][0] - strength[i - 1][0]; + if (gap > 0 && gap < span) span = gap; + } + strengthSpan = Number.isFinite(span) ? span : 5; + }, pipes: () => pipeList, chainages: () => pipeList.map((pipe) => Math.round(pipe.chainage_m * 100) / 100), selected: () => selectedIndex, @@ -211,6 +229,28 @@ export function createPipeEditor(onChange: () => void): PipeEditor { } return true; }, + drawStrength(context, view) { + if (strengthPeak <= 0 || route.length < 2) return; + context.save(); + context.lineCap = "butt"; + strength.forEach(([chainage, area]) => { + const from = chainageToXY(chainage); + const to = chainageToXY(Math.min(totalChainage, chainage + strengthSpan)); + if (!from || !to) return; + const start = metricToScreen(view, from.x, from.y); + const end = metricToScreen(view, to.x, to.y); + if (!start || !end) return; + // 강도는 편차가 커서(계곡 한 점에 수십 배 집중) 제곱근으로 눌러 표시한다. + const intensity = Math.sqrt(area / strengthPeak); + context.beginPath(); + context.moveTo(start.x, start.y); + context.lineTo(end.x, end.y); + context.lineWidth = 3 + 9 * intensity; + context.strokeStyle = `rgba(37, 99, 235, ${(0.15 + 0.5 * intensity).toFixed(3)})`; + context.stroke(); + }); + context.restore(); + }, draw(context, view, colorOf) { pipeList.forEach((pipe, position) => { const xy = chainageToXY(pipe.chainage_m); diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py b/B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Drainage_Watershed.py similarity index 100% rename from B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py rename to B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Drainage_Watershed.py diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Assemble.py b/B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Watershed_Assemble.py similarity index 100% rename from B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Assemble.py rename to B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Watershed_Assemble.py diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Subdivide.py b/B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Watershed_Subdivide.py similarity index 100% rename from B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Subdivide.py rename to B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Watershed_Subdivide.py diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py b/B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Watershed_Trace.py similarity index 100% rename from B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py rename to B05_wf2_Route/_legacy_watershed/B05_wf2_Route_Engine_Watershed_Trace.py diff --git a/B05_wf2_Route/_legacy_watershed/README.md b/B05_wf2_Route/_legacy_watershed/README.md new file mode 100644 index 00000000..7eab5969 --- /dev/null +++ b/B05_wf2_Route/_legacy_watershed/README.md @@ -0,0 +1,16 @@ +# _legacy_watershed (보관용, 실행 경로 아님) + +2026-07-31 배수유역 전면 재설계로 폐기된 **등고선 아크 추적 + 능선 행진** 방식 엔진 4종이다. +능선/계곡 분리가 안정적이지 않아 격자 흐름(D8 + 상류 BFS) 방식으로 교체되었다. + +| 파일 | 폐기 당시 역할 | +|---|---| +| `B05_wf2_Route_Engine_Drainage_Watershed.py` | 유역 산정 오케스트레이터 (`build_watershed_basins`) | +| `B05_wf2_Route_Engine_Watershed_Trace.py` | 등고선 아크 인덱싱·분수계 행진 | +| `B05_wf2_Route_Engine_Watershed_Assemble.py` | 아크+능선+도로선 폐합 폴리곤 조립 | +| `B05_wf2_Route_Engine_Watershed_Subdivide.py` | 메인 유역 내부 세부유역 분할 | + +**주의** +- 내용은 이동 당시 그대로이며 수정하지 않는다. 서로를 `B05_wf2_Route.B05_wf2_Route_Engine_Watershed_*` + 경로로 import하므로 이 폴더에서는 그대로 실행되지 않는다(의도된 상태 — 참고용 보관). +- 현행 엔진: `B05_wf2_Route_Engine_Watershed_Grid.py` / `_Flow.py` / `_Basin.py`. diff --git a/config/config_system.py b/config/config_system.py index 9a5bd69f..8710eb23 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -232,6 +232,49 @@ SKELETON_MAIN_RIDGE_ACC_THRESHOLD_CELLS = int( SKELETON_NODE_SPACING_M = float(os.getenv("SKELETON_NODE_SPACING_M", "10.0")) +# ───────────────────────────────────────────────────────────────────────── +# 5-3-1. 배수유역 격자 해석 파라미터 (B05 WF2 — 2026-07-31 전면 재설계) +# +# 도엽 등고선 TIN 보간 → 웅덩이 채움 → D8 물 방향 → 도로에서 상류 BFS(포인터 더블링) +# 순서로 유역을 정한다. 능선을 따로 탐지하지 않는다 — 도로로 물이 도달하는지 여부가 +# 유일한 판정 기준이며, 그 경계가 곧 능선이다. +# 라이다 DEM은 노선 주변만 커버해 유역 산정에 부족하므로 쓰지 않는다(사용자 지시). +# ───────────────────────────────────────────────────────────────────────── +# 해석 격자 한 변(m). 작을수록 정밀하나 셀 수가 제곱으로 늘어난다. +DRAINAGE_GRID_SIZE_M = float(os.getenv("DRAINAGE_GRID_SIZE_M", "1.0")) +# 1차 배수유역 반경(m). 정리된 세류선과 노선을 이 반경으로 버퍼해 초기 해석 범위를 잡는다. +DRAINAGE_INITIAL_RADIUS_M = float(os.getenv("DRAINAGE_INITIAL_RADIUS_M", "300.0")) +# 활성 셀이 격자 최외곽에 닿았을 때 한 번에 넓히는 폭(m). +DRAINAGE_EXPAND_STEP_M = float(os.getenv("DRAINAGE_EXPAND_STEP_M", "200.0")) +# 확장 반복 상한. 경계 링이 전부 비활성이 되면 그 전에 스스로 멈춘다(안전핀). +DRAINAGE_MAX_EXPAND_ROUNDS = int(os.getenv("DRAINAGE_MAX_EXPAND_ROUNDS", "6")) +# 격자 셀 수 상한. 초과하면 셀 크기를 자동으로 키워 맞춘다(메모리 보호). +DRAINAGE_MAX_GRID_CELLS = int(os.getenv("DRAINAGE_MAX_GRID_CELLS", "16000000")) +# 도로 폭(m). 이 폭으로 노선을 격자에 구워 D8 흐름이 도로를 대각선으로 건너뛰지 못하게 한다. +DRAINAGE_ROAD_WIDTH_M = float(os.getenv("DRAINAGE_ROAD_WIDTH_M", "4.0")) +# 계획선 최저점보다 이만큼 아래인 등고선은 상류 기여가 불가능하므로 보간에서 제외한다. +DRAINAGE_CONTOUR_MARGIN_M = float(os.getenv("DRAINAGE_CONTOUR_MARGIN_M", "10.0")) +# 이보다 짧은 등고선 파편은 노이즈로 보고 버린다. 봉우리 폐합 등고선은 이 값 이상이면 남는다. +DRAINAGE_CONTOUR_MIN_LENGTH_M = float(os.getenv("DRAINAGE_CONTOUR_MIN_LENGTH_M", "20.0")) +# 등고선 정점 재샘플 간격(m). 조밀할수록 TIN이 정확하나 Delaunay 비용이 커진다. +DRAINAGE_CONTOUR_RESAMPLE_M = float(os.getenv("DRAINAGE_CONTOUR_RESAMPLE_M", "5.0")) +# 평탄면 해소용 미세 경사(m/셀). 채움 후 흐름 방향이 없는 셀에 출구 쪽 경사를 만들어 준다. +DRAINAGE_FLAT_EPSILON_M = float(os.getenv("DRAINAGE_FLAT_EPSILON_M", "0.001")) +# 관 매설 최대 간격(m). 이 간격을 넘으면 흐름 강도가 가장 큰 지점에 관을 보충한다. +DRAINAGE_PIPE_MAX_SPACING_M = float(os.getenv("DRAINAGE_PIPE_MAX_SPACING_M", "300.0")) +# 관끼리 이보다 가까우면 같은 계곡으로 보고 하나로 합친다. +DRAINAGE_PIPE_MIN_SPACING_M = float(os.getenv("DRAINAGE_PIPE_MIN_SPACING_M", "20.0")) +# 측구 흐름(도로 셀 → 담당 관) 판정용 종단 계획선 샘플 간격(m). +DRAINAGE_DITCH_SAMPLE_M = float(os.getenv("DRAINAGE_DITCH_SAMPLE_M", "1.0")) +# 유역 폴리곤 단순화 허용오차(m). 격자 계단 경계를 매끄럽게 줄여 응답 크기를 낮춘다. +DRAINAGE_POLYGON_SIMPLIFY_M = float(os.getenv("DRAINAGE_POLYGON_SIMPLIFY_M", "2.0")) +# 이 면적(㎡) 미만의 유역 조각은 버린다(격자 노이즈 제거). +DRAINAGE_MIN_BASIN_AREA_M2 = float(os.getenv("DRAINAGE_MIN_BASIN_AREA_M2", "100.0")) +# 격자 해석 결과 캐시 파일명. 프로젝트 저장소의 B05_wf2_Route/drainage/ 아래에 놓인다. +DRAINAGE_CACHE_DIRNAME = "drainage" +DRAINAGE_CACHE_FILENAME = "watershed_grid.npz" + + # ───────────────────────────────────────────────────────────────────────── # 5-4. 종횡단 생성 파라미터 (B06 WF3) # ───────────────────────────────────────────────────────────────────────── diff --git a/pyproject.toml b/pyproject.toml index e3b1d9df..831e3b93 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,8 @@ [tool.ruff] target-version = "py313" line-length = 100 +# 폐기 엔진 보관 폴더는 이동 당시 원본 그대로 두기로 했으므로 린트/포맷 대상에서 뺀다. +extend-exclude = ["B05_wf2_Route/_legacy_watershed"] [tool.ruff.lint] select = ["E", "F", "I"]