Aislo/resources/tester/test_sheet_surface.py

"""도엽등고선 3D 서피스(B04_PreProcess_Engine_SheetSurface) 단위 검증.

합성 등고선(정사각 링 3단, EPSG:5186 → WGS84 역투영)으로 dtm_sheet.npz를 만들고,
종·횡단이 쓰는 build_surface_sampler("sheet", "dtm")로 왕복 조회가 되는지 본다.
"""

import json
import sys
from pathlib import Path

import numpy as np
import pytest
from pyproj import Transformer

PROJECT_ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(PROJECT_ROOT))

from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import (  # noqa: E402
    build_sheet_surface_model,
)
from common_util.common_util_surface_sampler import build_surface_sampler  # noqa: E402

EPSG = 5186
# 사업지 부근 임의 원점 (EPSG:5186 중부원점 좌표계, m)
X0, Y0 = 200000.0, 500000.0


def _square_ring(cx: float, cy: float, half: float) -> list[list[float]]:
    """정사각 폐합 링 (반시계)."""
    return [
        [cx - half, cy - half],
        [cx + half, cy - half],
        [cx + half, cy + half],
        [cx - half, cy + half],
        [cx - half, cy - half],
    ]


@pytest.fixture()
def sheet_project(tmp_path: Path) -> dict:
    """processed/도엽_등고선.geojson(WGS84)을 갖춘 가짜 프로젝트 저장소."""
    processed_dir = tmp_path / "B04_PreProcess" / "processed"
    models_dir = tmp_path / "B04_PreProcess" / "models"
    processed_dir.mkdir(parents=True)
    models_dir.mkdir(parents=True)

    to_wgs84 = Transformer.from_crs(f"EPSG:{EPSG}", "EPSG:4326", always_xy=True)
    features = []
    # 바깥 100m → 안쪽 120m: 중심으로 갈수록 높은 언덕
    for half, elev in ((450.0, 100.0), (300.0, 110.0), (150.0, 120.0)):
        ring = [list(to_wgs84.transform(x, y)) for x, y in _square_ring(X0, Y0, half)]
        features.append(
            {
                "type": "Feature",
                "properties": {"등고수치": elev},
                "geometry": {"type": "LineString", "coordinates": ring},
            }
        )
    (processed_dir / "도엽_등고선.geojson").write_text(
        json.dumps({"type": "FeatureCollection", "features": features}, ensure_ascii=False),
        encoding="utf-8",
    )
    return {"root": tmp_path, "processed": processed_dir, "models": models_dir}


def test_sheet_surface_npz_and_sampler_roundtrip(sheet_project: dict) -> None:
    # 노선: 중심을 지나는 100m 직선
    route_xy = np.array([[X0 - 50.0, Y0], [X0 + 50.0, Y0]], dtype=np.float64)
    models = build_sheet_surface_model(
        sheet_project["root"],
        sheet_project["processed"],
        sheet_project["models"],
        route_xy,
        EPSG,
        methods=["laplace"],
    )
    assert len(models) == 1
    model = models[0]
    assert model["model_type"] == "dtm"
    assert model["generation_params"]["source_filter"] == "sheet_laplace"

    npz_path = sheet_project["models"] / "dtm_sheet_laplace.npz"
    assert npz_path.is_file()
    with np.load(npz_path) as data:
        for key in ("x", "y", "z", "valid_mask"):
            assert key in data, key
        x, y = data["x"], data["y"]
        # DtmGridSampler 규약: 오름차순 축, z=(len(y), len(x))
        assert x[0] < x[-1] and y[0] < y[-1]
        assert data["z"].shape == (len(y), len(x))

    # 프리뷰 glb
    assert (sheet_project["models"] / "dtm_sheet_laplace_preview.glb").stat().st_size > 0

    # 종·횡단 공용 sampler로 왕복 — 등고선 안쪽은 유효 표고.
    # 최고 등고선(120) 안쪽은 마루 캡(경사 연장)으로 최대 +간격-0.5m까지 올라갈 수 있다.
    sampler = build_surface_sampler(sheet_project["models"], "sheet_laplace", "dtm", smooth=False)
    inner = np.array([[X0, Y0 + 200.0], [X0 + 100.0, Y0 - 100.0]], dtype=np.float64)
    z, valid = sampler.sample_xy(inner)
    assert valid.all()
    assert np.all((z >= 100.0 - 1e-6) & (z <= 124.5 + 1e-3))

    # 바깥(등고선 밖 + margin 밖)은 무효
    outer = np.array([[X0 + 10000.0, Y0]], dtype=np.float64)
    _, outer_valid = sampler.sample_xy(outer)
    assert not outer_valid.any()


