diff --git a/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts b/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts index fd71c833..2c36b24d 100644 --- a/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts +++ b/B05_wf2_Route/B05_wf2_Route_Api_Fetch.ts @@ -330,14 +330,22 @@ export interface DrainagePrimaryRegion { /** 셀별 흐름 방향과 도로 도달 여부. 등고선이 없어 판정을 못하면 null. */ flow: { encoding: "base64-uint8"; + /** 방위 분해능(32). 코드 0 = 화면 오른쪽, 시계방향 증가. */ + azimuth_steps: number; + /** 제자리(더 낮은 이웃 없음)를 뜻하는 코드. */ + sink_code: number; + /** 표고가 없어 판정 못한 셀 코드. */ + invalid_code: number; cells: number; reaches_road: number; no_road: number; /** 격자에는 있으나 등고선 TIN 밖이라 표고가 없어 판정 못한 셀. */ unanalyzed: number; + /** 확정된 상류 세류망을 따라 흐름을 강제로 새긴 셀 수. */ + burned: number; outer_seeds: number; interior_seeds: number; - /** 셀당 1바이트. 하위 4비트=3×3 방향코드(4=제자리, 15=무효), 0x80=도로 도달. + /** 셀당 1바이트. 하위 6비트=32방위 코드(32=제자리, 33=무효), 0x80=도로 도달. * 순서는 grid.row_spans를 행 → 구간 → 열 오름차순으로 훑은 순서와 같다. */ data: string; } | null; diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py index c15e9db1..3a03ff5b 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Basin.py @@ -38,20 +38,19 @@ from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import ( from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Flow import ( FlowClassification, RoadRaster, - border_contact, + burn_stream_flow, classify_flow, + expand_until_closed, largest_ring, outer_boundary, polygonize_labels, rasterize_road, - trace_flow, ) from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import ( GridSpec, TerrainGrid, build_contour_cloud, build_terrain_grid, - expand_grid_spec, route_elevation_floor, ) from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Stream import ( @@ -254,7 +253,10 @@ def preview_stages( started = time.perf_counter() terrain = build_terrain_grid(spec, cloud, region.cell_mask) road = rasterize_road(spec, route_line, DRAINAGE_ROAD_WIDTH_M) - flow = classify_flow(terrain, road) + # 확정된 상류 세류망을 따라 흐름을 새긴다 — 세류선 위 셀과 그리로 흘러드는 셀은 + # 반드시 도로에 도달해야 한다(TIN 보간면은 실제 물골을 재현하지 못한다). + terrain, burned = burn_stream_flow(terrain, road, region.split.upstream) + flow = classify_flow(terrain, road, burned) logger.info("배수유역: 흐름 판정 %.1fs (셀 %d개)", time.perf_counter() - started, spec.size) return StagePreview(region=region, terrain=terrain, road=road, flow=flow) @@ -291,40 +293,18 @@ def _solve_grid( logger.warning("배수유역: 1차 영역 안에 등고선이 없어 격자 해석을 건너뜁니다.") return None - terrain = road = flow = None - for round_index in range(DRAINAGE_MAX_EXPAND_ROUNDS + 1): - started = time.perf_counter() - terrain = build_terrain_grid(spec, cloud) - road = rasterize_road(spec, route_line, DRAINAGE_ROAD_WIDTH_M) - flow = trace_flow(terrain, road) - # 격자 크기(config DRAINAGE_GRID_SIZE_M)를 조정할 근거가 되도록 회차별 소요를 남긴다. - logger.info( - "배수유역: %d회차 해석 %.1fs (셀 %d개)", - round_index + 1, - time.perf_counter() - started, - spec.size, - ) - contact = border_contact(flow.active) - if not any(contact.values()): - break - if round_index == DRAINAGE_MAX_EXPAND_ROUNDS: - logger.warning( - "배수유역: 확장 상한(%d회)에 도달 — 경계 %s가 아직 활성입니다.", - DRAINAGE_MAX_EXPAND_ROUNDS, - [side for side, touched in contact.items() if touched], - ) - break - widened = expand_grid_spec(spec, contact, DRAINAGE_EXPAND_STEP_M) - if widened == spec: - break - logger.info( - "배수유역: 경계 %s 활성 — %.0fm 확장", - [side for side, touched in contact.items() if touched], - DRAINAGE_EXPAND_STEP_M, - ) - spec = widened - - assert terrain is not None and road is not None and flow is not None + # 확장 루프는 `Watershed_Flow.expand_until_closed()`로 분리했다(2026-07-31 사용자 지시). + # 단계 검증 미리보기(`preview_stages`)는 이 경로를 타지 않는다 — 확장 자체가 아직 검증 대상. + started = time.perf_counter() + expansion = expand_until_closed(spec, cloud, route_line) + spec, terrain, road, flow = expansion.spec, expansion.terrain, expansion.road, expansion.flow + logger.info( + "배수유역: 격자 해석 %.1fs (확장 %d회, %s, 셀 %d개)", + time.perf_counter() - started, + expansion.rounds, + "닫힘" if expansion.closed else "미닫힘", + spec.size, + ) solution = _GridSolution( spec=spec, elevation=terrain.elevation.reshape(-1), diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py index 7edd0366..dbd3340a 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Flow.py @@ -16,6 +16,7 @@ from __future__ import annotations import logging from dataclasses import dataclass +from typing import Any import numpy as np from rasterio.features import rasterize, shapes @@ -24,12 +25,18 @@ from shapely.geometry import LineString, MultiPolygon, Polygon, shape from shapely.ops import unary_union from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import ( + ContourCloud, GridSpec, TerrainGrid, - direction_codes, + build_cell_mask, + build_terrain_grid, + descent_azimuth, + expand_grid_spec, grid_transform, ) from config.config_system import ( + DRAINAGE_EXPAND_STEP_M, + DRAINAGE_MAX_EXPAND_ROUNDS, DRAINAGE_MIN_BASIN_AREA_M2, DRAINAGE_POLYGON_SIMPLIFY_M, DRAINAGE_ROAD_WIDTH_M, @@ -171,15 +178,103 @@ def trace_flow(terrain: TerrainGrid, road: RoadRaster) -> FlowResult: ) +def burn_stream_flow( + terrain: TerrainGrid, road: RoadRaster, streams: list[LineString] +) -> tuple[TerrainGrid, np.ndarray]: + """확정된 상류 세류망을 따라 격자 흐름 방향을 강제로 새긴다. + + 등고선 TIN 보간면은 실제 물골(thalweg)을 그대로 재현하지 못한다. 그래서 세류선 위 + 셀인데도 D8이 옆 사면으로 흘려보내 도로에 닿지 못하는 일이 생긴다. 세류선은 이미 + "도로를 건너 하류로 빠지는 물길"로 확정된 자료이므로, 그 위 셀의 흐름 방향은 추정할 + 것이 아니라 **그대로 따라야 한다**(2026-07-31 사용자 지시). + + 세류선 위 셀은 물길 방향의 다음 셀을 수신 셀로 삼는다. 그러면 세류선 셀은 물론, + 세류선으로 흘러드는 사면 셀까지 전부 도로에 도달한다. 도로 셀은 흡수점이므로 건드리지 않는다. + + 돌려주는 값: (흐름이 새겨진 지형, 새긴 셀 마스크). + """ + spec = terrain.spec + receiver = terrain.receiver.copy() + step_length = terrain.step_length.copy() + valid = terrain.valid.reshape(-1) + road_cells = road.mask.reshape(-1) + burned = np.zeros(spec.size, dtype=bool) + + tails: list[int] = [] + for line in streams: + chain = _line_cell_chain(spec, line) + if len(chain) < 2: + continue + for current, following in zip(chain, chain[1:]): + if not valid[current] or not valid[following] or road_cells[current]: + continue + if not _is_neighbour(spec, current, following): + continue + receiver[current] = following + step_length[current] = _cell_distance(spec, current, following) + burned[current] = True + tails.append(chain[-1]) + + # 하류 끝이 도로에도, 다른 세류 갈래에도 닿지 않은 갈래만 진짜 문제다. + # 지류가 본류 중간에 합류하는 경우는 끝점이 본류 셀이므로 정상이다. + detached = sum(1 for tail in tails if not road_cells[tail] and not burned[tail]) + if detached: + logger.warning( + "배수유역: 세류망 %d갈래의 하류 끝이 도로·다른 세류 어디에도 닿지 않습니다.", + detached, + ) + logger.info("배수유역: 세류망 흐름 새김 %d셀 (세류 %d갈래)", int(burned.sum()), len(streams)) + return ( + TerrainGrid( + spec=spec, + elevation=terrain.elevation, + valid=terrain.valid, + receiver=receiver, + step_length=step_length, + ), + burned, + ) + + +def _line_cell_chain(spec: GridSpec, line: LineString) -> list[int]: + """선을 따라 지나가는 셀을 순서대로 뽑는다(연속 중복 제거).""" + step = max(spec.cell_m / 2.0, 0.1) + positions = np.arange(0.0, line.length + step, step) + chain: list[int] = [] + for position in positions: + point = line.interpolate(float(position)) + col = int((point.x - spec.x_min) // spec.cell_m) + row = int((spec.y_max - point.y) // spec.cell_m) + if not (0 <= row < spec.n_rows and 0 <= col < spec.n_cols): + continue + index = row * spec.n_cols + col + if not chain or chain[-1] != index: + chain.append(index) + return chain + + +def _is_neighbour(spec: GridSpec, first: int, second: int) -> bool: + row_delta = abs(first // spec.n_cols - second // spec.n_cols) + col_delta = abs(first % spec.n_cols - second % spec.n_cols) + return max(row_delta, col_delta) == 1 + + +def _cell_distance(spec: GridSpec, first: int, second: int) -> float: + row_delta = abs(first // spec.n_cols - second // spec.n_cols) + col_delta = abs(first % spec.n_cols - second % spec.n_cols) + return spec.cell_m * float(np.hypot(row_delta, col_delta)) + + @dataclass class FlowClassification: """셀별 흐름 방향과 도로 도달 여부 — 확장 없이 현재 격자만 본 결과.""" - direction: np.ndarray # (R*C,) int8 — 3×3 방향 코드(4=제자리), 무효 셀은 −1 + direction: np.ndarray # (R*C,) int16 — 32방위 코드(32=제자리, 33=무효) reaches_road: np.ndarray # (R*C,) bool — 물길을 따라가면 도로에 닿는가 analyzed: np.ndarray # (R*C,) bool — 실제로 판정한 셀 outer_seeds: int # 최외곽에서 출발해 판정한 셀 수 interior_seeds: int # 최외곽 추적에 안 걸려 따로 출발시킨 내부 셀 수 + burned: np.ndarray | None = None # (R*C,) bool — 세류망을 따라 흐름을 새긴 셀 def outermost_cells(domain: np.ndarray) -> np.ndarray: @@ -196,33 +291,54 @@ def outermost_cells(domain: np.ndarray) -> np.ndarray: return exposed & domain -def classify_flow(terrain: TerrainGrid, road: RoadRaster) -> FlowClassification: +def classify_flow( + terrain: TerrainGrid, road: RoadRaster, burned: np.ndarray | None = None +) -> FlowClassification: """최외곽 셀부터 물길을 따라가며 도로 도달 여부를 판정한다. - 최외곽 셀에서 출발해 D8 수신 셀을 한 칸씩 따라가고, 경로가 끝나면 그 결과를 경로 전체에 - 되돌려 적는다. **한 번 판정한 셀은 다시 분석하지 않는다** — 다른 경로가 그 셀을 만나면 - 거기서 즉시 결론을 가져온다. 최외곽 추적으로 안 닿은 내부 셀은 그다음에 따로 출발시킨다 + 물이 흐르는 순서 그대로 따라간다 — 셀이 각자 도로를 바라보는 게 아니라, 한 셀에서 + 출발해 화살표가 가리키는 다음 셀, 그 셀이 가리키는 그다음 셀로 사슬처럼 이어 간다 (2026-07-31 사용자 지시). - 경로가 도로 셀에 닿으면 경로 전체가 도로 도달(적색), 싱크에서 멈추거나 해석 영역 밖으로 - 나가면 미도달(파랑)이다. 격자 확장은 하지 않는다 — 다음 단계다. + 사슬이 **도로 셀 또는 확정된 세류망 셀**에 닿으면 그 사슬 전체가 적색이다. 세류망은 + 이미 "도로를 건너 하류로 빠지는 물길"로 확정된 자료이므로, 물이 세류에 합류한 시점에 + 도로 도달이 결정된다. 사슬이 싱크에서 멈추거나 해석 영역 밖으로 나가면 전체가 미도달 + (파랑 채움 + 백색 화살표)이다. + + **한 번 판정한 셀은 다시 분석하지 않는다** — 사슬이 이미 색이 정해진 셀을 만나면 그 + 셀의 색을 그대로 물려받고 끝낸다. 최외곽 추적으로 안 닿은 내부 셀은 그다음에 따로 + 출발시킨다. + + 화살표 방위는 D8(8방위)이 아니라 지형 최급강하 32방위를 쓴다. `burned`(세류망을 따라 + 흐름을 새긴 셀)는 확정된 물길 방향을 그대로 쓴다. """ spec = terrain.spec valid = terrain.valid.reshape(-1) receiver = terrain.receiver - is_road = road.mask.reshape(-1) & valid + # 도로 셀과 세류망 셀 둘 다 "여기 닿으면 적색"인 종결점이다. + absorbing = road.mask.reshape(-1) & valid + if burned is not None: + absorbing = absorbing | (burned & valid) - direction = np.where(valid, direction_codes(spec, receiver), -1).astype(np.int8) + direction = descent_azimuth(spec, terrain.elevation, terrain.valid, receiver, burned) reaches = np.zeros(spec.size, dtype=bool) # 0=미방문, 1=경로에 올라 있음, 2=판정 완료 state = np.zeros(spec.size, dtype=np.int8) outer = np.flatnonzero(outermost_cells(terrain.valid)) remaining = np.flatnonzero(valid) - outer_done = _walk_from(outer, receiver, valid, is_road, state, reaches) - interior_done = _walk_from(remaining, receiver, valid, is_road, state, reaches) + outer_done = _walk_from(outer, receiver, valid, absorbing, state, reaches) + interior_done = _walk_from(remaining, receiver, valid, absorbing, state, reaches) analyzed = state == 2 + if burned is not None: + stranded = int((burned & analyzed & ~reaches).sum()) + if stranded: + logger.warning( + "배수유역: 세류망 새김 셀 %d개가 도로에 닿지 않습니다 — 하류 끝이 도로 셀과 " + "이어지지 않았는지 확인이 필요합니다.", + stranded, + ) logger.info( "배수유역: 흐름 판정 %d셀 (최외곽 출발 %d / 내부 보충 %d) — 도로 도달 %d, 미도달 %d", int(analyzed.sum()), @@ -237,6 +353,7 @@ def classify_flow(terrain: TerrainGrid, road: RoadRaster) -> FlowClassification: analyzed=analyzed, outer_seeds=outer_done, interior_seeds=interior_done, + burned=burned, ) @@ -244,11 +361,15 @@ def _walk_from( starts: np.ndarray, receiver: np.ndarray, valid: np.ndarray, - is_road: np.ndarray, + absorbing: np.ndarray, state: np.ndarray, reaches: np.ndarray, ) -> int: - """출발 셀 목록에서 물길을 따라가며 판정한다. 새로 판정한 셀 수를 돌려준다.""" + """출발 셀에서 물길 사슬을 따라가며 판정하고, 사슬 전체에 같은 색을 적는다. + + `absorbing`은 도로 셀과 확정된 세류망 셀 — 여기 닿으면 사슬 전체가 도로 도달이다. + 새로 판정한 셀 수를 돌려준다. + """ resolved = 0 path: list[int] = [] for start in starts.tolist(): @@ -258,15 +379,15 @@ def _walk_from( node = start while True: if state[node] == 2: - verdict = bool(reaches[node]) + verdict = bool(reaches[node]) # 이미 색이 정해진 셀 — 그 색을 물려받는다 break if state[node] == 1: # 방어: 채움·평탄해소 후에는 순환이 없어야 한다 verdict = False break state[node] = 1 path.append(node) - if is_road[node]: - verdict = True + if absorbing[node]: + verdict = True # 도로 또는 세류망에 합류 — 여기서 하류로 빠진다 break following = int(receiver[node]) if following == node or not valid[following]: @@ -280,6 +401,77 @@ def _walk_from( return resolved +# ── 격자 확장 (별도 단계) ─────────────────────────────────────────────────── + + +@dataclass +class ExpansionResult: + """확장 루프 결과. 확장을 쓰지 않는 경로에서는 이 모듈을 부르지 않는다.""" + + spec: GridSpec + terrain: TerrainGrid + road: RoadRaster + flow: FlowResult + rounds: int # 실제로 넓힌 횟수 (0 = 처음부터 닫혀 있었음) + closed: bool # 경계 링이 전부 비활성이 되어 스스로 멈췄는가 + + +def expand_until_closed( + spec: GridSpec, + cloud: ContourCloud, + route_line: LineString, + region_area: Any = None, + road_width_m: float = DRAINAGE_ROAD_WIDTH_M, + max_rounds: int = DRAINAGE_MAX_EXPAND_ROUNDS, + step_m: float = DRAINAGE_EXPAND_STEP_M, +) -> ExpansionResult: + """활성 셀이 격자 최외곽에 닿은 방향으로만 넓히며 유역이 닫힐 때까지 반복한다. + + **본 계산 경로에서만 쓰는 별도 단계다**(2026-07-31 사용자 지시로 분리). 단계 검증 + 미리보기는 확장 없이 현재 격자만 본다 — 확장 로직 자체가 아직 검증 대상이기 때문이다. + + 종료 조건은 반경 상한이 아니라 **경계 링 전체가 비활성**이 되는 것이다. 비활성 셀이 + 최외곽에 띠로 완성되면 그 바깥은 볼 필요가 없다. `max_rounds`는 무한 반복 방지용이다. + + `region_area`를 주면 확장된 격자에서도 그 영역에 걸치는 셀만 해석 대상으로 삼는다. + """ + terrain = road = flow = None + rounds = 0 + closed = False + for attempt in range(max_rounds + 1): + domain = build_cell_mask(spec, region_area) if region_area is not None else None + terrain = build_terrain_grid(spec, cloud, domain) + road = rasterize_road(spec, route_line, road_width_m) + flow = trace_flow(terrain, road) + contact = border_contact(flow.active) + if not any(contact.values()): + closed = True + break + if attempt == max_rounds: + logger.warning( + "배수유역: 확장 상한(%d회) 도달 — 경계 %s가 아직 활성입니다.", + max_rounds, + [side for side, touched in contact.items() if touched], + ) + break + widened = expand_grid_spec(spec, contact, step_m) + if widened == spec: + break + logger.info( + "배수유역: 경계 %s 활성 — %.0fm 확장 (%d회차)", + [side for side, touched in contact.items() if touched], + step_m, + attempt + 1, + ) + spec = widened + rounds += 1 + + assert terrain is not None and road is not None and flow is not None + return ExpansionResult( + spec=spec, terrain=terrain, road=road, flow=flow, rounds=rounds, closed=closed + ) + + def border_contact(active: np.ndarray) -> dict[str, bool]: """활성 셀이 격자 최외곽에 닿은 방향. 전부 False면 유역이 능선 안에서 닫힌 것이다.""" return { diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py index b4e72984..b6498d83 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Grid.py @@ -495,16 +495,58 @@ def build_terrain_grid( ) -def direction_codes(spec: GridSpec, receiver: np.ndarray) -> np.ndarray: - """수신 셀 인덱스를 3×3 방향 코드로 바꾼다. +# 화살표 방위 분해능. 0 = 화면상 오른쪽(+열), 시계방향으로 증가(행이 아래로 증가하므로). +AZIMUTH_STEPS = 32 +# 방위 코드 특수값. +AZIMUTH_SINK = AZIMUTH_STEPS # 32 = 제자리(더 낮은 이웃 없음) +AZIMUTH_INVALID = AZIMUTH_STEPS + 1 # 33 = 표고 없음(해석 불가) - 코드 = (행증분 + 1) * 3 + (열증분 + 1) → 0~8. 4는 제자리(싱크)를 뜻한다. - 화면이 셀마다 화살표를 그릴 수 있도록 방향만 뽑아낸 표현이다. + +def descent_azimuth( + spec: GridSpec, + surface: np.ndarray, + valid: np.ndarray, + receiver: np.ndarray, + forced: np.ndarray | None = None, +) -> np.ndarray: + """셀별 물 흐름 방위를 32방위 코드로 낸다. + + D8은 연결(도로 도달 판정)에는 충분하지만 화면에 8방위밖에 못 그린다. 실제 지표수는 + 지형 최급강하 방향으로 흐르고 그 방향은 연속값이므로, **표시는 지표면 기울기에서 뽑은 + 연속 방위를 32단계로 양자화**해 보여 준다(2026-07-31 사용자 지시). + + `forced`(세류망을 따라 흐름을 새긴 셀)는 기울기 대신 실제 수신 셀 방향을 쓴다 — 그 + 셀들은 지형 추정이 아니라 확정된 물길을 따르기 때문이다. 기울기가 0에 가까운 셀도 + 수신 셀 방향으로 대체한다. """ + rows, cols = spec.n_rows, spec.n_cols + filled = np.where(valid, surface, np.nan) + # np.gradient는 NaN이 번지므로 무효 셀을 주변 유효값으로 임시 대체한 뒤 기울기를 잡는다. + working = np.where(np.isfinite(filled), filled, np.nanmean(filled) if valid.any() else 0.0) + grad_row, grad_col = np.gradient(working.astype(np.float64), spec.cell_m) + # 내리막 방향 = 기울기 반대. 행은 아래로 증가하므로 화면 좌표와 부호가 같다. + move_row = -grad_row + move_col = -grad_col + magnitude = np.hypot(move_row, move_col) + index = np.arange(receiver.size, dtype=np.int64) - row_delta = receiver // spec.n_cols - index // spec.n_cols - col_delta = receiver % spec.n_cols - index % spec.n_cols - return ((row_delta + 1) * 3 + (col_delta + 1)).astype(np.int8) + receiver_row = (receiver // cols - index // cols).reshape(rows, cols).astype(np.float64) + receiver_col = (receiver % cols - index % cols).reshape(rows, cols).astype(np.float64) + use_receiver = magnitude < 1e-9 + if forced is not None: + use_receiver |= forced.reshape(rows, cols) + move_row = np.where(use_receiver, receiver_row, move_row) + move_col = np.where(use_receiver, receiver_col, move_col) + + angle = np.arctan2(move_row, move_col) + code = np.rint(angle / (2.0 * math.pi / AZIMUTH_STEPS)).astype(np.int16) % AZIMUTH_STEPS + # 수신 셀이 자기 자신이거나 이동량이 없는 셀은 방향이 없다. + is_sink = (receiver.reshape(rows, cols) == index.reshape(rows, cols)) | ( + (np.abs(move_row) < 1e-12) & (np.abs(move_col) < 1e-12) + ) + code = np.where(is_sink, AZIMUTH_SINK, code) + code = np.where(valid, code, AZIMUTH_INVALID) + return code.reshape(-1).astype(np.int16) def route_elevation_floor(route_z_values: list[float]) -> float | None: diff --git a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py index cdbedbb6..0d337bd0 100644 --- a/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py +++ b/B05_wf2_Route/B05_wf2_Route_Engine_Watershed_Stream.py @@ -37,7 +37,11 @@ logger = logging.getLogger(__name__) @dataclass class StreamSplit: - """세류망을 도로 교차점에서 상·하류로 가른 결과. 검증 화면이 그대로 그린다.""" + """세류망을 도로 교차점에서 상·하류로 가른 결과. 검증 화면이 그대로 그린다. + + `upstream`은 **물이 흐르는 방향(상류 → 하류)으로 정렬**돼 있다. 마지막 좌표가 도로에 + 가까운 끝이다. 격자 흐름에 세류 방향을 새겨 넣을 때 이 순서를 그대로 쓴다. + """ upstream: list[LineString] = field(default_factory=list) # 채택 — 1차 영역의 기준 downstream: list[LineString] = field(default_factory=list) # 도로 아래로 이어진 망 @@ -87,20 +91,24 @@ def split_streams_at_road( node_edges.setdefault(tail, []).append(index) sampler = ElevationSampler(cloud) - upper_seeds: set[int] = set() - lower_seeds: set[int] = set() + # 씨앗 조각 → 그 조각의 하류쪽 끝점(= 도로 교차 노드). 이 값이 물 흐름 방향의 기준이 된다. + upper_seeds: dict[int, tuple[float, float]] = {} + lower_seeds: dict[int, tuple[float, float]] = {} 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) + heights = sampler.at(np.array(touching, dtype=np.float64)) + crossing_node = touching[int(np.argmin(heights))] + if _mean_elevation(piece, sampler) > float(np.min(heights)): + upper_seeds[index] = crossing_node else: - lower_seeds.add(index) + lower_seeds[index] = crossing_node - upstream = _spread_network(upper_seeds, ends, node_edges, crossing_nodes) - downstream = _spread_network(lower_seeds, ends, node_edges, crossing_nodes) - upstream + upstream_flow = _spread_network(upper_seeds, ends, node_edges, crossing_nodes) + downstream_flow = _spread_network(lower_seeds, ends, node_edges, crossing_nodes) + upstream = set(upstream_flow) + downstream = set(downstream_flow) - upstream logger.info( "배수유역: 세류 조각 %d개 → 상류망 %d개 채택 / 하류망 %d개 · 미연결 %d개 제외", len(pieces), @@ -109,12 +117,26 @@ def split_streams_at_road( len(pieces) - len(upstream) - len(downstream), ) return StreamSplit( - upstream=[pieces[index] for index in sorted(upstream)], + # 상류망은 물 흐름 방향(상류 → 하류)으로 뒤집어 둔다 — 격자 흐름 새김에 그대로 쓴다. + upstream=[ + _oriented(pieces[index], upstream_flow[index], ends[index]) + for index in sorted(upstream) + ], downstream=[pieces[index] for index in sorted(downstream)], no_contact=len(pieces) - len(upstream) - len(downstream), ) +def _oriented( + piece: LineString, + downstream_node: tuple[float, float], + piece_ends: tuple[tuple[float, float], tuple[float, float]], +) -> LineString: + """조각을 하류쪽 끝이 마지막 좌표가 되도록 정렬한다.""" + head, _tail = piece_ends + return LineString(list(piece.coords)[::-1]) if head == downstream_node else piece + + def _cut_network_at_road( lines: list[LineString], route_line: LineString ) -> tuple[list[LineString], set[tuple[float, float]]]: @@ -151,13 +173,18 @@ def _cut_network_at_road( def _spread_network( - seeds: set[int], + seeds: dict[int, tuple[float, float]], 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) +) -> dict[int, tuple[float, float]]: + """씨앗 조각에서 끝점을 타고 퍼진다. 도로 교차 노드는 통과하지 않는다. + + 조각마다 **어느 끝점을 통해 도달했는지**를 함께 기록한다. 그 끝점이 도로에 더 가까운 + 쪽이므로 곧 그 조각의 하류 방향이다 — 세류망 전체의 물 흐름 방향이 이 한 번의 확산으로 + 같이 정해진다. + """ + downstream = dict(seeds) queue = list(seeds) while queue: index = queue.pop() @@ -165,10 +192,11 @@ def _spread_network( 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 + if neighbour in downstream: + continue + downstream[neighbour] = node + queue.append(neighbour) + return downstream def _node_key(x: float, y: float) -> tuple[float, float]: diff --git a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py index fb316b9c..d3862598 100644 --- a/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py +++ b/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py @@ -30,7 +30,12 @@ from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Basin import ( preview_stages, ) from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Export import write_grid_arrays, write_stage -from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import mask_row_spans +from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import ( + AZIMUTH_INVALID, + AZIMUTH_SINK, + AZIMUTH_STEPS, + mask_row_spans, +) from B05_wf2_Route.B05_wf2_Route_Repository import ( get_latest_route, get_route_points, @@ -282,19 +287,26 @@ def _write_stage_arrays(stored_path: str, preview: Any, region: Any, spec: Any) flow = preview.flow if flow is None: return + arrays = { + "direction": flow.direction.reshape(spec.n_rows, spec.n_cols), + "reaches_road": flow.reaches_road.reshape(spec.n_rows, spec.n_cols), + "analyzed": flow.analyzed.reshape(spec.n_rows, spec.n_cols), + } + if flow.burned is not None: + arrays["burned"] = flow.burned.reshape(spec.n_rows, spec.n_cols) write_grid_arrays( stored_path, "flow_direction", spec, + arrays, { - "direction": flow.direction.reshape(spec.n_rows, spec.n_cols), - "reaches_road": flow.reaches_road.reshape(spec.n_rows, spec.n_cols), - "analyzed": flow.analyzed.reshape(spec.n_rows, spec.n_cols), - }, - { + "azimuth_steps": AZIMUTH_STEPS, + "sink_code": AZIMUTH_SINK, + "invalid_code": AZIMUTH_INVALID, "analyzed": int(flow.analyzed.sum()), "reaches_road": int((flow.reaches_road & flow.analyzed).sum()), "no_road": int((~flow.reaches_road & flow.analyzed).sum()), + "burned": 0 if flow.burned is None else int(flow.burned.sum()), "outer_seeds": flow.outer_seeds, "interior_seeds": flow.interior_seeds, }, @@ -305,26 +317,31 @@ def _flow_payload(preview: Any, region: Any) -> dict[str, Any] | None: """셀별 흐름 방향·도로 도달 여부를 바이트 배열로 압축한다. 셀이 수십만 개라 JSON 객체로는 못 보낸다. 셀 하나당 1바이트로 줄이고 base64로 싣는다: - 하위 4비트 = 3×3 방향 코드(0~8, 4=제자리), 15 = 무효(표고 없음) - 최상위 비트(0x80) = 도로 도달(적색). 꺼져 있으면 미도달(파랑). + 하위 6비트(0x3F) = 32방위 코드(0~31, 0=화면 오른쪽·시계방향), 32=제자리, 33=표고 없음 + 최상위 비트(0x80) = 도로 도달(적색). 꺼져 있으면 미도달(백색 화살표). 바이트 순서는 `grid.row_spans`를 행 → 구간 → 열 오름차순으로 훑은 순서와 같다. """ flow = preview.flow if flow is None or region.cell_mask is None: return None order = np.flatnonzero(region.cell_mask.reshape(-1)) - codes = flow.direction[order] analyzed = flow.analyzed[order] reaches = flow.reaches_road[order] - packed = np.where(codes < 0, 15, codes).astype(np.uint8) + packed = np.clip(flow.direction[order], 0, AZIMUTH_INVALID).astype(np.uint8) packed |= np.where(reaches, 0x80, 0).astype(np.uint8) + burned = flow.burned return { "encoding": "base64-uint8", + "azimuth_steps": AZIMUTH_STEPS, + "sink_code": AZIMUTH_SINK, + "invalid_code": AZIMUTH_INVALID, "cells": int(order.size), "reaches_road": int((reaches & analyzed).sum()), "no_road": int((~reaches & analyzed).sum()), - # 격자에는 있으나 등고선 TIN 밖이라 표고가 없어 판정 못한 셀(방향 코드 15). + # 격자에는 있으나 등고선 TIN 밖이라 표고가 없어 판정 못한 셀. "unanalyzed": int((~analyzed).sum()), + # 확정된 상류 세류망을 따라 흐름을 강제로 새긴 셀 수. + "burned": 0 if burned is None else int(burned[order].sum()), "outer_seeds": flow.outer_seeds, "interior_seeds": flow.interior_seeds, "data": base64.b64encode(packed.tobytes()).decode("ascii"), diff --git a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts index b5afcb29..ff267960 100644 --- a/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts +++ b/B05_wf2_Route/B05_wf2_Route_UI_Drainage_Panel.ts @@ -53,11 +53,11 @@ const COLLAPSED_KEY = "b05-route-drainage-collapsed"; // 해석 격자 셀 선 — 등고선·세류 위에 얹으므로 흰색으로 둔다(2026-07-31 사용자 지시). const GRID_LINE_COLOR = "rgba(255, 255, 255, 0.55)"; -// 흐름 판정 색 — 도로로 물이 오는 셀은 적색, 오지 않는 셀은 파랑. +// 흐름 판정 색 — 도로로 물이 오는 셀은 적색, 오지 않는 셀은 파랑 채움 + 백색 화살표. const FLOW_TO_ROAD_FILL = "rgba(220, 38, 38, 0.28)"; const FLOW_TO_ROAD_LINE = "rgba(153, 27, 27, 0.95)"; const FLOW_AWAY_FILL = "rgba(37, 99, 235, 0.22)"; -const FLOW_AWAY_LINE = "rgba(30, 64, 175, 0.9)"; +const FLOW_AWAY_LINE = "rgba(255, 255, 255, 0.95)"; /** 등고선 TIN 밖이라 표고가 없어 판정하지 못한 셀 — 미도달(파랑)과 구분한다. */ const FLOW_UNKNOWN_FILL = "rgba(120, 113, 108, 0.18)"; /** 셀이 이보다 작으면 화살표가 뭉개져 읽히지 않으므로 채움색만 남긴다(px). */ @@ -357,6 +357,9 @@ export function createDrainagePanel(): DrainagePanel { } return; } + const sink = region.flow?.sink_code ?? 32; + const invalid = region.flow?.invalid_code ?? 33; + const steps = region.flow?.azimuth_steps ?? 32; for (let offset = 0; offset < count; offset += 1) { drawFlowCell( context, @@ -366,13 +369,14 @@ export function createDrainagePanel(): DrainagePanel { cellW, cellH, cellPx, + { sink, invalid, steps }, ); } }); context.restore(); } - /** 셀 하나 — 도달 여부로 칠하고, 충분히 크면 흐름 방향 화살표를 얹는다. */ + /** 셀 하나 — 도달 여부로 칠하고, 충분히 크면 32방위 흐름 화살표를 얹는다. */ function drawFlowCell( context: CanvasRenderingContext2D, code: number, @@ -381,11 +385,12 @@ export function createDrainagePanel(): DrainagePanel { cellW: number, cellH: number, cellPx: number, + codes: { sink: number; invalid: number; steps: number }, ): void { - const direction = code & 0x0f; + const azimuth = code & 0x3f; const reaches = (code & 0x80) !== 0; - // 코드 15 = 등고선 TIN 밖이라 표고가 없어 판정 못한 셀. 파랑(미도달)과 구분한다. - const unanalyzed = direction === 15; + // 표고가 없어 판정 못한 셀 — 미도달(파랑 채움)과 구분해야 오독이 없다. + const unanalyzed = azimuth === codes.invalid; context.fillStyle = unanalyzed ? FLOW_UNKNOWN_FILL : reaches @@ -397,11 +402,11 @@ export function createDrainagePanel(): DrainagePanel { context.lineWidth = 0.5; context.strokeRect(x, y, cellW, cellH); } - if (cellPx < ARROW_MIN_PX || direction === 15) return; + if (cellPx < ARROW_MIN_PX || unanalyzed) return; const stroke = reaches ? FLOW_TO_ROAD_LINE : FLOW_AWAY_LINE; const midX = x + cellW / 2; const midY = y + cellH / 2; - if (direction === 4) { + if (azimuth === codes.sink) { // 제자리(싱크) — 방향이 없으므로 점으로 표시한다. context.fillStyle = stroke; context.beginPath(); @@ -409,27 +414,26 @@ export function createDrainagePanel(): DrainagePanel { context.fill(); return; } - const colDelta = (direction % 3) - 1; - const rowDelta = Math.floor(direction / 3) - 1; - const length = Math.hypot(colDelta, rowDelta) || 1; - const reach = (cellPx * 0.38) / length; - const tipX = midX + colDelta * reach; - const tipY = midY + rowDelta * reach; + // 코드 0 = 화면 오른쪽(+x), 시계방향(캔버스 y는 아래가 +). + const angle = (azimuth * 2 * Math.PI) / codes.steps; + const unitX = Math.cos(angle); + const unitY = Math.sin(angle); + const reach = cellPx * 0.38; + const tipX = midX + unitX * reach; + const tipY = midY + unitY * reach; context.strokeStyle = stroke; context.lineWidth = Math.max(0.6, cellPx * 0.09); context.beginPath(); - context.moveTo(midX - colDelta * reach, midY - rowDelta * reach); + context.moveTo(midX - unitX * reach, midY - unitY * reach); context.lineTo(tipX, tipY); context.stroke(); // 촉 — 진행 방향 기준 좌우로 짧게 접는다. - const head = cellPx * 0.16; - const unitX = (colDelta / length) * head; - const unitY = (rowDelta / length) * head; + const head = cellPx * 0.18; context.beginPath(); context.moveTo(tipX, tipY); - context.lineTo(tipX - unitX - unitY * 0.7, tipY - unitY + unitX * 0.7); + context.lineTo(tipX - (unitX + unitY * 0.7) * head, tipY - (unitY - unitX * 0.7) * head); context.moveTo(tipX, tipY); - context.lineTo(tipX - unitX + unitY * 0.7, tipY - unitY - unitX * 0.7); + context.lineTo(tipX - (unitX - unitY * 0.7) * head, tipY - (unitY + unitX * 0.7) * head); context.stroke(); } @@ -624,11 +628,15 @@ export function createDrainagePanel(): DrainagePanel { region.flow && region.flow.unanalyzed > 0 ? ` / 표고없음 ${region.flow.unanalyzed.toLocaleString()}(회)` : ""; + const burned = + region.flow && region.flow.burned > 0 + ? ` · 세류망 새김 ${region.flow.burned.toLocaleString()}셀` + : ""; const flow = region.flow ? ` · 흐름 도로도달 ${region.flow.reaches_road.toLocaleString()}(적) / ` + `미도달 ${region.flow.no_road.toLocaleString()}(청)${unknown}, ` + `최외곽 출발 ${region.flow.outer_seeds.toLocaleString()} + ` + - `내부 보충 ${region.flow.interior_seeds.toLocaleString()}` + `내부 보충 ${region.flow.interior_seeds.toLocaleString()}${burned}` : " · 흐름 판정 없음"; return ( `1차 영역(반경 ${region.radius_m}m): 상류망 ${region.upstream_lines.length}조각 채택 / ` +