Merge remote-tracking branch 'origin/feat/las-free-sheet-surface' into feat/B07-cad-block-library

B04_PreProcess/B04_PreProcess_Engine.py

@@ -223,11 +223,64 @@ def run_surface_analysis(
             time.monotonic() - step_started,
         )
 
-    # 3-2. VWorld 지도 및 국가 GIS 벡터 다운로드 (기존 산출물이 있으면 스킵)
+    # 3-2·3-3. VWorld 지도·국가 GIS 벡터·수치지형도 도엽 (공용 블록 — 도엽 서피스도 사용)
     _report(90, "download_maps", "VWorld 지도 및 GIS 벡터 데이터 다운로드 중")
+    las_bounds_dict = {
+        "x": [float(bounds[0, 0]), float(bounds[0, 1])],
+        "y": [float(bounds[1, 0]), float(bounds[1, 1])],
+        "z": [float(bounds[2, 0]), float(bounds[2, 1])],
+    }
+    download_geodata(
+        project_root, processed_dir, las_bounds_dict, las_path.parent, rebuild, report=_report
+    )
+
+    # 3-4. 도엽등고선 3D 서피스 — LAS가 있어도 참고용으로 같이 만들어 영구저장한다
+    # (2026-08-30 사용자 확정). 실패해도 분석은 계속한다.
+    _report(94, "surface_model", "도엽등고선 3D 서피스 생성 중")
+    sheet_models: list[dict[str, Any]] = []
     try:
