diff --git a/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py b/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py index 10d11763..fb95b752 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py +++ b/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py @@ -17,16 +17,17 @@ import numpy as np logger = logging.getLogger(__name__) +# TIN(도엽선)이 물어 오는 격자 밖 여유 — 테두리가 삼각망 밖으로 나가지 않을 만큼만. +SHEET_TIN_CLIP_MARGIN_M = 100.0 + # 화면 버튼에 쓰는 이름 — 키는 surface_models.generation_params.source_filter 접미사다. SHEET_METHOD_LABELS: dict[str, str] = { - "distance": "거리비례", + "tin_sheet": "TIN(도엽선)", + "tin": "TIN 격자", "biharmonic": "TPS(박판)", "anudem": "ANUDEM형", "multires": "다중해상도", "laplace": "라플라스", - "tin": "TIN 선형", - "clough": "TIN 곡면", - "idw": "IDW", } @@ -50,7 +51,7 @@ def _grid_points(spec: Any) -> tuple[np.ndarray, np.ndarray]: return grid_x, grid_y -# ── ① 거리 비례 ────────────────────────────────────────────────────────────── +# ── ① 거리 비례 (버튼에서는 뺐지만 다른 방식의 초기추정으로 계속 쓴다) ─────── def build_distance(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: """가장 가까운 서로 다른 표고 두 라인 사이를 거리 비례로 나눈다. @@ -213,22 +214,18 @@ def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float return solve_min_curvature(burned, build_distance(spec, burned, features, cell_m)) -# ── ④ TIN 선형 / ⑤ TIN 곡면 ───────────────────────────────────────────────── -def _triangulated( - spec: Any, burned: np.ndarray, cell_m: float, smooth: bool, max_points: int = 120_000 +# ── ④ TIN ─────────────────────────────────────────────────────────────────── +def _triangulate( + spec: Any, shape_: tuple[int, int], points: np.ndarray, values: np.ndarray ) -> np.ndarray: - from scipy.interpolate import CloughTocher2DInterpolator, LinearNDInterpolator + """정점 구름을 Delaunay 삼각망 선형 보간해 격자로 편다. 삼각망 밖은 NaN.""" + from scipy.interpolate import LinearNDInterpolator - points, values = _contour_vertices(burned, spec) if len(points) < 3: - return np.full(burned.shape, np.nan, dtype=np.float32) - if len(points) > max_points: # 삼각망 비용은 정점 수에 비례한다 - step = int(np.ceil(len(points) / max_points)) - points, values = points[::step], values[::step] - factory = CloughTocher2DInterpolator if smooth else LinearNDInterpolator - interpolator = factory(points, values) + return np.full(shape_, np.nan, dtype=np.float32) + interpolator = LinearNDInterpolator(points, values) grid_x, grid_y = _grid_points(spec) - surface = np.empty(burned.shape, dtype=np.float32) + surface = np.empty(shape_, dtype=np.float32) chunk = max(1, int(4_000_000 // max(spec.n_cols, 1))) for start in range(0, spec.n_rows, chunk): stop = min(start + chunk, spec.n_rows) @@ -239,44 +236,79 @@ def _triangulated( def build_tin(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: - """등고선 정점 Delaunay 삼각망 선형 보간 — 가장 흔한 고전 방식. + """격자에 구운 등고 라인 셀을 Delaunay 삼각망 선형 보간한다. - 같은 표고 정점 3개로 이루어진 평탄 삼각형이 굴곡부·마루에 계단을 만든다. 비교 - 기준으로 남긴다. + 정점이 셀 중심에 맞춰져 있어 1m 계단이 삼각망에 그대로 실린다. 같은 표고 정점 + 3개로 이루어진 평탄 삼각형이 굴곡부·마루에 계단을 만든다. 비교 기준으로 남긴다. """ - return _triangulated(spec, burned, cell_m, smooth=False) - - -def build_clough(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: - """같은 삼각망을 C1 곡면(Clough–Tocher)으로 채운다 — TIN 계단을 곡면으로 눌러 준다.""" - return _triangulated(spec, burned, cell_m, smooth=True) - - -# ── ⑥ IDW ─────────────────────────────────────────────────────────────────── -def build_idw(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: - """등고 라인 셀을 표본으로 한 역거리가중(k=12, 거듭제곱 2). - - 구현이 단순해 널리 쓰이지만 표본이 곧 등고선이라 라인 주변에 표고가 뭉치는 - '황소눈' 결함이 잘 드러난다. 비교 기준. - """ - from scipy.spatial import cKDTree - points, values = _contour_vertices(burned, spec) - if len(points) < 2: + if len(points) > 120_000: # 삼각망 비용은 정점 수에 비례한다 + step = int(np.ceil(len(points) / 120_000)) + points, values = points[::step], values[::step] + return _triangulate(spec, burned.shape, points, values) + + +def build_tin_sheet(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: + """**원본 5m 도엽 등고선 정점**을 그대로 이은 고전 TIN (2026-08-30 사용자 지시). + + `tin`은 격자에 구운 라인 셀(=1m 계단으로 뭉개진 정점)을 쓰지만, 이쪽은 벡터 + 등고선의 정점을 좌표 그대로 쓴다 — 도면 등고선을 삼각망으로 잇는 측량 관행 그대로다. + 격자 밖 등고선은 물지 않는다(삼각망 비용만 커지고 결과는 같다). 다만 격자 테두리가 + 삼각망 밖으로 나가지 않도록 여유를 두고 자른다. + + 길이 필터는 두지 않는다 — 짧은 봉우리 폐합 등고선을 버리면 마루가 통째로 평평해진다. + """ + from shapely.geometry import shape as to_shape + + from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( + ELEVATION_KEYS, + iter_linestrings, + ) + + xs, ys = spec.cell_centers_x(), spec.cell_centers_y() + margin = max(SHEET_TIN_CLIP_MARGIN_M, cell_m * 2.0) + x_lo, x_hi = xs[0] - margin, xs[-1] + margin + y_lo, y_hi = ys[-1] - margin, ys[0] + margin + + coords: list[np.ndarray] = [] + levels: list[np.ndarray] = [] + for feature in features or []: + properties = feature.get("properties") or {} + elevation = next( + (float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None), + None, + ) + geometry = feature.get("geometry") + if elevation is None or not geometry: + continue + try: + parsed = to_shape(geometry) + except Exception: # noqa: BLE001 — 손상된 피처는 건너뛴다 + continue + for line in iter_linestrings(parsed): + point = np.asarray(line.coords, dtype=np.float64)[:, :2] + inside = ( + (point[:, 0] >= x_lo) + & (point[:, 0] <= x_hi) + & (point[:, 1] >= y_lo) + & (point[:, 1] <= y_hi) + ) + if not inside.any(): + continue + coords.append(point[inside]) + levels.append(np.full(int(inside.sum()), elevation, dtype=np.float64)) + + if not coords: + logger.warning("도엽 서피스: TIN(도엽선)에 쓸 등고선 정점이 없습니다.") return np.full(burned.shape, np.nan, dtype=np.float32) - tree = cKDTree(points) - grid_x, grid_y = _grid_points(spec) - query = np.column_stack([grid_x.ravel(), grid_y.ravel()]) - neighbours = min(12, len(points)) - distance, index = tree.query(query, k=neighbours, workers=-1) - if neighbours == 1: - distance = distance[:, None] - index = index[:, None] - weight = 1.0 / np.maximum(distance, 1e-6) ** 2 - surface = (weight * values[index]).sum(axis=1) / weight.sum(axis=1) - exact = distance[:, 0] < 1e-6 - surface[exact] = values[index[exact, 0]] - return surface.reshape(burned.shape).astype(np.float32) + + points = np.vstack(coords) + values = np.concatenate(levels) + # 도엽 이음매에서 같은 정점이 겹쳐 들어온다 — Qhull 비용만 늘어 미리 접는다. + _, unique = np.unique(np.round(points, 3), axis=0, return_index=True) + points, values = points[unique], values[unique] + logger.info("도엽 서피스: TIN(도엽선) 정점 %d개", len(points)) + return _triangulate(spec, burned.shape, points, values) # ── ⑦ ANUDEM형 (구조선 + 배수 강제) ───────────────────────────────────────── @@ -472,12 +504,10 @@ def build_multires(spec: Any, burned: np.ndarray, features: Any, cell_m: float) SHEET_METHOD_BUILDERS: dict[str, Callable[[Any, np.ndarray, Any, float], np.ndarray]] = { - "distance": build_distance, + "tin_sheet": build_tin_sheet, + "tin": build_tin, "biharmonic": build_biharmonic, "anudem": build_anudem, "multires": build_multires, "laplace": build_laplace, - "tin": build_tin, - "clough": build_clough, - "idw": build_idw, } diff --git a/B04_PreProcess/B04_PreProcess_UI_Page.ts b/B04_PreProcess/B04_PreProcess_UI_Page.ts index b6458a90..ae01a030 100644 --- a/B04_PreProcess/B04_PreProcess_UI_Page.ts +++ b/B04_PreProcess/B04_PreProcess_UI_Page.ts @@ -45,14 +45,12 @@ const DEFAULT_METHOD = "dtm"; const ROUTE_STAGE = ROUTES.B05_PROFILE; // 도엽 서피스 보간 방식 버튼 순서 — 백엔드 SHEET_SURFACE_METHODS와 같은 차례로 둔다. const SHEET_METHOD_ORDER = [ - "distance", + "tin_sheet", + "tin", "biharmonic", "anudem", "multires", "laplace", - "tin", - "clough", - "idw", ]; function L(key: keyof typeof ui_locales): string { diff --git a/config/config_system.py b/config/config_system.py index 06f61284..1b630579 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -194,12 +194,12 @@ SHEET_SURFACE_METHODS = [ method.strip() for method in os.getenv( "SHEET_SURFACE_METHODS", - "distance,biharmonic,anudem,multires,laplace,tin,clough,idw", + "tin_sheet,tin,biharmonic,anudem,multires,laplace", ).split(",") if method.strip() ] # 확정에 쓸 기본 방식 — 비교 후 사용자가 정한다. -SHEET_SURFACE_DEFAULT_METHOD = os.getenv("SHEET_SURFACE_DEFAULT_METHOD", "distance") +SHEET_SURFACE_DEFAULT_METHOD = os.getenv("SHEET_SURFACE_DEFAULT_METHOD", "multires") # 일반 사용자 WF1 자동 확정 기본값 SURFACE_CONFIRM_DEFAULT_FILTER = os.getenv("SURFACE_CONFIRM_DEFAULT_FILTER", "csf")