Aislo/D03_PlanAnalysis/D03_PlanAnalysis_Terrain.py

"""D 예정노선 지형 받기 — 지점 기준 수치지형도 · 위성 · 연도 구간별 DEM (요구 345 · PLAN 530-1).

지형은 D 프로젝트 storage 안 `B04_PreProcess/` 배치 그대로(지형 코드 주인 = B04 · 고칠 곳 0):

- 도엽 = 모든 지점 + 구간 표본(`SAMPLE_STEP_M`) 각 3×3 → `ensure_sheets`(없는 것만 · prune 없음)
- 병합 = 지난번 병합 도엽 ∪ 이번 도엽(옛 도엽 유지 · 열쇠 같으면 건너뜀)
- 위성 = 기준 도엽 도곽 ∪ 옛 위성 범위(덮으면 건너뜀)
- DEM = 연도 구간마다 `models_<연도>/`(긴 노선 메모리 · 비용면 한도) · 범위 = 구간 상자 + 탐색 복도
- 작업 좌표계 = 지점 epsg · `B03_FileInput/input/prj/plan.prj` 한 장(B04 · B05 좌표계 창구가 읽음)

계산 자리 = ③ 서버 단독. DB 행(`surface_models`)은 라우터가 씀.
"""

from __future__ import annotations

import json
import logging
import math
from pathlib import Path
from typing import Any, Callable

import numpy as np

from B04_PreProcess.B04_PreProcess_Engine_MapSheet import (
    neighbors_for_points,
    sheet5k_to_bounds,
    sheets_for_points,
)
from B04_PreProcess.B04_PreProcess_Engine_SheetParser import parse_and_merge_sheets
from B04_PreProcess.B04_PreProcess_Engine_SheetStore import (
    ensure_sheets,
    get_project_map_sheets_dir,
    get_sheet_path,
)
from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import build_sheet_surface_model
from B04_PreProcess.B04_PreProcess_Engine_VWorld import download_vworld_satellite_map
from common_util.common_util_json import atomic_write_json
from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD, SHEET_SURFACE_MARGIN_M

logger = logging.getLogger(__name__)

# 구간 표본 간격 — 도엽(≈2.8 × 2.3 km) 3×3 묶음이 끊기지 않게
SAMPLE_STEP_M = 1000.0
# 탐색 복도 반폭(D 만 · B 기본 100 m) — 실측(530-2): 노선이 복도 끝에 닿지 않는 폭
PLAN_CORRIDOR_M = 300.0
# DEM 이 복도 밖으로 더 덮는 폭(복도 끝 칸 기울기 · 횡단 표본 ±20 m)
DEM_EXTRA_M = 50.0
STAGE = Path("B04_PreProcess")
PRJ_NAME = "plan.prj"

Progress = Callable[[int, str], None]


def models_dir_of(project_root: Path, year: int) -> Path:
    """연도 구간 DEM 자리."""
    return project_root / STAGE / f"models_{int(year)}"


def segment_samples(
    start: tuple[float, float], end: tuple[float, float]
) -> list[tuple[float, float]]:
    """구간 위 표본(끝점 둘 포함 · `SAMPLE_STEP_M` 간격 이하)."""
    length = math.hypot(end[0] - start[0], end[1] - start[1])
    count = max(1, math.ceil(length / SAMPLE_STEP_M))
    return [
        (start[0] + (end[0] - start[0]) * i / count, start[1] + (end[1] - start[1]) * i / count)
        for i in range(count + 1)
    ]


def reference_points(segments: list[dict[str, Any]], epsg: int) -> list[tuple[float, float]]:
    """도엽 기준 좌표(위도, 경도) — 모든 지점 + 구간 표본."""
    from pyproj import Transformer

    to_wgs = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True)
    points: list[tuple[float, float]] = []
    for seg in segments:
        for x, y in segment_samples(seg["start"], seg["end"]):
            lon, lat = to_wgs.transform(x, y)
            points.append((lat, lon))
    return points


def write_plan_prj(project_root: Path, epsg: int) -> Path:
    """작업 좌표계 .prj — 내용이 같으면 그대로."""
    from pyproj import CRS

    path = project_root / "B03_FileInput" / "input" / "prj" / PRJ_NAME
    text = CRS.from_epsg(int(epsg)).to_wkt("WKT1_ESRI")
    if not path.is_file() or path.read_text(encoding="utf-8", errors="ignore") != text:
        path.parent.mkdir(parents=True, exist_ok=True)
        path.write_text(text, encoding="utf-8")
    return path


def dem_outline(start: tuple[float, float], end: tuple[float, float]) -> np.ndarray:
    """DEM 상자 꼭짓점 둘 — B 여유(`SHEET_SURFACE_MARGIN_M`)를 더하면 복도 + `DEM_EXTRA_M`."""
    pad = max(0.0, PLAN_CORRIDOR_M + DEM_EXTRA_M - SHEET_SURFACE_MARGIN_M)
    xs, ys = (start[0], end[0]), (start[1], end[1])
    return np.array(
        [[min(xs) - pad, min(ys) - pad], [max(xs) + pad, max(ys) + pad]], dtype=np.float64
    )


