eomsangdon/Aislo
1
0
Fork 0
Code Issues Pull Requests 2 Actions Packages Projects Releases Wiki Activity

auto: 2026-07-29 18:40 (EOMSANGDON-HOME)

Browse Source
This commit is contained in:
eomsangdon
2026-07-29 18:40:45 +09:00
parent 368764190f
commit 067de3a25f
2 changed files with 137 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py
+125 -27
View File
@@ -19,7 +19,7 @@ from dataclasses import dataclass, field
from typing import Any
from shapely.geometry import LineString, Point, Polygon
from shapely.ops import split, substring, unary_union
from shapely.ops import nearest_points, split, substring, unary_union
from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import (
StructureCandidate,
@@ -45,6 +45,8 @@ MIN_BASIN_AREA_M2 = 100.0
STREAM_COVER_BUFFER_M = 20.0
# 도로 절단용 노선 양끝 연장 길이(m). 노선 끝 너머로 새는 하류측 조각까지 잘라낸다.
ROAD_EXTEND_M = 500.0
# 폐합 등고선 기하 탐색: 분할선 경로와 등고선의 근접 허용 거리(m).
CLOSING_SEARCH_M = 30.0
# 산측 판정: 도로 접선의 좌우 수직 오프셋(m)과 표고 조회 반경(m).
UPHILL_PROBE_OFFSET_M = 40.0
UPHILL_PROBE_RADIUS_M = 60.0
@@ -124,11 +126,16 @@ def _road_segment_coords(
def _contour_arc(
line: Any, p_from: Point, p_to: Point, road_line: LineString
line: Any,
p_from: Point,
p_to: Point,
road_line: LineString,
network_union: Any = None,
) -> list[tuple[float, float]]:
"""등고선에서 두 분할선 접점 사이 아크(상측 경계)를 뽑는다.
폐합 등고선은 두 방향 아크가 생기므로 도로와 교차하지 않는(=산측) 쪽을 고른다.
폐합 등고선은 두 방향 아크가 생기므로, 세류 상류망을 가로지르지 않고(계곡을
자르지 않고) 도로와도 교차하지 않는(=산측) 쪽을 고른다.
"""
t1, t2 = sorted((line.project(p_from), line.project(p_to)))
arcs = []
@@ -145,36 +152,93 @@ def _contour_arc(
return []
scored = []
for arc in arcs:
crosses = arc.crosses(road_line)
crosses_stream = bool(network_union is not None and arc.crosses(network_union))
crosses_road = 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]
scored.append((crosses_stream, crosses_road, -road_line.distance(midpoint), arc))
scored.sort(key=lambda item: (item[0], item[1], item[2]))
arc = scored[0][3]
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)."""
def _junction(
left: list[DividerStep],
right: list[DividerStep],
min_z: float | None = None,
) -> tuple[int, int, int] | None:
"""두 분할선이 같은 등고선 지오메트리를 밟은 폐합 지점(좌 idx, 우 idx, geom idx).
min_z(세류 상류망 최고 표고)가 있으면 그 **이상인 가장 낮은** 공통 등고선을 고른다 —
"세류로 영역을 지정한 뒤 가까운 등고선으로 바로 올려치면 안 된다"(2026-07-29 사용자
지시). 계곡 발원부를 넘긴 첫 등고선이 유역 상측 경계가 된다. 없으면 최고 공통 등고선.
"""
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
matches: list[tuple[float, int, int, int]] = []
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)
matches.append((step.z, left_position, position, step.geom_index))
if not matches:
return None
if min_z is not None:
above = [match for match in matches if match[0] >= min_z]
if above:
best = min(above)
return best[1], best[2], best[3]
best = max(matches)
return best[1], best[2], best[3]
def _closing_contour(
left: list[DividerStep],
right: list[DividerStep],
contour_index: ContourIndex,
min_z: float,
) -> tuple[int, int, int, Point, Point] | None:
"""두 분할선 경로에 모두 근접한 등고선 중 min_z 이상 최저를 찾는다.
분할선이 같은 스텝에서 같은 지오메트리를 밟지 못해도(도엽 분할 등) 계곡 발원부
위를 지나는 폐합 등고선을 기하적으로 찾아낸다.
반환: (좌 절단 idx, 우 절단 idx, 등고선 geom idx, 좌 접점, 우 접점).
"""
if len(left) < 2 or len(right) < 2:
return None
left_line = LineString([step.point for step in left])
right_line = LineString([step.point for step in right])
shared = set(contour_index.query(left_line.buffer(CLOSING_SEARCH_M))) & set(
contour_index.query(right_line.buffer(CLOSING_SEARCH_M))
)
best: tuple[float, int] | None = None
for index in shared:
z = contour_index.zs[index]
if z < min_z:
continue
geom = contour_index.geoms[index]
if (
geom.distance(left_line) > CLOSING_SEARCH_M
or geom.distance(right_line) > CLOSING_SEARCH_M
):
continue
if best is None or z < best[0]:
best = (z, index)
if best is None:
return None
return best[1], best[2], best[3]
geom = contour_index.geoms[best[1]]
left_touch = nearest_points(geom, left_line)[0]
right_touch = nearest_points(geom, right_line)[0]
left_position = min(range(len(left)), key=lambda i: left[i].point.distance(left_touch))
right_position = min(range(len(right)), key=lambda i: right[i].point.distance(right_touch))
return left_position, right_position, best[1], left_touch, right_touch
def _assemble_polygon(
@@ -185,22 +249,41 @@ def _assemble_polygon(
right: list[DividerStep],
contour_index: ContourIndex,
road_line: LineString,
network_union: Any = None,
valley_top_z: float | None = None,
) -> Polygon | None:
"""도로 구간 + 우측 분할선 + 상측 등고선 아크 + 좌측 분할선으로 폴리곤을 폐합한다."""
"""도로 구간 + 우측 분할선 + 상측 등고선 아크 + 좌측 분할선으로 폴리곤을 폐합한다.
세류 유역(valley_top_z 지정)은 발원부 위를 지나는 폐합 등고선을 기하 탐색으로
먼저 찾고, 실패 시 같은 스텝 매칭(_junction)으로 폐합한다.
"""
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, right_used, arc = left, right, []
closure = (
_closing_contour(left, right, contour_index, valley_top_z)
if valley_top_z is not None
else None
)
if closure is not None:
left_position, right_position, geom_index, left_touch, right_touch = closure
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,
contour_index.geoms[geom_index], right_touch, left_touch, road_line, network_union
)
else:
left_used, right_used, arc = left, right, []
junction = _junction(left, right, min_z=valley_top_z)
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,
network_union,
)
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:]))
@@ -230,13 +313,15 @@ def _extended_road_line(vertices: list[Any]) -> LineString:
def _clip_to_uphill(
polygon: Polygon,
road_line: LineString,
extended_road: LineString,
uphill_sign: int,
keep_geom: Any = None,
) -> Polygon | None:
"""폴리곤을 도로선으로 절단해 산측 조각만 남긴다 — 도로가 유역의 한쪽 경계.
도로 하류측(성토부 아래)은 유역에 포함하지 않는다(2026-07-29 사용자 지시).
keep_geom(세류 상류망)이 걸린 조각은 부호와 무관하게 유지한다 — 상류망은 정의상
산측 배수인데, 노선 끝을 감아 도는 계곡은 연장 도로선 기준 부호가 뒤집힐 수 있다.
"""
try:
pieces = split(polygon, extended_road)
@@ -246,7 +331,11 @@ def _clip_to_uphill(
piece
for piece in getattr(pieces, "geoms", [pieces])
if piece.geom_type == "Polygon"
and side_sign(road_line, piece.representative_point()) == uphill_sign
and (
side_sign(extended_road, piece.representative_point()) == uphill_sign
# 교차점(도로 위 시작점) 접촉만으로는 부족 — 상류망이 실제로 지나가야 한다.
or (keep_geom is not None and piece.intersection(keep_geom).length > 5.0)
)
]
if not kept:
return None
@@ -328,6 +417,13 @@ def build_watershed_basins(
network, flow_length = trace_upstream_network(
outlet, stream_features, road_line, signs[position]
)
# 계곡 발원부(상류망 최고 표고) — 유역 상측 폐합 등고선은 이보다 높아야 한다.
network_union = unary_union(network) if network else None
valley_top_z = (
contour_index.max_elevation_within(network_union.buffer(10.0))
if network_union is not None
else None
)
polygon = _assemble_polygon(
vertices,
divides[position],
@@ -336,18 +432,20 @@ def build_watershed_basins(
dividers[position + 1],
contour_index,
road_line,
network_union,
valley_top_z,
)
if polygon is None:
continue
if network:
if network_union is not None:
# 상류망이 폴리곤 밖으로 뻗은 경우까지 유역이 덮도록 보정한다.
covered = polygon.union(unary_union(network).buffer(STREAM_COVER_BUFFER_M)).buffer(0)
covered = polygon.union(network_union.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
# 마지막에 도로선으로 절단해 하류측을 제거한다 — 도로가 유역의 한쪽 경계.
polygon = _clip_to_uphill(polygon, road_line, extended_road, signs[position])
polygon = _clip_to_uphill(polygon, extended_road, signs[position], network_union)
if polygon is None or polygon.area < MIN_BASIN_AREA_M2:
continue
B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Trace.py
+12 -6
View File
@@ -216,13 +216,20 @@ def trace_upstream_network(
used: set[int] = set()
total = 0.0
# ① 뿌리: 교차점을 지나는 세류의 산측 조각.
# ① 뿌리: 교차점을 지나는 세류의 산측 조각. 도엽 세류는 교차점 부근에서 별도
# 피처로 조각나 있는 경우가 많아, 근접(STREAM_JOIN_TOL_M) 조각도 뿌리로 받는다.
for index, line in enumerate(lines):
if line.distance(crossing) > 1.0:
distance = line.distance(crossing)
if distance > STREAM_JOIN_TOL_M:
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:
used.add(index)
if distance <= 1.0:
t = line.project(crossing)
pieces = [substring(line, 0.0, t), substring(line, t, line.length)]
else:
pieces = [line]
for piece in pieces:
if piece.geom_type != "LineString" or piece.length < 1.0:
continue
oriented = _oriented_from(piece, crossing)
probe = oriented.interpolate(min(10.0, oriented.length))
@@ -233,7 +240,6 @@ def trace_upstream_network(
continue
network.append((clipped, 0.0))
total += clipped.length
used.add(index)
if not network:
return [], 0.0