Aislo/common_util/common_util_crs.py

# common_util_crs.py
# PRJ(WKT) 좌표계 판별 — 정본은 WKT, EPSG는 라벨 (2026-08-31 확정).
#
# 배경: 원청 자료의 PRJ는 EPSG AUTHORITY가 없거나(UTM-K ESRI WKT), TOWGS84 7모수가
# EPSG DB 정의와 달라 `pyproj.CRS.to_epsg()`가 None을 주는 경우가 실측됐다.
# 그래서 ① 변환기는 파일의 WKT 그대로 만들고(지역 보정 보존), ② EPSG 코드는
# 메타데이터·로그용 라벨로만 판별한다. 라벨 실패는 차단 사유가 아니다.

from __future__ import annotations

import logging
import re
from pathlib import Path

from pyproj import CRS

logger = logging.getLogger(__name__)

# 한국에서 실제로 들어오는 투영 좌표계 후보 — 파라미터 지문 대조(사다리 ③)에만 쓴다.
# 사다리 ①(pyproj DB 대조)은 이 목록과 무관하게 세계 모든 EPSG를 판별한다.
KOREA_CANDIDATE_EPSG: tuple[int, ...] = (
    5185,
    5186,
    5187,
    5188,  # Korea 2000 서·중부·동부·동해 (GRS80, FN 600k)
    5179,
    5178,  # UTM-K (GRS80 / Bessel)
    5173,
    5174,
    5175,
    5176,
    5177,  # Korean 1985 구좌표계 (Bessel, FN 500k)
    32651,
    32652,  # WGS84 UTM 51N·52N
    4326,
)

_AUTHORITY_PATTERN = re.compile(r'AUTHORITY\["EPSG",\s*"(\d+)"\]')


def strip_bound(crs: CRS) -> CRS:
    """TOWGS84가 붙어 BoundCRS로 감싸진 경우 원래 좌표계를 꺼낸다.

    감싸진 채로는 `coordinate_operation`이 Helmert 7모수를 가리켜 투영
    파라미터가 안 보이고, `is_projected`도 False라 수평 성분 탐색이 실패한다.
    """
    if crs.is_bound and crs.source_crs is not None:
        return crs.source_crs
    return crs


def horizontal_crs(crs: CRS) -> CRS:
    """COMPD_CS(수직 CRS 결합)·BoundCRS를 벗겨 수평 성분만 꺼낸다."""
    crs = strip_bound(crs)
    if crs.is_compound:
        for component in crs.sub_crs_list:
            component = strip_bound(component)
            if component.is_projected or component.is_geographic:
                return component
    return crs


def _projection_fingerprint(crs: CRS) -> tuple | None:
    """투영 지문 — (중앙자오선, 원점위도, 축척, FE, FN, 타원체 장반경).

    TOWGS84는 뺀다: 같은 좌표계라도 파일마다 지역 보정 모수가 다를 수 있고,
    여기서 만드는 것은 변환기가 아니라 라벨이다.
    """
    crs = horizontal_crs(crs)
    if crs.is_geographic:
        return ("geographic", round(crs.ellipsoid.semi_major_metre, 3))
    operation = crs.coordinate_operation
    if operation is None:
        return None
    params = {p.name.lower(): p.value for p in operation.params}

    def pick(*needles: str) -> float | None:
        for needle in needles:
            for name, value in params.items():
                if needle in name:
                    return value
        return None

    return (
        round(pick("longitude of natural origin", "central meridian") or 0.0, 6),
        round(pick("latitude of natural origin") or 0.0, 6),
        round(pick("scale factor") or 1.0, 6),
        round(pick("false easting") or 0.0, 3),
        round(pick("false northing") or 0.0, 3),
        round(crs.ellipsoid.semi_major_metre, 3),
    )


