From a77ee317a58afefd6c474145a9de1dd398d3ac2f Mon Sep 17 00:00:00 2001 From: umsangdon Date: Fri, 31 Jul 2026 17:57:13 +0900 Subject: [PATCH] =?UTF-8?q?refactor(B05):=20=EC=84=B8=EB=A5=98=EB=A7=9D?= =?UTF-8?q?=C2=B71=EC=B0=A8=EC=98=81=EC=97=AD=EC=9D=84=20Watershed=5FStrea?= =?UTF-8?q?m=20=EB=AA=A8=EB=93=88=EB=A1=9C=20=EB=B6=84=EB=A6=AC=20(700?= =?UTF-8?q?=EC=A4=84=20=EC=A0=9C=ED=95=9C)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Grid.py 가 708줄로 제한을 넘어 관심사 단위로 갈랐다. - Watershed_Stream.py (292줄): 세류망 노딩·도로 절단·연결망 확산, ElevationSampler, PrimaryRegion, build_primary_region - Watershed_Grid.py (427줄): 등고선 구름, 격자 규격, TIN 보간, 웅덩이 채움, 평탄면 해소, D8 - _iter_linestrings -> iter_linestrings, _spec_from_bounds -> grid_spec_from_bounds 로 공개(모듈 간 재사용). 의존 방향은 Stream -> Grid 단방향. - 죽은 코드 제거: select_upstream_streams, _point_elevation 회귀 확인: 합성 유역 179,919㎡ 동일, 연결망 판정 동일. Co-Authored-By: Claude Fable 5 --- .../B05_wf2_Route_Engine_Watershed_Basin.py | 6 +- .../B05_wf2_Route_Engine_Watershed_Grid.py | 295 +----------------- .../B05_wf2_Route_Engine_Watershed_Stream.py | 292 +++++++++++++++++ 3 files changed, 303 insertions(+), 290 deletions(-) create mode 100644 B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py index d98a9fea..452ac168 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py @@ -45,13 +45,15 @@ from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Flow import ( ) from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import ( GridSpec, - PrimaryRegion, build_contour_cloud, - build_primary_region, build_terrain_grid, expand_grid_spec, route_elevation_floor, ) +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Stream import ( + PrimaryRegion, + build_primary_region, +) from config.config_system import ( DRAINAGE_DITCH_SAMPLE_M, DRAINAGE_EXPAND_STEP_M, diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py index 7f30333e..4d49bc2f 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py @@ -19,16 +19,14 @@ from __future__ import annotations import logging import math -from dataclasses import dataclass, field +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, MultiPolygon, Polygon, shape -from shapely.ops import substring, unary_union +from shapely.geometry import LineString, shape from skimage.morphology import reconstruction from config.config_system import ( @@ -37,7 +35,6 @@ from config.config_system import ( DRAINAGE_CONTOUR_MIN_LENGTH_M, DRAINAGE_CONTOUR_RESAMPLE_M, DRAINAGE_FLAT_EPSILON_M, - DRAINAGE_GRID_SIZE_M, DRAINAGE_MAX_GRID_CELLS, ) @@ -131,13 +128,13 @@ def _feature_elevation(properties: dict[str, Any]) -> float | None: return None -def _iter_linestrings(geometry: Any) -> list[LineString]: +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)) + lines.extend(iter_linestrings(part)) return lines return [] @@ -179,7 +176,7 @@ def build_contour_cloud( if clip_bounds is not None and _outside_bounds(parsed.bounds, clip_bounds): dropped_outside += 1 continue - for line in _iter_linestrings(parsed): + for line in iter_linestrings(parsed): if line.length < DRAINAGE_CONTOUR_MIN_LENGTH_M: dropped_short += 1 continue @@ -215,285 +212,7 @@ def _outside_bounds( return bounds[2] < clip[0] or bounds[0] > clip[2] or bounds[3] < clip[1] or bounds[1] > clip[3] -# ── ② 세류선 상류측만 남기기 ──────────────────────────────────────────────── - - -def select_upstream_streams( - route_line: LineString, - stream_features: list[dict[str, Any]], - cloud: ContourCloud, -) -> list[LineString]: - """도로 교차점 기준 상류측으로 이어진 세류 연결망 전체를 돌려준다.""" - return split_streams_at_road(route_line, stream_features, cloud).upstream - - -@dataclass -class StreamSplit: - """세류망을 도로 교차점에서 상·하류로 가른 결과. 검증 화면이 그대로 그린다.""" - - upstream: list[LineString] = field(default_factory=list) # 채택 — 1차 영역의 기준 - downstream: list[LineString] = field(default_factory=list) # 도로 아래로 이어진 망 - no_contact: int = 0 # 어느 쪽에도 이어지지 않아 제외한 조각 수 - - -def split_streams_at_road( - route_line: LineString, - stream_features: list[dict[str, Any]], - cloud: ContourCloud, -) -> StreamSplit: - """세류망을 도로에서 끊고, 교차점 상류측으로 **이어진 망 전체**를 채택한다. - - 도엽 하천중심선은 피처가 잘게 쪼개져 있고 지류가 본류 중간에 T자로 붙는다. 그래서 - 피처 단위로 보면 상류망이 통째로 빠진다. 순서를 이렇게 잡는다: - - ① 세류선끼리 `unary_union`으로 노딩 — 중간에서 만나는 지류도 연결로 인식된다 - ② 도로 교차점에서 한 번 더 잘라 상·하류 조각을 물리적으로 분리한다 - ③ 끝점 그래프를 만들고, **도로 교차 노드는 통과하지 못하게** 막는다 - ④ 도로에 접한 조각을 교차점 표고와 비교해 상·하류 씨앗으로 정한다 - ⑤ 씨앗에서 퍼뜨려 이어진 망 전체를 채택 — 도로를 넘어가지 못하므로 상·하류가 섞이지 않는다 - """ - lines: list[LineString] = [] - for feature in stream_features: - geometry = feature.get("geometry") - if not geometry: - continue - try: - parsed = shape(geometry) - except Exception: # noqa: BLE001 - continue - lines.extend(line for line in _iter_linestrings(parsed) if line.length > 0) - if not lines: - return StreamSplit() - - pieces, crossing_nodes = _cut_network_at_road(lines, route_line) - if not pieces: - return StreamSplit() - - node_edges: dict[tuple[float, float], list[int]] = {} - ends: list[tuple[tuple[float, float], tuple[float, float]]] = [] - for index, piece in enumerate(pieces): - head = _node_key(*piece.coords[0]) - tail = _node_key(*piece.coords[-1]) - ends.append((head, tail)) - node_edges.setdefault(head, []).append(index) - node_edges.setdefault(tail, []).append(index) - - sampler = ElevationSampler(cloud) - upper_seeds: set[int] = set() - lower_seeds: set[int] = set() - for index, piece in enumerate(pieces): - touching = [node for node in ends[index] if node in crossing_nodes] - if not touching: - continue - crossing_z = float(np.min(sampler.at(np.array(touching, dtype=np.float64)))) - if _mean_elevation(piece, sampler) > crossing_z: - upper_seeds.add(index) - else: - lower_seeds.add(index) - - upstream = _spread_network(upper_seeds, ends, node_edges, crossing_nodes) - downstream = _spread_network(lower_seeds, ends, node_edges, crossing_nodes) - upstream - logger.info( - "배수유역: 세류 조각 %d개 → 상류망 %d개 채택 / 하류망 %d개 · 미연결 %d개 제외", - len(pieces), - len(upstream), - len(downstream), - len(pieces) - len(upstream) - len(downstream), - ) - return StreamSplit( - upstream=[pieces[index] for index in sorted(upstream)], - downstream=[pieces[index] for index in sorted(downstream)], - no_contact=len(pieces) - len(upstream) - len(downstream), - ) - - -def _cut_network_at_road( - lines: list[LineString], route_line: LineString -) -> tuple[list[LineString], set[tuple[float, float]]]: - """세류망을 노딩한 뒤 도로 교차점에서 자르고, 그 교차 노드를 함께 돌려준다.""" - noded = unary_union(lines) - pieces: list[LineString] = [] - crossing_nodes: set[tuple[float, float]] = set() - for piece in _iter_linestrings(noded): - if not piece.intersects(route_line): - pieces.append(piece) - continue - hits = _intersection_points(piece.intersection(route_line)) - positions = sorted( - { - position - for position in (piece.project(point) for point in hits) - if 0.0 < position < piece.length - } - ) - for point in hits: - crossing_nodes.add(_node_key(point.x, point.y)) - if not positions: - # 끝점이 도로에 닿은 경우 — 자를 필요는 없고 그 끝점이 곧 교차 노드다. - pieces.append(piece) - continue - bounds = [0.0, *positions, piece.length] - for start, end in zip(bounds, bounds[1:]): - if end - start <= 0: - continue - cut = _substring(piece, start, end) - if cut is not None: - pieces.append(cut) - return pieces, crossing_nodes - - -def _spread_network( - seeds: set[int], - ends: list[tuple[tuple[float, float], tuple[float, float]]], - node_edges: dict[tuple[float, float], list[int]], - blocked: set[tuple[float, float]], -) -> set[int]: - """씨앗 조각에서 끝점을 타고 퍼진다. 도로 교차 노드는 통과하지 않는다.""" - reached = set(seeds) - queue = list(seeds) - while queue: - index = queue.pop() - for node in ends[index]: - if node in blocked: - continue - for neighbour in node_edges.get(node, ()): - if neighbour not in reached: - reached.add(neighbour) - queue.append(neighbour) - return reached - - -def _node_key(x: float, y: float) -> tuple[float, float]: - """끝점 일치 판정용 좌표 키. 노딩 후에도 부동소수 오차가 남아 mm로 반올림한다.""" - return (round(float(x), 3), round(float(y), 3)) - - -class ElevationSampler: - """등고선 구름에서 임의 지점 표고를 읽는다 — 상·하류 판정 전용. - - 최근접 등고선 정점만 쓰면 오차가 등고선 간격(주곡선 5m)만큼 나서, 계곡 교차점 표고가 - 실제보다 한 등고선 위로 잡히고 상류 조각이 통째로 하류로 오판된다. TIN 선형보간을 - 1차로 쓰고, TIN 밖(볼록껍질 외부)만 최근접 정점으로 메운다. - """ - - def __init__(self, cloud: ContourCloud) -> None: - self._z = cloud.z - if cloud.is_empty: - self._interpolator = None - self._tree = None - return - self._interpolator = LinearNDInterpolator(cloud.xy, cloud.z) - self._tree = cKDTree(cloud.xy) - - def at(self, xy: np.ndarray) -> np.ndarray: - """(N, 2) 좌표의 표고 (N,).""" - if self._interpolator is None or self._tree is None: - return np.zeros(xy.shape[0]) - values = np.asarray(self._interpolator(xy), dtype=np.float64) - missing = ~np.isfinite(values) - if missing.any(): - _, indices = self._tree.query(xy[missing]) - values[missing] = self._z[indices] - return values - - -def _point_elevation(point: Any, tree: cKDTree | None, cloud: ContourCloud) -> float: - """가장 가까운 등고선 정점의 표고. TIN이 없는 단계의 근사 표고다.""" - if tree is None: - return 0.0 - _, index = tree.query([[point.x, point.y]]) - return float(cloud.z[index[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, sampler: ElevationSampler) -> float: - """선을 10m 간격으로 훑은 평균 표고.""" - 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]) - return float(np.mean(sampler.at(points))) - - -# ── ③ 격자 범위 ───────────────────────────────────────────────────────────── - - -@dataclass -class PrimaryRegion: - """1차 배수유역 — 격자 범위의 근거. 검증 화면이 이 내용을 그대로 그린다.""" - - split: StreamSplit - # 상류 세류망을 반경 버퍼해 합친 영역. 노선은 버퍼하지 않는다. - area: Polygon | MultiPolygon | None - spec: GridSpec - radius_m: float - # 1차 영역 밖으로 나간 노선 길이(m). 그 구간 사면은 해석에서 빠진다는 경고 지표. - road_outside_m: float = 0.0 - - -def build_primary_region( - route_line: LineString, - stream_features: list[dict[str, Any]], - cloud: ContourCloud, - radius_m: float, - cell_m: float = DRAINAGE_GRID_SIZE_M, -) -> PrimaryRegion: - """**상류 세류망 + 계획 노선**을 반경 버퍼해 합친 범위 = 1차 배수유역, bbox = 해석 격자. - - 노선 버퍼는 세류 교차가 없는 구간의 도로도 격자 안에 들어오게 한다 — 그래야 그 구간 - 사면이 유역으로 잡힌다. 상류 세류망 선정이 정확해진 뒤 다시 포함했다(2026-07-31 사용자 지시). - - 노선이 이 영역 밖으로 나가는 길이는 따로 재서 남긴다 — 그 구간은 도로 셀이 격자에 - 없어 유역이 잡히지 않으므로 반경을 올릴지 판단하는 근거가 된다. - """ - split = split_streams_at_road(route_line, stream_features, cloud) - geometries = [route_line.buffer(radius_m)] - geometries.extend(line.buffer(radius_m) for line in split.upstream) - if not split.upstream: - logger.warning("배수유역: 상류 세류망이 없어 노선 버퍼만으로 1차 영역을 잡습니다.") - area = unary_union(geometries) - x_min, y_min, x_max, y_max = area.bounds - spec = _spec_from_bounds(x_min, y_min, x_max, y_max, cell_m) - outside = route_line.difference(area) - road_outside_m = float(outside.length) if not outside.is_empty else 0.0 - logger.info( - "배수유역: 1차 영역 %.0f㎡, 범위 %.0fm × %.0fm, 격자 %d×%d (%.2fm), 노선 이탈 %.0fm/%.0fm", - area.area, - x_max - x_min, - y_max - y_min, - spec.n_rows, - spec.n_cols, - spec.cell_m, - road_outside_m, - route_line.length, - ) - return PrimaryRegion( - split=split, area=area, spec=spec, radius_m=radius_m, road_outside_m=road_outside_m - ) - - -def _spec_from_bounds( +def grid_spec_from_bounds( x_min: float, y_min: float, x_max: float, y_max: float, cell_m: float ) -> GridSpec: """격자 크기는 절대 자동으로 바꾸지 않는다 — config 값이 그대로 쓰인다(사용자 지시). @@ -528,7 +247,7 @@ def expand_grid_spec(spec: GridSpec, sides: dict[str, bool], step_m: float) -> G 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) + return grid_spec_from_bounds(x_min, y_min, x_max, y_max, spec.cell_m) # ── ④ TIN 보간 ────────────────────────────────────────────────────────────── diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py new file mode 100644 index 00000000..cb629dea --- /dev/null +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py @@ -0,0 +1,292 @@ +"""세류망 상·하류 분리와 1차 배수유역 산정. + +도엽 하천중심선은 피처가 잘게 쪼개져 있고 지류가 본류 중간에 T자로 붙는다. 피처 단위로 +자르면 상류망이 통째로 빠지므로, 노딩 → 도로 절단 → 끝점 그래프 확산으로 **이어진 망 +전체**를 잡는다(2026-07-31 사용자 지시). + +여기서 정해진 1차 배수유역의 bbox가 곧 격자 해석 범위가 된다. +표고 해석·격자 생성은 `B05_wf2_Route_Engine_Watershed_Grid.py`가 맡는다. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass, field +from typing import Any + +import numpy as np +from scipy.interpolate import LinearNDInterpolator +from scipy.spatial import cKDTree +from shapely.geometry import LineString, MultiPolygon, Polygon, shape +from shapely.ops import substring, unary_union + +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import ( + ContourCloud, + GridSpec, + grid_spec_from_bounds, + iter_linestrings, +) +from config.config_system import DRAINAGE_GRID_SIZE_M + +logger = logging.getLogger(__name__) + + +# ── 세류망 상·하류 분리 ───────────────────────────────────────────────────── + + +@dataclass +class StreamSplit: + """세류망을 도로 교차점에서 상·하류로 가른 결과. 검증 화면이 그대로 그린다.""" + + upstream: list[LineString] = field(default_factory=list) # 채택 — 1차 영역의 기준 + downstream: list[LineString] = field(default_factory=list) # 도로 아래로 이어진 망 + no_contact: int = 0 # 어느 쪽에도 이어지지 않아 제외한 조각 수 + + +def split_streams_at_road( + route_line: LineString, + stream_features: list[dict[str, Any]], + cloud: ContourCloud, +) -> StreamSplit: + """세류망을 도로에서 끊고, 교차점 상류측으로 **이어진 망 전체**를 채택한다. + + 도엽 하천중심선은 피처가 잘게 쪼개져 있고 지류가 본류 중간에 T자로 붙는다. 그래서 + 피처 단위로 보면 상류망이 통째로 빠진다. 순서를 이렇게 잡는다: + + ① 세류선끼리 `unary_union`으로 노딩 — 중간에서 만나는 지류도 연결로 인식된다 + ② 도로 교차점에서 한 번 더 잘라 상·하류 조각을 물리적으로 분리한다 + ③ 끝점 그래프를 만들고, **도로 교차 노드는 통과하지 못하게** 막는다 + ④ 도로에 접한 조각을 교차점 표고와 비교해 상·하류 씨앗으로 정한다 + ⑤ 씨앗에서 퍼뜨려 이어진 망 전체를 채택 — 도로를 넘어가지 못하므로 상·하류가 섞이지 않는다 + """ + lines: list[LineString] = [] + for feature in stream_features: + geometry = feature.get("geometry") + if not geometry: + continue + try: + parsed = shape(geometry) + except Exception: # noqa: BLE001 + continue + lines.extend(line for line in iter_linestrings(parsed) if line.length > 0) + if not lines: + return StreamSplit() + + pieces, crossing_nodes = _cut_network_at_road(lines, route_line) + if not pieces: + return StreamSplit() + + node_edges: dict[tuple[float, float], list[int]] = {} + ends: list[tuple[tuple[float, float], tuple[float, float]]] = [] + for index, piece in enumerate(pieces): + head = _node_key(*piece.coords[0]) + tail = _node_key(*piece.coords[-1]) + ends.append((head, tail)) + node_edges.setdefault(head, []).append(index) + node_edges.setdefault(tail, []).append(index) + + sampler = ElevationSampler(cloud) + upper_seeds: set[int] = set() + lower_seeds: set[int] = set() + for index, piece in enumerate(pieces): + touching = [node for node in ends[index] if node in crossing_nodes] + if not touching: + continue + crossing_z = float(np.min(sampler.at(np.array(touching, dtype=np.float64)))) + if _mean_elevation(piece, sampler) > crossing_z: + upper_seeds.add(index) + else: + lower_seeds.add(index) + + upstream = _spread_network(upper_seeds, ends, node_edges, crossing_nodes) + downstream = _spread_network(lower_seeds, ends, node_edges, crossing_nodes) - upstream + logger.info( + "배수유역: 세류 조각 %d개 → 상류망 %d개 채택 / 하류망 %d개 · 미연결 %d개 제외", + len(pieces), + len(upstream), + len(downstream), + len(pieces) - len(upstream) - len(downstream), + ) + return StreamSplit( + upstream=[pieces[index] for index in sorted(upstream)], + downstream=[pieces[index] for index in sorted(downstream)], + no_contact=len(pieces) - len(upstream) - len(downstream), + ) + + +def _cut_network_at_road( + lines: list[LineString], route_line: LineString +) -> tuple[list[LineString], set[tuple[float, float]]]: + """세류망을 노딩한 뒤 도로 교차점에서 자르고, 그 교차 노드를 함께 돌려준다.""" + noded = unary_union(lines) + pieces: list[LineString] = [] + crossing_nodes: set[tuple[float, float]] = set() + for piece in iter_linestrings(noded): + if not piece.intersects(route_line): + pieces.append(piece) + continue + hits = _intersection_points(piece.intersection(route_line)) + positions = sorted( + { + position + for position in (piece.project(point) for point in hits) + if 0.0 < position < piece.length + } + ) + for point in hits: + crossing_nodes.add(_node_key(point.x, point.y)) + if not positions: + # 끝점이 도로에 닿은 경우 — 자를 필요는 없고 그 끝점이 곧 교차 노드다. + pieces.append(piece) + continue + bounds = [0.0, *positions, piece.length] + for start, end in zip(bounds, bounds[1:]): + if end - start <= 0: + continue + cut = _substring(piece, start, end) + if cut is not None: + pieces.append(cut) + return pieces, crossing_nodes + + +def _spread_network( + seeds: set[int], + ends: list[tuple[tuple[float, float], tuple[float, float]]], + node_edges: dict[tuple[float, float], list[int]], + blocked: set[tuple[float, float]], +) -> set[int]: + """씨앗 조각에서 끝점을 타고 퍼진다. 도로 교차 노드는 통과하지 않는다.""" + reached = set(seeds) + queue = list(seeds) + while queue: + index = queue.pop() + for node in ends[index]: + if node in blocked: + continue + for neighbour in node_edges.get(node, ()): + if neighbour not in reached: + reached.add(neighbour) + queue.append(neighbour) + return reached + + +def _node_key(x: float, y: float) -> tuple[float, float]: + """끝점 일치 판정용 좌표 키. 노딩 후에도 부동소수 오차가 남아 mm로 반올림한다.""" + return (round(float(x), 3), round(float(y), 3)) + + +class ElevationSampler: + """등고선 구름에서 임의 지점 표고를 읽는다 — 상·하류 판정 전용. + + 최근접 등고선 정점만 쓰면 오차가 등고선 간격(주곡선 5m)만큼 나서, 계곡 교차점 표고가 + 실제보다 한 등고선 위로 잡히고 상류 조각이 통째로 하류로 오판된다. TIN 선형보간을 + 1차로 쓰고, TIN 밖(볼록껍질 외부)만 최근접 정점으로 메운다. + """ + + def __init__(self, cloud: ContourCloud) -> None: + self._z = cloud.z + if cloud.is_empty: + self._interpolator = None + self._tree = None + return + self._interpolator = LinearNDInterpolator(cloud.xy, cloud.z) + self._tree = cKDTree(cloud.xy) + + def at(self, xy: np.ndarray) -> np.ndarray: + """(N, 2) 좌표의 표고 (N,).""" + if self._interpolator is None or self._tree is None: + return np.zeros(xy.shape[0]) + values = np.asarray(self._interpolator(xy), dtype=np.float64) + missing = ~np.isfinite(values) + if missing.any(): + _, indices = self._tree.query(xy[missing]) + values[missing] = self._z[indices] + return values + + +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, sampler: ElevationSampler) -> float: + """선을 10m 간격으로 훑은 평균 표고.""" + 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]) + return float(np.mean(sampler.at(points))) + + +# ── 1차 배수유역 ──────────────────────────────────────────────────────────── + + +@dataclass +class PrimaryRegion: + """1차 배수유역 — 격자 범위의 근거. 검증 화면이 이 내용을 그대로 그린다.""" + + split: StreamSplit + # 상류 세류망을 반경 버퍼해 합친 영역. 노선은 버퍼하지 않는다. + area: Polygon | MultiPolygon | None + spec: GridSpec + radius_m: float + # 1차 영역 밖으로 나간 노선 길이(m). 그 구간 사면은 해석에서 빠진다는 경고 지표. + road_outside_m: float = 0.0 + + +def build_primary_region( + route_line: LineString, + stream_features: list[dict[str, Any]], + cloud: ContourCloud, + radius_m: float, + cell_m: float = DRAINAGE_GRID_SIZE_M, +) -> PrimaryRegion: + """**상류 세류망 + 계획 노선**을 반경 버퍼해 합친 범위 = 1차 배수유역, bbox = 해석 격자. + + 노선 버퍼는 세류 교차가 없는 구간의 도로도 격자 안에 들어오게 한다 — 그래야 그 구간 + 사면이 유역으로 잡힌다. 상류 세류망 선정이 정확해진 뒤 다시 포함했다(2026-07-31 사용자 지시). + + 노선이 이 영역 밖으로 나가는 길이는 따로 재서 남긴다 — 그 구간은 도로 셀이 격자에 + 없어 유역이 잡히지 않으므로 반경을 올릴지 판단하는 근거가 된다. + """ + split = split_streams_at_road(route_line, stream_features, cloud) + geometries = [route_line.buffer(radius_m)] + geometries.extend(line.buffer(radius_m) for line in split.upstream) + if not split.upstream: + logger.warning("배수유역: 상류 세류망이 없어 노선 버퍼만으로 1차 영역을 잡습니다.") + area = unary_union(geometries) + x_min, y_min, x_max, y_max = area.bounds + spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, cell_m) + outside = route_line.difference(area) + road_outside_m = float(outside.length) if not outside.is_empty else 0.0 + logger.info( + "배수유역: 1차 영역 %.0f㎡, 범위 %.0fm × %.0fm, 격자 %d×%d (%.2fm), 노선 이탈 %.0fm/%.0fm", + area.area, + x_max - x_min, + y_max - y_min, + spec.n_rows, + spec.n_cols, + spec.cell_m, + road_outside_m, + route_line.length, + ) + return PrimaryRegion( + split=split, area=area, spec=spec, radius_m=radius_m, road_outside_m=road_outside_m + )