316 lines
10 KiB
Python
316 lines
10 KiB
Python
"""B05 지형 스켈레톤(주/지 능선·계곡) 추출.
|
|
|
|
수문학적 정의: D8 흐름누적이 임계값 이상인 셀=계곡, DEM 반전 시 능선.
|
|
누적값 크기의 2차 임계값으로 주/지를 나눈다. whitebox 우선, 실패 시 numpy
|
|
D8 폴백. 산출 polyline은 비용면과 동일 좌표계이며 B05_wf2_Route/route에 캐시된다.
|
|
"""
|
|
|
|
import json
|
|
import math
|
|
import tempfile
|
|
import time
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
|
|
from B05_wf2_Route.B05_wf2_Route_Engine_Solver import (
|
|
_MODELS_SUBDIR,
|
|
_ROUTE_CACHE_SUBDIR,
|
|
_cost_surface_signature,
|
|
_load_or_build_cost_surface,
|
|
)
|
|
from config.config_system import (
|
|
SKELETON_MAIN_RIDGE_ACC_THRESHOLD_CELLS,
|
|
SKELETON_MAIN_VALLEY_ACC_THRESHOLD_CELLS,
|
|
SKELETON_RIDGE_ACC_THRESHOLD_CELLS,
|
|
SKELETON_VALLEY_ACC_THRESHOLD_CELLS,
|
|
)
|
|
|
|
SKELETON_CLASSES = ("main_ridge", "minor_ridge", "main_valley", "minor_valley")
|
|
|
|
_D8_OFFSETS = [
|
|
(-1, -1),
|
|
(-1, 0),
|
|
(-1, 1),
|
|
(0, -1),
|
|
(0, 1),
|
|
(1, -1),
|
|
(1, 0),
|
|
(1, 1),
|
|
]
|
|
|
|
|
|
def _default_thresholds() -> dict[str, float]:
|
|
return {
|
|
"valley_acc": float(SKELETON_VALLEY_ACC_THRESHOLD_CELLS),
|
|
"main_valley_acc": float(SKELETON_MAIN_VALLEY_ACC_THRESHOLD_CELLS),
|
|
"ridge_acc": float(SKELETON_RIDGE_ACC_THRESHOLD_CELLS),
|
|
"main_ridge_acc": float(SKELETON_MAIN_RIDGE_ACC_THRESHOLD_CELLS),
|
|
}
|
|
|
|
|
|
def d8_flow_accumulation_numpy(z_grid: np.ndarray, valid_mask: np.ndarray) -> np.ndarray:
|
|
"""numpy 기반 D8 흐름누적 (whitebox 폴백)."""
|
|
rows, cols = z_grid.shape
|
|
n = rows * cols
|
|
z_flat = z_grid.ravel()
|
|
valid_flat = valid_mask.ravel()
|
|
|
|
receiver = np.full(n, -1, dtype=np.int64)
|
|
for idx in range(n):
|
|
if not valid_flat[idx]:
|
|
continue
|
|
r, c = divmod(idx, cols)
|
|
zc = z_flat[idx]
|
|
best_slope = 0.0
|
|
best = -1
|
|
for dr, dc in _D8_OFFSETS:
|
|
nr, nc = r + dr, c + dc
|
|
if not (0 <= nr < rows and 0 <= nc < cols):
|
|
continue
|
|
nidx = nr * cols + nc
|
|
if not valid_flat[nidx]:
|
|
continue
|
|
drop = zc - z_flat[nidx]
|
|
if drop <= 0:
|
|
continue
|
|
dist = math.sqrt(2.0) if (dr != 0 and dc != 0) else 1.0
|
|
slope = drop / dist
|
|
if slope > best_slope:
|
|
best_slope = slope
|
|
best = nidx
|
|
receiver[idx] = best
|
|
|
|
acc = np.where(valid_flat, 1.0, 0.0)
|
|
order = np.argsort(-z_flat, kind="stable")
|
|
for idx in order:
|
|
recv = receiver[idx]
|
|
if recv >= 0 and valid_flat[idx]:
|
|
acc[recv] += acc[idx]
|
|
return acc.reshape(rows, cols)
|
|
|
|
|
|
def _d8_flow_accumulation_whitebox(
|
|
z_grid: np.ndarray,
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
valid_mask: np.ndarray,
|
|
grid_res: float,
|
|
) -> np.ndarray | None:
|
|
"""WhiteboxTools로 D8 흐름누적을 계산한다 (실패 시 None → numpy 폴백)."""
|
|
try:
|
|
import rasterio
|
|
from rasterio.transform import from_origin
|
|
from whitebox import WhiteboxTools
|
|
except Exception:
|
|
return None
|
|
|
|
nodata = -9999.0
|
|
rows, cols = z_grid.shape
|
|
z_out = np.where(valid_mask, z_grid, nodata).astype(np.float32)[::-1, :]
|
|
transform = from_origin(
|
|
float(x_coords[0]) - grid_res / 2.0,
|
|
float(y_coords[-1]) + grid_res / 2.0,
|
|
grid_res,
|
|
grid_res,
|
|
)
|
|
try:
|
|
with tempfile.TemporaryDirectory(prefix="wbt_skel_") as tmp:
|
|
tmp_path = Path(tmp)
|
|
dem_tif = tmp_path / "dem.tif"
|
|
acc_tif = tmp_path / "acc.tif"
|
|
with rasterio.open(
|
|
dem_tif,
|
|
"w",
|
|
driver="GTiff",
|
|
height=rows,
|
|
width=cols,
|
|
count=1,
|
|
dtype="float32",
|
|
nodata=nodata,
|
|
transform=transform,
|
|
) as dst:
|
|
dst.write(z_out, 1)
|
|
|
|
wbt = WhiteboxTools()
|
|
wbt.set_verbose_mode(False)
|
|
wbt.set_working_dir(str(tmp_path))
|
|
if wbt.fill_depressions("dem.tif", "filled.tif") != 0:
|
|
raise RuntimeError("fill_depressions 실패")
|
|
if wbt.d8_flow_accumulation("filled.tif", "acc.tif", out_type="cells") != 0:
|
|
raise RuntimeError("d8_flow_accumulation 실패")
|
|
|
|
with rasterio.open(acc_tif) as src:
|
|
acc = src.read(1).astype(np.float64)[::-1, :]
|
|
return np.where(np.isfinite(acc) & (acc > 0) & valid_mask, acc, 0.0)
|
|
except Exception:
|
|
return None
|
|
|
|
|
|
def _flow_accumulation(
|
|
z_grid: np.ndarray,
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
valid_mask: np.ndarray,
|
|
grid_res: float,
|
|
) -> np.ndarray:
|
|
acc = _d8_flow_accumulation_whitebox(z_grid, x_coords, y_coords, valid_mask, grid_res)
|
|
if acc is None:
|
|
acc = d8_flow_accumulation_numpy(z_grid, valid_mask)
|
|
return acc
|
|
|
|
|
|
def _trace_polylines(mask: np.ndarray) -> list[list[tuple[int, int]]]:
|
|
"""1픽셀 폭 스켈레톤 마스크를 (r, c) polyline 목록으로 변환한다."""
|
|
pixels = set(zip(*np.nonzero(mask)))
|
|
if not pixels:
|
|
return []
|
|
|
|
def neighbors(p):
|
|
r, c = p
|
|
return [(r + dr, c + dc) for dr, dc in _D8_OFFSETS if (r + dr, c + dc) in pixels]
|
|
|
|
degree = {p: len(neighbors(p)) for p in pixels}
|
|
seeds = [p for p in pixels if degree[p] != 2]
|
|
visited_edges: set = set()
|
|
polylines: list[list[tuple[int, int]]] = []
|
|
|
|
def edge_key(a, b):
|
|
return (a, b) if a <= b else (b, a)
|
|
|
|
def walk(start, nxt):
|
|
path = [start, nxt]
|
|
visited_edges.add(edge_key(start, nxt))
|
|
prev, curr = start, nxt
|
|
while degree[curr] == 2:
|
|
candidates = [q for q in neighbors(curr) if q != prev]
|
|
if not candidates:
|
|
break
|
|
q = candidates[0]
|
|
if edge_key(curr, q) in visited_edges:
|
|
break
|
|
visited_edges.add(edge_key(curr, q))
|
|
path.append(q)
|
|
prev, curr = curr, q
|
|
return path
|
|
|
|
for seed in seeds:
|
|
for nb in neighbors(seed):
|
|
if edge_key(seed, nb) not in visited_edges:
|
|
polylines.append(walk(seed, nb))
|
|
|
|
for p in pixels:
|
|
if degree[p] == 2:
|
|
for nb in neighbors(p):
|
|
if edge_key(p, nb) not in visited_edges:
|
|
polylines.append(walk(p, nb))
|
|
|
|
return [pl for pl in polylines if len(pl) >= 2]
|
|
|
|
|
|
def _mask_to_polylines(
|
|
mask: np.ndarray, x_coords: np.ndarray, y_coords: np.ndarray, z_grid: np.ndarray
|
|
) -> list[dict[str, Any]]:
|
|
"""셀 마스크를 세선화한 뒤 모델좌표 polyline 목록으로 변환한다."""
|
|
if not mask.any():
|
|
return []
|
|
try:
|
|
from skimage.morphology import skeletonize
|
|
|
|
skel = skeletonize(mask)
|
|
except Exception:
|
|
skel = mask
|
|
out = []
|
|
for pixel_path in _trace_polylines(skel):
|
|
poly = [
|
|
[float(x_coords[c]), float(y_coords[r]), float(z_grid[r, c])] for r, c in pixel_path
|
|
]
|
|
out.append({"polyline": poly})
|
|
return out
|
|
|
|
|
|
def extract_skeleton_from_grid(
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
z_grid: np.ndarray,
|
|
valid_mask: np.ndarray,
|
|
grid_res: float,
|
|
thresholds: dict[str, float] | None = None,
|
|
use_whitebox: bool = True,
|
|
) -> dict[str, list[dict[str, Any]]]:
|
|
"""격자에서 주/지 능선·계곡 polyline을 추출한다 (캐시 없음, 테스트용 공개 API)."""
|
|
th = thresholds or _default_thresholds()
|
|
|
|
if use_whitebox:
|
|
acc_valley = _flow_accumulation(z_grid, x_coords, y_coords, valid_mask, grid_res)
|
|
acc_ridge = _flow_accumulation(-z_grid, x_coords, y_coords, valid_mask, grid_res)
|
|
else:
|
|
acc_valley = d8_flow_accumulation_numpy(z_grid, valid_mask)
|
|
acc_ridge = d8_flow_accumulation_numpy(-z_grid, valid_mask)
|
|
|
|
valley_mask = valid_mask & (acc_valley >= th["valley_acc"])
|
|
main_valley_mask = valley_mask & (acc_valley >= th["main_valley_acc"])
|
|
minor_valley_mask = valley_mask & ~main_valley_mask
|
|
|
|
ridge_mask = valid_mask & ~valley_mask & (acc_ridge >= th["ridge_acc"])
|
|
main_ridge_mask = ridge_mask & (acc_ridge >= th["main_ridge_acc"])
|
|
minor_ridge_mask = ridge_mask & ~main_ridge_mask
|
|
|
|
return {
|
|
"main_ridge": _mask_to_polylines(main_ridge_mask, x_coords, y_coords, z_grid),
|
|
"minor_ridge": _mask_to_polylines(minor_ridge_mask, x_coords, y_coords, z_grid),
|
|
"main_valley": _mask_to_polylines(main_valley_mask, x_coords, y_coords, z_grid),
|
|
"minor_valley": _mask_to_polylines(minor_valley_mask, x_coords, y_coords, z_grid),
|
|
}
|
|
|
|
|
|
def _skeleton_signature(models_dir: Path, filter_key: str, method: str, smooth: bool) -> str:
|
|
"""소스 모델 + 격자 해상도 + 분류 임계값 서명."""
|
|
th = _default_thresholds()
|
|
th_part = "|".join(f"{k}={v}" for k, v in sorted(th.items()))
|
|
return _cost_surface_signature(models_dir, filter_key, method, smooth) + "|" + th_part
|
|
|
|
|
|
def load_or_build_skeleton(
|
|
project_root: Path, filter_key: str, method: str, smooth: bool
|
|
) -> dict[str, Any]:
|
|
"""스켈레톤을 캐시에서 로드하거나 새로 계산해 저장한다."""
|
|
project_root = Path(project_root)
|
|
models_dir = project_root / _MODELS_SUBDIR
|
|
cache_dir = project_root / _ROUTE_CACHE_SUBDIR
|
|
suffix = "_smooth" if smooth else ""
|
|
cache_path = cache_dir / f"terrain_skeleton_{filter_key}_{method}{suffix}.json"
|
|
signature = _skeleton_signature(models_dir, filter_key, method, smooth)
|
|
|
|
if cache_path.exists():
|
|
try:
|
|
with open(cache_path, encoding="utf-8") as f:
|
|
cached = json.load(f)
|
|
if cached.get("signature") == signature:
|
|
return cached
|
|
except Exception:
|
|
pass
|
|
|
|
_t0 = time.time()
|
|
(x_coords, y_coords, z_grid, valid_mask, _dz_dx, _dz_dy, grid_res) = (
|
|
_load_or_build_cost_surface(project_root, models_dir, filter_key, method, smooth)
|
|
)
|
|
z_grid = np.asarray(z_grid, dtype=np.float64)
|
|
valid_mask = np.asarray(valid_mask, dtype=bool)
|
|
|
|
skeleton = extract_skeleton_from_grid(
|
|
np.asarray(x_coords), np.asarray(y_coords), z_grid, valid_mask, float(grid_res)
|
|
)
|
|
result: dict[str, Any] = dict(skeleton)
|
|
result["grid_res"] = float(grid_res)
|
|
result["signature"] = signature
|
|
|
|
try:
|
|
cache_dir.mkdir(parents=True, exist_ok=True)
|
|
with open(cache_path, "w", encoding="utf-8") as f:
|
|
json.dump(result, f, ensure_ascii=False)
|
|
except Exception:
|
|
pass
|
|
return result
|