auto: 2026-07-29 18:07 (EOMSANGDON-HOME)
This commit is contained in:
@@ -228,12 +228,18 @@ def _fill_spacing(
|
||||
start_m: float,
|
||||
end_m: float,
|
||||
) -> list[StructureCandidate]:
|
||||
"""[start, end] 구간이 300m를 넘으면 절토부 지점에 보충 측점을 만든다."""
|
||||
"""[start, end] 구간이 300m를 넘으면 보충 측점을 만든다.
|
||||
|
||||
종단도상 상대적으로 물이 모일 것으로 예상되는 지점(절토부 내 종단 저점)을
|
||||
우선 배치한다(2026-07-29 사용자 지시). 저점이 없으면 목표 인근 절토부로 대체한다.
|
||||
"""
|
||||
added: list[StructureCandidate] = []
|
||||
cursor = start_m
|
||||
while end_m - cursor > MAX_STRUCTURE_SPACING_M:
|
||||
target = cursor + MAX_STRUCTURE_SPACING_M
|
||||
placed = _nearest_uphill(vertices, target, end_m)
|
||||
placed = _gather_low_point(vertices, cursor, target, end_m)
|
||||
if placed is None:
|
||||
placed = _nearest_uphill(vertices, target, end_m)
|
||||
if placed is None:
|
||||
break
|
||||
x, y, _ = _interpolate_vertex(vertices, placed)
|
||||
@@ -242,6 +248,39 @@ def _fill_spacing(
|
||||
return added
|
||||
|
||||
|
||||
def _gather_low_point(
|
||||
vertices: list[RouteVertex],
|
||||
cursor_m: float,
|
||||
target_m: float,
|
||||
limit_m: float,
|
||||
step_m: float = 10.0,
|
||||
) -> float | None:
|
||||
"""탐색창 [cursor+150, target] 안 절토부의 종단 국소 저점(사그) 중 가장 낮은 지점.
|
||||
|
||||
창 하한을 간격의 절반으로 두어 보충 측점이 과밀하게 몰리지 않게 하고,
|
||||
국소 저점만 인정해 일정 오르막에서는 None(300m 규칙 폴백)을 돌려준다.
|
||||
"""
|
||||
window_start = cursor_m + MAX_STRUCTURE_SPACING_M / 2.0
|
||||
probes: list[float] = []
|
||||
probe = window_start - step_m
|
||||
while probe <= target_m + step_m:
|
||||
probes.append(probe)
|
||||
probe += step_m
|
||||
heights = [_interpolate_vertex(vertices, position)[2] for position in probes]
|
||||
best: tuple[float, float] | None = None # (계획고 z, 누가거리)
|
||||
for i in range(1, len(probes) - 1):
|
||||
position = probes[i]
|
||||
if position >= limit_m or position > target_m or position < window_start:
|
||||
continue
|
||||
# 국소 저점(양쪽이 같거나 높음) = 물이 모여 더 못 흐르는 지점. 앞쪽이 오르막인
|
||||
# 조건을 내포하므로 별도의 절토부(is_uphill_at) 판정은 두지 않는다.
|
||||
if heights[i] > heights[i - 1] or heights[i] > heights[i + 1]:
|
||||
continue
|
||||
if best is None or heights[i] < best[0]:
|
||||
best = (heights[i], position)
|
||||
return best[1] if best else None
|
||||
|
||||
|
||||
def _nearest_uphill(
|
||||
vertices: list[RouteVertex],
|
||||
target_m: float,
|
||||
|
||||
@@ -1,56 +1,50 @@
|
||||
"""배수유역 도로(노선) 기준 산정 엔진 — 측구 흐름 모델.
|
||||
"""배수유역 산정 엔진 — 세류 기반 등고선 기하 직접 분석 (2026-07-29 합의).
|
||||
|
||||
목적: 도로로 오는 물의 양 산정. 참고 도면(2026-07-28 사용자 제공)과 같이 **도로 산측
|
||||
사면 전체를 관(구조물 측점) 개수만큼 빈틈없이 분할**한다:
|
||||
DEM 보간·D8 전역 흐름분석을 쓰지 않는다. 계산 순서(사용자 정의 7단계):
|
||||
① 도로(노선)가 유역의 하측 경계 ② 도로를 가로지르는 세류 교차점에서 출발
|
||||
③ 세류 상류망을 추적하고 연관 등고선만 분석해 메인 유역 선정(하류 무의미)
|
||||
④ 세류 교차점 = 관매설 지점 ⑤ 300m 초과 구간은 종단 저점에 보충(제안 엔진 담당)
|
||||
⑥ 관 사이 물갈림 고개에서 오르는 분할선(능선 근사)으로 유역을 세분화하고
|
||||
번호·면적·표고차·유하장을 산출 ⑦ 관 추가·경로 변경 시 재호출로 재분석.
|
||||
|
||||
- 도엽 등고선·표고점을 격자 DEM으로 보간 (3D 라이다는 산 전체를 계측하지 않아 미사용).
|
||||
- 노선을 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로 도로 셀마다
|
||||
담당 관을 배정 — "도로에 닿은 물은 측구를 타고 내리막의 첫 관으로 들어간다".
|
||||
- 사면 각 셀은 D8 흐름으로 내려가 처음 닿는 도로 셀의 관을 물려받는다. 도로를 만나지
|
||||
못하는 셀(도로 하측 사면, 능선 너머)은 자동 제외.
|
||||
- 함몰 보정은 Whitebox `fill_depressions`(Engine_Skeleton과 동일 패턴), 실패 시 원본 진행.
|
||||
- 유역 경계의 산측이 분수령(능선)·지능선이고 하측이 도로선이다. 프론트가 능선 파선 표시.
|
||||
유역 폴리곤 = 도로 구간(하측) + 좌우 분할선 + 최상위 공통 등고선 아크(상측)로 폐합.
|
||||
등고선은 STRtree에서 필요한 것만 꺼내므로 분석량이 유역 크기에 비례한다(도엽 매수 무관).
|
||||
"""
|
||||
|
||||
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 shapely.geometry import LineString, Point, Polygon
|
||||
from shapely.ops import substring, unary_union
|
||||
|
||||
from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import (
|
||||
StructureCandidate,
|
||||
_interpolate_vertex,
|
||||
estimate_pipe_diameter_mm,
|
||||
)
|
||||
from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Trace import (
|
||||
LOCAL_MAX_STEPS,
|
||||
LOCAL_SEARCH_RADIUS_M,
|
||||
STREAM_JOIN_TOL_M,
|
||||
ContourIndex,
|
||||
DividerStep,
|
||||
_explode_lines,
|
||||
trace_divider,
|
||||
trace_upstream_network,
|
||||
)
|
||||
|
||||
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
|
||||
# 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)),
|
||||
)
|
||||
# 유효 유역 최소 면적(m²)과 상류망 커버 보정 버퍼(m).
|
||||
MIN_BASIN_AREA_M2 = 100.0
|
||||
STREAM_COVER_BUFFER_M = 20.0
|
||||
# 산측 판정: 도로 접선의 좌우 수직 오프셋(m)과 표고 조회 반경(m).
|
||||
UPHILL_PROBE_OFFSET_M = 40.0
|
||||
UPHILL_PROBE_RADIUS_M = 60.0
|
||||
|
||||
|
||||
@dataclass
|
||||
@@ -61,7 +55,7 @@ class WatershedBasin:
|
||||
chainage_m: float
|
||||
outlet_x: float
|
||||
outlet_y: float
|
||||
# 유역 경계(사업지 좌표계 m). 외곽선이 곧 분수령(능선).
|
||||
# 유역 경계(사업지 좌표계 m). 외곽선의 산측이 곧 분수령(능선), 하측이 도로선.
|
||||
boundary_xy: list[list[float]] = field(default_factory=list)
|
||||
area_m2: float = 0.0
|
||||
relief_m: float = 0.0
|
||||
@@ -69,164 +63,13 @@ class WatershedBasin:
|
||||
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 _divide_chainages(vertices: list[Any], ordered: list[StructureCandidate]) -> list[float]:
|
||||
"""유역 분할 누가거리 목록(양끝 포함, 관 개수+1개).
|
||||
|
||||
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,
|
||||
anchors_x: list[float],
|
||||
anchors_y: list[float],
|
||||
) -> tuple[np.ndarray, np.ndarray, np.ndarray, float] | None:
|
||||
"""표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도).
|
||||
|
||||
범위는 앵커(노선 전체 정점) bbox + 여유폭 — 유역이 노선 전 연장을 덮어야 하므로
|
||||
측점 bbox가 아니라 노선 bbox를 쓴다.
|
||||
인접한 두 관 사이 **종단 계획선의 최고점(물갈림 고개)**이 분할점이다 —
|
||||
"도로에 닿은 물은 측구를 타고 내리막의 첫 관으로 들어간다".
|
||||
"""
|
||||
if len(samples) < 10 or not anchors_x:
|
||||
return None
|
||||
min_x = min(anchors_x) - BBOX_MARGIN_M
|
||||
max_x = max(anchors_x) + BBOX_MARGIN_M
|
||||
min_y = min(anchors_y) - BBOX_MARGIN_M
|
||||
max_y = max(anchors_y) + 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 _ditch_intervals(
|
||||
vertices: list[Any],
|
||||
ordered: list[StructureCandidate],
|
||||
) -> list[tuple[float, float, int]]:
|
||||
"""측구(도로변 배수로) 흐름 모델의 담당 구간을 만든다.
|
||||
|
||||
도로에 닿은 물은 종단 내리막 방향으로 흘러 첫 번째 관으로 들어간다. 따라서 인접한
|
||||
두 관 사이 **종단 계획선의 최고점(물갈림 고개)**이 유역 분할점이 된다.
|
||||
반환: (구간 시작 chainage, 구간 끝 chainage, 관 라벨) 목록.
|
||||
"""
|
||||
total_length = vertices[-1].chainage_m if vertices else 0.0
|
||||
divides: list[float] = []
|
||||
divides = [vertices[0].chainage_m]
|
||||
for left, right in zip(ordered, ordered[1:]):
|
||||
window = [
|
||||
vertex for vertex in vertices if left.chainage_m < vertex.chainage_m < right.chainage_m
|
||||
@@ -235,179 +78,137 @@ def _ditch_intervals(
|
||||
divides.append(max(window, key=lambda vertex: vertex.z).chainage_m)
|
||||
else:
|
||||
divides.append((left.chainage_m + right.chainage_m) / 2.0)
|
||||
intervals: list[tuple[float, float, int]] = []
|
||||
start = 0.0
|
||||
for label, divide in enumerate(divides, start=1):
|
||||
intervals.append((start, divide, label))
|
||||
start = divide
|
||||
intervals.append((start, total_length + 1.0, len(ordered)))
|
||||
return intervals
|
||||
divides.append(vertices[-1].chainage_m)
|
||||
return divides
|
||||
|
||||
|
||||
def _rasterize_road(
|
||||
vertices: list[Any],
|
||||
intervals: list[tuple[float, float, int]],
|
||||
x_coords: np.ndarray,
|
||||
y_coords: np.ndarray,
|
||||
resolution: float,
|
||||
shape_rc: tuple[int, int],
|
||||
) -> dict[tuple[int, int], int]:
|
||||
"""노선 폴리라인을 격자에 새겨 도로 셀마다 담당 관 라벨을 붙인다.
|
||||
|
||||
대각 누수(도로가 셀 사이 대각으로 지나가 물이 새는 것)를 막으려고 도로 셀 주변
|
||||
3×3을 함께 같은 라벨로 칠한다.
|
||||
"""
|
||||
rows, cols = shape_rc
|
||||
|
||||
def _label_of(chainage: float) -> int:
|
||||
for start, end, label in intervals:
|
||||
if start <= chainage < end:
|
||||
return label
|
||||
return intervals[-1][2] if intervals else 0
|
||||
|
||||
road: dict[tuple[int, int], int] = {}
|
||||
step = resolution / 2.0
|
||||
for previous, current in zip(vertices, vertices[1:]):
|
||||
span = math.dist((previous.x, previous.y), (current.x, current.y))
|
||||
count = max(1, int(span / step))
|
||||
for i in range(count + 1):
|
||||
ratio = i / count
|
||||
x = previous.x + (current.x - previous.x) * ratio
|
||||
y = previous.y + (current.y - previous.y) * ratio
|
||||
chainage = previous.chainage_m + (current.chainage_m - previous.chainage_m) * ratio
|
||||
col = int(round((x - float(x_coords[0])) / resolution))
|
||||
row = int(round((y - float(y_coords[0])) / resolution))
|
||||
label = _label_of(chainage)
|
||||
for dr in (-1, 0, 1):
|
||||
for dc in (-1, 0, 1):
|
||||
r, c = row + dr, col + dc
|
||||
if 0 <= r < rows and 0 <= c < cols and (r, c) not in road:
|
||||
road[(r, c)] = label
|
||||
return road
|
||||
|
||||
|
||||
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,
|
||||
def _uphill_sign_at(vertices: list[Any], chainage_m: float, contour_index: ContourIndex) -> int:
|
||||
"""해당 측점의 산측이 도로 진행방향 기준 좌(+1)인지 우(-1)인지. 불명이면 0."""
|
||||
x, y, _ = _interpolate_vertex(vertices, chainage_m)
|
||||
back = _interpolate_vertex(vertices, max(0.0, chainage_m - 10.0))
|
||||
forward = _interpolate_vertex(vertices, chainage_m + 10.0)
|
||||
dx, dy = forward[0] - back[0], forward[1] - back[1]
|
||||
norm = math.hypot(dx, dy)
|
||||
if norm < 1e-6:
|
||||
return 0
|
||||
dx, dy = dx / norm, dy / norm
|
||||
# 좌측 법선 (-dy, dx) 방향 오프셋이 side_sign +1에 대응한다.
|
||||
left_z = contour_index.nearest_elevation(
|
||||
Point(x - dy * UPHILL_PROBE_OFFSET_M, y + dx * UPHILL_PROBE_OFFSET_M),
|
||||
UPHILL_PROBE_RADIUS_M,
|
||||
)
|
||||
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:
|
||||
right_z = contour_index.nearest_elevation(
|
||||
Point(x + dy * UPHILL_PROBE_OFFSET_M, y - dx * UPHILL_PROBE_OFFSET_M),
|
||||
UPHILL_PROBE_RADIUS_M,
|
||||
)
|
||||
if left_z is None or right_z is None or left_z == right_z:
|
||||
return 0
|
||||
return 1 if left_z > right_z else -1
|
||||
|
||||
|
||||
def _road_segment_coords(
|
||||
vertices: list[Any], start_m: float, end_m: float
|
||||
) -> list[tuple[float, float]]:
|
||||
"""분할점 사이 도로 구간의 평면 좌표열(유역 폴리곤의 하측 경계)."""
|
||||
sx, sy, _ = _interpolate_vertex(vertices, start_m)
|
||||
ex, ey, _ = _interpolate_vertex(vertices, end_m)
|
||||
coords = [(sx, sy)]
|
||||
coords.extend(
|
||||
(vertex.x, vertex.y) for vertex in vertices if start_m < vertex.chainage_m < end_m
|
||||
)
|
||||
coords.append((ex, ey))
|
||||
return coords
|
||||
|
||||
|
||||
def _contour_arc(
|
||||
line: Any, p_from: Point, p_to: Point, road_line: LineString
|
||||
) -> list[tuple[float, float]]:
|
||||
"""등고선에서 두 분할선 접점 사이 아크(상측 경계)를 뽑는다.
|
||||
|
||||
폐합 등고선은 두 방향 아크가 생기므로 도로와 교차하지 않는(=산측) 쪽을 고른다.
|
||||
"""
|
||||
t1, t2 = sorted((line.project(p_from), line.project(p_to)))
|
||||
arcs = []
|
||||
inner = substring(line, t1, t2)
|
||||
if inner.geom_type == "LineString" and len(inner.coords) >= 2:
|
||||
arcs.append(inner)
|
||||
if getattr(line, "is_closed", False):
|
||||
head = substring(line, t2, line.length)
|
||||
tail = substring(line, 0.0, t1)
|
||||
coords = list(head.coords) + list(tail.coords)[1:]
|
||||
if len(coords) >= 2:
|
||||
arcs.append(LineString(coords))
|
||||
if not arcs:
|
||||
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]
|
||||
scored = []
|
||||
for arc in arcs:
|
||||
crosses = arc.crosses(road_line)
|
||||
midpoint = arc.interpolate(0.5, normalized=True)
|
||||
scored.append((crosses, -road_line.distance(midpoint), arc))
|
||||
scored.sort(key=lambda item: (item[0], item[1]))
|
||||
arc = scored[0][2]
|
||||
coords = list(arc.coords)
|
||||
if Point(coords[0]).distance(p_from) > Point(coords[-1]).distance(p_from):
|
||||
coords.reverse()
|
||||
return [(float(x), float(y)) for x, y in coords]
|
||||
|
||||
|
||||
def _junction(left: list[DividerStep], right: list[DividerStep]) -> tuple[int, int, int] | None:
|
||||
"""두 분할선이 같은 등고선 지오메트리를 밟은 최고 표고 지점(좌 idx, 우 idx, geom idx)."""
|
||||
left_keys = {
|
||||
(step.z, step.geom_index): position
|
||||
for position, step in enumerate(left)
|
||||
if step.geom_index >= 0
|
||||
}
|
||||
best: tuple[float, int, int, int] | None = None
|
||||
for position, step in enumerate(right):
|
||||
if step.geom_index < 0:
|
||||
continue
|
||||
left_position = left_keys.get((step.z, step.geom_index))
|
||||
if left_position is None:
|
||||
continue
|
||||
if best is None or step.z > best[0]:
|
||||
best = (step.z, left_position, position, step.geom_index)
|
||||
if best is None:
|
||||
return None
|
||||
return best[1], best[2], best[3]
|
||||
|
||||
|
||||
def _assemble_polygon(
|
||||
vertices: list[Any],
|
||||
start_m: float,
|
||||
end_m: float,
|
||||
left: list[DividerStep],
|
||||
right: list[DividerStep],
|
||||
contour_index: ContourIndex,
|
||||
road_line: LineString,
|
||||
) -> Polygon | None:
|
||||
"""도로 구간 + 우측 분할선 + 상측 등고선 아크 + 좌측 분할선으로 폴리곤을 폐합한다."""
|
||||
ring = _road_segment_coords(vertices, start_m, end_m)
|
||||
junction = _junction(left, right)
|
||||
if junction is not None:
|
||||
left_position, right_position, geom_index = junction
|
||||
left_used = left[: left_position + 1]
|
||||
right_used = right[: right_position + 1]
|
||||
arc = _contour_arc(
|
||||
contour_index.geoms[geom_index],
|
||||
right_used[-1].point,
|
||||
left_used[-1].point,
|
||||
road_line,
|
||||
)
|
||||
else:
|
||||
left_used, right_used, arc = left, right, []
|
||||
ring.extend((step.point.x, step.point.y) for step in right_used[1:])
|
||||
ring.extend(arc)
|
||||
ring.extend((step.point.x, step.point.y) for step in reversed(left_used[1:]))
|
||||
if len(ring) < 4:
|
||||
return None
|
||||
polygon = Polygon(ring).buffer(0)
|
||||
if polygon.geom_type == "MultiPolygon":
|
||||
polygon = max(polygon.geoms, key=lambda part: part.area)
|
||||
if polygon.is_empty or polygon.geom_type != "Polygon":
|
||||
return None
|
||||
return polygon
|
||||
|
||||
|
||||
def build_watershed_basins(
|
||||
@@ -416,64 +217,116 @@ def build_watershed_basins(
|
||||
contour_features: list[dict[str, Any]],
|
||||
spot_features: list[dict[str, Any]],
|
||||
elevation_keys: tuple[str, ...],
|
||||
stream_features: list[dict[str, Any]] | None = None,
|
||||
) -> list[WatershedBasin]:
|
||||
"""도로(노선) 기준 배수유역을 산정한다 — 측구 흐름 모델.
|
||||
"""관 지점 배치를 기준으로 메인 배수유역을 세분화해 산정한다.
|
||||
|
||||
"도로에 닿은 물은 측구를 타고 종단 내리막 방향으로 흘러 첫 번째 관으로 들어간다"를
|
||||
전제로, 도로 산측 사면 전체를 관 개수만큼 빈틈없이 분할한다(참고 도면과 동일 개념):
|
||||
① 노선 전체를 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로
|
||||
각 도로 셀에 담당 관 라벨을 붙인다.
|
||||
② 사면 각 셀은 D8 흐름을 따라 내려가 처음 닿는 도로 셀의 관 라벨을 물려받는다.
|
||||
도로를 만나지 못하고 격자 밖으로 빠지는 셀(도로 하측 성토부 등)은 제외된다.
|
||||
반환된 boundary_xy 외곽선의 산측이 곧 분수령(능선)이며, 하측은 도로선을 따른다.
|
||||
세류 교차 관("stream")은 상류망을 추적해 넓은 한계로, 세류 없는 관은 도로 상측
|
||||
첫 능선까지 소범위 한계로 분할선을 올린다(작은 유역, 영역 선정 주의 — 사용자 지시).
|
||||
번호는 노선 시점에 가까운 순.
|
||||
"""
|
||||
if not candidates or len(vertices) < 2:
|
||||
return []
|
||||
samples = _collect_samples(contour_features, spot_features, elevation_keys)
|
||||
built = _build_dem(
|
||||
samples,
|
||||
[vertex.x for vertex in vertices],
|
||||
[vertex.y for vertex in vertices],
|
||||
)
|
||||
if built is None:
|
||||
logger.warning("DEM 보간 실패 — 표본 %d개", len(samples))
|
||||
contour_index = ContourIndex(contour_features, elevation_keys)
|
||||
if contour_index.tree is None:
|
||||
logger.warning("표고 속성이 있는 등고선이 없어 유역을 산정하지 못했습니다.")
|
||||
return []
|
||||
dem, x_coords, y_coords, resolution = built
|
||||
dem = _fill_depressions(dem, resolution)
|
||||
pointer = _d8_pointer(dem)
|
||||
spot_index = ContourIndex(spot_features, elevation_keys)
|
||||
road_line = LineString([(vertex.x, vertex.y) for vertex in vertices])
|
||||
|
||||
ordered = sorted(candidates, key=lambda item: item.chainage_m)
|
||||
intervals = _ditch_intervals(vertices, ordered)
|
||||
outlets = _rasterize_road(vertices, intervals, x_coords, y_coords, resolution, dem.shape)
|
||||
if not outlets:
|
||||
return []
|
||||
labels = _label_basins(pointer, outlets)
|
||||
lengths = _flow_lengths(pointer, labels, resolution)
|
||||
divides = _divide_chainages(vertices, ordered)
|
||||
signs = [_uphill_sign_at(vertices, item.chainage_m, contour_index) for item in ordered]
|
||||
majority = 1 if sum(signs) >= 0 else -1
|
||||
signs = [sign or majority for sign in signs]
|
||||
# 사용자 확정("confirmed") 측점도 세류에 닿아 있으면 세류 유역으로 취급한다.
|
||||
stream_lines = _explode_lines(stream_features) if stream_features else []
|
||||
is_stream = [
|
||||
item.reason == "stream"
|
||||
or any(line.distance(Point(item.x, item.y)) <= STREAM_JOIN_TOL_M for line in stream_lines)
|
||||
for item in ordered
|
||||
]
|
||||
|
||||
# 분할선은 인접 유역과 공유하므로 분할점마다 1회만 추적한다.
|
||||
dividers: list[list[DividerStep]] = []
|
||||
for position, chainage in enumerate(divides):
|
||||
neighbor_streams = []
|
||||
if position > 0:
|
||||
neighbor_streams.append(is_stream[position - 1])
|
||||
if position < len(ordered):
|
||||
neighbor_streams.append(is_stream[position])
|
||||
wide = any(neighbor_streams)
|
||||
sign = signs[position - 1] if position > 0 else signs[0]
|
||||
x, y, _ = _interpolate_vertex(vertices, chainage)
|
||||
if wide:
|
||||
steps = trace_divider(Point(x, y), contour_index, road_line, sign)
|
||||
else:
|
||||
steps = trace_divider(
|
||||
Point(x, y),
|
||||
contour_index,
|
||||
road_line,
|
||||
sign,
|
||||
radius_m=LOCAL_SEARCH_RADIUS_M,
|
||||
max_steps=LOCAL_MAX_STEPS,
|
||||
)
|
||||
dividers.append(steps)
|
||||
|
||||
basins: list[WatershedBasin] = []
|
||||
for label, candidate in enumerate(ordered, start=1):
|
||||
mask = labels == label
|
||||
cell_count = int(mask.sum())
|
||||
if cell_count < 4:
|
||||
for position, candidate in enumerate(ordered):
|
||||
outlet = Point(candidate.x, candidate.y)
|
||||
network: list[Any] = []
|
||||
flow_length = 0.0
|
||||
if is_stream[position] and stream_features:
|
||||
network, flow_length = trace_upstream_network(
|
||||
outlet, stream_features, road_line, signs[position]
|
||||
)
|
||||
polygon = _assemble_polygon(
|
||||
vertices,
|
||||
divides[position],
|
||||
divides[position + 1],
|
||||
dividers[position],
|
||||
dividers[position + 1],
|
||||
contour_index,
|
||||
road_line,
|
||||
)
|
||||
if polygon is None:
|
||||
continue
|
||||
boundary = _vectorize_basin(labels, label, x_coords, y_coords, resolution)
|
||||
if len(boundary) < 4:
|
||||
if network:
|
||||
# 상류망이 폴리곤 밖으로 뻗은 경우까지 유역이 덮도록 보정한다.
|
||||
covered = polygon.union(unary_union(network).buffer(STREAM_COVER_BUFFER_M)).buffer(0)
|
||||
if covered.geom_type == "MultiPolygon":
|
||||
covered = max(covered.geoms, key=lambda part: part.area)
|
||||
if covered.geom_type == "Polygon" and not covered.is_empty:
|
||||
polygon = covered
|
||||
if polygon.area < MIN_BASIN_AREA_M2:
|
||||
continue
|
||||
# 측점(관) 위치의 표고를 기준으로 낙차를 계산한다.
|
||||
col = int(round((candidate.x - float(x_coords[0])) / resolution))
|
||||
row = int(round((candidate.y - float(y_coords[0])) / resolution))
|
||||
in_grid = 0 <= row < dem.shape[0] and 0 <= col < dem.shape[1]
|
||||
outlet_z = float(dem[row, col]) if in_grid else float(dem[mask].min())
|
||||
|
||||
outlet_z = contour_index.nearest_elevation(outlet, UPHILL_PROBE_RADIUS_M)
|
||||
if outlet_z is None:
|
||||
outlet_z = _interpolate_vertex(vertices, candidate.chainage_m)[2]
|
||||
top_z = max(
|
||||
contour_index.max_elevation_within(polygon) or outlet_z,
|
||||
spot_index.max_elevation_within(polygon) or outlet_z,
|
||||
)
|
||||
boundary_line = polygon.simplify(5.0, preserve_topology=True)
|
||||
if boundary_line.is_empty or boundary_line.geom_type != "Polygon":
|
||||
boundary_line = polygon
|
||||
boundary = [[float(x), float(y)] for x, y in boundary_line.exterior.coords]
|
||||
if flow_length <= 0.0:
|
||||
flow_length = max(
|
||||
(math.dist((candidate.x, candidate.y), point) for point in boundary),
|
||||
default=0.0,
|
||||
)
|
||||
|
||||
basin = WatershedBasin(
|
||||
index=label,
|
||||
index=len(basins) + 1,
|
||||
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),
|
||||
area_m2=float(polygon.area),
|
||||
relief_m=max(0.0, float(top_z) - float(outlet_z)),
|
||||
flow_length_m=float(flow_length),
|
||||
)
|
||||
basin.pipe_diameter_mm = estimate_pipe_diameter_mm(
|
||||
basin.area_m2, basin.relief_m, basin.flow_length_m
|
||||
|
||||
@@ -0,0 +1,271 @@
|
||||
"""배수유역 추적 유틸 — 등고선 공간 인덱스·세류 상류망·분수령(능선) 추적.
|
||||
|
||||
등고선 기하 직접 분석(2026-07-29 합의)의 하위 도구 모음. DEM 보간 없이:
|
||||
- `ContourIndex`: 등고선(선)·표고점(점)을 STRtree에 1회 적재하고 필요한 것만 꺼낸다.
|
||||
- `trace_upstream_network`: 세류 교차점에서 도로 산측 상류망만 추적한다(하류 무시).
|
||||
- `trace_divider`: 물갈림 지점에서 상향 등고선을 한 겹씩 따라 오르는 유역 분할선(능선 근사).
|
||||
|
||||
전체 등고선을 순회하는 연산을 두지 않아 분석량이 유역 크기에 비례한다(도엽 매수와 무관).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import Any
|
||||
|
||||
from shapely.geometry import LineString, Point, shape
|
||||
from shapely.ops import nearest_points, substring
|
||||
from shapely.strtree import STRtree
|
||||
|
||||
# 분할선(능선 근사) 추적: 다음 상위 등고선을 찾는 탐색 반경(m)과 최대 단계 수.
|
||||
DIVIDER_SEARCH_RADIUS_M = 120.0
|
||||
DIVIDER_MAX_STEPS = 60
|
||||
# 세류 없는 소규모 유역: 도로 상측 첫 능선까지만 오르도록 좁힌 한계(영역 선정 주의).
|
||||
LOCAL_SEARCH_RADIUS_M = 80.0
|
||||
LOCAL_MAX_STEPS = 12
|
||||
# 세류 연결 판정 이격(m)과 상류망 총연장 상한(m).
|
||||
STREAM_JOIN_TOL_M = 15.0
|
||||
MAX_UPSTREAM_TOTAL_M = 5000.0
|
||||
|
||||
|
||||
@dataclass
|
||||
class DividerStep:
|
||||
"""분할선의 한 단계 — 어느 등고선(geom_index)의 어느 지점을 밟았는지."""
|
||||
|
||||
point: Point
|
||||
z: float
|
||||
geom_index: int
|
||||
|
||||
|
||||
class ContourIndex:
|
||||
"""표고 속성이 있는 등고선·표고점 피처의 STRtree 래퍼."""
|
||||
|
||||
def __init__(self, features: list[dict[str, Any]], elevation_keys: tuple[str, ...]) -> None:
|
||||
self.geoms: list[Any] = []
|
||||
self.zs: list[float] = []
|
||||
for feature in features:
|
||||
elevation = _feature_elevation(feature, elevation_keys)
|
||||
if elevation is None:
|
||||
continue
|
||||
geometry = feature.get("geometry") or {}
|
||||
try:
|
||||
geom = shape(geometry)
|
||||
except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다
|
||||
continue
|
||||
if geom.is_empty:
|
||||
continue
|
||||
parts = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom]
|
||||
for part in parts:
|
||||
self.geoms.append(part)
|
||||
self.zs.append(elevation)
|
||||
self.tree = STRtree(self.geoms) if self.geoms else None
|
||||
|
||||
def query(self, geometry: Any) -> list[int]:
|
||||
"""geometry 근방(bbox 교차) 피처의 인덱스만 돌려준다."""
|
||||
if self.tree is None:
|
||||
return []
|
||||
return [int(i) for i in self.tree.query(geometry)]
|
||||
|
||||
def nearest_elevation(self, point: Point, radius_m: float) -> float | None:
|
||||
"""point에서 radius 안 가장 가까운 피처의 표고. 없으면 None."""
|
||||
best_z: float | None = None
|
||||
best_distance = radius_m
|
||||
for index in self.query(point.buffer(radius_m)):
|
||||
distance = self.geoms[index].distance(point)
|
||||
if distance <= best_distance:
|
||||
best_distance = distance
|
||||
best_z = self.zs[index]
|
||||
return best_z
|
||||
|
||||
def max_elevation_within(self, polygon: Any) -> float | None:
|
||||
"""polygon과 실제로 교차하는 피처들의 최고 표고."""
|
||||
best: float | None = None
|
||||
for index in self.query(polygon):
|
||||
if not polygon.intersects(self.geoms[index]):
|
||||
continue
|
||||
if best is None or self.zs[index] > best:
|
||||
best = self.zs[index]
|
||||
return best
|
||||
|
||||
|
||||
def _feature_elevation(feature: dict[str, Any], elevation_keys: tuple[str, ...]) -> float | None:
|
||||
properties = feature.get("properties") or {}
|
||||
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 side_sign(road_line: LineString, point: Point) -> int:
|
||||
"""도로선 기준 point가 어느 쪽인지(+1/-1, 선상이면 0). 국소 접선과의 외적 부호."""
|
||||
t = road_line.project(point)
|
||||
a = road_line.interpolate(max(0.0, t - 5.0))
|
||||
b = road_line.interpolate(min(road_line.length, t + 5.0))
|
||||
cross = (b.x - a.x) * (point.y - a.y) - (b.y - a.y) * (point.x - a.x)
|
||||
if cross > 0:
|
||||
return 1
|
||||
if cross < 0:
|
||||
return -1
|
||||
return 0
|
||||
|
||||
|
||||
def trace_divider(
|
||||
start: Point,
|
||||
contour_index: ContourIndex,
|
||||
road_line: LineString,
|
||||
uphill_sign: int,
|
||||
radius_m: float = DIVIDER_SEARCH_RADIUS_M,
|
||||
max_steps: int = DIVIDER_MAX_STEPS,
|
||||
) -> list[DividerStep]:
|
||||
"""물갈림 지점에서 상향 등고선을 한 겹씩 밟아 오르는 분할선을 만든다.
|
||||
|
||||
각 단계에서 반경 안의 "현재보다 높은 등고선 중 가장 낮은 것"의 최근접점으로 이동한다.
|
||||
도로 산측(uphill_sign)을 벗어나거나 도로 쪽으로 되돌아가는 이동은 막는다.
|
||||
더 높은 등고선이 반경 안에 없으면 능선(분수령)에 닿은 것으로 보고 멈춘다.
|
||||
"""
|
||||
z = contour_index.nearest_elevation(start, radius_m)
|
||||
if z is None:
|
||||
return []
|
||||
steps = [DividerStep(point=start, z=z, geom_index=-1)]
|
||||
current = start
|
||||
road_distance = road_line.distance(start)
|
||||
for _ in range(max_steps):
|
||||
best: tuple[float, float, Point, int] | None = None
|
||||
for index in contour_index.query(current.buffer(radius_m)):
|
||||
candidate_z = contour_index.zs[index]
|
||||
if candidate_z <= z + 0.01:
|
||||
continue
|
||||
if best is not None and candidate_z > best[0]:
|
||||
continue
|
||||
point = nearest_points(contour_index.geoms[index], current)[0]
|
||||
distance = current.distance(point)
|
||||
if distance > radius_m:
|
||||
continue
|
||||
# 도로 반대편·도로 방향 후퇴 금지 — 분할선은 산측으로만 오른다.
|
||||
if side_sign(road_line, point) == -uphill_sign:
|
||||
continue
|
||||
if road_line.distance(point) + 1.0 < road_distance:
|
||||
continue
|
||||
if (
|
||||
best is None
|
||||
or candidate_z < best[0]
|
||||
or (candidate_z == best[0] and distance < best[1])
|
||||
):
|
||||
best = (candidate_z, distance, point, index)
|
||||
if best is None:
|
||||
break
|
||||
z, _, current, geom_index = best
|
||||
road_distance = max(road_distance, road_line.distance(current))
|
||||
steps.append(DividerStep(point=current, z=z, geom_index=geom_index))
|
||||
return steps
|
||||
|
||||
|
||||
def _explode_lines(stream_features: list[dict[str, Any]]) -> list[LineString]:
|
||||
lines: list[LineString] = []
|
||||
for feature in stream_features:
|
||||
geometry = feature.get("geometry") or {}
|
||||
try:
|
||||
geom = shape(geometry)
|
||||
except Exception: # noqa: BLE001
|
||||
continue
|
||||
if geom.is_empty:
|
||||
continue
|
||||
parts = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom]
|
||||
lines.extend(part for part in parts if part.geom_type == "LineString")
|
||||
return lines
|
||||
|
||||
|
||||
def _oriented_from(line: LineString, origin: Point) -> LineString:
|
||||
"""origin에 가까운 끝이 시작점이 되도록 방향을 맞춘다."""
|
||||
if Point(line.coords[0]).distance(origin) <= Point(line.coords[-1]).distance(origin):
|
||||
return line
|
||||
return LineString(list(line.coords)[::-1])
|
||||
|
||||
|
||||
def _clip_uphill(line: LineString, road_line: LineString, origin: Point) -> LineString | None:
|
||||
"""도로를 다시 가로지르면 교차 지점에서 잘라 origin 쪽 조각만 남긴다."""
|
||||
if not line.crosses(road_line):
|
||||
return line
|
||||
t = line.project(nearest_points(line.intersection(road_line), origin)[0])
|
||||
piece = substring(line, 0.0, t) if line.project(origin) < t else substring(line, t, line.length)
|
||||
if piece.geom_type != "LineString" or piece.length < 1.0:
|
||||
return None
|
||||
return piece
|
||||
|
||||
|
||||
def trace_upstream_network(
|
||||
crossing: Point,
|
||||
stream_features: list[dict[str, Any]],
|
||||
road_line: LineString,
|
||||
uphill_sign: int,
|
||||
) -> tuple[list[LineString], float]:
|
||||
"""세류 교차점에서 도로 산측으로 뻗는 상류망을 추적한다.
|
||||
|
||||
① 교차한 세류를 교차점에서 잘라 산측 조각을 뿌리로 삼는다.
|
||||
② 끝점이 기존 망에 근접(STREAM_JOIN_TOL_M)한 세류를 반복 편입한다(분기 포함).
|
||||
도로를 다시 가로지르는 조각은 절단하고, 총연장 상한을 두어 폭주를 막는다.
|
||||
반환: (상류망 폴리라인 목록, 최장 유하 경로 길이 m).
|
||||
"""
|
||||
lines = _explode_lines(stream_features)
|
||||
network: list[tuple[LineString, float]] = [] # (폴리라인, 뿌리에서 시작점까지 누적거리)
|
||||
used: set[int] = set()
|
||||
total = 0.0
|
||||
|
||||
# ① 뿌리: 교차점을 지나는 세류의 산측 조각.
|
||||
for index, line in enumerate(lines):
|
||||
if line.distance(crossing) > 1.0:
|
||||
continue
|
||||
t = line.project(crossing)
|
||||
for piece in (substring(line, 0.0, t), substring(line, t, line.length)):
|
||||
if piece.geom_type != "LineString" or piece.length < STREAM_JOIN_TOL_M:
|
||||
continue
|
||||
oriented = _oriented_from(piece, crossing)
|
||||
probe = oriented.interpolate(min(10.0, oriented.length))
|
||||
if side_sign(road_line, probe) != uphill_sign:
|
||||
continue
|
||||
clipped = _clip_uphill(oriented, road_line, crossing)
|
||||
if clipped is None:
|
||||
continue
|
||||
network.append((clipped, 0.0))
|
||||
total += clipped.length
|
||||
used.add(index)
|
||||
if not network:
|
||||
return [], 0.0
|
||||
|
||||
# ② 편입 반복: 끝점이 망에 닿는 세류를 상류로 붙인다.
|
||||
grew = True
|
||||
while grew and total < MAX_UPSTREAM_TOTAL_M:
|
||||
grew = False
|
||||
for index, line in enumerate(lines):
|
||||
if index in used:
|
||||
continue
|
||||
attach: tuple[float, Point, LineString, float] | None = None
|
||||
for endpoint in (Point(line.coords[0]), Point(line.coords[-1])):
|
||||
for parent, parent_cum in network:
|
||||
distance = parent.distance(endpoint)
|
||||
if distance > STREAM_JOIN_TOL_M:
|
||||
continue
|
||||
cum = parent_cum + parent.project(endpoint)
|
||||
if attach is None or distance < attach[0]:
|
||||
attach = (distance, endpoint, parent, cum)
|
||||
if attach is None:
|
||||
continue
|
||||
used.add(index)
|
||||
_, endpoint, _, cum = attach
|
||||
if side_sign(road_line, line.interpolate(0.5, normalized=True)) == -uphill_sign:
|
||||
continue
|
||||
oriented = _oriented_from(line, endpoint)
|
||||
clipped = _clip_uphill(oriented, road_line, endpoint)
|
||||
if clipped is None:
|
||||
continue
|
||||
network.append((clipped, cum))
|
||||
total += clipped.length
|
||||
grew = True
|
||||
|
||||
flow_length = max((cum + line.length for line, cum in network), default=0.0)
|
||||
return [line for line, _ in network], flow_length
|
||||
@@ -194,7 +194,12 @@ async def post_drainage_basins(
|
||||
candidates = propose_structure_stations(vertices, prepared["streams"])
|
||||
|
||||
basins = build_watershed_basins(
|
||||
vertices, candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS
|
||||
vertices,
|
||||
candidates,
|
||||
prepared["contours"],
|
||||
prepared["spots"],
|
||||
_ELEVATION_KEYS,
|
||||
stream_features=prepared["streams"],
|
||||
)
|
||||
to_lonlat = prepared["to_lonlat"]
|
||||
return {
|
||||
|
||||
Reference in New Issue
Block a user