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