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