def identify_epsg(crs: CRS, wkt_text: str | None = None) -> int | None:
    """EPSG 라벨 판별 사다리. 못 고르면 None — 기본값으로 때우지 않는다.

    ① pyproj DB 대조(전 세계 EPSG) → ② 원문 WKT 마지막 AUTHORITY 태그
    → ③ 투영 파라미터 지문을 한국 후보와 대조.
    """
    flat = horizontal_crs(crs)
    # WKT1 재직렬화는 축 순서 등 부차 정보가 빠져 기본 신뢰도(70)로는 DB 대조가
    # 실패한다. 라벨 용도이므로 느슨한 대조(25)를 쓴다 — always_xy 변환에는 무해.
    epsg = flat.to_epsg(min_confidence=25)
    if epsg is not None:
        return epsg

    if wkt_text:
        codes = _AUTHORITY_PATTERN.findall(wkt_text)
        # WKT1은 마지막 AUTHORITY가 PROJCS 자신의 코드다 (앞의 것들은 타원체·측지계 등).
        for code_text in reversed(codes):
            code = int(code_text)
            if code in KOREA_CANDIDATE_EPSG:
                return code

    fingerprint = _projection_fingerprint(crs)
    if fingerprint is None:
        return None
    for code in KOREA_CANDIDATE_EPSG:
        try:
            if _projection_fingerprint(CRS.from_epsg(code)) == fingerprint:
                return code
        except Exception:  # pragma: no cover — pyproj DB 이상 시 다음 후보로
            continue
    return None


_TERRAIN_DATA_GLOBS = ("las/*", "laz/*", "tif/*", "tfw/*")


def find_project_prj(project_root: Path) -> Path | None:
    """프로젝트 **작업 좌표계**를 정하는 PRJ를 고른다.

    자료가 둘 이상의 좌표계로 들어온다(실측 2026-08-31 — 노선 shapefile은 UTM-K,
    지형은 동부원점 Bessel). 서피스 격자가 모델좌표의 주인이므로 **지형 PRJ**를 쓴다.
    노선 PRJ는 shapefile 세트 폴더(`input/shp/`)에 있어 여기서 섞이지 않는다.

    `input/prj/`에도 PRJ가 여럿 쌓인다 — 자료를 다시 올려도 파일명이 다르면 옛 PRJ가
    남기 때문이다. 이름 정렬로 고르면 잔재를 집는다(실측 2026-08-31: 옛 `result.prj`
    (5187)가 새 `용화.prj`(5176)보다 앞서 뽑혀 노선이 Y로 100,000m 어긋났고 B05
    경로 계산이 400으로 실패했다). 그래서 **지금 쓰는 지형 자료(LAS·TIF·TFW)와
    basename이 같은 PRJ**를 먼저 찾고, 못 찾으면 가장 최근 것을 쓴다.
    """
    prj_dir = project_root / "B03_FileInput" / "input" / "prj"
    candidates = sorted(prj_dir.glob("*.prj")) if prj_dir.is_dir() else []
    if candidates:
        if len(candidates) == 1:
            return candidates[0]
        input_root = project_root / "B03_FileInput" / "input"
        terrain_stems = {
            path.stem
            for pattern in _TERRAIN_DATA_GLOBS
            for path in input_root.glob(pattern)
            if path.is_file()
        }
        paired = [path for path in candidates if path.stem in terrain_stems]
        pool = paired or candidates
        chosen = max(pool, key=lambda path: path.stat().st_mtime)
        if not paired:
            logger.warning(
                "지형 자료와 짝이 되는 PRJ를 찾지 못해 가장 최근 PRJ를 씁니다: %s", chosen.name
            )
        return chosen
    remainder = sorted(project_root.glob("B03_FileInput/**/*.prj"))
    # 노선 세트 폴더의 PRJ는 노선 좌표계라 프로젝트 좌표계가 될 수 없다.
    remainder = [path for path in remainder if path.parent.name != "shp"] or remainder
    return remainder[0] if remainder else None


