diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py index b4f0c4e5c..32324686f 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py @@ -1,14 +1,18 @@ -"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 요구 443 · 486). +"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 616 · 요구 443 · 497). 0차곡 = 뚜렷한 유로가 없는 곡두의 집수지형(`노선_배치.md:31` · 교본 `3. 임도망의 배치 절차.md:156`) -→ 사업 DEM(평활 안 한 `models_plan/dtm_sheet_<방법>.npz`)을 2 m 로 솎아 D8 흐름누적 · +→ 공용 지형 창고 2 m 격자에서 D8 흐름누적 · 위 집수 넓이가 `min_m2` 이상 · `max_m2` 미만인 칸 = 0차곡(그보다 크면 1차 계류) → 세선화 · 선 뽑기(B05 `_mask_to_polylines`) · 경위도 선 + 가닥 · 길이 합. 임계값 = 지식DB 수치 없음 → 정함(브레인) 0.1 ha · 1 ha · 바꿈 = 요청 `min_m2` · `max_m2`. -파생 저장본(열쇠 = 격자 지문 · 임계 · 칸 · 쓴 도엽 누적 지문). -흐름누적 = 공용 지형 창고 도엽마다(도엽 + 여유 200 m · 창고에 있는 이웃 도엽만) 한 번 셈 → -`storage/_shared/terrain/<도엽>/zero_acc__<칸>m.npz`(열쇠 = 쓴 도엽 지형 열쇠) · -임계 바꿈 = 저장 누적에서 칸 고르기 · 선 뽑기만 · 사업 DEM 이 창고 격자가 아니면(옛 받기) 통째 셈. +흐름누적 = 창고 도엽마다(도엽 + 여유 200 m · 창고에 있는 이웃 도엽만) 한 번 셈 → +`storage/_shared/terrain/<도엽>/zero_acc__<칸>m.npz`(열쇠 = 쓴 도엽 지형 열쇠). +기본 임계 선 = 도엽마다 미리 뽑아 창고에(`zero_lines__<칸>m.npz` · 도엽 + 여유 200 m 에서 뽑아 +끝점(작은 쪽)이 그 도엽 안인 선만 · 단순화 끝) → 요청 = 도엽 선 이어 붙여 받은 범위로 자름(616-2). +다른 임계 = 받은 범위 도엽 누적 이어 붙여 칸 고르기 · 선 뽑기(616-1 · 사업 DEM 격자 안 씀). +다른 임계만 파생 저장본(열쇠 = 받은 범위 · 임계 · 칸 · 쓴 도엽 누적 지문). +무거운 셈(창고에 없는 지형 · 누적 · 선 · 다른 임계 첫 셈) = +서버 본체 밖 일꾼 하나(616-3 · `heavy_pool`). 계산 자리 = ③ 서버 단독. """ @@ -16,7 +20,14 @@ from __future__ import annotations import json +import logging import math +import os +import sys +import threading +import time +import types +from concurrent.futures import ProcessPoolExecutor from pathlib import Path from typing import Any @@ -33,13 +44,17 @@ from common_util.common_util_derived_store import ( from common_util.common_util_json import atomic_write_json from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD +logger = logging.getLogger(__name__) + ZERO_ORDER_KIND = "plan_zero_order" -ZERO_ORDER_VERSION = 2 +ZERO_ORDER_VERSION = 3 # 창고 격자 · 도엽 선(616) +ACC_VERSION = 2 # 창고 누적 열쇠 판 — 올리면 전국 누적을 다시 셈 +LINES_VERSION = 1 # 창고 도엽 선 열쇠 판 ZERO_ORDER_MIN_M2 = 1_000.0 # 집수 시작(0.1 ha) ZERO_ORDER_MAX_M2 = 10_000.0 # 이 넓이부터 1차 계류(1 ha) -ZERO_ORDER_CELL_M = 2.0 # 흐름누적 칸(1 m 격자는 솎음) +ZERO_ORDER_CELL_M = 2.0 # 흐름누적 칸 ZERO_ORDER_MIN_LEN_M = 10.0 # 세선화 갈림목 토막 거름 -ZERO_ORDER_MARGIN_M = 200.0 # 도엽 누적 여유(통째 셈과 0차곡 칸 차 0 · 분석 486 3-1) +ZERO_ORDER_MARGIN_M = 200.0 # 도엽 누적 · 선 여유(통째 셈과 0차곡 칸 차 0 · 분석 486 3-1) ACC_CAP = 65535 # 저장 누적 상한(칸 수 · 2 m 면 26 ha · 임계 1 ha 보다 넉넉) SIMPLIFY_M = 1.0 @@ -159,6 +174,13 @@ def zero_order_lines( return mask_lines(d8_accumulation(z, valid), z, valid, x, y, cell_m, min_m2, max_m2) +class NeedsCompute(Exception): + """`compute=False` 인데 창고에 없거나 낡아 셈이 필요함 — 부르는 쪽이 일꾼으로 넘김.""" + + +Memo = dict[tuple[str, int, float], Any] + + def _read_json(path: Path) -> Any: try: return json.loads(path.read_text(encoding="utf-8")) @@ -166,29 +188,62 @@ def _read_json(path: Path) -> Any: return None -def _acc_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: +def _store_names(prefix: str, sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: folder = terrain_store.TERRAIN_ROOT / str(sheet_no) - stem = f"zero_acc_{int(epsg)}_{cell_m:g}m" + stem = f"{prefix}_{int(epsg)}_{cell_m:g}m" return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock" -def _dem_key(sheet_no: str, epsg: int, cell_m: float) -> dict[str, Any] | None: - """창고 지형 열쇠(셈 시간 뺌) — 안 셌거나 원천이 바뀌었으면 None.""" +def _acc_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: + return _store_names("zero_acc", sheet_no, epsg, cell_m) + + +def _lines_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]: + return _store_names("zero_lines", sheet_no, epsg, cell_m) + + +def _dem_key(sheet_no: str, epsg: int, cell_m: float, memo: Memo | None = None) -> dict | None: + """창고 지형 열쇠(셈 시간 뺌) — 안 셌거나 원천이 바뀌었으면 None · `memo` = 한 요청 안.""" + slot = (str(sheet_no), int(epsg), float(cell_m)) + if memo is not None and slot in memo: + return memo[slot] stored = terrain_store._stored_ok(sheet_no, epsg, cell_m) - return None if stored is None else {k: v for k, v in stored.items() if k != "seconds"} + key = None if stored is None else {k: v for k, v in stored.items() if k != "seconds"} + if memo is not None: + memo[slot] = key + return key -def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL_M) -> Path | None: - """도엽 하나의 흐름누적(창고) — 도엽 + 여유 200 m 지형(창고에 있는 도엽만 · 없는 이웃 = - 출구)으로 셈 · 도엽 칸만 저장(uint16 · `ACC_CAP` 로 자름) · 열쇠 = 쓴 도엽 지형 열쇠 · - 지형이 없으면 None.""" - own = _dem_key(sheet_no, epsg, cell_m) - dem_npz = terrain_store._names(sheet_no, epsg, cell_m)[0] - if own is None or own.get("empty") or not dem_npz.is_file(): - return None - with np.load(dem_npz) as data: +def _has_terrain( + sheet_no: str, epsg: int, cell_m: float, compute: bool = True, memo: Memo | None = None +) -> bool: + """창고 지형이 있음(없으면 셈) — 빈 도엽 · 받기 실패 = False.""" + dem = _dem_key(sheet_no, epsg, cell_m, memo) + if dem is None: + if not compute: + raise NeedsCompute(sheet_no) + try: + terrain_store.ensure_sheet_terrain(sheet_no, epsg, cell_m) + except RuntimeError as error: + logger.warning("0차곡 지형 없음 %s: %s", sheet_no, error) + return False + if memo is not None: + memo.pop((str(sheet_no), int(epsg), float(cell_m)), None) + dem = _dem_key(sheet_no, epsg, cell_m, memo) + return dem is not None and not dem.get("empty") + + +def _sheet_frame(path: Path, cell_m: float) -> tuple[int, int, int, int, tuple[float, ...]]: + """창고 도엽 파일(지형 · 누적) → (ix0, iy0, 행, 열, 여유 200 m 상자) — 큰 칸은 머리만 읽음.""" + import zipfile + + with np.load(path) as data: ix0, iy0 = int(data["ix0"]), int(data["iy0"]) - rows, cols = data["z"].shape + name = "z.npy" if "z" in data.files else "acc.npy" + with zipfile.ZipFile(path) as archive, archive.open(name) as head: + version = np.lib.format.read_magic(head) + read = np.lib.format.read_array_header_1_0 + rows, cols = (read if version == (1, 0) else np.lib.format.read_array_header_2_0)(head)[0] pad = math.ceil(ZERO_ORDER_MARGIN_M / cell_m) box = ( (ix0 - pad) * cell_m, @@ -196,15 +251,36 @@ def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL (ix0 + cols + pad) * cell_m, (iy0 + rows + pad) * cell_m, ) + return ix0, iy0, rows, cols, box + + +def sheet_accumulation( + sheet_no: str, + epsg: int, + cell_m: float = ZERO_ORDER_CELL_M, + compute: bool = True, + memo: Memo | None = None, +) -> Path | None: + """도엽 하나의 흐름누적(창고) — 도엽 + 여유 200 m 지형(창고에 있는 도엽만 · 없는 이웃 = + 출구)으로 셈 · 도엽 칸만 저장(uint16 · `ACC_CAP` 로 자름) · 열쇠 = 쓴 도엽 지형 열쇠 · + 지형이 없으면 None.""" + own = _dem_key(sheet_no, epsg, cell_m, memo) + dem_npz = terrain_store._names(sheet_no, epsg, cell_m)[0] + if own is None or own.get("empty") or not dem_npz.is_file(): + return None + ix0, iy0, rows, cols, box = _sheet_frame(dem_npz, cell_m) + pad = math.ceil(ZERO_ORDER_MARGIN_M / cell_m) used = {} for no in terrain_store.sheets_for_box(box, epsg): - dem = _dem_key(no, epsg, cell_m) + dem = _dem_key(no, epsg, cell_m, memo) if dem is not None and not dem.get("empty"): used[no] = dem - key = {"version": ZERO_ORDER_VERSION, "margin_m": ZERO_ORDER_MARGIN_M, "dem": used} + key = {"version": ACC_VERSION, "margin_m": ZERO_ORDER_MARGIN_M, "dem": used} npz, key_path, lock_path = _acc_names(sheet_no, epsg, cell_m) if npz.is_file() and _read_json(key_path) == key: return npz + if not compute: + raise NeedsCompute(sheet_no) with terrain_store._Lock(lock_path): if npz.is_file() and _read_json(key_path) == key: # 기다리는 동안 남이 셌을 수 있음 return npz @@ -222,30 +298,10 @@ def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL return npz -def store_accumulation_paths( - x: np.ndarray, y: np.ndarray, epsg: int, cell_m: float -) -> list[Path] | None: - """사업 DEM 격자(창고 격자 위)를 덮는 도엽 누적 — 창고 격자 아님 · 지형 안 센 도엽 = None.""" - fx, fy = x[0] / cell_m - 0.5, y[0] / cell_m - 0.5 - if cell_m != ZERO_ORDER_CELL_M or abs(fx - round(fx)) > 1e-6 or abs(fy - round(fy)) > 1e-6: - return None - half = cell_m / 2 - box = (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half) - paths: list[Path] = [] - for no in terrain_store.sheets_for_box(box, epsg): - dem = _dem_key(no, epsg, cell_m) - if dem is None: - return None - path = None if dem.get("empty") else sheet_accumulation(no, epsg, cell_m) - if path is not None: - paths.append(path) - return paths - - def paste_accumulation( paths: list[Path], x: np.ndarray, y: np.ndarray, cell_m: float ) -> np.ndarray: - """도엽 누적을 사업 격자(y 오름차순)에 이어 붙임 — 없는 칸 0.""" + """도엽 누적을 창고 격자(y 오름차순)에 이어 붙임 — 없는 칸 0 · 겹친 칸은 뒤 파일.""" ix0, iy0 = round(x[0] / cell_m - 0.5), round(y[0] / cell_m - 0.5) acc = np.zeros((len(y), len(x))) for path in paths: @@ -259,69 +315,251 @@ def paste_accumulation( return acc +def _box_lines( + acc_paths: list[Path], + epsg: int, + box: tuple[float, ...], + cell_m: float, + min_m2: float, + max_m2: float, +) -> list[Any]: + """누적 이어 붙여 상자 안 창고 격자에서 선 뽑기 → 단순화한 shapely 선(작업 좌표계 m).""" + import shapely + + sheets = [p.parent.name for p in acc_paths] + grid = terrain_store.load_terrain(epsg, sheets=sheets, box=box, cell_m=cell_m) + x, y = grid["x"], grid["y"] + if not len(x) or not len(y): + return [] + acc = paste_accumulation(acc_paths, x, y, cell_m) + lines = mask_lines(acc, grid["z"], grid["valid_mask"], x, y, cell_m, min_m2, max_m2) + if not lines: + return [] + index = np.repeat(np.arange(len(lines)), [len(line) for line in lines]) + geoms = shapely.linestrings(np.concatenate(lines), indices=index) + return list(shapely.simplify(geoms, SIMPLIFY_M)) + + +def _owned_lines(lines: list[Any], sheet_no: str, epsg: int) -> list[Any]: + """끝점 둘 중 작은 쪽(x, y 차례)이 그 도엽(경위도) 안인 선 — 이웃과 겹쳐 뽑힌 선은 한 곳만.""" + import shapely + from pyproj import Transformer + + from B04_PreProcess.B04_PreProcess_Engine_MapSheet import sheet5k_to_bounds + + if not lines: + return [] + first = shapely.get_coordinates(shapely.get_point(lines, 0)) + last = shapely.get_coordinates(shapely.get_point(lines, -1)) + swap = (last[:, 0] < first[:, 0]) | ((last[:, 0] == first[:, 0]) & (last[:, 1] < first[:, 1])) + ends = np.where(swap[:, None], last, first) + lon, lat = Transformer.from_crs( + terrain_store._crs(epsg), "EPSG:4326", always_xy=True + ).transform(ends[:, 0], ends[:, 1]) + lon0, lat0, lon1, lat1 = sheet5k_to_bounds(sheet_no) + keep = (lon >= lon0) & (lon < lon1) & (lat >= lat0) & (lat < lat1) + return [g for g, k in zip(lines, keep) if k] + + +def sheet_lines( + sheet_no: str, + epsg: int, + cell_m: float = ZERO_ORDER_CELL_M, + compute: bool = True, + memo: Memo | None = None, +) -> Path | None: + """도엽 하나의 기본 임계 0차곡 선(창고 · 616-2) — 도엽 + 여유 200 m 에서 뽑아(창고에 있는 이웃 + 누적만) 끝점이 그 도엽 안인 선만 · 단순화 끝 · 열쇠 = 쓴 누적 지문 · 지형 없으면 None.""" + import shapely + + own = sheet_accumulation(sheet_no, epsg, cell_m, compute, memo) + if own is None: + return None + box = _sheet_frame(own, cell_m)[4] + used = {} + for no in terrain_store.sheets_for_box(box, epsg): + path = own if no == sheet_no else _acc_names(no, epsg, cell_m)[0] + if (stamp := file_stamp(path)) is not None: + used[no] = stamp + key = { + "version": LINES_VERSION, + "min_m2": ZERO_ORDER_MIN_M2, + "max_m2": ZERO_ORDER_MAX_M2, + "min_len_m": ZERO_ORDER_MIN_LEN_M, + "simplify_m": SIMPLIFY_M, + "margin_m": ZERO_ORDER_MARGIN_M, + "acc": used, + } + npz, key_path, lock_path = _lines_names(sheet_no, epsg, cell_m) + if npz.is_file() and _read_json(key_path) == key: + return npz + if not compute: + raise NeedsCompute(sheet_no) + with terrain_store._Lock(lock_path): + if npz.is_file() and _read_json(key_path) == key: # 기다리는 동안 남이 셌을 수 있음 + return npz + paths = [_acc_names(no, epsg, cell_m)[0] for no in sorted(used)] + lines = _box_lines(paths, epsg, box, cell_m, ZERO_ORDER_MIN_M2, ZERO_ORDER_MAX_M2) + mine = _owned_lines(lines, sheet_no, epsg) + coords, index = shapely.get_coordinates(mine, return_index=True) + key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게 + atomic_write_npz(npz, coords=coords, index=index.astype(np.int32)) + atomic_write_json(key_path, key) + return npz + + +def _read_lines(path: Path) -> list[Any]: + import shapely + + with np.load(path) as data: + coords, index = data["coords"], data["index"] + return list(shapely.linestrings(coords, indices=index)) if len(coords) else [] + + +def _project_box(npz: Path) -> tuple[float, float, float, float]: + """사업 DEM 범위(칸 모서리) — 창고 격자에서 이 상자 안 칸을 씀.""" + with np.load(npz) as dem: + x, y = np.asarray(dem["x"], float), np.asarray(dem["y"], float) + cell = float(dem["resolution"][0]) if "resolution" in dem.files else float(x[1] - x[0]) + half = cell / 2 + return (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half) + + +def _collection(shapes: list[Any], epsg: int) -> dict[str, Any]: + """작업 좌표계 선 → 경위도 GeoJSON(좌표 변환 한 번).""" + import shapely + from pyproj import Transformer + + features: list[dict[str, Any]] = [] + if shapes: + coords, index = shapely.get_coordinates(shapes, return_index=True) + lon, lat = Transformer.from_crs( + terrain_store._crs(epsg), "EPSG:4326", always_xy=True + ).transform(coords[:, 0], coords[:, 1]) + points = np.column_stack([lon, lat]) + for part in np.split(points, np.flatnonzero(np.diff(index)) + 1): + geometry = {"type": "LineString", "coordinates": part.tolist()} + features.append({"type": "Feature", "geometry": geometry, "properties": {}}) + return {"type": "FeatureCollection", "features": features} + + def zero_order_layer( project_root: Path, models_dir: Path, epsg: int, min_m2: float = ZERO_ORDER_MIN_M2, max_m2: float = ZERO_ORDER_MAX_M2, -) -> dict[str, Any]: - """0차곡 겹 — 지문(격자 파일 · 임계 · 칸 · 도엽 누적) 같으면 저장본 · 경위도 선 + 가닥 · - 길이 합 m.""" - from shapely.geometry import LineString, mapping - from shapely.ops import transform - - from D03_PlanAnalysis.D03_PlanAnalysis_Layers import _transforms - + compute: bool = True, +) -> dict[str, Any] | None: + """0차곡 겹 — 받은 범위의 창고 도엽 선(기본 임계) 또는 누적(다른 임계 · 지문 같으면 저장본) · + 경위도 선 + 가닥 · 길이 합 m. `compute=False` = 셈이 필요하면 None(일꾼으로 넘김).""" npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz" - stamp = file_stamp(npz) - if stamp is None: + if file_stamp(npz) is None: raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.") - with np.load(npz) as dem: - x, y = np.asarray(dem["x"], float), np.asarray(dem["y"], float) - cell = float(dem["resolution"][0]) if "resolution" in dem.files else float(x[1] - x[0]) - step = max(1, round(ZERO_ORDER_CELL_M / cell)) - x, y, cell_m = x[::step], y[::step], cell * step - paths = store_accumulation_paths(x, y, epsg, cell_m) - key = derived_key( - ZERO_ORDER_KIND, - ZERO_ORDER_VERSION, - npz.name, - stamp, - int(epsg), - float(min_m2), - float(max_m2), - ZERO_ORDER_CELL_M, - ZERO_ORDER_MIN_LEN_M, - SIMPLIFY_M, - None if paths is None else [f"{p.parent.name}:{file_stamp(p)}" for p in paths], - ) - saved = read_derived(project_root, ZERO_ORDER_KIND, key) - if saved is not None: - return saved["value"] - with np.load(npz) as dem: - z = np.asarray(dem["z"], float)[::step, ::step] - valid = np.asarray(dem["valid_mask"], bool)[::step, ::step] - if paths is None: - lines = zero_order_lines(z, valid, x, y, cell_m, min_m2, max_m2) - else: - acc = paste_accumulation(paths, x, y, cell_m) - lines = mask_lines(acc, z, valid, x, y, cell_m, min_m2, max_m2) - _, to_ll = _transforms(epsg) - shapes = [LineString(line).simplify(SIMPLIFY_M) for line in lines] - value = { - "min_m2": float(min_m2), - "max_m2": float(max_m2), - "cell_m": cell_m, - "count": len(shapes), - "length_m": round(sum(s.length for s in shapes), 1), - "collection": { - "type": "FeatureCollection", - "features": [ - {"type": "Feature", "geometry": mapping(transform(to_ll, s)), "properties": {}} - for s in shapes - ], - }, - } + box = _project_box(npz) + cell_m = ZERO_ORDER_CELL_M + default = (float(min_m2), float(max_m2)) == (ZERO_ORDER_MIN_M2, ZERO_ORDER_MAX_M2) + memo: Memo = {} + try: + paths: dict[str, Path] = {} + for no in terrain_store.sheets_for_box(box, epsg): + if _has_terrain(no, epsg, cell_m, compute, memo): + if (got := sheet_accumulation(no, epsg, cell_m, compute, memo)) is not None: + paths[no] = got + if default: # 누적을 모두 갖춘 뒤 선(이웃 누적이 선 열쇠에 듦) + for no in list(paths): + paths[no] = sheet_lines( + no, epsg, cell_m, compute, memo + ) # 창고 도엽 선 = 저장본 · 이어 붙여 자르기만(파생 저장본 안 씀) + return _layer_value(_clipped_lines(list(paths.values()), box), epsg, min_m2, max_m2) + key = derived_key( + ZERO_ORDER_KIND, + ZERO_ORDER_VERSION, + [round(v, 3) for v in box], + int(epsg), + float(min_m2), + float(max_m2), + cell_m, + ZERO_ORDER_MIN_LEN_M, + SIMPLIFY_M, + [f"{p.name}:{no}:{file_stamp(p)}" for no, p in sorted(paths.items())], + ) + saved = read_derived(project_root, ZERO_ORDER_KIND, key) + if saved is not None: + return saved["value"] + if not compute: + raise NeedsCompute("임계") + except NeedsCompute: + return None + shapes = _box_lines([paths[no] for no in sorted(paths)], epsg, box, cell_m, min_m2, max_m2) + value = _layer_value(shapes, epsg, min_m2, max_m2) write_derived(project_root, ZERO_ORDER_KIND, key, value) return value + + +def _clipped_lines(paths: list[Path], box: tuple[float, ...]) -> list[Any]: + """도엽 선 이어 붙여 받은 범위로 자름 — 잘려 짧아진 토막(`ZERO_ORDER_MIN_LEN_M` 미만) 거름.""" + import shapely + + lines = [g for p in paths for g in _read_lines(p)] + if not lines: + return [] + parts = shapely.get_parts(shapely.clip_by_rect(lines, *box)) + keep = (shapely.get_type_id(parts) == 1) & (shapely.length(parts) >= ZERO_ORDER_MIN_LEN_M) + return list(parts[keep]) + + +def _layer_value(shapes: list[Any], epsg: int, min_m2: float, max_m2: float) -> dict[str, Any]: + import shapely + + return { + "min_m2": float(min_m2), + "max_m2": float(max_m2), + "cell_m": ZERO_ORDER_CELL_M, + "count": len(shapes), + "length_m": round(float(shapely.length(shapes).sum()), 1) if shapes else 0.0, + "collection": _collection(shapes, epsg), + } + + +# ── 무거운 셈 일꾼(616-3) ───────────────────────────────────────────── + +_POOL: ProcessPoolExecutor | None = None + + +def _watch_parent(pid: int) -> None: + """일꾼 시작 — 서버 프로세스가 끝나면(다시 켬 · 강제 종료) 따라 끝남(고아 일꾼 막음).""" + import psutil + + parent = psutil.Process(pid) + + def loop() -> None: + while parent.is_running(): + time.sleep(2.0) + os._exit(0) + + threading.Thread(target=loop, daemon=True).start() + + +def _warm() -> None: + """일꾼 띄울 때 무거운 모듈을 미리 읽음(첫 셈이 읽기를 안 기다리게).""" + import shapely # noqa: F401 + import skimage.morphology # noqa: F401 + + from B05_Profile import B05_Profile_Engine_Skeleton # noqa: F401 + + +def heavy_pool() -> ProcessPoolExecutor: + """서버 본체 밖 일꾼 하나 — 처음 부를 때 띄워 둠 · 죽었으면 새로(이벤트 루프 실에서만 부름).""" + global _POOL + if _POOL is None or getattr(_POOL, "_broken", False): + _POOL = ProcessPoolExecutor(1, initializer=_watch_parent, initargs=(os.getpid(),)) + # ponytail: spawn 일꾼이 main.py(라우터 전부 · 약 10 s)를 다시 읽지 않게 띄우는 순간만 + # 빈 __main__ 으로 바꿈 — 일꾼은 이 모듈만 읽음 · 파이썬 판이 바뀌어 안 먹히면 느려질 뿐 + main = sys.modules["__main__"] + sys.modules["__main__"] = types.ModuleType("__main__") + try: + _POOL.submit(_warm) + finally: + sys.modules["__main__"] = main + return _POOL diff --git a/resources/tester/test_584_1_plan_geo_lines.py b/resources/tester/test_584_1_plan_geo_lines.py index 529380914..bb9ae2cf1 100644 --- a/resources/tester/test_584_1_plan_geo_lines.py +++ b/resources/tester/test_584_1_plan_geo_lines.py @@ -182,31 +182,6 @@ def test_zero_order_lines_head_only() -> None: assert sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= zo.ZERO_ORDER_MIN_LEN_M -def test_zero_order_layer_saved_and_threshold_key(tmp_path, monkeypatch) -> None: - n = 60 - models = tmp_path / "models" - models.mkdir() - np.savez( - models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz", - x=230000.0 + np.arange(n) * 2.0, - y=450000.0 + np.arange(n) * 2.0, - z=_valley(n), - valid_mask=np.ones((n, n), bool), - resolution=np.array([2.0, 2.0]), - ) - calls = [] - real = zo.zero_order_lines - monkeypatch.setattr(zo, "zero_order_lines", lambda *a, **k: calls.append(1) or real(*a, **k)) - first = zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 2000.0) - assert first["count"] >= 1 and first["length_m"] > 0 and first["cell_m"] == 2.0 - assert first["collection"]["features"][0]["geometry"]["type"] == "LineString" - again = zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 2000.0) - assert (again["count"], again["length_m"]) == (first["count"], first["length_m"]) - assert len(calls) == 1 # 두 번째 = 저장본 - zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 3000.0) # 임계 바뀜 = 다시 셈 - assert len(calls) == 2 - - def test_zero_order_no_dem(tmp_path) -> None: with pytest.raises(FileNotFoundError): zo.zero_order_layer(tmp_path, tmp_path, EPSG) diff --git a/resources/tester/test_608_plan_zero_order_store.py b/resources/tester/test_608_plan_zero_order_store.py index c399e4b0a..6f242bced 100644 --- a/resources/tester/test_608_plan_zero_order_store.py +++ b/resources/tester/test_608_plan_zero_order_store.py @@ -1,5 +1,5 @@ """D 0차곡 흐름누적 = 도엽별 공용 창고(PLAN 608-2) — 도엽 + 여유 200 m 조각 = 통째 셈 · -처음만 셈 · 이웃 도엽 지형이 새로 생기면 다시 · 임계 바꿈은 누적 그대로 · 창고 격자 아닌 DEM = 통째. +처음만 셈 · 이웃 도엽 지형이 새로 생기면 다시 · 임계 바꿈은 누적 그대로(사업 DEM 격자 무관 = 616). ./venv/Scripts/python.exe -m pytest resources/tester/test_608_plan_zero_order_store.py -q """ @@ -77,7 +77,7 @@ def test_tiles_equal_whole_and_saved(shared) -> None: def test_layer_uses_store_and_threshold_reuses(shared, tmp_path, monkeypatch) -> None: - """사업 DEM 이 창고 격자 위면 도엽 누적을 이어 붙임 · 임계 바꿈 = 누적 다시 안 셈.""" + """받은 범위 도엽 누적을 이어 붙임 · 임계 바꿈 = 누적 다시 안 셈.""" calls = shared monkeypatch.setattr(zo, "ZERO_ORDER_CELL_M", CELL) _put_sheet(WEST) @@ -100,5 +100,3 @@ def test_layer_uses_store_and_threshold_reuses(shared, tmp_path, monkeypatch) -> assert len(calls) == 2 # 도엽 둘 한 번씩 zo.zero_order_layer(tmp_path, models, EPSG, 1 * cells, 10 * cells) assert len(calls) == 2 # 임계 바꿈 = 저장 누적에서 칸 고르기만 - # 창고 격자에서 반 칸 어긋난 DEM = 통째 셈 - assert zo.store_accumulation_paths(x + CELL / 2, y, EPSG, CELL) is None diff --git a/resources/tester/test_616_zero_order_fast.py b/resources/tester/test_616_zero_order_fast.py new file mode 100644 index 000000000..3ecf7f68c --- /dev/null +++ b/resources/tester/test_616_zero_order_fast.py @@ -0,0 +1,189 @@ +"""PLAN 616 — 0차곡 첫 계산 빠르게: 창고 격자에서 셈(616-1) · 도엽 선 미리 뽑아 이어 붙임(616-2) · +무거운 셈 = 서버 밖 일꾼(616-3) · 급경사 값 그림 바이너리(616-4). + +./venv/Scripts/python.exe -m pytest resources/tester/test_616_zero_order_fast.py -q +""" + +from __future__ import annotations + +import base64 +import os +import sys +from pathlib import Path + +import numpy as np +import pytest +import shapely + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) + +from common_util import common_util_terrain_store as store # noqa: E402 +from common_util.common_util_atomic import atomic_write_npz # noqa: E402 +from common_util.common_util_json import atomic_write_json # noqa: E402 +from D03_PlanAnalysis import D03_PlanAnalysis_Layers as layers # noqa: E402 +from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zo # noqa: E402 +from tools import precompute_terrain as tool # noqa: E402 + +EPSG = 5187 +WEST, EAST = "36901074", "36901075" +CELL = 10.0 # 시험은 굵은 칸(여유 200 m = 20 칸 · 기본 임계 = 10 ~ 100 칸) + + +def _height(x, y): + """기운 면 + 골 · 마루 — 웅덩이 없음 · 작은 집수.""" + return 0.08 * x + 6.0 * np.sin(y / 38.0) + 2.0 * np.sin(x / 55.0) + + +def _put_sheet(no: str) -> None: + ix0, iy0, own = store._owned(no, EPSG, CELL) + rows, cols = own.shape + gx, gy = np.meshgrid((np.arange(cols) + ix0 + 0.5) * CELL, (np.arange(rows) + iy0 + 0.5) * CELL) + z = _height(gx - 300000.0, gy).astype(np.float32) + z[~own] = np.nan + npz, key_path, _ = store._names(no, EPSG, CELL) + atomic_write_npz(npz, z=z, ix0=np.int64(ix0), iy0=np.int64(iy0), cell_m=np.float64(CELL)) + atomic_write_json(key_path, {**store._key(EPSG, CELL, {}), "source": "gpkg", "seconds": 1}) + + +@pytest.fixture +def shared(tmp_path, monkeypatch): + """임시 창고 · 칸 10 m · d8 부름 셈.""" + monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path / "terrain") + monkeypatch.setattr(zo, "ZERO_ORDER_CELL_M", CELL) + calls = [] + real = zo.d8_accumulation + monkeypatch.setattr(zo, "d8_accumulation", lambda *a: calls.append(1) or real(*a)) + return calls + + +def _models(tmp_path: Path, shift: float = 0.0) -> Path: + """두 도엽 가운데를 덮는 사업 DEM(칸 = 시험 칸 · 값은 안 씀 · 범위만).""" + grid = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL) + x, y = grid["x"][30:-30] + shift, grid["y"][30:-30] + models = tmp_path / f"models{shift:g}" + models.mkdir() + np.savez( + models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz", + x=x, + y=y, + z=np.zeros((len(y), len(x))), + valid_mask=np.ones((len(y), len(x)), bool), + resolution=np.array([CELL], np.float32), + ) + return models + + +def test_sheet_lines_join_equals_box(shared, tmp_path) -> None: + """도엽 선(여유 200 m 에서 뽑아 끝점 도엽만) 이어 붙임 = 두 도엽 통째 상자에서 뽑은 선.""" + _put_sheet(WEST) + _put_sheet(EAST) + paths = [zo.sheet_lines(no, EPSG, CELL) for no in (WEST, EAST)] + joined = [g for p in paths for g in zo._read_lines(p)] + grid = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL) + half = CELL / 2 + box = (grid["x"][0] - half, grid["y"][0] - half, grid["x"][-1] + half, grid["y"][-1] + half) + accs = [zo._acc_names(no, EPSG, CELL)[0] for no in (WEST, EAST)] + whole = zo._box_lines(accs, EPSG, box, CELL, zo.ZERO_ORDER_MIN_M2, zo.ZERO_ORDER_MAX_M2) + assert len(whole) >= 20 + # 여유 200 m 보다 긴 선(시험 지형 = 동서로 긴 골)만 창 끝에서 잘려 다를 수 있음 — 1 % 안 + assert abs(len(joined) - len(whole)) <= 0.01 * len(whole) + total = sum(g.length for g in whole) + assert abs(sum(g.length for g in joined) - total) <= 0.01 * total + keys = {tuple(map(tuple, shapely.get_coordinates(g).round(3))) for g in joined} + assert len(keys) == len(joined) # 겹쳐 뽑힌 선은 한 도엽에만 + + +def test_layer_default_reads_lines_light(shared, tmp_path) -> None: + """기본 임계 = 창고 도엽 선 이어 붙여 자름 — 창고에 없으면 compute=False 가 None(일꾼 몫) · + 있으면 d8 · 선 뽑기 없이 바로.""" + calls = shared + _put_sheet(WEST) + _put_sheet(EAST) + models = _models(tmp_path) + assert zo.zero_order_layer(tmp_path, models, EPSG, compute=False) is None + assert calls == [] + first = zo.zero_order_layer(tmp_path, models, EPSG) + assert first["count"] >= 10 and first["cell_m"] == CELL + assert len(calls) == 2 # 도엽 둘 누적 한 번씩 + light = zo.zero_order_layer(tmp_path, models, EPSG, compute=False) + assert (light["count"], light["length_m"]) == (first["count"], first["length_m"]) + assert len(calls) == 2 + assert light["collection"]["features"][0]["geometry"]["type"] == "LineString" + # 받은 범위 밖 선은 잘림 — 모든 점이 범위 안 + with np.load(models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz") as dem: + x0, x1 = dem["x"][0] - CELL / 2, dem["x"][-1] + CELL / 2 + from pyproj import Transformer + + to_m = Transformer.from_crs("EPSG:4326", f"EPSG:{EPSG}", always_xy=True) + pts = np.array( + [p for f in light["collection"]["features"] for p in f["geometry"]["coordinates"]] + ) + mx, _ = to_m.transform(pts[:, 0], pts[:, 1]) + assert mx.min() >= x0 - 0.01 and mx.max() <= x1 + 0.01 + + +def test_layer_custom_threshold_saved_and_grid_free(shared, tmp_path) -> None: + """다른 임계 = 창고 누적에서 셈(사업 DEM 격자 안 씀 · 반 칸 어긋나도 창고) · 둘째 = 저장본.""" + calls = shared + _put_sheet(WEST) + _put_sheet(EAST) + models = _models(tmp_path) + cells = CELL * CELL + assert zo.zero_order_layer(tmp_path, models, EPSG, 5 * cells, 40 * cells, False) is None + first = zo.zero_order_layer(tmp_path, models, EPSG, 5 * cells, 40 * cells) + assert first["count"] >= 1 and len(calls) == 2 + again = zo.zero_order_layer(tmp_path, models, EPSG, 5 * cells, 40 * cells, False) + assert again == first # 저장본 + shifted = zo.zero_order_layer( + tmp_path, _models(tmp_path, CELL / 2), EPSG, 5 * cells, 40 * cells + ) + assert shifted["count"] >= 1 and len(calls) == 2 # 옛 격자 갈래 없음 = 누적 다시 안 셈 + + +def test_lines_recomputed_when_neighbor_acc_appears(shared) -> None: + _put_sheet(WEST) + _put_sheet(EAST) + path = zo.sheet_lines(WEST, EPSG, CELL) # 동쪽 누적 없음 = 빼고 셈 + stamp = path.stat().st_mtime_ns + assert zo.sheet_lines(WEST, EPSG, CELL, compute=False) == path + zo.sheet_accumulation(EAST, EPSG, CELL) + with pytest.raises(zo.NeedsCompute): + zo.sheet_lines(WEST, EPSG, CELL, compute=False) + assert zo.sheet_lines(WEST, EPSG, CELL).stat().st_mtime_ns != stamp + + +def test_acc_key_version_unchanged() -> None: + """누적 열쇠 판을 올리면 전국 미리 계산 누적을 통째 다시 셈 — 616 은 안 올림.""" + assert zo.ACC_VERSION == 2 + + +def test_precompute_zero_lines_step(shared, monkeypatch) -> None: + _put_sheet(WEST) + _put_sheet(EAST) + monkeypatch.setattr(tool, "CELL_M", CELL) + assert "zero_lines" in tool.STEPS and tool.STEPS.index("zero_lines") == 2 + _, _, info = tool.run_step("zero_lines", WEST) + assert info["ok"] and info["fresh"] and not info["empty"] + _, _, again = tool.run_step("zero_lines", WEST) + assert again["ok"] and not again["fresh"] + + +def test_raster_bytes_roundtrip() -> None: + png = b"\x89PNG\r\n\x1a\nfake" + value = {"image": "data:image/png;base64," + base64.b64encode(png).decode(), "counts": [1, 2]} + got, meta = layers.raster_bytes(value) + assert got == png and meta == {"counts": [1, 2]} + + +def test_heavy_pool_runs_outside() -> None: + """일꾼 = 다른 프로세스 · 띄울 때 __main__ 을 되돌려 놓음.""" + main = sys.modules["__main__"] + pool = zo.heavy_pool() + try: + assert sys.modules["__main__"] is main + assert pool.submit(os.getpid).result(timeout=120) != os.getpid() + assert zo.heavy_pool() is pool + finally: + pool.shutdown() + zo._POOL = None diff --git a/tools/precompute_terrain.py b/tools/precompute_terrain.py index b5a78e4a3..2e3b9a8b0 100644 --- a/tools/precompute_terrain.py +++ b/tools/precompute_terrain.py @@ -2,7 +2,8 @@ """전국 미리 계산 일괄 도구(PLAN 612-1 · 요구 493) — 공용 지형 창고에 도엽마다 미리 셈. 실행: ./venv/Scripts/python.exe tools/precompute_terrain.py --sheets 36901 --workers 4 - --steps dem,zero,e07 단계(차례 고정 = 지형 → 0차곡 누적 → E07 · 뒤에 tiles) + --steps dem,zero,zero_lines,e07 단계(차례 고정 = 지형 → 0차곡 누적 → 0차곡 선 → E07 · + 뒤에 tiles · zero_lines = 두 PC zero 단계가 다 끝난 뒤 · 이웃 누적 없으면 그 이웃 빼고 셈) --sheets 36901,36712089 도엽 번호 또는 앞자리(1:50,000 묶음) · 주면 높이 차 거름 안 씀 --box lon0,lat0,lon1,lat1 도엽 가운데가 이 경위도 상자 안 --min-relief 150 도엽 안 높이 차(m) 이상만(산지) · --limit N 앞 N 장만 @@ -14,8 +15,8 @@ - 도엽 목록 = 전국 등고선 파일 꼭짓점이 닿는 도엽(0.025°) — `_precompute_sheets.json`(전국 파일 지문이 같으면 다시 안 셈) · 차례 = 등고선 길이 큰 것(산지)부터 · 같으면 높이 차. - 좌표계 = 도엽 가운데 경도의 TM 띠(`tm_epsg_for_lon` · 5185~5187) · 격자 2 m. -- 셈 = 공용 함수만 부름(창고 지형 · D 0차곡 누적 · E07 도엽) — 저장 열쇠가 같으면 함수가 건너뜀 · - 쓰기 = 함수 쪽 원자 쓰기(끊겨도 반쯤 쓴 파일 없음). +- 셈 = 공용 함수만 부름(창고 지형 · D 0차곡 누적 · 선 · E07 도엽) — 저장 열쇠가 같으면 함수가 + 건너뜀 · 쓰기 = 함수 쪽 원자 쓰기(끊겨도 반쯤 쓴 파일 없음). - 이어 하기 = PC 별 `_precompute_progress_<컴퓨터 이름>.json`(단계 · 도엽마다 성공 · 시간 · 시놀로지가 두 PC 사이 옮김 → 한 파일은 한 PC 만 씀) · 옛 `_precompute_progress.json` 의 성공은 읽어서 이어 받음(안 씀) · 성공한 도엽은 다시 안 부름 · 실패는 다음에 다시 · 전국 파일이 바뀌면 @@ -59,7 +60,7 @@ from common_util import common_util_terrain_store as store # noqa: E402 from common_util.common_util_crs import tm_epsg_for_lon # noqa: E402 from common_util.common_util_json import atomic_write_json # noqa: E402 -STEPS = ("dem", "zero", "e07") +STEPS = ("dem", "zero", "zero_lines", "e07") CELL_M = 2.0 PROGRESS_NAME = "_precompute_progress.json" # 옛 한 벌 기록(읽기만) PC_NAME = platform.node() or "pc" @@ -70,7 +71,12 @@ INDEX_CHUNK = 50_000 # 가닥 묶음(약 0.5 GB 읽기) INDEX_STRIDE = 8 # 꼭짓점 솎음(도엽 2 km 에 비해 넉넉) SAVE_EVERY_S = 10.0 # 단계별 창고 파일 앞머리(크기 셈) -STEP_FILES = {"dem": ("dem_",), "zero": ("zero_acc_",), "e07": ("e07_", "e07si_")} +STEP_FILES = { + "dem": ("dem_",), + "zero": ("zero_acc_",), + "zero_lines": ("zero_lines_",), + "e07": ("e07_", "e07si_"), +} _ENVELOPE = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64} @@ -322,6 +328,13 @@ def run_step(step: str, job: str, part: Part | None = None) -> tuple[str, str, d before = _mtime(target) empty = zero.sheet_accumulation(sheet_no, epsg, CELL_M) is None fresh = _mtime(target) != before + elif step == "zero_lines": + from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zero + + target = zero._lines_names(sheet_no, epsg, CELL_M)[0] + before = _mtime(target) + empty = zero.sheet_lines(sheet_no, epsg, CELL_M) is None + fresh = _mtime(target) != before elif step == "e07": from E07_StabArea import E07_StabArea_Router as e07_router from E07_StabArea import E07_StabArea_Sheet as e07 @@ -558,7 +571,7 @@ def summary(sheets: list[str], progress: dict[str, Any]) -> dict[str, Any]: def main(argv: list[str] | None = None) -> int: parser = argparse.ArgumentParser(description="공용 지형 창고 전국 미리 계산(PLAN 612)") - parser.add_argument("--steps", default="dem,zero,e07") + parser.add_argument("--steps", default="dem,zero,zero_lines,e07") parser.add_argument("--sheets", default="") parser.add_argument("--box", default="") parser.add_argument("--min-relief", type=float, default=150.0)