-        # 입력 LAS와 같은 폴더의 .prj를 우선 사용 (PLAN E-1: 전체 재귀 glob 제거)
-        prj_candidates = sorted(las_path.parent.glob("*.prj")) or sorted(
+        from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import (
+            build_sheet_surface_from_route,
+        )
+
+        sheet_models = build_sheet_surface_from_route(project_root, processed_dir, models_dir)
+    except Exception as exc:
+        logger.warning("도엽등고선 서피스 생성 실패: %s", exc)
+
+    _report(95, "saving", "결과 저장 중")
+
+    return _collect_analysis_result(
+        project_root,
+        models_dir,
+        structured_path,
+        bounds_dict,
+        stats,
+        total_points,
+        ground_summary,
+        manifest,
+        sheet_models,
+        total_started,
+    )
+
+
+def download_geodata(
+    project_root: Path,
+    processed_dir: Path,
+    las_bounds_dict: dict[str, list[float]],
+    prj_search_dir: Path,
+    rebuild: bool,
+    *,
+    default_epsg: str = "EPSG:5186",
+    report: Any = None,
+) -> None:
+    """VWorld 지도·국가 GIS 벡터·수치지형도 도엽 확보 (공용 블록).
+
+    LAS 분석(run_surface_analysis)과 LAS 없는 도엽 서피스 분석이 같이 쓴다.
+    실패해도 예외를 밖으로 던지지 않는다 — 분석 본체를 막지 않는다.
+    """
+    try:
+        # 입력 파일과 같은 폴더의 .prj를 우선 사용 (PLAN E-1: 전체 재귀 glob 제거)
+        prj_candidates = sorted(prj_search_dir.glob("*.prj")) or sorted(
             project_root.glob("B03_FileInput/**/*.prj")
         )
         prj_path = prj_candidates[0] if prj_candidates else project_root / "result.prj"
@@ -243,12 +296,7 @@ def run_surface_analysis(
             get_epsg_from_prj,
         )
 
-        las_bounds_dict = {
-            "x": [float(bounds[0, 0]), float(bounds[0, 1])],
-            "y": [float(bounds[1, 0]), float(bounds[1, 1])],
-            "z": [float(bounds[2, 0]), float(bounds[2, 1])],
-        }
-        project_epsg = "EPSG:5186"
+        project_epsg = default_epsg
         if prj_path.exists():
             project_epsg = get_epsg_from_prj(prj_path.read_text(encoding="utf-8", errors="ignore"))
         # 국가 GIS 벡터는 라이다∪계획노선 범위로 받는다.
@@ -305,7 +353,8 @@ def run_surface_analysis(
         # 3-3. 1:5,000 수치지형도 도엽 확보 → 프로젝트 영구저장소
         # 기준은 계획노선 시점·종점 (같은 도엽이면 9매, 이웃 도엽에 걸치면 12매).
         # (실패해도 분석은 계속 — 폴백은 수동 다운로드 + 인제스트)
-        _report(92, "download_maps", "수치지형도 도엽 확보 중")
+        if report is not None:
+            report(92, "download_maps", "수치지형도 도엽 확보 중")
         try:
             from B04_PreProcess.B04_PreProcess_Engine_Extent import sheet_reference_points_wgs84
             from B04_PreProcess.B04_PreProcess_Engine_MapSheet import neighbors_for_points
@@ -356,8 +405,20 @@ def run_surface_analysis(
     except Exception as exc:
         logger.warning("B04 지도·GIS 다운로드 단계 실패: %s", exc)
 
-    _report(95, "saving", "결과 저장 중")
 
+def _collect_analysis_result(
+    project_root: Path,
+    models_dir: Path,
+    structured_path: Path,
+    bounds_dict: dict[str, float],
+    stats: dict[str, Any],
+    total_points: int,
+    ground_summary: dict[str, Any],
+    manifest: dict[str, Any],
+    sheet_models: list[dict[str, Any]],
+    total_started: float,
+) -> dict[str, Any]:
+    """manifest에서 모델 목록을 추려 분석 결과 dict를 조립한다."""
     processed = {
         "processed_file_path": _relative_to_project(project_root, structured_path),
         "converted_file_path": None,
@@ -403,6 +464,7 @@ def run_surface_analysis(
                     "layers": layers,
                 }
             )
+    models.extend(sheet_models)
 
     logger.info(
         "B04 WF1 분석 완료: 모델 %d개, 총 %.1fs", len(models), time.monotonic() - total_started

B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py

@@ -0,0 +1,513 @@
+"""도엽등고선 → 표고 격자 보간 방식 모음 (비교용).
+
+문헌(Hutchinson 1988/89 ANUDEM; Chaplot 2006; Arun 2013 등)은 지형 복잡도·자료 밀도에
+따라 우열이 갈리며 단일 최적해가 없다고 본다. 그래서 방식을 하나로 고르지 않고 여기
+모아 두고 B04 화면에서 바꿔 가며 보게 한다(2026-08-30 사용자 지시).
+
+각 builder는 `(spec, burned, features, cell_m) -> (R, C) float32` 격자를 돌려준다.
+`burned`는 등고 라인이 구워진 격자(라인 셀 = 표고, 그 외 NaN)다. 폐합 링 안쪽 처리와
+프리뷰·저장은 호출측(`_SheetSurface`)이 방식과 무관하게 똑같이 해 준다.
+"""
+
+import logging
+import warnings
+from typing import Any, Callable
+
+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] = {
+    "tin_sheet": "TIN(도엽선)",
+    "tin": "TIN 격자",
+    "biharmonic": "TPS(박판)",
+    "anudem": "ANUDEM형",
+    "multires": "다중해상도",
+    "laplace": "라플라스",
+}
+
+
+def _laplacian(values: np.ndarray) -> np.ndarray:
+    padded = np.pad(values, 1, mode="edge")
+    return (
+        padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:] - 4.0 * values
+    )
+
+
+def _contour_vertices(burned: np.ndarray, spec: Any) -> tuple[np.ndarray, np.ndarray]:
+    """등고 라인 셀을 (N,2) 세계좌표와 표고로 바꾼다."""
+    rows, cols = np.nonzero(np.isfinite(burned))
+    xs = spec.cell_centers_x()[cols]
+    ys = spec.cell_centers_y()[rows]
+    return np.column_stack([xs, ys]), burned[rows, cols].astype(np.float64)
+
+
+def _grid_points(spec: Any) -> tuple[np.ndarray, np.ndarray]:
+    grid_x, grid_y = np.meshgrid(spec.cell_centers_x(), spec.cell_centers_y())
+    return grid_x, grid_y
+
+
+# ── ① 거리 비례 (버튼에서는 뺐지만 다른 방식의 초기추정으로 계속 쓴다) ───────
+def build_distance(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
+    """가장 가까운 서로 다른 표고 두 라인 사이를 거리 비례로 나눈다.
+
+        z = (L1·d2 + L2·d1) / (d1 + d2)
+
+    지도 제작의 고전적 손보간을 그대로 옮긴 것이다. 원뿔·능선(z=r) 형상을 정확히
+    재현하고 계단이 생기지 않는다. 표고별 거리장을 돌며 가장 작은 두 값을 추적하므로
+    L1≠L2가 보장된다.
+    """
+    from scipy.ndimage import distance_transform_edt
+
+    levels = np.unique(burned[np.isfinite(burned)])
+    if len(levels) < 2:
+        return np.full(burned.shape, levels[0] if len(levels) else np.nan, dtype=np.float32)
+
+    infinity = np.float32(np.inf)
+    first_d = np.full(burned.shape, infinity, dtype=np.float32)
+    first_z = np.zeros(burned.shape, dtype=np.float32)
+    second_d = np.full(burned.shape, infinity, dtype=np.float32)
+    second_z = np.zeros(burned.shape, dtype=np.float32)
+    for level in levels:
+        distance = distance_transform_edt(burned != level, sampling=cell_m).astype(np.float32)
+        beats_first = distance < first_d
+        second_d = np.where(beats_first, first_d, second_d)
+        second_z = np.where(beats_first, first_z, second_z)
+        first_d = np.where(beats_first, distance, first_d)
+        first_z = np.where(beats_first, np.float32(level), first_z)
+        beats_second = ~beats_first & (distance < second_d)
+        second_d = np.where(beats_second, distance, second_d)
+        second_z = np.where(beats_second, np.float32(level), second_z)
+
+    total = first_d + second_d
+    usable = np.isfinite(second_d) & (total > 1e-9)
+    surface = first_z.astype(np.float32)
+    surface[usable] = (
+        (
+            first_z[usable].astype(np.float64) * second_d[usable].astype(np.float64)
+            + second_z[usable].astype(np.float64) * first_d[usable].astype(np.float64)
+        )
+        / total[usable].astype(np.float64)
+    ).astype(np.float32)
+    return surface
+
+
+def relax_laplace(surface: np.ndarray, fixed: np.ndarray, iterations: int) -> None:
+    """등고 라인을 고정한 채 이웃 평균으로 다듬는다 (in-place, red-black 순서)."""
+    if iterations <= 0:
+        return
+    free = ~fixed & np.isfinite(surface)
+    if not free.any():
+        return
+    rows, cols = np.indices(surface.shape)
+    red = free & (((rows + cols) & 1) == 0)
+    black = free & ~red
+    padded = np.zeros((surface.shape[0] + 2, surface.shape[1] + 2), dtype=np.float32)
+    for _ in range(iterations):
+        for colour in (red, black):
+            padded[1:-1, 1:-1] = surface
+            padded[0, 1:-1] = surface[0]
+            padded[-1, 1:-1] = surface[-1]
+            padded[1:-1, 0] = surface[:, 0]
+            padded[1:-1, -1] = surface[:, -1]
+            neighbours = (
+                padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:]
+            ) * np.float32(0.25)
+            surface[colour] = neighbours[colour]
+
+
+# ── ② 라플라스(조화) ─────────────────────────────────────────────────────────
+def build_laplace(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
+    """등고선을 경계값으로 두고 Δz=0을 푼다.
+
+    면이 매끈해지지만 z=r(원뿔·능선)은 조화함수가 아니라 마루가 눌린다. 비교 기준으로
+    남겨 둔다 — ANUDEM이 라플라스 대신 박판 스플라인을 쓰는 이유를 눈으로 보기 위함.
+    """
+    surface = build_distance(spec, burned, features, cell_m)
+    relax_laplace(surface, np.isfinite(burned), 400)
+    return surface
+
+
+# ── ③ 박판 스플라인(중조화) ──────────────────────────────────────────────────
+def solve_min_curvature(constrained: np.ndarray, guess: np.ndarray) -> np.ndarray:
+    """제약 셀을 고정하고 Δ²z=0(박판 스플라인)을 최소곡률 최소제곱으로 푼다.
+
+    ANUDEM/Topo to Raster가 쓰는 박판 스플라인과 같은 연산자다. 라플라스와 달리 z=r을
+    그대로 통과시켜 능선·마루가 눌리지 않고, 경사가 등고선 너머로 자연스럽게 이어진다.
+
+    `constrained`: 값이 고정된 셀(등고 라인, 필요하면 구조선 앵커) — 그 외는 NaN.
+    `guess`: 시작값이자 감쇠 기준(보통 거리 보간 결과).
+    """
+    from scipy.sparse.linalg import LinearOperator, lsmr
+
+    guess = guess.astype(np.float64)
+    burned = constrained
+    # 제약 셀 + **격자 테두리**를 고정한다. 최외곽 등고선 바깥이 통째로 자유면
+    # 1차함수가 Δ²의 영공간에 남아 해가 하나로 정해지지 않고 켤레기울기가 발산한다
+    # (2026-08-30 실측: |Δz| 2092m). 테두리는 거리 보간값으로 묶는다.
+    fixed = np.isfinite(burned)
+    fixed[0, :] = fixed[-1, :] = True
+    fixed[:, 0] = fixed[:, -1] = True
+    free = ~fixed & np.isfinite(guess)
+    if not free.any():
+        return guess.astype(np.float32)
+    index = np.flatnonzero(free.ravel())
+    values = np.where(np.isfinite(burned), burned, guess)
+    base = np.where(fixed, np.nan_to_num(values), 0.0)
+
+    # Δ²z=0을 정규방정식(CG)으로 풀면 조건수가 격자변 4제곱이라 발산한다(실측).
+    # 대신 **최소곡률** 최소제곱으로 세운다 — 자유 셀에 대해 ‖Δz‖를 최소화하며,
+    # 그 정상해가 곧 Δ²z=0이다. 조건수가 제곱으로 줄어 LSMR이 안정적으로 푼다
+    # (Briggs 1974의 최소곡률 격자화와 같은 목적함수).
+    # 최소곡률만으로는 제약(등고선)에서 먼 영역이 정해지지 않아 해가 폭주한다.
+    # ANUDEM의 거칠기 벌점과 같은 취지로 감쇠항을 붙여 거리 보간값에 묶어 둔다:
+    #     minimize ‖Δz‖² + λ‖z − 거리보간‖²
+    # λ가 작을수록 더 매끈하고 클수록 거리 보간에 가깝다.
+    total_cells = base.size
+    free_count = len(index)
+    damping = np.float64(np.sqrt(0.02))
+    anchor = guess.ravel()[index]
+
+    def forward(vector: np.ndarray) -> np.ndarray:
+        # 반드시 **선형**이어야 한다 — 고정 셀 기여(base)를 여기서 더하면 아핀이 되어
+        # LSMR의 전제가 깨지고 해가 폭주한다. base 몫은 우변으로만 넘긴다.
+        scattered = np.zeros_like(base)
+        scattered.ravel()[index] = vector
+        return np.concatenate([_laplacian(scattered).ravel(), damping * vector])
+
+    def adjoint(vector: np.ndarray) -> np.ndarray:
+        curvature = _laplacian(vector[:total_cells].reshape(base.shape)).ravel()[index]
+        return curvature + damping * vector[total_cells:]
+
+    rhs = np.concatenate([-_laplacian(base).ravel(), damping * anchor])
+    linear = LinearOperator(
+        (total_cells + free_count, free_count),
+        matvec=forward,
+        rmatvec=adjoint,
+        dtype=np.float64,
+    )
+    result = lsmr(linear, rhs, x0=anchor, maxiter=400, atol=1e-8, btol=1e-8)
+    solution, info = result[0], result[1]
+    surface = base.copy()
+    surface.ravel()[index] = solution
+    surface[fixed] = values[fixed]
+    # 안전장치 — 발산하면 조용히 틀린 지형을 넘기지 말고 거리 보간으로 되돌린다.
+    drift = float(np.nanmax(np.abs(surface - guess)))
+    span = float(np.nanmax(guess) - np.nanmin(guess))
+    if not np.isfinite(drift) or drift > max(span, 1.0):
+        logger.warning(
+            "도엽 서피스(TPS): 해가 발산해(최대 %.1fm) 거리 보간으로 되돌립니다 (info=%s).",
+            drift,
+            info,
+        )
+        return guess.astype(np.float32)
+    logger.info("도엽 서피스(TPS): 최소제곱 info=%s, 최대 변화 %.2fm", info, drift)
+    return surface.astype(np.float32)
+
+
+def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
+    """등고선만 제약으로 둔 박판 스플라인."""
+    return solve_min_curvature(burned, build_distance(spec, burned, features, cell_m))
+
+
+# ── ④ TIN ───────────────────────────────────────────────────────────────────
+def _triangulate(
+    spec: Any, shape_: tuple[int, int], points: np.ndarray, values: np.ndarray
+) -> np.ndarray:
+    """정점 구름을 Delaunay 삼각망 선형 보간해 격자로 편다. 삼각망 밖은 NaN."""
+    from scipy.interpolate import LinearNDInterpolator
+
+    if len(points) < 3:
+        return np.full(shape_, np.nan, dtype=np.float32)
+    interpolator = LinearNDInterpolator(points, values)
+    grid_x, grid_y = _grid_points(spec)
+    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)
+        surface[start:stop] = interpolator(grid_x[start:stop], grid_y[start:stop]).astype(
+            np.float32
+        )
+    return surface
+
+
+def build_tin(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
+    """격자에 구운 등고 라인 셀을 Delaunay 삼각망 선형 보간한다.
+
+    정점이 셀 중심에 맞춰져 있어 1m 계단이 삼각망에 그대로 실린다. 같은 표고 정점
+    3개로 이루어진 평탄 삼각형이 굴곡부·마루에 계단을 만든다. 비교 기준으로 남긴다.
+    """
+    points, values = _contour_vertices(burned, spec)
+    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)
+
+    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형 (구조선 + 배수 강제) ─────────────────────────────────────────
+def _contour_corner_anchors(
+    spec: Any, burned: np.ndarray, guess: np.ndarray, cell_m: float
+) -> tuple[np.ndarray, np.ndarray] | None:
+    """등고선의 국소 최대 곡률점(코너)에서 능선·계곡 구조선 앵커를 만든다.
+
+    ANUDEM은 등고선 자체의 곡률에서 능선·계곡망을 먼저 뽑아 흐름 구조를 세운다
+    (Hutchinson 1988/89). 여기서도 같은 순서를 따른다.
+
+      ① 라인마다 정점 곡률을 재 국소 최대점(V자 꼭짓점)을 고른다
+      ② 굽은 안쪽이 더 높으면 **계곡**(등고선 V가 상류를 가리킴), 낮으면 **능선**
+      ③ 같은 종류의 코너를 이웃 표고끼리 이어 그 사이를 선형 보간해 앵커로 심는다
+
+    앵커는 박판 해의 제약으로 들어가 계곡 바닥이 이어져 내려가고 능선 마루가 선다.
+    """
+    from scipy.ndimage import label
+    from scipy.spatial import cKDTree
+
+    levels = np.unique(burned[np.isfinite(burned)])
+    if len(levels) < 2:
+        return None
+    interval = float(np.diff(levels).min())
+    xs = spec.cell_centers_x()
+    ys = spec.cell_centers_y()
+
+    corners: list[tuple[float, float, float, int]] = []  # x, y, level, +1 계곡 / -1 능선
+    span = 6  # 곡률을 재는 정점 간격(px) — 짧으면 노이즈, 길면 꼭짓점을 놓친다
+    for level in levels:
+        labelled, count = label(burned == level)
+        for component_id in range(1, count + 1):
+            line_rows, line_cols = np.nonzero(labelled == component_id)
+            if len(line_rows) < 3 * span:
+                continue
+            # 라인 셀을 한 줄로 세운다 — 좌표 정렬로 근사한다(정밀 추적은 과하다).
+            order = np.argsort(line_cols + line_rows * 1e-3)
+            path = np.column_stack([line_cols[order], line_rows[order]]).astype(np.float64)
+            before = np.roll(path, span, axis=0)
+            after = np.roll(path, -span, axis=0)
+            first = path - before
+            second = after - path
+            first_len = np.hypot(first[:, 0], first[:, 1])
+            second_len = np.hypot(second[:, 0], second[:, 1])
+            valid = (first_len > 1e-6) & (second_len > 1e-6)
+            cosine = np.ones(len(path))
+            cosine[valid] = (first[valid] * second[valid]).sum(axis=1) / (
+                first_len[valid] * second_len[valid]
+            )
+            sharp = np.flatnonzero(valid & (cosine < 0.3))  # 70도 이상 꺾인 자리
+            for i in sharp[:: max(1, span)]:
+                # 굽은 안쪽 방향 = 두 변 단위벡터 합의 반대
+                inward = -(first[i] / first_len[i] + second[i] / second_len[i])
+                norm = float(np.hypot(inward[0], inward[1]))
+                if norm < 1e-6:
+                    continue
+                probe = path[i] + inward / norm * 6.0
+                probe_col = int(round(probe[0]))
+                probe_row = int(round(probe[1]))
+                if not (0 <= probe_row < guess.shape[0] and 0 <= probe_col < guess.shape[1]):
+                    continue
+                inside = float(guess[probe_row, probe_col])
+                if not np.isfinite(inside) or abs(inside - level) < interval * 0.15:
+                    continue
+                corners.append(
+                    (
+                        float(xs[int(path[i, 0])]),
+                        float(ys[int(path[i, 1])]),
+                        float(level),
+                        1 if inside > level else -1,
+                    )
+                )
+    if len(corners) < 4:
+        logger.info("도엽 서피스(ANUDEM형): 등고선 코너가 부족해 구조선을 건너뜁니다.")
+        return None
+
+    array = np.asarray(corners, dtype=np.float64)
+    anchor_xy: list[np.ndarray] = []
+    anchor_z: list[np.ndarray] = []
+    reach = interval * 20.0  # 이보다 먼 코너는 같은 구조선으로 보지 않는다
+    for level in levels[:-1]:
+        upper = level + interval
+        lower_set = array[np.abs(array[:, 2] - level) < 1e-6]
+        upper_set = array[np.abs(array[:, 2] - upper) < 1e-6]
+        if not len(lower_set) or not len(upper_set):
+            continue
+        tree = cKDTree(upper_set[:, :2])
+        distance, index = tree.query(lower_set[:, :2], k=1)
+        for i in range(len(lower_set)):
+            j = int(index[i])
+            if distance[i] > reach or lower_set[i, 3] != upper_set[j, 3]:
+                continue
+            start = lower_set[i, :2]
+            end = upper_set[j, :2]
+            steps = max(2, int(distance[i] / max(cell_m, 1e-6) / 4))
+            fraction = np.linspace(0.0, 1.0, steps + 1)[1:-1]
+            if not len(fraction):
+                continue
+            anchor_xy.append(start + (end - start) * fraction[:, None])
+            anchor_z.append(level + interval * fraction)
+    if not anchor_xy:
+        return None
+    return np.vstack(anchor_xy), np.concatenate(anchor_z)
+
+
+def _enforce_drainage(surface: np.ndarray, epsilon: float = 0.01) -> int:
+    """가짜 웅덩이를 메운다 — ANUDEM의 배수 강제와 같은 목적.
+
+    등고선만으로 만든 면에는 흐름이 끊기는 웅덩이가 남는다. 형태학적 재구성(erosion)
+    으로 채우되 완전 평탄해지지 않게 아주 작은 값을 얹는다. 채운 셀 수를 반환한다.
+    """
+    from skimage.morphology import reconstruction
+
+    if not np.isfinite(surface).all():
+        return 0
+    seed = np.full(surface.shape, float(surface.max()), dtype=np.float64)
+    seed[0, :] = surface[0, :]
+    seed[-1, :] = surface[-1, :]
+    seed[:, 0] = surface[:, 0]
+    seed[:, -1] = surface[:, -1]
+    filled = reconstruction(seed, surface.astype(np.float64), method="erosion")
+    raised = filled > surface + 1e-6
+    if not raised.any():
+        return 0
+    # 그냥 채우면 웅덩이가 통째로 평탄해져 흐름 방향이 없어진다. 채운 영역 안쪽으로
+    # 갈수록 아주 조금 높아지게 해서 물이 가장자리(넘침점)로 빠져나가게 둔다
+    # (Garbrecht·Martz의 평탄면 해소를 간단히 옮긴 것 — 표고 변화는 cm 단위다).
+    from scipy.ndimage import distance_transform_edt
+
+    inner = distance_transform_edt(raised)
+    surface[raised] = (filled[raised] + epsilon * inner[raised]).astype(surface.dtype)
+    return int(raised.sum())
+
+
+def build_anudem(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
+    """ANUDEM형 — 등고선 곡률에서 능선·계곡 구조선을 뽑아 제약에 더하고, 박판으로 풀고,
+    가짜 웅덩이를 메운다. Topo to Raster가 밟는 세 단계를 그대로 옮긴 것이다."""
+    guess = build_distance(spec, burned, features, cell_m).astype(np.float64)
+    constrained = burned.astype(np.float64).copy()
+    anchors = _contour_corner_anchors(spec, burned, guess, cell_m)
+    if anchors is not None:
+        xy, z = anchors
+        row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy())
+        inside = (row >= 0) & (col >= 0)
+        row, col, z = row[inside], col[inside], z[inside]
+        free = ~np.isfinite(constrained[row, col])
+        constrained[row[free], col[free]] = z[free]
+        logger.info("도엽 서피스(ANUDEM형): 구조선 앵커 %d셀", int(free.sum()))
+    surface = solve_min_curvature(constrained, guess).astype(np.float32)
+    logger.info("도엽 서피스(ANUDEM형): 가짜 웅덩이 %d셀 메움", _enforce_drainage(surface))
+    return surface
+
+
+# ── ⑧ 다중해상도 (coarse → fine) ────────────────────────────────────────────
+def build_multires(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
+    """성긴 격자에서 풀고 점차 세밀화한다 — ANUDEM의 다중해상도 전략.
+
+    전체 형상은 성긴 격자에서 싸게 잡고, 세밀한 격자에서는 등고선 근처만 다듬는다.
+    한 해상도에서만 풀 때보다 넓은 밴드가 고르게 퍼지고 값싸게 수렴한다.
+    """
+    from scipy.ndimage import zoom
+
+    surface: np.ndarray | None = None
+    for factor in (8, 4, 2, 1):
+        if factor == 1:
+            coarse = burned
+        else:
+            # 성긴 격자의 제약 — 블록 안 등고 라인의 평균 표고를 대표로 쓴다.
+            rows = burned.shape[0] // factor * factor
+            cols = burned.shape[1] // factor * factor
+            blocks = burned[:rows, :cols].reshape(rows // factor, factor, cols // factor, factor)
+            # 라인이 하나도 없는 블록은 NaN이 정상이라 경고를 삼킨다.
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", RuntimeWarning)
+                coarse = np.nanmean(blocks, axis=(1, 3)).astype(np.float32)
+        level = build_distance(spec, coarse, features, cell_m * factor)
+        if surface is not None:
+            # 앞 단계 해를 지금 해상도로 올려 절반씩 섞는다 — 성긴 단계의 넓은 추세를
+            # 이어받되 이번 해상도의 등고선 정보를 덮지 않는다.
+            scale = (level.shape[0] / surface.shape[0], level.shape[1] / surface.shape[1])
+            upscaled = zoom(surface, scale, order=1)
+            free = ~np.isfinite(coarse)
+            level[free] = (level[free] + upscaled[free]) * 0.5
+        relax_laplace(level, np.isfinite(coarse), 8)
+        surface = level
+    assert surface is not None
+    if surface.shape != burned.shape:  # 블록 자르기로 남은 가장자리 보정
+        scale = (burned.shape[0] / surface.shape[0], burned.shape[1] / surface.shape[1])
+        surface = zoom(surface, scale, order=1)
+    line = np.isfinite(burned)
+    surface[line] = burned[line]
+    return surface.astype(np.float32)
+
+
+SHEET_METHOD_BUILDERS: dict[str, Callable[[Any, np.ndarray, Any, float], np.ndarray]] = {
+    "tin_sheet": build_tin_sheet,
+    "tin": build_tin,
+    "biharmonic": build_biharmonic,
+    "anudem": build_anudem,
+    "multires": build_multires,
+    "laplace": build_laplace,
+}

B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py

@@ -0,0 +1,550 @@
+"""도엽등고선 3D 서피스 — 1:5,000 수치지형도 등고선으로 DTM 격자를 만든다.
+
+LAS 없는 설계(2026-08-30 사용자 확정)의 지형 원천이자, LAS가 있어도 참고용으로
+같이 만들어 두는 서피스다. 산출 형식은 LAS 파이프라인의 DTM과 완전히 같게 맞춘다
+(`dtm_sheet.npz`: x/y/z/valid_mask) — 종·횡단·배수 세부설계가 쓰는
+`build_surface_sampler(models_dir, "sheet", "dtm", smooth=False)`가 무수정으로 돈다.
+
+절취 범위는 노선 XY bbox + `SHEET_SURFACE_MARGIN_M`(300m) 직사각형(사용자 확정),
+격자는 `SHEET_SURFACE_GRID_M`(1m).
+
+순서는 **2D 먼저, 메시는 맨 마지막**이다(2026-08-30 사용자 지시):
+  ① 등고 라인을 격자에 굽고 ② 등고선 사이 거리 비례 보간으로 표고 격자를 만든 뒤
+  ③ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장하고 ④ 라인을 고정한 채
+  완화(라플라스)해 등고 간격을 고르게 한다.
+  격자에서 뽑는 1m 등고선이 곧 2D 보간선이며, 메시(glb)는 그 격자의 표현일 뿐이다.
+"""
+
+import json
+import logging
+import time
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+from pyproj import Transformer
+
+from B04_PreProcess.B04_PreProcess_Engine_ModelContext import (
+    atomic_npz,
+    clip_and_compact_mesh,
+    grid_faces,
+    grid_vertices,
+    write_glb,
+)
+from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import (
+    SHEET_METHOD_BUILDERS,
+    SHEET_METHOD_LABELS,
+)
+from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import grid_spec_from_bounds
+from config.config_system import (
+    SHEET_SURFACE_GRID_M,
+    SHEET_SURFACE_MARGIN_M,
+    SHEET_SURFACE_METHODS,
+    SURFACE_MAX_PREVIEW_VERTICES,
+    SURFACE_SMOOTHING_DTM_PREVIEW_RESOLUTION_M,
+    SURFACE_SMOOTHING_DTM_SIGMA_M,
+    SURFACE_SMOOTHING_DTM_SPLINE_SMOOTH,
+)
+
+logger = logging.getLogger(__name__)
+
+# 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값)
+_CONTOUR_FILE = "도엽_등고선.geojson"
+
+# 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다.
+SHEET_SOURCE_FILTER = "sheet"
+
+
+def _load_features_metric(
+    processed_dir: Path, epsg: int, filename: str = _CONTOUR_FILE
+) -> list[dict[str, Any]]:
+    """병합 도엽 레이어(WGS84)를 읽어 사업지 CRS(m)로 재투영한다."""
+    path = processed_dir / filename
+    if not path.exists():
+        logger.warning("도엽 서피스: 도엽 레이어 파일이 없습니다: %s", path)
+        return []
+    try:
+        with path.open("r", encoding="utf-8") as file:
+            data = json.load(file)
+    except (OSError, json.JSONDecodeError):
+        logger.warning("도엽 서피스: 등고선 GeoJSON을 읽지 못했습니다: %s", path)
+        return []
+    features = data.get("features")
+    if not isinstance(features, list):
+        return []
+    transformer = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True)
+
+    def _map(coords: Any) -> Any:
+        if not isinstance(coords, list):
+            return coords
+        if coords and isinstance(coords[0], (int, float)):
+            x, y = transformer.transform(coords[0], coords[1])
+            return [x, y, *coords[2:]]
+        return [_map(item) for item in coords]
+
+    converted: list[dict[str, Any]] = []
+    for feature in features:
+        geometry = feature.get("geometry") or {}
+        coordinates = _map(geometry.get("coordinates"))
+        if coordinates is None:
+            continue
+        converted.append(
+            {
+                "type": "Feature",
+                "properties": feature.get("properties") or {},
+                "geometry": {"type": geometry.get("type"), "coordinates": coordinates},
+            }
+        )
+    return converted
+
+
+def _preview_mesh(
+    x: np.ndarray, y: np.ndarray, z: np.ndarray, valid: np.ndarray
+) -> tuple[np.ndarray, np.ndarray]:
+    """프리뷰용 정점·면 — 정점 수가 상한을 넘으면 격자를 성기게 딴다.
+
+    격자를 그대로 잇는다. 한때 NURBS 곡면을 걸었으나(2026-08-30) DTM 스무딩이
+    들어오면서 곡면 적합이 이중으로 걸려 되돌렸다 — 스무딩은 표고 정본(npz)에서
+    한 번만 한다.
+    """
+    stride = 1
+    while (len(x) // stride + 1) * (len(y) // stride + 1) > SURFACE_MAX_PREVIEW_VERTICES:
+        stride += 1
+    px, py = x[::stride], y[::stride]
+    pz, pv = z[::stride, ::stride], valid[::stride, ::stride]
+    vertices = grid_vertices(px, py, pz.astype(np.float64))
+    faces = grid_faces(len(py), len(px))
+    return clip_and_compact_mesh(vertices, faces, pv.reshape(-1))
+
+
+def _write_smoothed(
+    models_dir: Path,
+    stem: str,
+    x: np.ndarray,
+    y: np.ndarray,
+    z: np.ndarray,
+    valid: np.ndarray,
+    bounds: np.ndarray,
+) -> None:
+    """`{stem}_smooth.npz`·`_smooth_preview.glb`를 만든다 — LAS DTM 스무딩과 같은 절차.
+
+    `B04_PreProcess_Engine_Smooth.smooth_dtm()`은 `TerrainContext`(라이다 발자국)를
+    받으므로 그대로 못 쓴다. 그래서 계수는 **config 값을 그대로** 두고 같은 두 단계만
+    옮긴다(2026-08-30 사용자 지시 — 계수 변경 금지):
+
+      ① 무효 영역이 번지지 않는 정규화 가우시안 (`smoothing_dtm_sigma_meters`)
+      ② C² 바이큐빅 B-spline 재평가 (`kx=ky=3`, `s=smoothing_dtm_spline_smooth`)를
+         `smoothing_dtm_preview_resolution_meters` 격자에서
+
+    화면 스무딩 토글과 확정 스냅샷이 이 파일을 찾으므로 이름 규칙을 지켜야 한다.
+    """
+    from scipy.interpolate import RectBivariateSpline
+
+    from B04_PreProcess.B04_PreProcess_Engine_Smooth import _masked_gaussian_filter
+
+    cell_m = float(x[1] - x[0]) if len(x) > 1 else SHEET_SURFACE_GRID_M
+    sigma_pixels = SURFACE_SMOOTHING_DTM_SIGMA_M / cell_m if cell_m > 0 else 0.0
+    # 결측이 하나라도 있으면 스플라인 결과가 통째로 NaN이 된다(TIN·TIN 곡면은 볼록껍질
+    # 밖이 결측이다). 최근접 표고로 메워 적합하고 아래에서 원래 마스크로 되돌린다.
+    filled = z.astype(np.float64)
+    if not valid.all():
+        from scipy.ndimage import distance_transform_edt
+
+        _, (near_row, near_col) = distance_transform_edt(~valid, return_indices=True)
+        filled = filled[near_row, near_col]
+    z_pre = _masked_gaussian_filter(filled, valid, sigma_pixels)
+    try:
+        spline = RectBivariateSpline(y, x, z_pre, kx=3, ky=3, s=SURFACE_SMOOTHING_DTM_SPLINE_SMOOTH)
+    except Exception as exc:  # noqa: BLE001 — 스무딩 실패가 원본 산출을 막으면 안 된다
+        logger.warning("도엽 서피스(%s): 스무딩 스플라인 실패(%s) — 건너뜁니다.", stem, exc)
+        return
+
+    # 재평가 격자는 config의 프리뷰 해상도를 쓰되 원본보다 성기게 잡지 않는다 —
+    # 이 npz는 화면용이자 **스무딩 확정 시 종·횡단이 샘플링하는 표고 정본**이라,
+    # 정점 상한(화면 사정)으로 해상도를 깎으면 설계 정밀도가 같이 깎인다.
+    # 메시 정점 수는 _preview_mesh가 알아서 성기게 딴다.
+    step = max(SURFACE_SMOOTHING_DTM_PREVIEW_RESOLUTION_M, cell_m)
+    sx = np.arange(x[0], x[-1] + step * 0.5, step)
+    sy = np.arange(y[0], y[-1] + step * 0.5, step)
+    sz = np.asarray(spline(sy, sx), dtype=np.float32)
+    # 원본 유효 마스크를 최근접으로 옮겨 무효 영역을 그대로 지킨다.
+    col = np.clip(np.searchsorted(x, sx) - 1, 0, len(x) - 1)
+    row = np.clip(np.searchsorted(y, sy) - 1, 0, len(y) - 1)
+    svalid = valid[np.ix_(row, col)]
+    sz[~svalid] = np.nan
+
+    atomic_npz(
+        models_dir / f"{stem}_smooth.npz",
+        x=sx,
+        y=sy,
+        z=sz,
+        valid_mask=svalid,
+        bounds=bounds,
+        resolution=np.array([step], np.float32),
+    )
+    vertices, faces = _preview_mesh(sx, sy, np.nan_to_num(sz, nan=float(bounds[2, 0])), svalid)
+    write_glb(models_dir / f"{stem}_smooth_preview.glb", vertices, faces, bounds)
+
+
+def _rasterize_contour_levels(spec: Any, features: list[dict[str, Any]]) -> np.ndarray:
+    """등고 라인을 격자에 굽는다 — 셀 = 그 위를 지나는 라인의 표고, 그 외 NaN.
+
+    배수유역의 `rasterize_contours()`와 두 가지가 다르다(둘 다 서피스 품질 때문이다):
+      · `all_touched=False` — 스치는 셀까지 칠하면 라인이 2px 두께가 되고, 그 폭만큼
+        정확히 등고 표고인 **평탄 띠**가 생겨 사이 1m 등고선 간격이 찌그러진다
+        (2026-08-30 사용자 지적: 보간선이 등간격이 아님).
+      · 길이 필터 없음 — 봉우리 폐합 링 같은 짧은 등고선을 버리면 그 일대가 통째로
+        평평해진다. 배수유역은 노이즈를 버려야 하지만 지형면은 다 있어야 한다.
+
+    같은 셀을 두 표고가 지나면 낮은 쪽을 남긴다(배수유역과 같은 규칙).
+    """
+    from rasterio.features import rasterize
+    from shapely.geometry import shape
+
+    from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
+        ELEVATION_KEYS,
+        grid_transform,
+        iter_linestrings,
+    )
+
+    by_level: dict[float, list[Any]] = {}
+    for feature in features:
+        geometry = feature.get("geometry")
+        if not geometry:
+            continue
+        properties = feature.get("properties") or {}
+        elevation = next(
+            (float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None),
+            None,
+        )
+        if elevation is None:
+            continue
+        try:
+            parsed = shape(geometry)
+        except Exception:  # noqa: BLE001
+            continue
+        for line in iter_linestrings(parsed):
+            by_level.setdefault(elevation, []).append(line)
+
+    burned = np.full((spec.n_rows, spec.n_cols), np.nan, dtype=np.float32)
+    transform = grid_transform(spec)
+    for elevation in sorted(by_level, reverse=True):
+        stamp = rasterize(
+            [(line, 1) for line in by_level[elevation]],
+            out_shape=(spec.n_rows, spec.n_cols),
+            transform=transform,
+            fill=0,
+            dtype="uint8",
+            all_touched=False,
+        ).astype(bool)
+        burned[stamp] = elevation
+    logger.info(
+        "도엽 서피스: 등고 라인 %d단을 격자에 굽어 %d셀",
+        len(by_level),
+        int(np.isfinite(burned).sum()),
+    )
+    return burned
+
+
+def _resolve_enclosed_interiors(
+    burned: np.ndarray, present: list[float], surface: np.ndarray, cell_m: float, interval_m: float
+) -> np.ndarray:
+    """폐합 등고선 안쪽(봉우리·웅덩이)을 바깥 사면 경사로 연장한다 (in-place).
+
+    거리 보간은 "가장 가까운 서로 다른 두 라인 사이"를 채우므로, 마지막 등고선
+    안쪽에는 더 높은 라인이 없어 아래쪽 라인 쪽으로 끌려 **분화구처럼 파인다**.
+    등고선이 'ㅜ'에서 점점 짧아지다 사라지는 마루가 바로 이 자리다(2026-08-30
+    사용자 지적).
+
+    폐합 라인 내부에 다른 표고 제약이 하나도 없으면 그 안이 마루(바깥이 낮을 때)
+    또는 웅덩이(바깥이 높을 때)다. 바깥 사면의 국소 경사를 안쪽으로 연장하되
+    ±(간격−0.5m)로 제한한다 — 다음 등고선이 없다는 사실과 모순되지 않는 범위다.
+    처리한 영역의 마스크를 반환한다 — 이어지는 완화에서 함께 고정해야 한다.
+    """
+    from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label
+
+    constrained = np.isfinite(burned)
+    band_cells = 8
+    limit = max(interval_m - 0.5, 0.5)
+    resolved = 0
+    handled = np.zeros(burned.shape, dtype=bool)
+    for level in present:
+        mask = burned == level
+        interior = binary_fill_holes(mask) & ~mask
+        if not interior.any():
+            continue
+        components, count = label(interior)
+        for component_id in range(1, count + 1):
+            component = components == component_id
+            if (constrained & component).any():
+                continue  # 안에 다른 제약이 있으면 마루가 아니다(보통의 감싸는 링)
+            ring = binary_dilation(binary_fill_holes(mask) | mask) & ~component & ~mask
+            ring &= np.isfinite(surface)
+            if not ring.any():
+                continue
+            outside_mean = float(surface[ring].mean())
+            direction = 1.0 if outside_mean < level else -1.0
+            inner = distance_transform_edt(component, sampling=cell_m)
+            # 바깥 사면 경사 — 라인 밖 band_cells 이내 유효 셀의 (낙차 / 거리) 평균.
+            outer_distance = distance_transform_edt(~(mask | component), sampling=cell_m)
+            band = (outer_distance > 0) & (outer_distance <= band_cells * cell_m)
+            band &= np.isfinite(surface) & ~component & ~mask
+            if band.any():
+                slope = float(
+                    np.mean(np.abs(level - surface[band].astype(np.float64)) / outer_distance[band])
+                )
+            else:
+                slope = 0.0
+            if slope > 1e-3:
+                offset = np.minimum(slope * inner[component], limit)
+            else:
+                peak = float(inner.max())
+                offset = (limit / 2.0) * (inner[component] / peak) if peak > 0 else 0.0
+            surface[component] = level + direction * offset
+            handled |= component
+            resolved += 1
+    if resolved:
+        logger.info("도엽 서피스: 폐합 등고선 내부 %d곳을 사면 경사로 연장", resolved)
+    return handled
+
+
+def _write_method_model(
+    project_root: Path,
+    models_dir: Path,
+    spec: Any,
+    surface: np.ndarray,
+    method_key: str,
+) -> dict[str, Any] | None:
+    """방식 하나의 표고 격자를 npz·프리뷰 glb로 저장하고 등록용 dict를 만든다."""
+    # DtmGridSampler 규약에 맞춰 y 오름차순으로 뒤집어 저장한다.
+    x_coords = spec.cell_centers_x()
+    y_coords = spec.cell_centers_y()[::-1]
+    z_grid = surface[::-1, :].astype(np.float32)
+    valid_grid = np.isfinite(z_grid)
+    if not valid_grid.any():
+        logger.warning("도엽 서피스(%s): 유효 표고 셀이 없습니다.", method_key)
+        return None
+
+    source_filter = f"{SHEET_SOURCE_FILTER}_{method_key}"
+    stem = f"dtm_{source_filter}"
+    model_path = models_dir / f"{stem}.npz"
+    preview_path = models_dir / f"{stem}_preview.glb"
+    finite_z = z_grid[valid_grid]
+    # bounds를 npz에 같이 넣는다 — 등고선 API가 이 값을 화면 원점으로 쓴다. 없으면
+    # LAS structured.npz로 폴백해 메시(glb) 원점과 어긋난다(LAS 없는 설계는 아예 실패).
+    bounds = np.array(
+        [
+            [x_coords[0], x_coords[-1]],
+            [y_coords[0], y_coords[-1]],
+            [float(finite_z.min()), float(finite_z.max())],
+        ]
+    )
+    atomic_npz(
+        model_path,
+        x=x_coords,
+        y=y_coords,
+        z=z_grid,
+        valid_mask=valid_grid,
+        bounds=bounds,
+        resolution=np.array([SHEET_SURFACE_GRID_M], np.float32),
+    )
+    vertices, faces = _preview_mesh(x_coords, y_coords, z_grid, valid_grid)
+    write_glb(preview_path, vertices, faces, bounds)
+    _write_smoothed(models_dir, stem, x_coords, y_coords, z_grid, valid_grid, bounds)
+    return {
+        "model_type": "dtm",
+        "source_filter": source_filter,
+        "representation": "regular_grid",
+        "model_file_path": str(model_path.relative_to(project_root)).replace("\\", "/"),
+        "resolution_m": SHEET_SURFACE_GRID_M,
+        "generation_params": {
+            "source_filter": source_filter,
+            "representation": "regular_grid",
+            "source": "map_sheet_contours",
+            "interpolation": method_key,
+            "interpolation_label": SHEET_METHOD_LABELS.get(method_key, method_key),
+            "margin_m": SHEET_SURFACE_MARGIN_M,
+        },
+        "layers": [
+            {
+                "layer_name": f"{stem}_preview",
+                "geometry_type": "MESH",
+                "file_path": str(preview_path.relative_to(project_root)).replace("\\", "/"),
+                "file_format": "glb",
+            }
+        ],
+    }
+
+
+def build_sheet_surface_model(
+    project_root: Path,
+    processed_dir: Path,
+    models_dir: Path,
+    route_xy: np.ndarray,
+    epsg: int,
+    methods: list[str] | None = None,
+) -> list[dict[str, Any]]:
+    """도엽등고선으로 방식별 DTM npz·프리뷰 glb를 만들고 등록용 dict 목록을 돌려준다.
+
+    방식을 하나로 고르지 않고 전부 만들어 두는 이유: 문헌상 지형에 따라 우열이 갈려
+    화면에서 바꿔 보며 정해야 한다(2026-08-30 사용자 지시). 실패하면 빈 목록 —
+    호출측은 분석을 계속한다(도엽 미확보 지역 폴백).
+
+    `route_xy`: (N, 2) 노선 정점 XY(사업지 CRS, m).
+    """
+    started = time.monotonic()
+    features = _load_features_metric(processed_dir, epsg)
+    if not features:
+        return []
+
+    x_min = float(np.min(route_xy[:, 0])) - SHEET_SURFACE_MARGIN_M
+    x_max = float(np.max(route_xy[:, 0])) + SHEET_SURFACE_MARGIN_M
+    y_min = float(np.min(route_xy[:, 1])) - SHEET_SURFACE_MARGIN_M
+    y_max = float(np.max(route_xy[:, 1])) + SHEET_SURFACE_MARGIN_M
+    spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M)
+
+    # ① 2D — 등고 라인을 격자에 굽는다(라인 셀 = 표고, 그 외 NaN).
+    burned = _rasterize_contour_levels(spec, features)
+    present = sorted(np.unique(burned[np.isfinite(burned)]).tolist())
+    if len(present) < 2:
+        logger.warning("도엽 서피스: 절취 범위 안에 등고선이 부족합니다.")
+        return []
+    # 등고 간격(m) — 마루 연장 상한의 근거. 레벨이 하나뿐이면 5m(주곡선) 폴백.
+    interval_m = float(np.diff(np.array(present)).min()) if len(present) > 1 else 5.0
+
+    selected = methods or list(SHEET_SURFACE_METHODS)
+    models: list[dict[str, Any]] = []
+    for method_key in selected:
+        builder = SHEET_METHOD_BUILDERS.get(method_key)
+        if builder is None:
+            logger.warning("도엽 서피스: 알 수 없는 보간 방식 %s — 건너뜁니다.", method_key)
+            continue
+        step_started = time.monotonic()
+        try:
+            # ② 2D 보간 — 여기서 나온 격자에서 1m 등고선을 뽑으므로 화면 등고선이 곧
+            #    2D 보간선이다. 메시(glb)는 그 격자의 표현일 뿐이다(사용자 지시).
+            surface = builder(spec, burned, features, spec.cell_m)
+            # ③ 폐합 등고선 안쪽(마루·웅덩이)은 방식과 무관하게 같은 규칙으로 채운다.
+            _resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m)
+        except Exception as exc:  # noqa: BLE001 — 한 방식이 죽어도 나머지는 만든다
+            logger.warning("도엽 서피스(%s) 생성 실패: %s", method_key, exc)
+            continue
+        model = _write_method_model(project_root, models_dir, spec, surface, method_key)
+        if model is not None:
+            models.append(model)
+            logger.info("도엽 서피스(%s) 완료 (%.1fs)", method_key, time.monotonic() - step_started)
+
+    logger.info(
+        "도엽 서피스 생성 완료: %d×%d 격자, 등고 %d단, 방식 %d개 (%.1fs)",
+        spec.n_rows,
+        spec.n_cols,
+        len(present),
+        len(models),
+        time.monotonic() - started,
+    )
+    return models
+
+
+def build_sheet_surface_from_route(
+    project_root: Path, processed_dir: Path, models_dir: Path
+) -> list[dict[str, Any]]:
+    """B03 업로드 계획노선 CSV를 찾아 방식별 도엽 서피스를 만든다. 없으면 빈 목록."""
+    from common_util.common_util_route_geometry import (
+        find_planned_route_file,
+        read_planned_route_csv,
+    )
+
+    route_file = find_planned_route_file(project_root / "B03_FileInput" / "input")
+    if route_file is None:
+        logger.warning("도엽 서피스: 계획 노선 파일이 없습니다.")
+        return []
+    planned = read_planned_route_csv(route_file)
+    if planned is None or len(planned.vertices) < 2:
+        logger.warning("도엽 서피스: 계획 노선 파일을 읽지 못했습니다: %s", route_file.name)
+        return []
+    route_xy = np.array([(v.x, v.y) for v in planned.vertices], dtype=np.float64)
+    return build_sheet_surface_model(
+        project_root, processed_dir, models_dir, route_xy, planned.epsg or 5186
+    )
+
+
+def run_sheet_surface_analysis(
+    project_root: Path,
+    route_csv_path: Path,
+    *,
+    on_progress: Any = None,
+) -> dict[str, Any]:
+    """LAS 없는 WF1 — 도엽 확보 후 도엽등고선 서피스만으로 분석 결과를 만든다.
+
+    반환 형식은 `run_surface_analysis()`와 같다(save_surface_analysis_to_db 호환).
+    """
+    from common_util.common_util_route_geometry import read_planned_route_csv
+
+    def _report(percent: int, stage: str, message: str) -> None:
+        if on_progress is not None:
+            on_progress(percent, stage, message)
+
+    stage_root = project_root / "B04_PreProcess"
+    processed_dir = stage_root / "processed"
+    models_dir = stage_root / "models"
+    processed_dir.mkdir(parents=True, exist_ok=True)
+    models_dir.mkdir(parents=True, exist_ok=True)
+
+    planned = read_planned_route_csv(route_csv_path)
+    if planned is None or len(planned.vertices) < 2:
+        raise ValueError(f"계획 노선 파일을 읽지 못했습니다: {route_csv_path.name}")
+    epsg = planned.epsg or 5186
+    route_xy = np.array([(v.x, v.y) for v in planned.vertices], dtype=np.float64)
+
+    bounds_dict = {
+        "x": [float(route_xy[:, 0].min()), float(route_xy[:, 0].max())],
+        "y": [float(route_xy[:, 1].min()), float(route_xy[:, 1].max())],
+        "z": [
+            float(min(v.z for v in planned.vertices)),
+            float(max(v.z for v in planned.vertices)),
+        ],
+    }
+
+    # VWorld 지도·GIS 벡터·도엽 확보 — LAS 경로와 같은 공용 블록 (지연 import로 순환 회피)
+    _report(30, "download_maps", "VWorld 지도 및 수치지형도 도엽 확보 중")
+    from B04_PreProcess.B04_PreProcess_Engine import download_geodata
+
+    download_geodata(
+        project_root,
+        processed_dir,
+        bounds_dict,
+        route_csv_path.parent,
+        rebuild=False,
+        default_epsg=f"EPSG:{epsg}",
+        report=_report,
+    )
+
+    _report(70, "surface_model", "도엽등고선 3D 서피스 생성 중")
+    models = build_sheet_surface_model(project_root, processed_dir, models_dir, route_xy, epsg)
+    if not models:
+        raise ValueError("도엽등고선으로 지표면을 만들지 못했습니다 — 도엽 확보를 확인하세요.")
+
+    _report(95, "saving", "결과 저장 중")
+    return {
+        "processed": {
+            "processed_file_path": str(
+                (processed_dir / _CONTOUR_FILE).relative_to(project_root)
+            ).replace("\\", "/"),
+            "converted_file_path": None,
+            "point_count": int(len(route_xy)),
+            "bounds": {
+                "x_min": bounds_dict["x"][0],
+                "x_max": bounds_dict["x"][1],
+                "y_min": bounds_dict["y"][0],
+                "y_max": bounds_dict["y"][1],
+            },
+            "statistics": {
+                "min_z": bounds_dict["z"][0],
+                "max_z": bounds_dict["z"][1],
+                "mean_z": None,
+            },
+        },
+        "ground_summary": {},
+        "manifest": {"status": "sheet_only"},
+        "models": models,
+    }

B04_PreProcess/B04_PreProcess_Router.py

@@ -432,8 +432,26 @@ async def get_confirmed_surface(project_id: UUID) -> SurfaceConfirmedResponse |
                     "z_max": float(bounds[2, 1]),
                 }
 