def crs_input_from_prj(prj_text: str) -> str | None:
    """PRJ 텍스트를 `Transformer.from_crs` 입력 문자열로 정규화한다.

    pyproj가 DB 정의와 일치를 확정하면 `"EPSG:n"`, 아니면 **원문 WKT 그대로**
    반환한다 — 파일에 박힌 TOWGS84 지역 보정을 보존하기 위해서다.
    WKT 자체가 불량이면 None (차단은 이때만).
    """
    text = (prj_text or "").strip()
    if not text:
        return None
    try:
        crs = CRS.from_wkt(text)
    except Exception:
        try:
            crs = CRS.from_user_input(text)
        except Exception:
            return None
    horizontal_raw = crs
    if crs.is_compound:
        for component in crs.sub_crs_list:
            flattened = strip_bound(component)
            if flattened.is_projected or flattened.is_geographic:
                horizontal_raw = component
                break
    if horizontal_raw.is_bound:
        # 수평 좌표계에 지역 보정(TOWGS84)이 박혀 있다 — EPSG 코드로 갈아타면
        # 그 보정이 사라지므로 라벨과 무관하게 원문 WKT로 변환한다.
        # (수직 성분의 bound — KNGeoid 지오이드 — 는 2D 변환에 무관하므로 무시.)
        label = identify_epsg(crs, text)
        logger.info(
            "PRJ에 사용자 TOWGS84 보정 — 파일 WKT로 변환 (라벨 %s)",
            f"EPSG:{label}" if label else "미상",
        )
        return text
    epsg = horizontal_crs(crs).to_epsg(min_confidence=25)
    if epsg is not None:
        return f"EPSG:{epsg}"
    label = identify_epsg(crs, text)
    if label is not None:
        logger.info("PRJ 좌표계 라벨 EPSG:%d — 변환은 파일 WKT로 수행(보정 모수 보존)", label)
    else:
        logger.warning("PRJ 좌표계의 EPSG 라벨을 찾지 못함 — 파일 WKT로 변환: %s", crs.name)
    return text


def project_prj_crs(project_root: Path) -> str | None:
    """프로젝트 **작업 좌표계** — 지형 PRJ 원문. PRJ가 없으면 None.

    작업 좌표계는 지형 PRJ 로 확정(2026-09-03 사용자 결정, 종전 동작 유지). 서피스 격자가
    모델좌표의 주인이고 노선은 `load_design_route()` 가 이 좌표계로 옮겨 오기 때문이다.
    노선 shapefile 은 다른 좌표계로 들어온다(실측 2026-08-31 — 노선 UTM-K / 지형 동부원점).

    기본값을 섞지 않는다 — 사다리(`resolve_project_crs`)가 다음 근거로 내려갈 수 있어야 한다.
    """
    prj_path = find_project_prj(project_root)
    if prj_path is None:
        return None
    from B04_PreProcess.B04_PreProcess_Engine_VWorld import get_epsg_from_prj

    return get_epsg_from_prj(prj_path.read_text(encoding="utf-8", errors="ignore"))


def resolve_project_crs(
    project_root: Path,
    *,
    route_crs_input: str | None = None,
    file_label_epsg: int | None = None,
    db_epsg: int | None = None,
) -> str:
    """작업 좌표계를 정하는 **단일 창구** — pyproj 입력 문자열(`EPSG:n` 또는 WKT 원문).

    사다리(위가 이김):
      ① `load_design_route()` 가 돌려준 노선 `crs_input` — 이미 작업 좌표계로 맞춰 둔 값
      ② `file_label_epsg` — **원본 파일 좌표를 그대로 해석할 때만** 넘긴다(노선 CSV 의
         `crs_epsg` 열). 그 자리에서는 파일이 스스로 밝힌 좌표계가 유일한 근거다.
      ③ 지형 PRJ(`project_prj_crs`) — 작업 좌표계 정본
      ④ DB `surface_models.crs_epsg`
      ⑤ 중부원점 `EPSG:5186` — 최후 폴백

    **변환이 끝난 좌표를 다룰 때는 ②를 넘기지 않는다.** 그 라벨은 실좌표계와 다른 사례가
    실측됐다(2026-09-01 용화 — 라벨 5179 / 실제 5176). 이 창구를 두기 전에는 같은 사다리가
    네 곳에 서로 다른 모양으로 흩어져 있었다.
    """
    if route_crs_input:
        return route_crs_input
    if file_label_epsg:
        return f"EPSG:{file_label_epsg}"
    prj_crs = project_prj_crs(project_root)
    if prj_crs:
        return prj_crs
    return f"EPSG:{db_epsg or 5186}"