diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Analyze.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Analyze.py index 9d386804c..72b5cb1ed 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Analyze.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Analyze.py @@ -96,7 +96,9 @@ def resolve_primary_region( """ floor = route_elevation_floor([vertex.z for vertex in vertices]) near_bounds = route_line.buffer(DRAINAGE_INITIAL_RADIUS_M * 2.0).bounds + cloud_started = time.perf_counter() cloud = build_contour_cloud(contour_features, floor, near_bounds) + logger.info("[계측] 등고선 점 구름 %.2fs", time.perf_counter() - cloud_started) if cloud.is_empty: logger.warning("배수유역: 노선 주변에 등고선이 없어 1차 영역을 정할 수 없습니다.") return None diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py index 5afc21cdb..9c52ab4d0 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py @@ -365,12 +365,12 @@ _DESCENT_KEEP = 2 def descent_input_key( contour_key: str, contour_features: list[dict[str, Any]], - elevation_floor_m: float | None, cell_m: float, ) -> str: - """방향장 입력 지문 — 도엽 등고선 지문 + 남긴 등고 단 + 셀 + 계산 규칙 값. + """방향장 입력 지문 — 도엽 등고선 지문 + 등고 단 + 셀 + 계산 규칙 값. - 남긴 단은 속성 표고만 본다 — 같은 도엽 · 같은 단 집합이면 남는 선도 같다(길이 거름은 결정적). + 노선과 무관 — 노선 최저점 하한을 두지 않아 프로젝트마다 한 종류(520-6). + 단은 속성 표고만 본다 — 같은 도엽 · 같은 단 집합이면 선도 같다(길이 거름은 결정적). """ levels = sorted( { @@ -378,7 +378,6 @@ def descent_input_key( for feature in contour_features if feature.get("geometry") and (elevation := _feature_elevation(feature.get("properties") or {})) is not None - and (elevation_floor_m is None or elevation >= elevation_floor_m) } ) text = json.dumps( diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py index ebcc2c8fc..8865fb76f 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Expand.py @@ -14,7 +14,11 @@ 하강 방향장은 처음에 **넉넉한 큰 격자**로 한 번 짓고, 회차마다 그 격자에서 잘라 쓴다 (`DESCENT_SPARE_M` · 342-3). 격자가 커질 때마다 다시 짓던 것(8회 · 30 s)을 한 번으로 줄인다. 저장소(`descent_store` = 폴더 · 도엽 등고선 지문)를 주면 지은 것을 남기고, 덮는 저장본이 있으면 -읽기만 한다(347-1). +읽기만 한다(347-1). 방향장은 노선 최저점 하한 없이 짓는다 — 노선이 바뀌어도 같은 저장본(520-6). + +확장은 **큰 방향장 위 한 번 셈**으로 낸다(520-6) — 큰 격자 전체에서 도로 도달을 한 번 판정하고 +「도로로 오는 셀 + 확장 폭 띠 ∪ 1차 영역」으로 잘라 돌려준다. 도로로 오는 셀이 큰 격자 +가장자리에 닿으면 바깥에서 더 들어올 수 있으므로 위 회차 확장(넓혀 다시 짓기)으로 간다. """ from __future__ import annotations @@ -27,6 +31,7 @@ from pathlib import Path from typing import Any import numpy as np +from scipy.ndimage import distance_transform_cdt from shapely.geometry import LineString from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import ( @@ -91,7 +96,6 @@ def analyze_domain( contour_features: list[dict[str, Any]], route_line: LineString, upstream_streams: list[LineString], - elevation_floor_m: float | None = None, descent: ContourDescent | None = None, ) -> GridAnalysis | None: """주어진 격자·해석 영역에 대해 등고선 하강 방향과 도로 도달 색을 한 번 계산한다. @@ -105,7 +109,7 @@ def analyze_domain( marks = [("시작", time.perf_counter())] fitted = crop_descent(descent, spec) if descent is not None else None if fitted is None: - descent = build_contour_descent(spec, contour_features, None, elevation_floor_m) + descent = build_contour_descent(spec, contour_features) else: descent = fitted marks.append( @@ -149,7 +153,10 @@ def expand_by_red_boundary( descent_store: tuple[Path, str] | None = None, should_stop: Callable[[], bool] | None = None, ) -> RedExpansion | None: - """최외곽 적색 셀 주변으로 넓히며, 새로 추가한 셀에 적색이 없을 때까지 반복한다.""" + """큰 방향장 위 한 번 셈 — 가장자리에 닿으면 최외곽 적색 셀 주변으로 넓히며 반복한다. + + `elevation_floor_m` 은 받기만 한다 — 방향장은 노선과 무관(520-6). + """ _expand_started = time.perf_counter() band_cells = max(1, int(round(band_m / spec.cell_m))) # 1차 영역의 bbox는 영역에 딱 붙어 있어 첫 회차부터 격자를 넓혀야 한다. 미리 여유를 @@ -162,7 +169,7 @@ def expand_by_red_boundary( folder, contour_key = descent_store store = ( folder, - descent_input_key(contour_key, contour_features, elevation_floor_m, spec.cell_m), + descent_input_key(contour_key, contour_features, spec.cell_m), ) def source_for(grid: GridSpec, mask: np.ndarray) -> ContourDescent: @@ -176,17 +183,21 @@ def expand_by_red_boundary( saved.spec.n_cols, ) return saved - built = build_contour_descent( - _pad_spec(grid, mask, spare_cells)[0], contour_features, None, elevation_floor_m - ) + built = build_contour_descent(_pad_spec(grid, mask, spare_cells)[0], contour_features) if store is not None: save_descent(*store, built) return built source = source_for(spec, domain) - analysis = analyze_domain( - spec, domain, contour_features, route_line, upstream_streams, elevation_floor_m, source + if should_stop is not None and should_stop(): + raise AnalysisStopped + once = _expand_once( + source, spec, domain, contour_features, route_line, upstream_streams, band_cells ) + if once is not None: + logger.info("[계측] 확장 한 번 셈 %.1fs", time.perf_counter() - _expand_started) + return once + analysis = analyze_domain(spec, domain, contour_features, route_line, upstream_streams, source) if analysis is None: return None @@ -235,7 +246,6 @@ def expand_by_red_boundary( contour_features, route_line, upstream_streams, - elevation_floor_m, descent=source, ) if widened_analysis is None: @@ -274,6 +284,75 @@ def expand_by_red_boundary( return RedExpansion(analysis=analysis, rounds=rounds, closed=closed, added_cells=added) +def _expand_once( + source: ContourDescent, + spec: GridSpec, + domain: np.ndarray, + contour_features: list[dict[str, Any]], + route_line: LineString, + upstream_streams: list[LineString], + band_cells: int, +) -> RedExpansion | None: + """큰 방향장 전체에서 도로 도달을 한 번 판정 → 「도로로 오는 셀 + 띠 ∪ 1차 영역」으로 자름. + + 도로로 오는 셀은 물길 아래쪽으로 닫혀 있다(도달 셀의 다음 셀도 도달) — 잘라도 그 물길이 + 영역 안에 다 있어 판정이 같다. 가장자리에 닿으면 None(회차 확장으로). + """ + started = time.perf_counter() + big = source.spec + whole = analyze_domain( + big, + np.ones((big.n_rows, big.n_cols), dtype=bool), + contour_features, + route_line, + upstream_streams, + descent=source, + ) + if whole is None: + return None + reaches = whole.flow.reaches_road.reshape(big.n_rows, big.n_cols) + if (outermost_cells(np.ones_like(reaches)) & reaches).any(): + logger.info("배수유역: 도로로 오는 셀이 큰 격자 가장자리에 닿아 넓혀 가며 다시 셉니다.") + return None + row0 = int(round((big.y_max - spec.y_max) / spec.cell_m)) + col0 = int(round((spec.x_min - big.x_min) / spec.cell_m)) + wanted = np.zeros_like(reaches) + wanted[row0 : row0 + spec.n_rows, col0 : col0 + spec.n_cols] = domain + if reaches.any(): + wanted |= distance_transform_cdt(~reaches, metric="chessboard") <= band_cells + rows = np.flatnonzero(wanted.any(axis=1)) + cols = np.flatnonzero(wanted.any(axis=0)) + top, left = int(rows[0]), int(cols[0]) + cut = GridSpec( + x_min=big.x_min + left * big.cell_m, + y_max=big.y_max - top * big.cell_m, + cell_m=big.cell_m, + n_rows=int(rows[-1]) - top + 1, + n_cols=int(cols[-1]) - left + 1, + ) + cut_domain = wanted[top : top + cut.n_rows, left : left + cut.n_cols] + analysis = analyze_domain( + cut, cut_domain, contour_features, route_line, upstream_streams, descent=source + ) + if analysis is None: + return None + added = int(cut_domain.sum()) - int(domain.sum()) + logger.info( + "[계측] 한 번 셈 %.1fs · 큰 격자 %d×%d · 도로로 오는 셀 %d · 자른 격자 %d×%d · 셀 %d (+%d)", + time.perf_counter() - started, + big.n_rows, + big.n_cols, + int(reaches.sum()), + cut.n_rows, + cut.n_cols, + int(cut_domain.sum()), + added, + ) + return RedExpansion( + analysis=analysis, rounds=1 if added > 0 else 0, closed=True, added_cells=added + ) + + def _dilate_by(mask: np.ndarray, steps: int) -> np.ndarray: """8이웃 팽창을 `steps`회 반복한다(정사각 커널이라 반경 = steps 셀).""" rows, cols = mask.shape diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Grid.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Grid.py index 57d9f62ff..2a8542f3e 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Grid.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Grid.py @@ -143,6 +143,10 @@ def iter_linestrings(geometry: Any) -> list[LineString]: return [] +# 점 구름 메모리 보관(520-6) — 열쇠 = 피처 목록(객체째 들고 있음) · 개수 · 자르는 범위 · 2 벌. +_CLOUD_CACHE: list[tuple[Any, int, Any, ContourCloud]] = [] + + def build_contour_cloud( contour_features: list[dict[str, Any]], elevation_floor_m: float | None = None, @@ -156,11 +160,30 @@ def build_contour_cloud( `clip_bounds`(x_min, y_min, x_max, y_max)를 주면 그 밖 등고선은 읽지 않는다. 도엽 전체 등고선을 다 물고 가면 TIN 삼각망 비용만 커지고 결과는 같다. + + 하한 없는 구름을 보관해 두고 하한은 점 표고로 거른다 — 피처를 거른 것과 점 · 차례가 같다. """ + for holder, count, bounds, cached in _CLOUD_CACHE: + if holder is contour_features and count == len(contour_features) and bounds == clip_bounds: + cloud = cached + break + else: + cloud = _contour_cloud(contour_features, clip_bounds) + _CLOUD_CACHE.append((contour_features, len(contour_features), clip_bounds, cloud)) + del _CLOUD_CACHE[:-2] + if elevation_floor_m is None: + return cloud + keep = cloud.z >= elevation_floor_m + return ContourCloud(cloud.xy[keep], cloud.z[keep]) + + +def _contour_cloud( + contour_features: list[dict[str, Any]], + clip_bounds: tuple[float, float, float, float] | None, +) -> ContourCloud: xs: list[np.ndarray] = [] ys: list[np.ndarray] = [] zs: list[np.ndarray] = [] - dropped_low = 0 dropped_short = 0 dropped_outside = 0 for feature in contour_features: @@ -170,9 +193,6 @@ def build_contour_cloud( elevation = _feature_elevation(feature.get("properties") or {}) if elevation is None: continue - if elevation_floor_m is not None and elevation < elevation_floor_m: - dropped_low += 1 - continue try: parsed = shape(geometry) except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다 @@ -192,8 +212,7 @@ def build_contour_cloud( zs.append(np.full(coords.shape[0], elevation, dtype=np.float64)) if not xs: logger.warning( - "배수유역: 사용할 등고선이 없습니다(저지대 %d, 파편 %d, 범위밖 %d 제외).", - dropped_low, + "배수유역: 사용할 등고선이 없습니다(파편 %d, 범위밖 %d 제외).", dropped_short, dropped_outside, ) @@ -201,9 +220,8 @@ def build_contour_cloud( xy = np.column_stack((np.concatenate(xs), np.concatenate(ys))) z = np.concatenate(zs) logger.info( - "배수유역: 등고선 정점 %d개 (저지대 %d, 파편 %d, 범위밖 %d 제외)", + "배수유역: 등고선 정점 %d개 (파편 %d, 범위밖 %d 제외)", xy.shape[0], - dropped_low, dropped_short, dropped_outside, ) diff --git a/resources/tester/test_347_1_descent_store.py b/resources/tester/test_347_1_descent_store.py index 932f5923d..efe344135 100644 --- a/resources/tester/test_347_1_descent_store.py +++ b/resources/tester/test_347_1_descent_store.py @@ -1,7 +1,7 @@ """347-1 흐름 방향 · 도엽 좌표 저장 — 저장본을 읽어도 결과 같음 · 두 번째는 짓지 않음. - 격자 원점 = 좌표 정수 m(셀 모서리 0,0) — 노선이 바뀌어도 같은 자리 셀은 같은 자리. -- 방향장 저장본: 입력 지문(도엽 지문 · 남긴 등고 단 · 셀) 같고 범위가 덮으면 읽기만 · 최근 2벌. +- 방향장 저장본: 입력 지문(도엽 지문 · 등고 단 · 셀) 같고 범위가 덮으면 읽기만 · 최근 2벌. - 도엽 좌표 저장본: 원본 자국 · 좌표계가 같으면 다시 안 바꿈. 실측(검증_sub3 · 차이 0칸 · 108 → 10.8 s 붐빔 · 340 기준 약 6 s)은 bench_342_watershed.py. """ @@ -46,13 +46,12 @@ def test_grid_origin_is_integer_metre(): assert spec.x_min == 10.0 and spec.y_max == -3.0 -def test_input_key_follows_contours_and_floor(): +def test_input_key_follows_contours(): contours = _contours() - base = descent_input_key("도엽-1", contours, None, 1.0) - assert base == descent_input_key("도엽-1", contours, None, 1.0) - assert base != descent_input_key("도엽-2", contours, None, 1.0) # 도엽 · 좌표계 바뀜 - assert base != descent_input_key("도엽-1", contours, 200.0, 1.0) # 남긴 등고 단 바뀜 - assert base == descent_input_key("도엽-1", contours, 100.5, 1.0) # 단 집합 그대로(전부 남음) + base = descent_input_key("도엽-1", contours, 1.0) + assert base == descent_input_key("도엽-1", contours, 1.0) + assert base != descent_input_key("도엽-2", contours, 1.0) # 도엽 · 좌표계 바뀜 + assert base != descent_input_key("도엽-1", contours[:-1], 1.0) # 등고 단 바뀜 def _run(store, monkeypatch, counter): @@ -88,7 +87,7 @@ def test_store_gives_same_result_and_reads_second_time(tmp_path, monkeypatch): assert np.array_equal(analysis.descent.receiver, plain.descent.receiver) saved = list(tmp_path.glob("descent_*.npz")) assert 1 <= len(saved) <= 2 - key = descent_input_key("도엽-1", _contours(), None, 1.0) + key = descent_input_key("도엽-1", _contours(), 1.0) assert next(iter_saved_descents(tmp_path, key)).spec.n_rows >= plain.spec.n_rows diff --git a/resources/tester/test_520_6_watershed_once.py b/resources/tester/test_520_6_watershed_once.py new file mode 100644 index 000000000..5b37810e2 --- /dev/null +++ b/resources/tester/test_520_6_watershed_once.py @@ -0,0 +1,88 @@ +"""520-6 배수유역 확장 한 번 셈 · 노선 무관 방향장 저장본. + +- 한 번 셈(큰 방향장 전체 판정 → 도로로 오는 셀 + 띠 ∪ 1차 영역으로 자름) = 회차 확장과 도로로 + 오는 셀 · 도로 셀 귀속이 셀마다 같음(합성 등고선 · 능선까지 비탈). +- 방향장 열쇠 · 짓기가 노선 최저점 하한과 무관 — 하한이 달라도 저장본 하나를 읽기만. +""" + +import numpy as np +from shapely.geometry import LineString + +from B04_PreProcess import B04_PreProcess_Engine_Watershed_Expand as expand +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import trace_flow +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec + + +def _ridge_contours(): + """도로 북쪽 비탈 — y = −60 능선까지 오르고 그 밖 · 도로 남쪽은 내려감(곧은 등고선).""" + features = [] + for y in range(-250, 0, 5): + level = 290.0 - abs(y + 60.0) + features.append( + { + "geometry": {"type": "LineString", "coordinates": [[0.0, y], [300.0, y]]}, + "properties": {"elevation": level}, + } + ) + return features + + +def _case(): + spec = GridSpec(x_min=110.0, y_max=-130.0, cell_m=1.0, n_rows=40, n_cols=80) + domain = np.zeros((40, 80), bool) + domain[5:35, 5:75] = True + route = LineString([(120.0, -150.0), (180.0, -150.0)]) + return spec, domain, route + + +def _run(monkeypatch, once, floor=None, store=None): + monkeypatch.setattr(expand, "DESCENT_SPARE_M", 60.0) + if not once: + monkeypatch.setattr(expand, "_expand_once", lambda *args: None) + spec, domain, route = _case() + return expand.expand_by_red_boundary( + spec, domain, _ridge_contours(), route, [], floor, band_m=10.0, descent_store=store + ) + + +def _world(analysis): + """도로로 오는 셀 → {(x, y): 도로 셀 누가거리}.""" + spec = analysis.spec + routing = trace_flow(analysis.terrain, analysis.road) + red = analysis.flow.reaches_road & analysis.flow.analyzed + chainage = analysis.road.chainage + out = {} + for index in np.flatnonzero(red): + slot = routing.road_slot[index] + x = spec.x_min + (index % spec.n_cols) * spec.cell_m + y = spec.y_max - (index // spec.n_cols) * spec.cell_m + out[(round(x, 3), round(y, 3))] = None if slot < 0 else round(float(chainage[slot]), 3) + return out + + +def test_once_matches_round_expansion(monkeypatch): + rounds = _run(monkeypatch, once=False) + monkeypatch.undo() + once = _run(monkeypatch, once=True) + assert rounds.rounds >= 2 and rounds.closed # 회차 확장이 실제로 여러 번 넓힘 + assert once.closed and once.rounds == 1 and once.added_cells > 0 + old, new = _world(rounds.analysis), _world(once.analysis) + assert len(old) > 3_000 + assert old == new + + +def test_descent_key_and_build_ignore_floor(tmp_path, monkeypatch): + built = [] + real = expand.build_contour_descent + + def counting(*args, **kwargs): + built.append(args[2:] + tuple(kwargs.values())) + return real(*args, **kwargs) + + monkeypatch.setattr(expand, "build_contour_descent", counting) + store = (tmp_path, "도엽-1") + low = _run(monkeypatch, once=True, floor=None, store=store) + high = _run(monkeypatch, once=True, floor=230.0, store=store) + assert len(built) == 1 and built[0] == () # 하한 없이 한 번만 지음 · 두 번째는 읽기만 + assert len(list(tmp_path.glob("descent_*.npz"))) == 1 + assert _world(low.analysis) == _world(high.analysis)