-        # 지도(2D) 초기 화면을 도로 기준으로 맞추기 위한 계획노선 범위(없으면 None).
+        # LAS 없이 설계한 프로젝트는 위 두 파일이 아예 없다(도엽등고선으로 만든
+        # 서피스가 정본). 그때는 확정 모델 격자의 bounds를 그대로 쓴다 — 없으면
+        # B05가 "지표면 범위 정보를 찾을 수 없습니다"로 3D를 못 띄운다(2026-08-30).
         project_root = processed_dir.parent.parent
+        if bounds_payload is None and confirmed and confirmed.get("model_file_path"):
+            model_path = project_root / str(confirmed["model_file_path"])
+            if model_path.is_file():
+                with np.load(model_path) as stored:
+                    if "bounds" in stored:
+                        bounds = np.asarray(stored["bounds"], dtype=np.float64)
+                        bounds_payload = {
+                            "x_min": float(bounds[0, 0]),
+                            "x_max": float(bounds[0, 1]),
+                            "y_min": float(bounds[1, 0]),
+                            "y_max": float(bounds[1, 1]),
+                            "z_min": float(bounds[2, 0]),
+                            "z_max": float(bounds[2, 1]),
+                        }
+
+        # 지도(2D) 초기 화면을 도로 기준으로 맞추기 위한 계획노선 범위(없으면 None).
         route_bounds = planned_route_bounds(project_root, project_epsg_from_prj(project_root))
 
         signature = "|".join(

B04_PreProcess/B04_PreProcess_Router_Basins.py

@@ -239,7 +239,9 @@ async def _resolve(
 
     requested = parse_pipe_points((payload or {}).get("points"))
     signature = route_signature(context.vertices)
-    stored = load_pipe_points(context.stored_path, signature) if use_stored else None
+    stored = (
+        load_pipe_points(context.stored_path, signature, context.vertices) if use_stored else None
+    )
     points = requested or stored or None
 
     result = await asyncio.to_thread(_build, context.stored_path, context, points)
@@ -305,7 +307,10 @@ async def put_pipe_points(
     context, detail, points, _ = resolved
 
     signature = route_signature(context.vertices)
-    saved = await asyncio.to_thread(save_pipe_points, context.stored_path, signature, points)
+    # 좌표를 같이 남긴다 — 다른 선(B05 최적 경로)으로 읽어도 그 자리에 되놓는다.
+    saved = await asyncio.to_thread(
+        save_pipe_points, context.stored_path, signature, points, context.vertices
+    )
     await asyncio.to_thread(
         save_detail_basins, context.stored_path, _basin_features(context, detail, points)
     )

B04_PreProcess/B04_PreProcess_UI_Page.ts

@@ -9,7 +9,10 @@ import {
 } from "@ui/ui_template_elements";
 import { createWorkflowLayout } from "@ui/ui_template_workflow_layout";
 import { fetchDashboardMe } from "../B01_Dashboard/B01_Dashboard_Api_Fetch";
-import { clearPreloadMark, purgeOtherProjects } from "../A00_Common/b_asset_cache";
+import {
+  clearPreloadMark,
+  purgeOtherProjects,
+} from "../A00_Common/b_asset_cache";
 import { clearRouteLatestCache } from "../B05_Profile/B05_Profile_Api_Fetch";
 import { workflowSteps } from "../A00_Common/b_page_scaffold";
 import {
@@ -40,6 +43,15 @@ const MODEL_METHODS = ["tin", "dtm", "nurbs", "implicit", "meshfree"] as const;
 const DEFAULT_FILTER = "csf";
 const DEFAULT_METHOD = "dtm";
 const ROUTE_STAGE = ROUTES.B05_PROFILE;
+// 도엽 서피스 보간 방식 버튼 순서 — 백엔드 SHEET_SURFACE_METHODS와 같은 차례로 둔다.
+const SHEET_METHOD_ORDER = [
+  "tin_sheet",
+  "tin",
+  "biharmonic",
+  "anudem",
+  "multires",
+  "laplace",
+];
 
 function L(key: keyof typeof ui_locales): string {
   return ui_locales[key][currentLanguageIndex];
@@ -105,11 +117,17 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
   if (guardedProjectId) {
     const user = await fetchDashboardMe();
     if (user.role !== "SYSTEM_ADMIN") {
-      const workflowState = await fetchWorkflowState(guardedProjectId).catch(() => undefined);
-      const surfaceStage = workflowState?.stages.find((stage) => stage.stage_no === 1);
+      const workflowState = await fetchWorkflowState(guardedProjectId).catch(
+        () => undefined,
+      );
+      const surfaceStage = workflowState?.stages.find(
+        (stage) => stage.stage_no === 1,
+      );
       goToWorkflowStage(
         guardedProjectId,
-        surfaceStage?.state === "COMPLETE" ? ROUTES.B05_PROFILE : ROUTES.B03_FILE_INPUT,
+        surfaceStage?.state === "COMPLETE"
+          ? ROUTES.B05_PROFILE
+          : ROUTES.B03_FILE_INPUT,
       );
       return;
     }
@@ -140,6 +158,58 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
   const viewer = createSurfacePointCloudViewer();
   const terrainViewer = createSurfaceTerrainViewer();
   const mapViewer = createSurfaceMapViewer();
+  // 도엽등고 3D 서피스 — 전처리에서 함께 생성되는 참고 서피스(LAS 없는 설계의 지형 원천,
+  // 2026-08-30). 모델 목록에 sheet/dtm이 있을 때만 별도 컨테이너로 보여준다.
+  const sheetViewer = createSurfaceTerrainViewer();
+  const sheetSection = document.createElement("section");
+  sheetSection.className = "b04-surface__sheet-section ui-sidebar-section";
+  const sheetTitle = document.createElement("h3");
+  sheetTitle.className = "b04-surface__panel-title";
+  sheetTitle.textContent = L("B04_Surface_SheetSurface");
+  // 보간 방식 전환 줄 — 어느 방식이 이 지형에 맞는지 눈으로 비교해 정한다
+  // (2026-08-30 사용자 지시). 버튼 목록은 실제 생성된 모델에서 만든다.
+  const sheetToolbar = document.createElement("div");
+  sheetToolbar.className = "b04-surface__sheet-toolbar";
+  const sheetMethodButtons = new Map<string, HTMLButtonElement>();
+  let sheetMethod = "";
+
+  function selectSheetMethod(method: string): void {
+    sheetMethod = method;
+    for (const [key, button] of sheetMethodButtons) {
+      button.classList.toggle("is-active", key === method);
+    }
+    const projectId = getProjectId();
+    if (!projectId) return;
+    sheetViewer.setSelection(`sheet_${method}`, "dtm");
+    sheetViewer.render(projectId, models);
+  }
+
+  // 라이다 지표면 겹쳐 보기 — 확정 필터의 DTM을 반투명으로 얹는다.
+  const lidarLabel = document.createElement("label");
+  lidarLabel.className = "toggle-label toggle-button b04-surface__sheet-lidar";
+  const lidarCheck = document.createElement("input");
+  lidarCheck.type = "checkbox";
+  lidarLabel.append(
+    lidarCheck,
+    document.createTextNode(` ${L("B04_Surface_SheetLidar")}`),
+  );
+  lidarCheck.addEventListener("change", () => {
+    void sheetViewer
+      .showOverlay(
+        lidarCheck.checked ? filterGroup.select.value : "",
+        "dtm",
+        terrainViewer.isSmoothingEnabled(),
+      )
+      .then((loaded) => {
+        if (lidarCheck.checked && !loaded) {
+          showToast(L("B04_Surface_SheetLidar_Missing"), "warning");
+          lidarCheck.checked = false;
+        }
+      });
+  });
+
+  sheetSection.append(sheetTitle, sheetToolbar, sheetViewer.root);
+  sheetSection.hidden = true;
 
   let syncingCamera = false;
   viewer.onCameraChange((state) => {
@@ -200,14 +270,20 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
   const actionRow = document.createElement("div");
   actionRow.className = "ui-sidebar-actions";
   actionRow.append(confirmButton, resetButton);
-  panel.append(inputGroup, analysisGroup, displayGroup, viewer.controlsGroup, actionRow);
+  panel.append(
+    inputGroup,
+    analysisGroup,
+    displayGroup,
+    viewer.controlsGroup,
+    actionRow,
+  );
 
   const viewers = document.createElement("div");
   viewers.className = "b04-surface__viewers";
   viewers.append(viewer.root, terrainViewer.root);
   const workspace = document.createElement("div");
   workspace.className = "b04-surface__workspace";
-  workspace.append(statusBox, viewers, mapViewer.root);
+  workspace.append(statusBox, viewers, sheetSection, mapViewer.root);
 
   let workflowState: WorkflowState | undefined;
   const layoutProjectId = localStorage.getItem(CURRENT_PROJECT_ID_KEY);
@@ -229,7 +305,8 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
     currentStage: workflowState?.current_stage,
     routes: WORKFLOW_STEP_ROUTES,
     onStepClick: (stepIndex) => {
-      if (layoutProjectId) goToWorkflowStage(layoutProjectId, WORKFLOW_STEP_ROUTES[stepIndex]);
+      if (layoutProjectId)
+        goToWorkflowStage(layoutProjectId, WORKFLOW_STEP_ROUTES[stepIndex]);
     },
   });
 
@@ -263,7 +340,11 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
       return;
     }
     const variant =
-      status.status === "completed" ? "success" : status.status === "failed" ? "danger" : "warning";
+      status.status === "completed"
+        ? "success"
+        : status.status === "failed"
+          ? "danger"
+          : "warning";
     statusBox.append(
       createTag(`${status.progress_percent}%`, variant),
       document.createTextNode(status.message),
@@ -280,7 +361,9 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
     inputInfo.append(
       buildInfoLine(
         "좌표계",
-        selectedInputFile.crs_epsg ? `EPSG:${selectedInputFile.crs_epsg}` : null,
+        selectedInputFile.crs_epsg
+          ? `EPSG:${selectedInputFile.crs_epsg}`
+          : null,
       ),
       buildInfoLine(
         "크기",
@@ -289,7 +372,10 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
           : `${selectedInputFile.file_size_mb.toFixed(2)} MB`,
       ),
       buildInfoLine("포인트 수", pointCloud?.point_count.toLocaleString()),
-      buildInfoLine("표시 포인트 수", pointCloud?.sampled_count.toLocaleString()),
+      buildInfoLine(
+        "표시 포인트 수",
+        pointCloud?.sampled_count.toLocaleString(),
+      ),
       buildInfoLine("높이 범위", heightRange),
     );
   }
@@ -328,7 +414,10 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
   function updateSelectedModel(): void {
     const projectId = getProjectId();
     if (!projectId) return;
-    terrainViewer.setSelection(filterGroup.select.value, methodGroup.select.value);
+    terrainViewer.setSelection(
+      filterGroup.select.value,
+      methodGroup.select.value,
+    );
     terrainViewer.render(projectId, models);
     confirmButton.disabled = !findSelectedModel();
   }
@@ -339,14 +428,20 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
     showLoadingOverlay();
     viewer.setLoading("포인트 데이터 로딩 중…");
     try {
-      pointCloud = await fetchSurfacePointCloud(projectId, filterGroup.select.value);
+      pointCloud = await fetchSurfacePointCloud(
+        projectId,
+        filterGroup.select.value,
+      );
       terrainViewer.setReferenceBounds(pointCloud.bounds);
       viewer.render(pointCloud);
       renderInputInfo();
     } catch (error) {
       pointCloud = null;
       viewer.render(null);
-      const detail = error instanceof Error ? error.message : "지면 포인트 조회에 실패했습니다.";
+      const detail =
+        error instanceof Error
+          ? error.message
+          : "지면 포인트 조회에 실패했습니다.";
       showToast(detail, "error");
     } finally {
       hideLoadingOverlay();
@@ -366,7 +461,8 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
     // 확정본과 같은 조합에서 시작해야 B05와 같은 파일을 보고, 보관함도 한 벌만 쓴다.
     // 확정 이력이 없을 때만 개발 기본값(csf·dtm)으로 둔다.
     if (confirmed.model_id) {
-      if (confirmed.source_filter) filterGroup.select.value = confirmed.source_filter;
+      if (confirmed.source_filter)
+        filterGroup.select.value = confirmed.source_filter;
       if (confirmed.method) methodGroup.select.value = confirmed.method;
       terrainViewer.setSmoothing(confirmed.smooth ?? false);
     }
@@ -376,7 +472,10 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
     renderStatus(status);
     viewer.setLoading("포인트 데이터 로딩 중…");
     try {
-      pointCloud = await fetchSurfacePointCloud(projectId, filterGroup.select.value);
+      pointCloud = await fetchSurfacePointCloud(
+        projectId,
+        filterGroup.select.value,
+      );
       terrainViewer.setReferenceBounds(pointCloud.bounds);
       viewer.render(pointCloud);
       // 지도(2D)는 계획노선 기준으로 연다 — 라이다 범위와 다루는 범위가 다르다.
@@ -388,6 +487,51 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
     }
     renderInputInfo();
     updateSelectedModel();
+
+    // 도엽등고 3D 서피스 — sheet_* 모델이 있으면 별도 컨테이너로 보여준다.
+    // 보간 방식마다 모델이 하나씩 있으므로 버튼으로 갈아 끼운다.
+    const sheetMethods = models
+      .filter(
+        (model) =>
+          model.model_type.toLowerCase() === "dtm" &&
+          getModelFilter(model).startsWith("sheet_"),
+      )
+      .map((model) => ({
+        key: getModelFilter(model).slice("sheet_".length),
+        label:
+          typeof model.generation_params?.interpolation_label === "string"
+            ? (model.generation_params.interpolation_label as string)
+            : getModelFilter(model).slice("sheet_".length),
+      }))
+      // 모델 목록은 최신순이라 버튼이 뒤섞인다 — 정의 순서로 고정한다.
+      .sort(
+        (a, b) =>
+          (SHEET_METHOD_ORDER.indexOf(a.key) + 1 || 99) -
+          (SHEET_METHOD_ORDER.indexOf(b.key) + 1 || 99),
+      );
+    sheetSection.hidden = sheetMethods.length === 0;
+    if (sheetMethods.length) {
+      sheetToolbar.replaceChildren();
+      sheetMethodButtons.clear();
+      for (const method of sheetMethods) {
+        const button = document.createElement("button");
+        button.type = "button";
+        button.className = "b04-surface__sheet-method";
+        button.textContent = method.label;
+        button.addEventListener("click", () => selectSheetMethod(method.key));
+        sheetToolbar.append(button);
+        sheetMethodButtons.set(method.key, button);
+      }
+      // 스무딩 드롭다운과 라이다 토글은 오른쪽 끝에 함께 둔다.
+      sheetViewer.smoothingField.classList.add("b04-surface__sheet-smoothing");
+      sheetToolbar.append(sheetViewer.smoothingField, lidarLabel);
+      sheetViewer.setSmoothing(true);
+      selectSheetMethod(
+        sheetMethods.some((method) => method.key === sheetMethod)
+          ? sheetMethod
+          : sheetMethods[0].key,
+      );
+    }
   }
 
   async function onB04_Surface_Confirm_Click(): Promise<void> {
@@ -414,14 +558,20 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
         L("B04_Surface_Confirm_Success")
           .replace("{filter}", filterGroup.select.value)
           .replace("{method}", methodGroup.select.value)
-          .replace("{smoothing}", terrainViewer.isSmoothingEnabled() ? "ON" : "OFF"),
+          .replace(
+            "{smoothing}",
+            terrainViewer.isSmoothingEnabled() ? "ON" : "OFF",
+          ),
         "success",
       );
       await loadProjectData(projectId);
       enableRouteStep(projectId);
       goToWorkflowStage(projectId, ROUTE_STAGE);
     } catch (error) {
-      const detail = error instanceof Error ? error.message : L("B04_Surface_Confirm_Failed");
+      const detail =
+        error instanceof Error
+          ? error.message
+          : L("B04_Surface_Confirm_Failed");
       showToast(`${L("B04_Surface_Confirm_Failed")} ${detail}`, "error");
     } finally {
       hideLoadingOverlay();
@@ -446,7 +596,8 @@ export async function renderB04Surface(root: HTMLElement): Promise<void> {
   }
 
   inputSelect.addEventListener("change", () => {
-    selectedInputFile = inputFiles.find((file) => String(file.id) === inputSelect.value) ?? null;
+    selectedInputFile =
+      inputFiles.find((file) => String(file.id) === inputSelect.value) ?? null;
     renderInputInfo();
   });
   filterGroup.select.addEventListener("change", () => {

B04_PreProcess/B04_PreProcess_UI_Style.css

@@ -793,3 +793,65 @@
     grid-template-columns: 1fr;
   }
 }
+
+/* 도엽등고 3D 서피스 컨테이너 — 2026-08-30 */
+.b04-surface__sheet-section {
+  margin: 0 var(--spacing-24) var(--spacing-16);
+  padding: var(--spacing-16);
+  box-sizing: border-box;
+}
+
+.b04-surface__sheet-section > .terrain-model-group {
+  margin-top: var(--spacing-12);
+}
+
+/* 도엽 서피스 보간 방식 전환 줄 — 2026-08-30 */
+.b04-surface__sheet-toolbar {
+  display: flex;
+  flex-wrap: wrap;
+  align-items: center;
+  gap: var(--spacing-8);
+  margin: var(--spacing-12) 0;
+}
+
+.b04-surface__sheet-method {
+  padding: 4px 10px;
+  font-size: var(--text-caption);
+  color: var(--color-text-body);
+  background: var(--color-surface);
+  border: 1px solid var(--color-border);
+  border-radius: var(--radius-cards);
+  cursor: pointer;
+}
+
+.b04-surface__sheet-method:hover {
+  border-color: var(--color-primary);
+}
+
+.b04-surface__sheet-method.is-active {
+  color: var(--color-on-primary, #fff);
+  background: var(--color-primary);
+  border-color: var(--color-primary);
+}
+
+.b04-surface__sheet-lidar {
+  margin-left: auto;
+}
+
+/* 도엽 서피스 스무딩 드롭다운 — 방식 버튼 줄 오른쪽 끝 */
+.b04-surface__sheet-smoothing {
+  margin-left: auto;
+  display: flex;
+  align-items: center;
+  gap: var(--spacing-8);
+  font-size: var(--text-caption);
+}
+
+.b04-surface__sheet-smoothing .b04-surface__select {
+  width: auto;
+  min-width: 96px;
+}
+
+.b04-surface__sheet-toolbar .b04-surface__sheet-lidar {
+  margin-left: 0;
+}

B04_PreProcess/B04_PreProcess_UI_TerrainViewer.ts

@@ -6,7 +6,10 @@ import { API_BASE_URL } from "@config/config_frontend";
 import { fetchCachedBytes, fetchCachedJson } from "../A00_Common/b_asset_cache";
 import { currentLanguageIndex, ui_locales } from "@ui/ui_template_locale";
 import { createProgressCircle } from "@ui/ui_template_progress";
-import type { SurfaceBounds, SurfaceModelSummary } from "./B04_PreProcess_Api_Fetch";
+import type {
+  SurfaceBounds,
+  SurfaceModelSummary,
+} from "./B04_PreProcess_Api_Fetch";
 import {
   bindCursorPivotControls,
   bindSurfaceViewerTheme,
@@ -18,6 +21,10 @@ import {
   type SurfaceCameraState,
 } from "./B04_PreProcess_UI_Camera";
 
+/** 화면에 띄우는 등고 라벨 상한 — 긴 등고선부터 채운다. 조각이 많은 지형에서
+ *  라벨이 수백 개가 되면 매 프레임 위치 재계산이 화면을 멈춰 세운다(2026-08-30). */
+const MAX_CONTOUR_LABELS = 40;
+
 function L(key: keyof typeof ui_locales): string {
   return ui_locales[key][currentLanguageIndex];
 }
@@ -31,6 +38,12 @@ export interface SurfaceTerrainViewer {
   render: (projectId: string, models: readonly SurfaceModelSummary[]) => void;
   setReferenceBounds: (bounds: SurfaceBounds) => void;
   setSelection: (sourceFilter: string, method: string) => void;
+  /** 다른 모델(예: 라이다 지표면)을 반투명으로 겹쳐 본다. 빈 문자열이면 걷어낸다. */
+  showOverlay: (
+    sourceFilter: string,
+    method: string,
+    smooth: boolean,
+  ) => Promise<boolean>;
   applyCameraState: (state: SurfaceCameraState) => void;
   onCameraChange: (listener: (state: SurfaceCameraState) => void) => void;
   onAxesVisibilityChange: (listener: (visible: boolean) => void) => void;
@@ -225,7 +238,12 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
     scene.background = new THREE.Color(color);
   });
 
-  const camera = new THREE.PerspectiveCamera(SURFACE_CAMERA_FOV, 1, 0.01, 100000);
+  const camera = new THREE.PerspectiveCamera(
+    SURFACE_CAMERA_FOV,
+    1,
+    0.01,
+    100000,
+  );
   const renderer = new THREE.WebGLRenderer({ canvas, antialias: true });
   renderer.setPixelRatio(Math.min(window.devicePixelRatio || 1, 2));
 
@@ -262,7 +280,8 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
 
   function disposeObject(obj: THREE.Object3D) {
     obj.traverse((child) => {
-      const renderable = child as THREE.Mesh | THREE.Points | THREE.LineSegments;
+      const renderable = child as
+        THREE.Mesh | THREE.Points | THREE.LineSegments;
       renderable.geometry?.dispose();
       const material = renderable.material;
       if (Array.isArray(material)) material.forEach((item) => item.dispose());
@@ -278,6 +297,76 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
     }
   }
 
+  // ── 겹쳐 보기 메시 ─────────────────────────────────────────────────────────
+  // 도엽등고 서피스 위에 라이다 지표면을 겹쳐 두 지형을 눈으로 대조한다
+  // (2026-08-30 사용자 지시). 본 메시와 카메라·좌표계를 공유하므로 같은 자리에 겹친다.
+  let overlayMesh: THREE.Object3D | null = null;
+  let overlayGeneration = 0;
+
+  function clearOverlay() {
+    if (overlayMesh) {
+      scene.remove(overlayMesh);
+      disposeObject(overlayMesh);
+      overlayMesh = null;
+    }
+  }
+
+  async function loadOverlay(
+    projectId: string,
+    models: readonly SurfaceModelSummary[],
+    sourceFilter: string,
+    method: string,
+    smooth: boolean,
+  ): Promise<boolean> {
+    const generation = ++overlayGeneration;
+    clearOverlay();
+    const match = models.find((model) => {
+      const configured = model.generation_params?.source_filter;
+      return (
+        model.model_type.toLowerCase() === method.toLowerCase() &&
+        typeof configured === "string" &&
+        configured.toLowerCase() === sourceFilter.toLowerCase()
+      );
+    });
+    if (!match) return false;
+    const url = `${API_BASE_URL}/projects/${projectId}/surface/models/${match.id}/preview?smooth=${smooth}`;
+    try {
+      const buffer = await fetchCachedBytes(projectId, url);
+      if (generation !== overlayGeneration) return false;
+      return await new Promise<boolean>((resolve) => {
+        new GLTFLoader().parse(
+          buffer,
+          "",
+          (gltf) => {
+            if (generation !== overlayGeneration) {
+              disposeObject(gltf.scene);
+              resolve(false);
+              return;
+            }
+            // 겹친 두 면을 구분하려고 반투명 단색으로 덮어씌운다.
+            gltf.scene.traverse((child) => {
+              if (child instanceof THREE.Mesh) {
+                child.material = new THREE.MeshStandardMaterial({
+                  color: 0x60a5fa,
+                  transparent: true,
+                  opacity: 0.45,
+                  side: THREE.DoubleSide,
+                  flatShading: false,
+                });
+              }
+            });
+            overlayMesh = gltf.scene;
+            scene.add(gltf.scene);
+            resolve(true);
+          },
+          () => resolve(false),
+        );
+      });
+    } catch {
+      return false;
+    }
+  }
+
   function clearContours() {
     while (contourGroup.children.length > 0) {
       const child = contourGroup.children[0];
@@ -303,8 +392,11 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
 
   const fitCamera = (object: THREE.Object3D) => {
     const { span } = getFitParams(object);
-    const aspect = viewerArea.clientWidth / Math.max(viewerArea.clientHeight, 1);
-    const distance = referenceBounds ? getTopFitDistance(referenceBounds, aspect) : span * 1.2;
+    const aspect =
+      viewerArea.clientWidth / Math.max(viewerArea.clientHeight, 1);
+    const distance = referenceBounds
+      ? getTopFitDistance(referenceBounds, aspect)
+      : span * 1.2;
     controls.target.set(0, 0, 0);
     // 정확히 수직이면 lookAt이 화면 방향을 못 정해 첫 드래그에 화면이 뒤집힌다.
     camera.position.set(0, distance, distance * TOP_VIEW_TILT);
@@ -366,12 +458,15 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
     // model_type is TIN / DTM / NURBS / Implicit / Meshfree (we match activeMethod)
     // model_file_path contains the activeFilter (e.g. csf, pmf, grid_min_z)
     const match = currentModelsList.find((m) => {
-      const typeMatches = m.model_type.toLowerCase() === activeMethod.toLowerCase();
+      const typeMatches =
+        m.model_type.toLowerCase() === activeMethod.toLowerCase();
       const configuredFilter = m.generation_params?.source_filter;
       const filterMatches =
         (typeof configuredFilter === "string" &&
           configuredFilter.toLowerCase() === activeFilter.toLowerCase()) ||
-        Boolean(m.model_file_path?.toLowerCase().includes(activeFilter.toLowerCase()));
+        Boolean(
+          m.model_file_path?.toLowerCase().includes(activeFilter.toLowerCase()),
+        );
       return typeMatches && filterMatches;
     });
 
@@ -382,7 +477,8 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
     }
 
     const modelId = match.id;
-    const isSmooth = (activeMethod === "tin" || activeMethod === "dtm") && smoothingOn();
+    const isSmooth =
+      (activeMethod === "tin" || activeMethod === "dtm") && smoothingOn();
     currentModelId = modelId;
     currentModelSmooth = isSmooth;
     const generation = ++loadGeneration;
@@ -429,7 +525,8 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
             gltf.scene.traverse((child) => {
               if (child instanceof THREE.Mesh) {
                 child.material.side = THREE.DoubleSide;
-                child.material.vertexColors = child.geometry.hasAttribute("color");
+                child.material.vertexColors =
+                  child.geometry.hasAttribute("color");
               }
             });
             gltf.scene.visible = surfCheck.checked;
@@ -442,7 +539,8 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
           },
           () => {
             if (generation !== loadGeneration) return;
-            statusSpan.textContent = "3D 메쉬 파일이 없거나 로드할 수 없습니다.";
+            statusSpan.textContent =
+              "3D 메쉬 파일이 없거나 로드할 수 없습니다.";
             showProgress(null, null);
           },
         );
@@ -498,6 +596,11 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
       // 주곡선·보조곡선 두 덩어리로 합쳐 객체 2개만 만든다(2026-08-01).
       const majorPoints: THREE.Vector3[] = [];
       const minorPoints: THREE.Vector3[] = [];
+      const labelCandidates: {
+        level: number;
+        position: THREE.Vector3;
+        length: number;
+      }[] = [];
 
       data.contours.forEach((c: any) => {
         if (c.level < minH) minH = c.level;
@@ -512,47 +615,63 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
           bucket.push(points[i], points[i + 1]);
         }
 
+        // 라벨은 여기서 만들지 않고 후보만 모은다 — 등고선이 잘게 쪼개지면 조각마다
+        // 라벨이 붙어 수백 개가 되고, 매 프레임 위치 재계산이 화면을 멈춰 세운다
+        // (2026-08-30 사용자 보고). 아래에서 긴 것부터 상한만큼만 만든다.
         if (isMajor && points.length > 4) {
-          const labelPos = points[Math.floor(points.length / 2)];
-          const labelDiv = document.createElement("div");
-          labelDiv.className = "contour-label";
-          labelDiv.innerText = `${Math.round(c.level)}m`;
-          labelDiv.style.position = "absolute";
-          labelDiv.style.background = "rgba(255, 255, 255, 0.85)";
-          labelDiv.style.border = "1px solid #d97706";
-          labelDiv.style.color = "#b45309";
-          labelDiv.style.padding = "1px 4px";
-          labelDiv.style.borderRadius = "3px";
-          labelDiv.style.fontSize = "9px";
-          labelDiv.style.fontWeight = "bold";
-          labelDiv.style.pointerEvents = "none";
-          labelDiv.style.zIndex = "5";
-          labelDiv.style.transform = "translate(-50%, -50%)";
-
-          (labelDiv as any).__updateLabelPos = () => {
-            if (!contourCheck.checked) {
-              labelDiv.style.display = "none";
-              return;
-            }
-            const proj = labelPos.clone().project(camera);
-            const x = (proj.x * 0.5 + 0.5) * viewerArea.clientWidth;
-            const y = (-(proj.y * 0.5) + 0.5) * viewerArea.clientHeight;
-
-            if (proj.z > 1) {
-              labelDiv.style.display = "none";
-            } else {
-              labelDiv.style.display = "block";
-              labelDiv.style.left = `${x}px`;
-              labelDiv.style.top = `${y}px`;
-            }
-          };
-
-          viewerArea.appendChild(labelDiv);
-          labelElements.push(labelDiv);
-          labelsDirty = true;
+          let length = 0;
+          for (let i = 0; i < points.length - 1; i++) {
+            length += points[i].distanceTo(points[i + 1]);
+          }
+          labelCandidates.push({
+            level: c.level,
+            position: points[Math.floor(points.length / 2)],
+            length,
+          });
         }
       });
 
+      labelCandidates.sort((a, b) => b.length - a.length);
+      for (const candidate of labelCandidates.slice(0, MAX_CONTOUR_LABELS)) {
+        const labelPos = candidate.position;
+        const labelDiv = document.createElement("div");
+        labelDiv.className = "contour-label";
+        labelDiv.innerText = `${Math.round(candidate.level)}m`;
+        labelDiv.style.position = "absolute";
+        labelDiv.style.background = "rgba(255, 255, 255, 0.85)";
+        labelDiv.style.border = "1px solid #d97706";
+        labelDiv.style.color = "#b45309";
+        labelDiv.style.padding = "1px 4px";
+        labelDiv.style.borderRadius = "3px";
+        labelDiv.style.fontSize = "9px";
+        labelDiv.style.fontWeight = "bold";
+        labelDiv.style.pointerEvents = "none";
+        labelDiv.style.zIndex = "5";
+        labelDiv.style.transform = "translate(-50%, -50%)";
+
+        (labelDiv as any).__updateLabelPos = () => {
+          if (!contourCheck.checked) {
+            labelDiv.style.display = "none";
+            return;
+          }
+          const proj = labelPos.clone().project(camera);
+          const x = (proj.x * 0.5 + 0.5) * viewerArea.clientWidth;
+          const y = (-(proj.y * 0.5) + 0.5) * viewerArea.clientHeight;
+
+          if (proj.z > 1) {
+            labelDiv.style.display = "none";
+          } else {
+            labelDiv.style.display = "block";
+            labelDiv.style.left = `${x}px`;
+            labelDiv.style.top = `${y}px`;
+          }
+        };
+
+        viewerArea.appendChild(labelDiv);
+        labelElements.push(labelDiv);
+        labelsDirty = true;
+      }
+
       // 합쳐 둔 점들을 주곡선·보조곡선 각각 한 덩어리로 올린다.
       [
         { points: minorPoints, color: 0xf59e0b },
@@ -632,18 +751,26 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
     if (terrainMesh && terrainMesh.visible) {
       scaleBar.hidden = false;
       const dist = camera.position.distanceTo(controls.target);
-      const metersPerPixel = targetPlaneMetersPerPixel(dist, viewerArea.clientHeight);
+      const metersPerPixel = targetPlaneMetersPerPixel(
+        dist,
+        viewerArea.clientHeight,
+      );
       const roughMeters = 100 * metersPerPixel;
       const prettyMeters = niceScaleDistance(roughMeters);
       scaleBar.style.width = `${prettyMeters / metersPerPixel}px`;
       scaleLabel.textContent =
-        prettyMeters >= 1000 ? `${(prettyMeters / 1000).toFixed(0)} km` : `${prettyMeters} m`;
+        prettyMeters >= 1000
+          ? `${(prettyMeters / 1000).toFixed(0)} km`
+          : `${prettyMeters} m`;
     } else {
       scaleBar.hidden = true;
     }
 
     // 등고 라벨 위치 — 화면이 실제로 움직였을 때만 다시 계산한다(매 프레임 재계산은 낭비).
-    if (labelsDirty || !cameraMatrixSnapshot.equals(camera.matrixWorldInverse)) {
+    if (
+      labelsDirty ||
+      !cameraMatrixSnapshot.equals(camera.matrixWorldInverse)
+    ) {
       labelsDirty = false;
       cameraMatrixSnapshot.copy(camera.matrixWorldInverse);
       labelElements.forEach((label) => {
@@ -685,7 +812,8 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
   intervalForm.addEventListener("submit", async (e) => {
     e.preventDefault();
     const interval = Number(intervalInput.value);
-    if (!Number.isFinite(interval) || interval < 0.5 || currentModelId === null) return;
+    if (!Number.isFinite(interval) || interval < 0.5 || currentModelId === null)
+      return;
     intervalSubmit.disabled = true;
     await loadSelectedContours(currentModelId, currentModelSmooth, true);
     intervalSubmit.disabled = false;
@@ -724,6 +852,19 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
       activeMethod = method;
       syncSmoothingSupport();
     },
+    showOverlay(sourceFilter, method, smooth) {
+      if (!sourceFilter) {
+        clearOverlay();
+        return Promise.resolve(false);
+      }
+      return loadOverlay(
+        currentProjectId,
+        currentModelsList,
+        sourceFilter,
+        method,
+        smooth,
+      );
+    },
     applyCameraState,
     onCameraChange(listener) {
       cameraListener = listener;
@@ -739,7 +880,8 @@ export function createSurfaceTerrainViewer(): SurfaceTerrainViewer {
       syncSmoothingSupport();
     },
     setContourInterval(interval) {
-      if (Number.isFinite(interval) && interval > 0) intervalInput.value = String(interval);
+      if (Number.isFinite(interval) && interval > 0)
+        intervalInput.value = String(interval);
     },
     getContourInterval() {
       return Number.parseFloat(intervalInput.value);