auto: 2026-07-28 19:53 (EOMSANGDON-HOME)
This commit is contained in:
@@ -521,3 +521,30 @@ export function drawFilledRing(
|
||||
context.fillStyle = "#1f2937";
|
||||
context.fillText(entry.label, centerX, centerY);
|
||||
}
|
||||
|
||||
/** 유역 경계(분수령=능선)를 능선 스타일(갈색 파선)로 강조해 그린다. */
|
||||
export function drawRidgeRing(
|
||||
context: CanvasRenderingContext2D,
|
||||
ring: ReadonlyArray<readonly [number, number]>,
|
||||
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();
|
||||
}
|
||||
|
||||
@@ -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))
|
||||
@@ -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
|
||||
@@ -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),
|
||||
|
||||
@@ -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);
|
||||
});
|
||||
}
|
||||
// 등고선을 얇게 깔고 세류·표고점을 그 위에, 노선을 맨 위에 둔다.
|
||||
|
||||
Reference in New Issue
Block a user