"""도엽등고선 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). 등고선→정점 구름→Delaunay TIN 보간은 배수유역
엔진(`_Watershed_Grid`)의 검증된 경로를 그대로 쓴다.
"""

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_Watershed_Grid import (
    build_contour_cloud,
    grid_spec_from_bounds,
    interpolate_elevation,
)
from config.config_system import (
    SHEET_SURFACE_GRID_M,
    SHEET_SURFACE_MARGIN_M,
    SHEET_SURFACE_MIDLINE_MAX_M,
    SURFACE_MAX_PREVIEW_VERTICES,
)

logger = logging.getLogger(__name__)

# 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값)
_CONTOUR_FILE = "도엽_등고선.geojson"

# 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다.
SHEET_SOURCE_FILTER = "sheet"


def _load_contour_features_metric(processed_dir: Path, epsg: int) -> list[dict[str, Any]]:
    """병합 도엽 등고선(WGS84)을 읽어 사업지 CRS(m)로 재투영한다."""
    path = processed_dir / _CONTOUR_FILE
    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]:
    """프리뷰용 정점·면 — 정점 수가 상한을 넘으면 격자를 성기게 딴다."""
    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 _densify_between_contours(spec: Any, features: list[dict[str, Any]], cloud: Any) -> Any:
    """인접 표고 등고선 쌍의 등거리 중간선을 중간 표고 정점으로 추가한다.

    등고선 정점만으로 Delaunay TIN을 치면 같은 표고 정점 3개짜리 평탄 삼각형이
    굴곡부·능선에서 계단(terrace)을 만든다. 두 등고선의 "사이점"(등거리선)에
    평균 표고 정점을 넣어 삼각형이 반드시 서로 다른 표고를 잇게 한다
    (2026-08-30 사용자 지적). 최고 등고선 안쪽 봉우리·같은 표고 사이 골짜기
    바닥은 짝이 없어 평탄으로 남는다 — 원천 데이터 한계.
    """
    from scipy.ndimage import distance_transform_edt

    from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours
    from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ContourCloud

    burned, levels = rasterize_contours(spec, features, None)
    present = [level for level in sorted(levels) if bool((burned == level).any())]
    if len(present) < 2:
        return cloud

    def _distance(level: float) -> np.ndarray:
        return distance_transform_edt(burned != level, sampling=spec.cell_m).astype(np.float32)

    # 1차: 셀별 최근접 등고선 레벨 — 등거리 조건은 멀리 있는 다른 표고 쌍에도 우연히
    # 성립하므로(예: 525~530 골짜기에 555/560 중간점이 박힘), "그 지점의 최근접
    # 등고선이 바로 그 쌍"일 때만 사이점으로 인정한다.
    best_distance = np.full(burned.shape, np.inf, dtype=np.float32)
    best_index = np.full(burned.shape, -1, dtype=np.int16)
    for index, level in enumerate(present):
        distance = _distance(level)
        take = distance < best_distance
        best_distance[take] = distance[take]
        best_index[take] = index

    xs = spec.cell_centers_x()
    ys = spec.cell_centers_y()
    extra_xy: list[np.ndarray] = []
    extra_z: list[np.ndarray] = []
    previous_distance = _distance(present[0])
    for index in range(1, len(present)):
        distance = _distance(present[index])
        near = (previous_distance < SHEET_SURFACE_MIDLINE_MAX_M) & (
            distance < SHEET_SURFACE_MIDLINE_MAX_M
        )
        pair_is_nearest = (best_index == index - 1) | (best_index == index)
        midline = near & pair_is_nearest & (np.abs(previous_distance - distance) <= spec.cell_m)
        rows, cols = np.nonzero(midline)
        if len(rows):
            # 등고선 정점 재샘플 밀도에 맞춰 성기게 담는다 — TIN 비용 절제.
            rows, cols = rows[::5], cols[::5]
            extra_xy.append(np.column_stack([xs[cols], ys[rows]]))
            extra_z.append(np.full(len(rows), (present[index - 1] + present[index]) / 2.0))
        previous_distance = distance
    if not extra_xy:
        return cloud
    added = int(sum(len(a) for a in extra_z))
    logger.info("도엽 서피스: 중간 보간선 정점 %d개 추가 (등고 %d단)", added, len(present))
    return ContourCloud(
        xy=np.vstack([cloud.xy, *extra_xy]),
        z=np.concatenate([cloud.z, *extra_z]),
    )