def _read_json(path: Path) -> dict[str, Any] | None:
    try:
        return json.loads(path.read_text(encoding="utf-8"))
    except (OSError, ValueError):
        return None


def _merge_sheets(processed: Path, store: Path, available: list[str]) -> list[str]:
    """병합 = 지난번 병합 도엽(zip 남은 것) ∪ 이번 도엽 · 열쇠 같으면 건너뜀."""
    key_path = processed / "sheet_merge.key.json"
    old = (_read_json(key_path) or {}).get("sheets") or []
    sheets = sorted({*available, *(s for s in old if get_sheet_path(store, s) is not None)})
    key = {"sheets": sheets}
    if sheets and (_read_json(key_path) != key or not any(processed.glob("도엽_*.geojson"))):
        parse_and_merge_sheets(store, sheets, processed)
        atomic_write_json(key_path, key)
    return sheets


def _satellite(processed: Path, prj: Path, points: list[tuple[float, float]], epsg: int) -> str:
    """위성 = 기준 도엽 도곽 ∪ 옛 범위 · 이미 덮으면 건너뜀. 돌려줌 = skip · done · failed."""
    from pyproj import Transformer

    to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True)
    xs: list[float] = []
    ys: list[float] = []
    for sheet in sheets_for_points(points):
        lon0, lat0, lon1, lat1 = sheet5k_to_bounds(sheet)
        for lon, lat in ((lon0, lat0), (lon0, lat1), (lon1, lat0), (lon1, lat1)):
            x, y = to_metric.transform(lon, lat)
            xs.append(x)
            ys.append(y)
    need = {"x": [min(xs), max(xs)], "y": [min(ys), max(ys)]}
    meta = _read_json(processed / "vworld_satellite_meta.json")
    try:
        old = [float(meta[k]) for k in ("x_min", "x_max", "y_min", "y_max")] if meta else None
    except (KeyError, TypeError, ValueError):
        old = None
    if old and old[0] <= need["x"][0] and old[1] >= need["x"][1]:
        if old[2] <= need["y"][0] and old[3] >= need["y"][1]:
            return "skip"
    if old:
        need = {
            "x": [min(old[0], need["x"][0]), max(old[1], need["x"][1])],
            "y": [min(old[2], need["y"][0]), max(old[3], need["y"][1])],
        }
    try:
        done = download_vworld_satellite_map(prj, need, processed, layer_name="Satellite")
    except Exception as exc:  # noqa: BLE001 — 위성은 배경 · 지형 받기를 막지 않음
        logger.warning("D 위성 지도 받기 실패: %s", exc)
        return "failed"
    return "done" if done else "failed"


def ensure_plan_terrain(
    project_root: Path,
    segments: list[dict[str, Any]],
    epsg: int,
    on_progress: Progress | None = None,
) -> dict[str, Any]:
    """지점 기준 지형 받기 — `segments` = [{year, start: (x, y), end: (x, y)}](작업 좌표계 m).

    돌려줌 = 도엽(받은 · 실패 · 병합) · 위성 상태 · 연도별 모델 등록용 dict 목록.
    """

    def report(percent: int, message: str) -> None:
        if on_progress is not None:
            on_progress(percent, message)

    if not segments:
        raise ValueError("연도 구간이 없습니다.")
    processed = project_root / STAGE / "processed"
    processed.mkdir(parents=True, exist_ok=True)
    prj = write_plan_prj(project_root, epsg)
    crs = f"EPSG:{int(epsg)}"

    report(5, "수치지형도 도엽 고르는 중")
    points = reference_points(segments, epsg)
    wanted = neighbors_for_points(points)
    store = get_project_map_sheets_dir(project_root)
    before = {s for s in wanted if get_sheet_path(store, s) is not None}
    report(10, f"수치지형도 도엽 {len(wanted) - len(before)}매 받는 중")
    got = ensure_sheets(store, wanted)
    if not got["available"]:
        raise ValueError(f"수치지형도 도엽을 하나도 받지 못했습니다: {got['failed']}")
    report(45, "도엽 레이어 병합 중")
    merged = _merge_sheets(processed, store, list(got["available"]))
    report(60, "위성 지도 받는 중")
    satellite = _satellite(processed, prj, points, epsg)

    years: dict[int, list[dict[str, Any]]] = {}
    for i, seg in enumerate(segments):
        year = int(seg["year"])
        report(65 + 30 * i // len(segments), f"{year}년 구간 지형면 만드는 중")
        models_dir = models_dir_of(project_root, year)
        models_dir.mkdir(parents=True, exist_ok=True)
        years[year] = build_sheet_surface_model(
            project_root,
            processed,
            models_dir,
            dem_outline(seg["start"], seg["end"]),
            crs,
            [SHEET_SURFACE_DEFAULT_METHOD],
        )
        if not years[year]:
            raise ValueError(f"{year}년 구간 지형면을 만들지 못했습니다 — 도엽 등고선 확인")
    report(100, "지형 받기 끝")
    return {
        "sheets": {
            "wanted": len(wanted),
            "downloaded": sorted(set(got["available"]) - before),
            "failed": got["failed"],
            "merged": merged,
        },
        "satellite": satellite,
        "years": years,
    }