"""도엽등고선 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,
    evaluate_nurbs_spline,
    fit_nurbs_spline,
    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_NURBS_CONTROL_POINTS_PER_AXIS,
    SURFACE_NURBS_DEGREE,
    SURFACE_NURBS_PATCH_SIZE_M,
)

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(B-spline) 곡면**으로 만든다.

    격자 삼각형을 그대로 쓰면 1m 셀 경계가 계단처럼 보인다. LAS 파이프라인이 쓰는
    것과 같은 곡면 적합(`fit_nurbs_spline`·`evaluate_nurbs_spline`)을 걸어 매끈한
    면을 얻는다(2026-08-30 사용자 지시). **표고 정본(npz)은 격자 그대로**이고 여기서
    바꾸는 것은 화면에 보이는 메시뿐이다 — 종·횡단은 격자를 샘플링한다.

    적합에 쓰는 제어 간격은 LAS NURBS와 같은 config 값(패치 크기 ÷ 축당 제어점 수)이며,
    제어 격자가 격자보다 성기지 않도록 최소 1셀로 묶는다.
    """
    stride = 1
    while (len(x) // stride + 1) * (len(y) // stride + 1) > SURFACE_MAX_PREVIEW_VERTICES:
        stride += 1
    px, py = x[::stride], y[::stride]
    pv = valid[::stride, ::stride]

    cell_m = float(x[1] - x[0]) if len(x) > 1 else SHEET_SURFACE_GRID_M
    control_m = max(
        SURFACE_NURBS_PATCH_SIZE_M / max(SURFACE_NURBS_CONTROL_POINTS_PER_AXIS - 1, 1), cell_m
    )
    control_step = max(1, int(round(control_m / max(cell_m, 1e-6))))
    # 곡면 적합은 결측을 못 받는다 — 무효 셀은 가장 가까운 유효 표고로 메우고,
    # 메시를 자를 때 원래 유효 마스크로 다시 도려낸다.
    filled = z.astype(np.float64)
    if not valid.all():
        from scipy.ndimage import distance_transform_edt

        _, (rows, cols) = distance_transform_edt(~valid, return_indices=True)
        filled = filled[rows, cols]
    cx, cy = x[::control_step], y[::control_step]
    cz = filled[::control_step, ::control_step]
    try:
        spline, z_range = fit_nurbs_spline(cx, cy, cz, SURFACE_NURBS_DEGREE)
        pz = evaluate_nurbs_spline(spline, py, px, z_range, "도엽 서피스")
    except Exception as exc:  # noqa: BLE001 — 적합 실패 시 격자 메시로 되돌린다
        logger.warning("도엽 서피스: NURBS 곡면 적합 실패(%s) — 격자 메시를 씁니다.", exc)
        pz = filled[::stride, ::stride]

    vertices = grid_vertices(px, py, np.asarray(pz, dtype=np.float64))
    faces = grid_faces(len(py), len(px))
    return clip_and_compact_mesh(vertices, faces, pv.reshape(-1))


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)
    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,
    }