def _cap_flat_summits(surface: np.ndarray, cell_m: float, interval_m: float) -> int:
    """마지막 등고선 안쪽 평탄 마루를 완만한 능선 돔으로 올린다 (in-place).

    다음 등고선이 없는 봉우리·능선 마루는 TIN이 마지막 등고 표고로 평탄하게 채운다
    ('ㅜ'가 갑자기 'ㅡ'로 변하는 단 — 2026-08-30 사용자 지적). 실제 마루는 그 표고와
    +간격 사이이므로, 경계에서 가장 먼 골격선이 +간격/2가 되도록 거리 비례로 올린다
    (기대값 캡, 오차 상한 ±간격/2 = 원천 데이터 불확실성 이내).
    적용 조건: 등표고 평탄 컴포넌트이고 바깥 유효 이웃이 전부 더 낮을 때(=마루).
    골짜기 바닥은 배수 방향을 모르므로 건드리지 않는다. 돋운 컴포넌트 수를 반환.
    """
    from scipy.ndimage import binary_dilation, distance_transform_edt, label

    finite = np.isfinite(surface)
    quantized = np.round(surface * 100.0)
    capped = 0
    # 평탄 후보: 상하좌우 이웃 중 같은 표고가 있는 셀만 모아 컴포넌트를 짠다.
    same_right = np.zeros_like(finite)
    same_right[:, :-1] = finite[:, :-1] & finite[:, 1:] & (quantized[:, :-1] == quantized[:, 1:])
    same_down = np.zeros_like(finite)
    same_down[:-1, :] = finite[:-1, :] & finite[1:, :] & (quantized[:-1, :] == quantized[1:, :])
    flat = same_right | same_down
    flat[:, 1:] |= same_right[:, :-1]
    flat[1:, :] |= same_down[:-1, :]

    # 서로 다른 표고의 평탄면이 맞닿아 있을 수 있으므로 표고값별로 컴포넌트를 짠다.
    # 값별 셀 수 25 미만은 컴포넌트도 25 미만이므로 라벨링 전에 거른다(속도).
    flat_values, flat_counts = np.unique(quantized[flat], return_counts=True)
    for value in flat_values[flat_counts >= 25]:
        value_labels, value_count = label(flat & (quantized == value))
        sizes = np.bincount(value_labels.reshape(-1))
        for component_id in range(1, value_count + 1):
            if sizes[component_id] < 25:  # 25㎡ 미만은 시각적으로 단이 아니다
                continue
            rows, cols = np.nonzero(value_labels == component_id)
            if (
                rows.min() == 0
                or cols.min() == 0
                or rows.max() == surface.shape[0] - 1
                or cols.max() == surface.shape[1] - 1
            ):
                continue  # 격자 가장자리에 닿으면 바깥을 모른다
            # 이후 연산은 컴포넌트 bbox(+1) 창에서만 — 전체 격자 반복을 피한다.
            window = (
                slice(rows.min() - 1, rows.max() + 2),
                slice(cols.min() - 1, cols.max() + 2),
            )
            component = value_labels[window] == component_id
            ring = binary_dilation(component) & ~component & finite[window]
            if not ring.any():
                continue
            level = float(surface[rows[0], cols[0]])
            if float(surface[window][ring].max()) >= level - 1e-3:
                continue  # 더 높은 이웃이 있으면 마루가 아니다(사면 벤치·골짜기)
            inner = distance_transform_edt(component, sampling=cell_m)
            peak = float(inner.max())
            if peak <= 0.0:
                continue
            patch = surface[window]
            patch[component] += (interval_m / 2.0) * (inner[component] / peak).astype(surface.dtype)
            capped += 1
    if capped:
        logger.info(
            "도엽 서피스: 평탄 마루 %d곳을 +%.1fm 골격 캡으로 돋움", capped, interval_m / 2.0
        )
    return capped