def test_sheet_surface_without_contours_returns_none(tmp_path: Path) -> None:
    processed_dir = tmp_path / "B04_PreProcess" / "processed"
    models_dir = tmp_path / "B04_PreProcess" / "models"
    processed_dir.mkdir(parents=True)
    models_dir.mkdir(parents=True)
    route_xy = np.array([[0.0, 0.0], [10.0, 0.0]])
    assert build_sheet_surface_model(tmp_path, processed_dir, models_dir, route_xy, EPSG) == []


def test_interpolate_between_contours_is_continuous() -> None:
    """등고선 사이 z가 거리 비례로 연속 변화 — 계단(같은 표고 평탄면)이 없다."""
    from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import build_distance

    burned = np.full((10, 101), np.nan, dtype=np.float32)
    burned[:, 0] = 100.0  # 좌측 벽 100m
    burned[:, 100] = 105.0  # 우측 벽 105m
    z = build_distance(None, burned, None, 1.0)
    middle = z[5]
    assert abs(float(middle[0]) - 100.0) < 1e-3
    assert abs(float(middle[100]) - 105.0) < 1e-3
    assert abs(float(middle[50]) - 102.5) < 0.05  # 중앙 = 정확히 절반
    # 단조 증가이고 계단(연속 동일값)이 없다
    steps = np.diff(middle)
    assert np.all(steps > 0)


def test_resolve_enclosed_interiors_raises_summit_lowers_pit() -> None:
    """폐합 등고선 안쪽 — 바깥이 낮으면 마루로 올리고, 높으면 웅덩이로 내린다."""
    from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import _resolve_enclosed_interiors

    burned = np.full((60, 60), np.nan, dtype=np.float32)
    burned[10:50, 10] = 100.0
    burned[10:50, 49] = 100.0
    burned[10, 10:50] = 100.0
    burned[49, 10:50] = 100.0  # 100m 폐합 사각 링
    surface = np.full((60, 60), 95.0, dtype=np.float32)  # 바깥이 낮은 지형
    surface[11:49, 11:49] = 97.0  # 보간이 만든 분화구
    handled = _resolve_enclosed_interiors(burned, [100.0], surface, 1.0, 5.0)
    assert handled.any()
    center = float(surface[29:31, 29:31].max())
    assert 100.0 < center <= 100.0 + 4.5 + 1e-3  # 상한 +간격-0.5m

    # 바깥이 높으면 웅덩이 — 안쪽이 내려간다
    surface2 = np.full((60, 60), 110.0, dtype=np.float32)
    surface2[11:49, 11:49] = 103.0
    assert _resolve_enclosed_interiors(burned, [100.0], surface2, 1.0, 5.0).any()
    assert float(surface2[30, 30]) < 100.0


def test_cone_reproduction_beats_over_relaxation() -> None:
    """원뿔(동심원 등고선) 재현 — 거리 보간이 정답에 가깝고, 과도한 완화는 악화된다.

    z=r 형상은 harmonic이 아니라 biharmonic이라, 라플라스 완화를 수렴시키면 마루가
    눌린다. ANUDEM이 thin plate spline을 쓰는 이유(Hutchinson 1988/89)를 지킨다.
    """
    from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import build_distance, relax_laplace

    n = 301
    centre = n // 2
    rows, cols = np.indices((n, n))
    radius = np.hypot(rows - centre, cols - centre)
    truth = 140.0 - radius * 0.1  # 반경 1m마다 0.1m 하강 = 1m 등고선 간격 10m
    burned = np.full((n, n), np.nan, dtype=np.float32)
    for level in (110.0, 115.0, 120.0, 125.0, 130.0, 135.0):
        burned[np.abs(radius - (140.0 - level) / 0.1) < 0.5] = level
    band = (radius > 55) & (radius < 295)

    base = build_distance(None, burned, None, 1.0)
    base_error = float(np.abs(base - truth)[band].mean())
    assert base_error < 0.5  # 거리 보간은 원뿔을 그대로 재현한다

    light = base.copy()
    relax_laplace(light, np.isfinite(burned), 5)
    assert float(np.abs(light - truth)[band].mean()) <= base_error + 0.02

    heavy = base.copy()
    relax_laplace(heavy, np.isfinite(burned), 200)
    # 수렴시키면 harmonic 해로 끌려가 오차가 눈에 띄게 커진다
    assert float(np.abs(heavy - truth)[band].mean()) > base_error