diff --git a/B04_wf1_Surface/B04_wf1_Surface_UI_MapRender.ts b/B04_wf1_Surface/B04_wf1_Surface_UI_MapRender.ts index 3892c42e..ad2a9efe 100644 --- a/B04_wf1_Surface/B04_wf1_Surface_UI_MapRender.ts +++ b/B04_wf1_Surface/B04_wf1_Surface_UI_MapRender.ts @@ -521,3 +521,30 @@ export function drawFilledRing( context.fillStyle = "#1f2937"; context.fillText(entry.label, centerX, centerY); } + +/** 유역 경계(분수령=능선)를 능선 스타일(갈색 파선)로 강조해 그린다. */ +export function drawRidgeRing( + context: CanvasRenderingContext2D, + ring: ReadonlyArray, + normalizer: Normalizer, + view: ViewState, +): void { + if (ring.length < 3) return; + const affine = affineOf(view); + context.beginPath(); + ring.forEach(([lon, lat], index) => { + const nx = (lon - normalizer.lonMin) / normalizer.lonRange; + const ny = 1 - (lat - normalizer.latMin) / normalizer.latRange; + const x = nx * affine.ax + affine.bx; + const y = ny * affine.ay + affine.by; + if (index === 0) context.moveTo(x, y); + else context.lineTo(x, y); + }); + context.closePath(); + context.save(); + context.strokeStyle = "#92400e"; + context.lineWidth = 1.8; + context.setLineDash([7, 4]); + context.stroke(); + context.restore(); +} diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Basin.py b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Basin.py deleted file mode 100644 index 854dd34f..00000000 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Basin.py +++ /dev/null @@ -1,257 +0,0 @@ -"""배수유역 경계 산정. - -구조물 측점(관 매설 지점)에서 산정상부까지 역추적해 밀폐된 유역 경계를 만든다. -지형 판단 근거는 도엽 등고선·세류선·표고점뿐이다(3D 미사용, 2026-07-28 사용자 지시). - -정상부 판정 규칙(사용자 지시): - 표고점 데이터는 산 정상부가 아닌 경우가 많다. 따라서 **등고선의 동심 폐합 패턴** - (안쪽으로 갈수록 표고가 높아지는 폐합 등고선의 최내곽)으로 정상부를 먼저 판단하고, - 표고점은 그 판정을 보조·검증하는 용도로만 쓴다. -""" - -from __future__ import annotations - -import logging -import math -from typing import Any - -from shapely.geometry import LineString, Point, Polygon, shape -from shapely.ops import unary_union - -from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( - MAX_BASIN_RADIUS_M, - DrainageBasin, - RouteVertex, - StructureCandidate, - estimate_pipe_diameter_mm, -) - -logger = logging.getLogger(__name__) - -# 등고선을 폐합으로 볼 때 허용하는 시종점 이격(m). 도엽 경계에서 잘린 선을 걸러낸다. -CLOSED_TOLERANCE_M = 1.0 -# 폐합 등고선이 정상부 후보가 되는 최대 둘레(m). 이보다 크면 산체 전체라 정상부로 보지 않는다. -MAX_SUMMIT_PERIMETER_M = 1200.0 - - -class ContourField: - """도엽 등고선 피처 모음을 표고 조회·정상부 판정에 쓸 수 있게 감싼 것.""" - - def __init__(self, features: list[dict[str, Any]], elevation_keys: tuple[str, ...]): - self.lines: list[tuple[LineString, float]] = [] - self.closed: list[tuple[Polygon, float]] = [] - for feature in features: - elevation = _read_elevation(feature.get("properties") or {}, elevation_keys) - if elevation is None: - continue - geometry = feature.get("geometry") - if not geometry: - continue - try: - geom = shape(geometry) - except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다 - continue - for line in _iter_lines(geom): - self.lines.append((line, elevation)) - polygon = _as_closed_polygon(line) - if polygon is not None: - self.closed.append((polygon, elevation)) - - def elevation_at(self, x: float, y: float, radius_m: float = 200.0) -> float | None: - """가장 가까운 등고선의 표고를 그 지점의 표고로 본다.""" - point = Point(x, y) - best: tuple[float, float] | None = None - for line, elevation in self.lines: - distance = line.distance(point) - if distance > radius_m: - continue - if best is None or distance < best[0]: - best = (distance, elevation) - return None if best is None else best[1] - - def find_summits(self, near: Point, radius_m: float) -> list[tuple[Polygon, float]]: - """주변의 정상부 후보를 찾는다. - - 동심 폐합 등고선 중 **자기보다 높은 폐합 등고선을 안에 품지 않은 것**이 최내곽, - 즉 정상부다. 도엽 경계에서 잘린 선과 산체 전체를 감싸는 큰 폐합은 제외한다. - """ - nearby = [ - (polygon, elevation) - for polygon, elevation in self.closed - if polygon.length <= MAX_SUMMIT_PERIMETER_M and polygon.distance(near) <= radius_m - ] - summits: list[tuple[Polygon, float]] = [] - for polygon, elevation in nearby: - has_higher_inside = any( - other_elevation > elevation and polygon.contains(other.representative_point()) - for other, other_elevation in nearby - if other is not polygon - ) - if not has_higher_inside: - summits.append((polygon, elevation)) - return summits - - -def _read_elevation(properties: dict[str, Any], keys: tuple[str, ...]) -> float | None: - for key in keys: - value = properties.get(key) - if value is None: - continue - try: - return float(value) - except (TypeError, ValueError): - continue - return None - - -def _iter_lines(geometry: Any) -> list[LineString]: - if geometry.geom_type == "LineString": - return [geometry] - if geometry.geom_type == "MultiLineString": - return list(geometry.geoms) - return [] - - -def _as_closed_polygon(line: LineString) -> Polygon | None: - coords = list(line.coords) - if len(coords) < 4: - return None - if math.dist(coords[0], coords[-1]) > CLOSED_TOLERANCE_M: - return None - try: - polygon = Polygon(coords) - except Exception: # noqa: BLE001 - return None - return polygon if polygon.is_valid and polygon.area > 0 else None - - -def build_basins( - vertices: list[RouteVertex], - candidates: list[StructureCandidate], - contours: ContourField, - stream_features: list[dict[str, Any]], - to_lonlat: Any, -) -> list[DrainageBasin]: - """확정된 구조물 측점별 배수유역을 만든다. - - 측점에 물을 보내는 세류 가지를 따라 위로 올라가 정상부 폐합 등고선까지 닿는 범위를 - 유역으로 본다. 정상부는 ContourField.find_summits가 등고선 폐합 패턴으로 판정한다. - 번호는 노선 시점에 가까운 순서(측점 누가거리 오름차순)로 1부터 매긴다. - """ - route_line = ( - LineString([(vertex.x, vertex.y) for vertex in vertices]) if len(vertices) > 1 else None - ) - streams = _stream_lines(stream_features) - basins: list[DrainageBasin] = [] - for index, candidate in enumerate( - sorted(candidates, key=lambda item: item.chainage_m), start=1 - ): - outlet = Point(candidate.x, candidate.y) - uphill = _uphill_streams(outlet, streams, contours) - summits = contours.find_summits(outlet, MAX_BASIN_RADIUS_M) - boundary = _basin_polygon(outlet, uphill, summits, route_line) - if boundary is None or boundary.is_empty: - continue - outlet_elevation = contours.elevation_at(candidate.x, candidate.y) or 0.0 - top_elevation = max((elevation for _, elevation in summits), default=outlet_elevation) - basin = DrainageBasin( - index=index, - chainage_m=candidate.chainage_m, - outlet_x=candidate.x, - outlet_y=candidate.y, - polygon_lonlat=[list(to_lonlat(x, y)) for x, y in boundary.exterior.coords], - area_m2=float(boundary.area), - relief_m=float(max(0.0, top_elevation - outlet_elevation)), - flow_length_m=_flow_length(outlet, uphill, boundary), - ) - basin.pipe_diameter_mm = estimate_pipe_diameter_mm( - basin.area_m2, basin.relief_m, basin.flow_length_m - ) - basins.append(basin) - return basins - - -def _stream_lines(features: list[dict[str, Any]]) -> list[LineString]: - lines: list[LineString] = [] - for feature in features: - geometry = feature.get("geometry") - if not geometry: - continue - try: - lines.extend(_iter_lines(shape(geometry))) - except Exception: # noqa: BLE001 - continue - return lines - - -def _uphill_streams( - outlet: Point, - streams: list[LineString], - contours: ContourField, - tolerance_m: float = 30.0, -) -> list[LineString]: - """측점에 연결된 세류 가지 중 위쪽(표고가 높아지는 방향)으로 뻗은 것만 모은다.""" - connected = [line for line in streams if line.distance(outlet) <= tolerance_m] - uphill: list[LineString] = [] - outlet_elevation = contours.elevation_at(outlet.x, outlet.y) - for line in connected: - far = _far_end(line, outlet) - far_elevation = contours.elevation_at(far.x, far.y) - if outlet_elevation is None or far_elevation is None or far_elevation >= outlet_elevation: - uphill.append(line) - return uphill - - -def _far_end(line: LineString, outlet: Point) -> Point: - start = Point(line.coords[0]) - end = Point(line.coords[-1]) - return end if start.distance(outlet) <= end.distance(outlet) else start - - -def _basin_polygon( - outlet: Point, - uphill: list[LineString], - summits: list[tuple[Polygon, float]], - route_line: LineString | None, -) -> Polygon | None: - """유역 경계를 만든다. - - 측점 + 상류 세류 + 정상부 폐합 등고선을 함께 감싸는 볼록 껍질을 1차 경계로 삼고, - 노선 아래쪽(성토부 방향)은 노선을 경계로 잘라낸다. 세류가 없으면 정상부까지의 - 반경 안에서 만들어지는 범위만 남는다. - """ - parts: list[Any] = [outlet.buffer(5.0)] - parts.extend(uphill) - parts.extend(polygon for polygon, _ in summits) - if len(parts) <= 1: - return None - hull = unary_union(parts).convex_hull - if hull.geom_type != "Polygon": - return None - if route_line is not None: - hull = _clip_downhill(hull, route_line, outlet) - return hull if hull is not None and hull.geom_type == "Polygon" else None - - -def _clip_downhill(hull: Polygon, route_line: LineString, outlet: Point) -> Polygon | None: - """노선을 경계로 유역을 잘라 산 쪽(상류) 조각만 남긴다.""" - try: - pieces = hull.difference(route_line.buffer(0.5)) - except Exception: # noqa: BLE001 - return hull - if pieces.is_empty: - return hull - parts = list(pieces.geoms) if pieces.geom_type == "MultiPolygon" else [pieces] - # 상류 조각 판별이 애매할 때를 대비해 면적이 가장 큰 조각을 채택한다. - best = max(parts, key=lambda part: part.area, default=None) - return best if best is not None and best.geom_type == "Polygon" else hull - - -def _flow_length(outlet: Point, uphill: list[LineString], boundary: Polygon) -> float: - """유하거리: 측점에서 유역 최상단까지의 물길 길이(m). - - 상류 세류가 있으면 그 물길 길이의 최댓값을, 없으면 유역 안 최원점까지의 직선거리를 쓴다. - """ - if uphill: - return float(max(line.length for line in uphill)) - return float(max((outlet.distance(Point(xy)) for xy in boundary.exterior.coords), default=0.0)) diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py new file mode 100644 index 00000000..bdcfe002 --- /dev/null +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py @@ -0,0 +1,449 @@ +"""배수유역 능선(분수령) 기반 산정 엔진. + +목적: 도로(관 매설 지점)로 모이는 물의 양을 알기 위한 유역 산정. 유역 경계는 반드시 +분수령(능선)을 따라야 하므로, 도엽 등고선·표고점을 격자 DEM으로 보간한 뒤 D8 흐름 +방향으로 "각 셀의 물이 어느 측점으로 흘러가는가"를 직접 추적한다. + +- 3D 라이다는 산 전체를 계측하지 않으므로 쓰지 않는다. 도엽 데이터만 사용(사용자 확정). +- 함몰 보정은 Whitebox `fill_depressions`를 쓰고(기존 Engine_Skeleton과 동일 패턴), + 실패하면 원본 DEM으로 진행한다(결과 저하 가능하나 계산은 지속). +- 유역 폴리곤 외곽선이 곧 분수령(능선)이며 프론트가 능선 스타일로 표시한다. +""" + +from __future__ import annotations + +import logging +import math +import tempfile +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any + +import numpy as np +from scipy.interpolate import griddata +from shapely.geometry import shape + +from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( + StructureCandidate, + estimate_pipe_diameter_mm, +) + +logger = logging.getLogger(__name__) + +# 격자 해상도(m)와 최대 격자 크기. 도엽 9매 범위라도 이 상한 안에서 해상도를 낮춰 계산한다. +GRID_RES_M = 10.0 +MAX_GRID_CELLS = 1_400_000 +# 유역 계산 범위: 측점 bbox + 여유폭(m). 주변 8도엽까지 확보되어 있어 넉넉히 잡는다. +BBOX_MARGIN_M = 1500.0 +# pour point 스냅 반경(m): 측점을 주변 흐름 누적 최대 셀로 옮겨 세류 격자 정합 오차를 흡수. +SNAP_RADIUS_M = 50.0 +# DEM 보간 표본 상한(속도 확보용 간축). 초과 시 균등 간격으로 추린다. +MAX_SAMPLE_POINTS = 250_000 + +# D8 이웃: (행 오프셋, 열 오프셋, 거리 계수) +_D8 = ( + (-1, -1, math.sqrt(2.0)), + (-1, 0, 1.0), + (-1, 1, math.sqrt(2.0)), + (0, -1, 1.0), + (0, 1, 1.0), + (1, -1, math.sqrt(2.0)), + (1, 0, 1.0), + (1, 1, math.sqrt(2.0)), +) + + +@dataclass +class WatershedBasin: + """능선 기반으로 산정된 배수유역 1개.""" + + index: int + chainage_m: float + outlet_x: float + outlet_y: float + # 유역 경계(사업지 좌표계 m). 외곽선이 곧 분수령(능선). + boundary_xy: list[list[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 + + +def _collect_samples( + contour_features: list[dict[str, Any]], + spot_features: list[dict[str, Any]], + elevation_keys: tuple[str, ...], +) -> np.ndarray: + """등고선 정점·표고점을 (x, y, z) 표본 배열로 모은다.""" + xs: list[float] = [] + ys: list[float] = [] + zs: list[float] = [] + + def _walk(coordinates: Any, elevation: float) -> None: + if not isinstance(coordinates, list) or not coordinates: + return + if isinstance(coordinates[0], (int, float)): + xs.append(float(coordinates[0])) + ys.append(float(coordinates[1])) + zs.append(elevation) + return + for item in coordinates: + _walk(item, elevation) + + for feature in [*contour_features, *spot_features]: + properties = feature.get("properties") or {} + elevation: float | None = None + for key in elevation_keys: + value = properties.get(key) + if value is None: + continue + try: + elevation = float(value) + break + except (TypeError, ValueError): + continue + if elevation is None: + continue + geometry = feature.get("geometry") or {} + _walk(geometry.get("coordinates"), elevation) + + if not xs: + return np.empty((0, 3)) + samples = np.column_stack([xs, ys, zs]) + if len(samples) > MAX_SAMPLE_POINTS: + step = len(samples) // MAX_SAMPLE_POINTS + 1 + samples = samples[::step] + return samples + + +def _build_dem( + samples: np.ndarray, + candidates: list[StructureCandidate], +) -> tuple[np.ndarray, np.ndarray, np.ndarray, float] | None: + """표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도).""" + if len(samples) < 10 or not candidates: + return None + min_x = min(candidate.x for candidate in candidates) - BBOX_MARGIN_M + max_x = max(candidate.x for candidate in candidates) + BBOX_MARGIN_M + min_y = min(candidate.y for candidate in candidates) - BBOX_MARGIN_M + max_y = max(candidate.y for candidate in candidates) + BBOX_MARGIN_M + # 표본 범위 밖으로는 나가지 않는다(외삽 방지). + min_x = max(min_x, float(samples[:, 0].min())) + max_x = min(max_x, float(samples[:, 0].max())) + min_y = max(min_y, float(samples[:, 1].min())) + max_y = min(max_y, float(samples[:, 1].max())) + if max_x - min_x < GRID_RES_M * 4 or max_y - min_y < GRID_RES_M * 4: + return None + + resolution = GRID_RES_M + while ((max_x - min_x) / resolution) * ((max_y - min_y) / resolution) > MAX_GRID_CELLS: + resolution *= 1.5 + x_coords = np.arange(min_x, max_x + resolution, resolution) + y_coords = np.arange(min_y, max_y + resolution, resolution) + grid_x, grid_y = np.meshgrid(x_coords, y_coords) + + points = samples[:, :2] + values = samples[:, 2] + dem = griddata(points, values, (grid_x, grid_y), method="linear") + # linear 보간 밖(볼록 껍질 바깥)은 nearest로 메워 유역 추적이 끊기지 않게 한다. + holes = ~np.isfinite(dem) + if holes.any(): + dem[holes] = griddata(points, values, (grid_x[holes], grid_y[holes]), method="nearest") + return dem.astype(np.float64), x_coords, y_coords, resolution + + +def _fill_depressions(dem: np.ndarray, resolution: float) -> np.ndarray: + """Whitebox로 함몰을 메운다. 실패하면 원본 그대로 진행한다.""" + try: + import rasterio + from rasterio.transform import from_origin + from whitebox import WhiteboxTools + except Exception: # noqa: BLE001 + return dem + rows, cols = dem.shape + try: + with tempfile.TemporaryDirectory(prefix="wbt_drain_") as tmp: + tmp_path = Path(tmp) + transform = from_origin(0.0, rows * resolution, resolution, resolution) + with rasterio.open( + tmp_path / "dem.tif", + "w", + driver="GTiff", + height=rows, + width=cols, + count=1, + dtype="float32", + nodata=-9999.0, + crs="EPSG:3857", + transform=transform, + ) as dst: + dst.write(dem.astype(np.float32)[::-1, :], 1) + wbt = WhiteboxTools() + wbt.set_verbose_mode(False) + wbt.set_working_dir(str(tmp_path)) + if wbt.fill_depressions("dem.tif", "filled.tif") != 0: + raise RuntimeError("fill_depressions 실패") + with rasterio.open(tmp_path / "filled.tif") as src: + filled = src.read(1).astype(np.float64)[::-1, :] + return np.where(np.isfinite(filled), filled, dem) + except Exception: # noqa: BLE001 + logger.warning("Whitebox 함몰 보정 실패 — 원본 DEM으로 진행") + return dem + + +def _d8_pointer(dem: np.ndarray) -> np.ndarray: + """각 셀의 최급강하 이웃 인덱스(0~7, 배수구 없으면 -1).""" + rows, cols = dem.shape + pointer = np.full((rows, cols), -1, dtype=np.int8) + best_drop = np.zeros((rows, cols), dtype=np.float64) + for direction, (dr, dc, distance) in enumerate(_D8): + shifted = np.full_like(dem, np.inf) + r_src = slice(max(0, -dr), rows - max(0, dr)) + c_src = slice(max(0, -dc), cols - max(0, dc)) + r_dst = slice(max(0, dr), rows - max(0, -dr)) + c_dst = slice(max(0, dc), cols - max(0, -dc)) + shifted[r_src, c_src] = dem[r_dst, c_dst] + drop = (dem - shifted) / distance + better = drop > best_drop + pointer[better] = direction + best_drop[better] = drop[better] + return pointer + + +def _flow_accumulation(pointer: np.ndarray) -> np.ndarray: + """D8 포인터 기반 흐름 누적(자기 자신 포함 셀 수). 위상 순서로 한 번에 계산.""" + rows, cols = pointer.shape + accumulation = np.ones((rows, cols), dtype=np.float64) + indegree = np.zeros((rows, cols), dtype=np.int32) + for direction, (dr, dc, _) in enumerate(_D8): + sources = np.argwhere(pointer == direction) + for r, c in sources: + nr, nc = r + dr, c + dc + if 0 <= nr < rows and 0 <= nc < cols: + indegree[nr, nc] += 1 + stack = [tuple(cell) for cell in np.argwhere(indegree == 0)] + while stack: + r, c = stack.pop() + direction = pointer[r, c] + if direction < 0: + continue + dr, dc, _ = _D8[direction] + nr, nc = r + dr, c + dc + if not (0 <= nr < rows and 0 <= nc < cols): + continue + accumulation[nr, nc] += accumulation[r, c] + indegree[nr, nc] -= 1 + if indegree[nr, nc] == 0: + stack.append((nr, nc)) + return accumulation + + +def _snap_outlet( + accumulation: np.ndarray, + row: int, + col: int, + radius_cells: int, +) -> tuple[int, int]: + """측점 주변 반경 안에서 흐름 누적이 가장 큰 셀로 옮긴다(물길 위로 스냅).""" + rows, cols = accumulation.shape + r0 = max(0, row - radius_cells) + r1 = min(rows, row + radius_cells + 1) + c0 = max(0, col - radius_cells) + c1 = min(cols, col + radius_cells + 1) + window = accumulation[r0:r1, c0:c1] + local = np.unravel_index(int(np.argmax(window)), window.shape) + return r0 + int(local[0]), c0 + int(local[1]) + + +def _label_basins( + pointer: np.ndarray, + outlets: dict[tuple[int, int], int], +) -> np.ndarray: + """각 셀이 흐름을 따라 처음 만나는 pour point의 라벨을 붙인다(경로 메모이제이션).""" + rows, cols = pointer.shape + labels = np.zeros((rows, cols), dtype=np.int32) # 0 = 미소속 + for (r, c), label in outlets.items(): + labels[r, c] = label + flat_pointer = pointer.ravel() + flat_labels = labels.ravel() + for start in range(flat_labels.size): + if flat_labels[start] != 0: + continue + path: list[int] = [] + current = start + label = 0 + while True: + if flat_labels[current] != 0: + label = flat_labels[current] + break + direction = flat_pointer[current] + if direction < 0: + label = -1 # 배수구 없음(격자 밖 유출) — 어떤 유역에도 속하지 않음 + break + path.append(current) + dr, dc, _ = _D8[direction] + r, c = divmod(current, cols) + nr, nc = r + dr, c + dc + if not (0 <= nr < rows and 0 <= nc < cols): + label = -1 + break + current = nr * cols + nc + for cell in path: + flat_labels[cell] = label + return labels + + +def _flow_lengths(pointer: np.ndarray, labels: np.ndarray, resolution: float) -> dict[int, float]: + """라벨별 최장 흐름 경로(셀→해당 pour point) 길이.""" + rows, cols = pointer.shape + distance = np.full((rows, cols), -1.0, dtype=np.float64) + # pour point 셀은 자기 라벨의 시작점이므로 거리 0. + longest: dict[int, float] = {} + flat_pointer = pointer.ravel() + flat_labels = labels.ravel() + flat_distance = distance.ravel() + + def _resolve(start: int) -> float: + chain: list[int] = [] + current = start + total = 0.0 + while True: + if flat_distance[current] >= 0: + total = flat_distance[current] + break + direction = flat_pointer[current] + if direction < 0: + total = 0.0 + break + r, c = divmod(current, cols) + dr, dc, factor = _D8[direction] + nr, nc = r + dr, c + dc + if not (0 <= nr < rows and 0 <= nc < cols): + total = 0.0 + break + next_cell = nr * cols + nc + # 다음 셀이 다른 라벨이면(=pour point 통과) 여기서 경로가 끝난 것으로 본다. + chain.append(current) + if flat_labels[next_cell] != flat_labels[current]: + total = 0.0 + break + current = next_cell + # 뒤에서부터 거리를 되채운다. + for cell in reversed(chain): + direction = flat_pointer[cell] + factor = _D8[direction][2] if direction >= 0 else 0.0 + total += factor * resolution + flat_distance[cell] = total + return total + + for start in range(flat_labels.size): + label = int(flat_labels[start]) + if label <= 0: + continue + length = _resolve(start) + if length > longest.get(label, 0.0): + longest[label] = length + return longest + + +def _vectorize_basin( + labels: np.ndarray, + label: int, + x_coords: np.ndarray, + y_coords: np.ndarray, + resolution: float, +) -> list[list[float]]: + """유역 셀 집합을 폴리곤 외곽 링(사업지 좌표계)으로 벡터화한다.""" + try: + from rasterio import features as rio_features + from rasterio.transform import from_origin + except Exception: # noqa: BLE001 + return [] + mask = (labels == label).astype(np.uint8) + if mask.sum() == 0: + return [] + transform = from_origin( + float(x_coords[0]) - resolution / 2.0, + float(y_coords[-1]) + resolution / 2.0, + resolution, + resolution, + ) + shapes = rio_features.shapes(mask[::-1, :], mask=mask[::-1, :] > 0, transform=transform) + polygons = [shape(geometry) for geometry, value in shapes if value == 1] + if not polygons: + return [] + merged = max(polygons, key=lambda polygon: polygon.area) + simplified = merged.simplify(resolution, preserve_topology=True) + if simplified.is_empty or simplified.geom_type != "Polygon": + simplified = merged + return [[float(x), float(y)] for x, y in simplified.exterior.coords] + + +def build_watershed_basins( + candidates: list[StructureCandidate], + contour_features: list[dict[str, Any]], + spot_features: list[dict[str, Any]], + elevation_keys: tuple[str, ...], +) -> list[WatershedBasin]: + """능선(분수령) 기반 배수유역을 산정한다. + + 반환된 boundary_xy 외곽선이 곧 분수령(능선)이다. 번호는 노선 시점에 가까운 순. + """ + if not candidates: + return [] + samples = _collect_samples(contour_features, spot_features, elevation_keys) + built = _build_dem(samples, candidates) + if built is None: + logger.warning("DEM 보간 실패 — 표본 %d개", len(samples)) + return [] + dem, x_coords, y_coords, resolution = built + dem = _fill_depressions(dem, resolution) + pointer = _d8_pointer(dem) + accumulation = _flow_accumulation(pointer) + + ordered = sorted(candidates, key=lambda item: item.chainage_m) + radius_cells = max(1, int(SNAP_RADIUS_M / resolution)) + outlets: dict[tuple[int, int], int] = {} + outlet_cells: dict[int, tuple[int, int]] = {} + for label, candidate in enumerate(ordered, start=1): + col = int(round((candidate.x - float(x_coords[0])) / resolution)) + row = int(round((candidate.y - float(y_coords[0])) / resolution)) + if not (0 <= row < dem.shape[0] and 0 <= col < dem.shape[1]): + continue + snapped = _snap_outlet(accumulation, row, col, radius_cells) + outlets[snapped] = label + outlet_cells[label] = snapped + + if not outlets: + return [] + labels = _label_basins(pointer, outlets) + lengths = _flow_lengths(pointer, labels, resolution) + + basins: list[WatershedBasin] = [] + for label, candidate in enumerate(ordered, start=1): + cell = outlet_cells.get(label) + if cell is None: + continue + mask = labels == label + cell_count = int(mask.sum()) + if cell_count < 4: + continue + boundary = _vectorize_basin(labels, label, x_coords, y_coords, resolution) + if len(boundary) < 4: + continue + outlet_z = float(dem[cell[0], cell[1]]) + basin = WatershedBasin( + index=label, + chainage_m=candidate.chainage_m, + outlet_x=candidate.x, + outlet_y=candidate.y, + boundary_xy=boundary, + area_m2=cell_count * resolution * resolution, + relief_m=max(0.0, float(dem[mask].max()) - outlet_z), + flow_length_m=lengths.get(label, 0.0), + ) + basin.pipe_diameter_mm = estimate_pipe_diameter_mm( + basin.area_m2, basin.relief_m, basin.flow_length_m + ) + basins.append(basin) + return basins diff --git a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py index 359fe649..3d137c3f 100644 --- a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py @@ -20,7 +20,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_Basin import ContourField, build_basins +from B05_wf2_Route.B05_wf2_Route_Engine_Drainage_Watershed import build_watershed_basins from B05_wf2_Route.B05_wf2_Route_Repository import ( get_latest_route, get_route_points, @@ -35,8 +35,9 @@ router = APIRouter(prefix="/api/projects", tags=["B05 Route Drainage"]) # 도엽 레이어 파일명 (B04 전처리 산출물과 동일 위치) _CONTOUR_FILE = "도엽_등고선.geojson" _STREAM_FILE = "도엽_하천중심선.geojson" -# 도엽 등고선의 표고 속성 키. gpkg 등고선(CTRLN_HG)도 함께 본다. -_ELEVATION_KEYS = ("등고수치", "CTRLN_HG", "elevation", "ELEV") +_SPOT_FILE = "도엽_표고점.geojson" +# 표고 속성 키: 도엽 등고선(등고수치)·gpkg 등고선(CTRLN_HG)·표고점(수치/표고) 통합. +_ELEVATION_KEYS = ("등고수치", "CTRLN_HG", "수치", "표고", "높이", "elevation", "ELEV") def _sheet_dir(stored_path: str) -> Path: @@ -144,11 +145,15 @@ 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, "to_lonlat": lambda x, y: to_lonlat_transformer.transform(x, y), } @@ -188,10 +193,10 @@ async def post_drainage_basins( else: candidates = propose_structure_stations(vertices, prepared["streams"]) - contours = ContourField(prepared["contours"], _ELEVATION_KEYS) - basins = build_basins( - vertices, candidates, contours, prepared["streams"], prepared["to_lonlat"] + basins = build_watershed_basins( + candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS ) + to_lonlat = prepared["to_lonlat"] return { "status": "success", "project_id": str(project_id), @@ -200,7 +205,8 @@ async def post_drainage_basins( { "index": basin.index, "chainage_m": round(basin.chainage_m, 2), - "polygon_lonlat": basin.polygon_lonlat, + # 유역 경계 외곽선 = 분수령(능선). 프론트가 파스텔 채움 + 능선 파선으로 표시한다. + "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), "flow_length_m": round(basin.flow_length_m, 1), 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 5cf7fd33..ac1f419a 100644 --- a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts +++ b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts @@ -10,6 +10,7 @@ import { createNormalizer, drawFilledRing, drawPreparedLayer, + drawRidgeRing, prepareLayer, prepareMetricPolyline, type GeoJsonCollection, @@ -116,6 +117,8 @@ export function createDrainagePanel(): DrainagePanel { let normalizer: Normalizer | null = null; let basins: DrainageBasin[] = []; let selectedBasin: number | null = null; + // 유역 경계 외곽선 = 분수령(능선). 사용자 지시로 기본 표시. + let showRidge = true; let scale = 1; let offsetX = 0; let offsetY = 0; @@ -144,6 +147,21 @@ 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})`; } @@ -181,6 +199,8 @@ export function createDrainagePanel(): DrainagePanel { ? color : color.replace(/0\.45\)$/, "0.18)"), ); + // 유역 경계 = 분수령이므로 그 외곽선을 능선 파선으로 강조한다. + if (showRidge) drawRidgeRing(context, basin.polygon_lonlat, normalizer!, view); }); } // 등고선을 얇게 깔고 세류·표고점을 그 위에, 노선을 맨 위에 둔다.