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

auto: 2026-07-28 20:03 (EOMSANGDON-HOME)

Browse Source
This commit is contained in:
eomsangdon
2026-07-28 20:03:32 +09:00
parent cb35e8c1cf
commit 121d20ef80
2 changed files with 114 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files
B05_wf2_Route/B05_wf2_Route_Engine_Drainage_Watershed.py
+113 -80
View File
@@ -1,13 +1,15 @@
"""배수유역 능선(분수령) 기반 산정 엔진.
"""배수유역 도로(노선) 기준 산정 엔진 — 측구 흐름 모델.
목적: 도로(관 매설 지점)로 모이는 물의 양을 알기 위한 유역 산정. 유역 경계는 반드시
분수령(능선)을 따라야 하므로, 도엽 등고선·표고점을 격자 DEM으로 보간한 뒤 D8 흐름
방향으로 "각 셀의 물이 어느 측점으로 흘러가는가"를 직접 추적한다.
목적: 도로로 오는 물의 양 산정. 참고 도면(2026-07-28 사용자 제공)과 같이 **도로 산측
사면 전체를 관(구조물 측점) 개수만큼 빈틈없이 분할**한다:
- 3D 라이다는 산 전체를 계측하지 않으므로 쓰지 않는다. 도엽 데이터만 사용(사용자 확정).
- 함몰 보정은 Whitebox `fill_depressions`를 쓰고(기존 Engine_Skeleton과 동일 패턴),
실패하면 원본 DEM으로 진행한다(결과 저하 가능하나 계산은 지속).
- 유역 폴리곤 외곽선이 곧 분수령(능선)이며 프론트가 능선 스타일로 표시한다.
- 도엽 등고선·표고점을 격자 DEM으로 보간 (3D 라이다는 산 전체를 계측하지 않아 미사용).
- 노선을 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로 도로 셀마다
담당 관을 배정 — "도로에 닿은 물은 측구를 타고 내리막의 첫 관으로 들어간다".
- 사면 각 셀은 D8 흐름으로 내려가 처음 닿는 도로 셀의 관을 물려받는다. 도로를 만나지
못하는 셀(도로 하측 사면, 능선 너머)은 자동 제외.
- 함몰 보정은 Whitebox `fill_depressions`(Engine_Skeleton과 동일 패턴), 실패 시 원본 진행.
- 유역 경계의 산측이 분수령(능선)·지능선이고 하측이 도로선이다. 프론트가 능선 파선 표시.
"""
from __future__ import annotations
@@ -33,10 +35,8 @@ logger = logging.getLogger(__name__)
# 격자 해상도(m)와 최대 격자 크기. 도엽 9매 범위라도 이 상한 안에서 해상도를 낮춰 계산한다.
GRID_RES_M = 10.0
MAX_GRID_CELLS = 1_400_000
# 유역 계산 범위: 측점 bbox + 여유폭(m). 주변 8도엽까지 확보되어 있어 넉넉히 잡는다.
# 유역 계산 범위: 노선 bbox + 여유폭(m). 주변 8도엽까지 확보되어 있어 넉넉히 잡는다.
BBOX_MARGIN_M = 1500.0
# pour point 스냅 반경(m): 측점을 주변 흐름 누적 최대 셀로 옮겨 세류 격자 정합 오차를 흡수.
SNAP_RADIUS_M = 50.0
# DEM 보간 표본 상한(속도 확보용 간축). 초과 시 균등 간격으로 추린다.
MAX_SAMPLE_POINTS = 250_000
@@ -118,15 +118,20 @@ def _collect_samples(
def _build_dem(
samples: np.ndarray,
candidates: list[StructureCandidate],
anchors_x: list[float],
anchors_y: list[float],
) -> tuple[np.ndarray, np.ndarray, np.ndarray, float] | None:
"""표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도)."""
if len(samples) < 10 or not candidates:
"""표본을 격자 DEM으로 보간한다. 반환: (dem, x좌표, y좌표, 해상도).
범위는 앵커(노선 전체 정점) bbox + 여유폭 — 유역이 노선 전 연장을 덮어야 하므로
측점 bbox가 아니라 노선 bbox를 쓴다.
"""
if len(samples) < 10 or not anchors_x:
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 = 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()))
@@ -210,49 +215,75 @@ def _d8_pointer(dem: np.ndarray) -> np.ndarray:
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 _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] = []
for left, right in zip(ordered, ordered[1:]):
window = [
vertex for vertex in vertices if left.chainage_m < vertex.chainage_m < right.chainage_m
]
if window:
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
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 _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(
@@ -380,40 +411,41 @@ def _vectorize_basin(
def build_watershed_basins(
vertices: list[Any],
candidates: list[StructureCandidate],
contour_features: list[dict[str, Any]],
spot_features: list[dict[str, Any]],
elevation_keys: tuple[str, ...],
) -> list[WatershedBasin]:
"""능선(분수령) 기반 배수유역을 산정한다.
"""도로(노선) 기준 배수유역을 산정한다 — 측구 흐름 모델.
반환된 boundary_xy 외곽선이 곧 분수령(능선)이다. 번호는 노선 시점에 가까운 순.
"도로에 닿은 물은 측구를 타고 종단 내리막 방향으로 흘러 첫 번째 관으로 들어간다"를
전제로, 도로 산측 사면 전체를 관 개수만큼 빈틈없이 분할한다(참고 도면과 동일 개념):
① 노선 전체를 도로 셀로 래스터화하고, 관 사이 종단 최고점(물갈림 고개)을 경계로
각 도로 셀에 담당 관 라벨을 붙인다.
② 사면 각 셀은 D8 흐름을 따라 내려가 처음 닿는 도로 셀의 관 라벨을 물려받는다.
도로를 만나지 못하고 격자 밖으로 빠지는 셀(도로 하측 성토부 등)은 제외된다.
반환된 boundary_xy 외곽선의 산측이 곧 분수령(능선)이며, 하측은 도로선을 따른다.
번호는 노선 시점에 가까운 순.
"""
if not candidates:
if not candidates or len(vertices) < 2:
return []
samples = _collect_samples(contour_features, spot_features, elevation_keys)
built = _build_dem(samples, candidates)
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))
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
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)
@@ -421,9 +453,6 @@ def build_watershed_basins(
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:
@@ -431,7 +460,11 @@ def build_watershed_basins(
boundary = _vectorize_basin(labels, label, x_coords, y_coords, resolution)
if len(boundary) < 4:
continue
outlet_z = float(dem[cell[0], cell[1]])
# 측점(관) 위치의 표고를 기준으로 낙차를 계산한다.
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())
basin = WatershedBasin(
index=label,
chainage_m=candidate.chainage_m,
B05_wf2_Route/B05_wf2_Route_Router_Drainage.py
+1 -1
View File
@@ -194,7 +194,7 @@ async def post_drainage_basins(
candidates = propose_structure_stations(vertices, prepared["streams"])
basins = build_watershed_basins(
candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS
vertices, candidates, prepared["contours"], prepared["spots"], _ELEVATION_KEYS
)
to_lonlat = prepared["to_lonlat"]
return {