def build_sheet_surface_model(
    project_root: Path,
    processed_dir: Path,
    models_dir: Path,
    route_xy: np.ndarray,
    epsg: int,
) -> dict[str, Any] | None:
    """도엽등고선으로 DTM npz·프리뷰 glb를 만들고 surface_models 등록용 dict를 돌려준다.

    `route_xy`: (N, 2) 노선 정점 XY(사업지 CRS, m). 실패하면 None — 호출측은
    분석을 계속한다(도엽 미확보 지역 폴백).
    """
    started = time.monotonic()
    features = _load_contour_features_metric(processed_dir, epsg)
    if not features:
        return None

    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

    # 저지대 제거(floor) 없이 전부 쓴다 — 종·횡단은 낮은 지반도 필요하다.
    cloud = build_contour_cloud(features, None, (x_min, y_min, x_max, y_max))
    if cloud.is_empty:
        logger.warning("도엽 서피스: 절취 범위 안에 등고선 정점이 없습니다.")
        return None

    spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M)
    # 등고 간격(m) — 마루 캡 크기의 근거. 레벨이 하나뿐이면 5m(1:5,000 주곡선) 폴백.
    contour_levels = np.unique(cloud.z)
    interval_m = float(np.diff(contour_levels).min()) if len(contour_levels) > 1 else 5.0
    # 계단(평탄 삼각형) 해소 — 인접 등고선 사이 중간 보간선을 정점으로 추가.
    cloud = _densify_between_contours(spec, features, cloud)
    surface = interpolate_elevation(spec, cloud)  # (R, C), 북→남 행 순서, 외부 NaN
    # 마지막 등고선 안쪽 평탄 마루를 완만한 돔으로 (2026-08-30 사용자 확정).
    _cap_flat_summits(surface, spec.cell_m, interval_m)

    # 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("도엽 서피스: 유효 표고 셀이 없습니다.")
        return None

    stem = f"dtm_{SHEET_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)

    logger.info(
        "도엽 서피스 생성 완료: %d×%d 격자, 정점 %d개 (%.1fs)",
        spec.n_rows,
        spec.n_cols,
        cloud.xy.shape[0],
        time.monotonic() - started,
    )
    return {
        "model_type": "dtm",
        "source_filter": SHEET_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": SHEET_SOURCE_FILTER,
            "representation": "regular_grid",
            "source": "map_sheet_contours",
            "margin_m": SHEET_SURFACE_MARGIN_M,
        },
        "layers": [
            {
                "layer_name": f"dtm_{SHEET_SOURCE_FILTER}_preview",
                "geometry_type": "MESH",
                "file_path": str(preview_path.relative_to(project_root)).replace("\\", "/"),
                "file_format": "glb",
            }
        ],
    }


def build_sheet_surface_from_route(
    project_root: Path, processed_dir: Path, models_dir: Path
) -> dict[str, Any] | None:
    """B03 업로드 계획노선 CSV를 찾아 도엽 서피스를 만든다. 없거나 실패하면 None."""
    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 None
    planned = read_planned_route_csv(route_file)
    if planned is None or len(planned.vertices) < 2:
        logger.warning("도엽 서피스: 계획 노선 파일을 읽지 못했습니다: %s", route_file.name)
        return None
    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 서피스 생성 중")
    model = build_sheet_surface_model(project_root, processed_dir, models_dir, route_xy, epsg)
    if model is None:
        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": [model],
    }