Aislo/B04_PreProcess/B04_PreProcess_Engine_VWorld.py

"""VWorld 위성 맵 및 타일 다운로드 엔진."""

from __future__ import annotations

import json
import logging
import math
import urllib.parse
import urllib.request
from io import BytesIO
from pathlib import Path

from PIL import Image
from pyproj import Transformer

logger = logging.getLogger(__name__)

# 설정 불러오기
try:
    from config import config_system

    VWORLD_API_KEY = getattr(
        config_system, "VWORLD_API_KEY", "3DBD7306-7DBD-38BB-B292-267C5ED7AC6B"
    )
    SURFACE_MAP_MAX_TILES_PER_SIDE = getattr(config_system, "SURFACE_MAP_MAX_TILES_PER_SIDE", 12)
    SURFACE_MAP_MAX_ZOOM = getattr(config_system, "SURFACE_MAP_MAX_ZOOM", 18)
    SURFACE_MAP_MIN_ZOOM = getattr(config_system, "SURFACE_MAP_MIN_ZOOM", 14)
except ImportError:
    VWORLD_API_KEY = "3DBD7306-7DBD-38BB-B292-267C5ED7AC6B"
    SURFACE_MAP_MAX_TILES_PER_SIDE = 12
    SURFACE_MAP_MAX_ZOOM = 18
    SURFACE_MAP_MIN_ZOOM = 14


def get_epsg_from_prj(prj_content: str) -> str:
    """PRJ 텍스트를 좌표 변환 입력 문자열로 정규화합니다.

    pyproj가 EPSG DB 정의와 일치를 확정하면 `"EPSG:n"`, 아니면 파일 WKT 그대로
    반환합니다(TOWGS84 지역 보정 보존) — 둘 다 `Transformer.from_crs` 입력으로
    동작합니다. WKT가 불량·공백일 때만 기존 기본값 중부원점을 유지합니다.

    (2026-08-31 수리: 예전 문자열 매칭은 WKT에 항상 있는 `false_easting` 때문에
    "EAST" 분기가 무조건 걸려 모든 PRJ를 EPSG:5187로 판정했다.)
    """
    from common_util.common_util_crs import crs_input_from_prj

    crs_input = crs_input_from_prj(prj_content)
    if crs_input is None:
        logger.warning("PRJ 좌표계 해석 실패 — 기본값 EPSG:5186 사용")
        return "EPSG:5186"
    return crs_input


def latlon_to_tile(lat: float, lon: float, zoom: int) -> tuple[int, int]:
    """위도, 경도를 OSM/VWorld 표준 지도 타일 좌표 (X, Y)로 변환합니다."""
    lat_rad = math.radians(lat)
    n = 2.0**zoom
    xtile = int((lon + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return xtile, ytile


def tile_to_latlon(xtile: int, ytile: int, zoom: int) -> tuple[float, float]:
    """OSM/VWorld 표준 타일 좌표 (X, Y)를 위도, 경도 좌표로 역변환합니다."""
    n = 2.0**zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return lat_deg, lon_deg


def download_vworld_satellite_map(
    prj_path: Path, bounds: dict, output_dir: Path, layer_name: str = "Satellite", ext: str = "jpeg"
) -> Path | None:
    """PRJ 좌표 정보와 지형 범위(bounds)를 바탕으로,
    해당 임도 산림 지역을 감싸는 VWorld 지도 타일을 다운로드하여 하나의 PNG 텍스처로 병합합니다.
    """
    output_dir.mkdir(parents=True, exist_ok=True)
    file_name = f"vworld_{layer_name.lower()}.png"
    map_png_path = output_dir / file_name
    meta_json_path = output_dir / f"vworld_{layer_name.lower()}_meta.json"

    # 1. PRJ 내용 분석
    if not prj_path.exists():
        return None
    prj_content = prj_path.read_text(encoding="utf-8", errors="ignore")
    src_epsg = get_epsg_from_prj(prj_content)

    # 2. 중심점 및 경계 좌표 위경도 변환
    # bounds: {'x': [x_min, x_max], 'y': [y_min, y_max], 'z': [z_min, z_max]}
    x_min, x_max = bounds["x"][0], bounds["x"][1]
    y_min, y_max = bounds["y"][0], bounds["y"][1]

    # 좌표 변환기 생성 (지정 좌표계 -> 경위도 EPSG:4326)
    try:
        transformer = Transformer.from_crs(src_epsg, "EPSG:4326", always_xy=True)
    except Exception:
        # 실패 시 강제로 Korea Central -> WGS84 생성
        transformer = Transformer.from_crs("EPSG:5186", "EPSG:4326", always_xy=True)

    lon_min, lat_min = transformer.transform(x_min, y_min)
    lon_max, lat_max = transformer.transform(x_max, y_max)

    # 3. 요청 범위를 그대로 덮는 타일 범위를 잡는다(주변으로 넓히지 않는다).
    #    한 변 타일 수가 한도를 넘으면 zoom을 한 단계씩 낮춘다 — 도엽 1매를 zoom 18로 받으면
    #    한 변이 19타일(4,864px)이라 파일이 지나치게 커진다(2026-08-01 사용자 지시).
    zoom = SURFACE_MAP_MAX_ZOOM
    while True:
        x1, y1 = latlon_to_tile(lat_max, lon_min, zoom)
        x2, y2 = latlon_to_tile(lat_min, lon_max, zoom)
        x_start, x_end = min(x1, x2), max(x1, x2)
        y_start, y_end = min(y1, y2), max(y1, y2)
        tile_w = x_end - x_start + 1
        tile_h = y_end - y_start + 1
        if zoom <= SURFACE_MAP_MIN_ZOOM:
            break
        if max(tile_w, tile_h) <= SURFACE_MAP_MAX_TILES_PER_SIDE:
            break
        zoom -= 1

    # 4. 개별 타일 다운로드 및 이미지 병합
    map_img = Image.new("RGBA", (tile_w * 256, tile_h * 256))
    headers = {"User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64)"}

    for ty in range(y_start, y_end + 1):
        for tx in range(x_start, x_end + 1):
            try:
                normalized_layer = layer_name.lower()
                if normalized_layer == "satellite":
                    sat_url = f"http://api.vworld.kr/req/wmts/1.0.0/{VWORLD_API_KEY}/Satellite/{zoom}/{ty}/{tx}.jpeg"
                    sat_req = urllib.request.Request(sat_url, headers=headers)
                    with urllib.request.urlopen(sat_req, timeout=5) as sat_res:
                        sat_data = sat_res.read()
                    base_img = Image.open(BytesIO(sat_data)).convert("RGBA")
                elif normalized_layer == "hybrid":
                    hy_url = f"http://api.vworld.kr/req/wmts/1.0.0/{VWORLD_API_KEY}/Hybrid/{zoom}/{ty}/{tx}.png"
                    hy_req = urllib.request.Request(hy_url, headers=headers)
                    with urllib.request.urlopen(hy_req, timeout=5) as hy_res:
                        hy_data = hy_res.read()
                    base_img = Image.open(BytesIO(hy_data)).convert("RGBA")
                elif normalized_layer == "white":
                    white_url = f"http://api.vworld.kr/req/wmts/1.0.0/{VWORLD_API_KEY}/white/{zoom}/{ty}/{tx}.png"
                    white_req = urllib.request.Request(white_url, headers=headers)
                    with urllib.request.urlopen(white_req, timeout=5) as white_res:
                        white_data = white_res.read()
                    base_img = Image.open(BytesIO(white_data)).convert("RGBA")
                else:
                    layer_url = f"http://api.vworld.kr/req/wmts/1.0.0/{VWORLD_API_KEY}/{layer_name}/{zoom}/{ty}/{tx}.{ext}"
                    layer_req = urllib.request.Request(layer_url, headers=headers)
                    with urllib.request.urlopen(layer_req, timeout=5) as layer_res:
                        layer_data = layer_res.read()
                    base_img = Image.open(BytesIO(layer_data)).convert("RGBA")

                # 병합 좌표 계산
                px = (tx - x_start) * 256
                py = (ty - y_start) * 256
                map_img.paste(base_img, (px, py))
            except Exception:
                # 다운로드 실패 시 회색 격자 텍스처로 대체
                fallback_tile = Image.new("RGBA", (256, 256), (200, 200, 200, 255))
                px = (tx - x_start) * 256
                py = (ty - y_start) * 256
                map_img.paste(fallback_tile, (px, py))

    # 5. 병합된 이미지의 실제 위경도 바운더리 계산
    lat_top, lon_left = tile_to_latlon(x_start, y_start, zoom)
    lat_bottom, lon_right = tile_to_latlon(x_end + 1, y_end + 1, zoom)

    # 위경도를 다시 원본 Local 3D 미터 좌표계로 변환
    try:
        rev_transformer = Transformer.from_crs("EPSG:4326", src_epsg, always_xy=True)
    except Exception:
        rev_transformer = Transformer.from_crs("EPSG:4326", "EPSG:5186", always_xy=True)

    img_x_min, img_y_max = rev_transformer.transform(lon_left, lat_top)
    img_x_max, img_y_min = rev_transformer.transform(lon_right, lat_bottom)

    # 6. 이미지 저장 및 메타데이터 작성
    map_img.save(map_png_path, "PNG")

    meta = {
        "x_min": img_x_min,
        "x_max": img_x_max,
        "y_min": img_y_min,
        "y_max": img_y_max,
        "width_meters": img_x_max - img_x_min,
        "height_meters": img_y_max - img_y_min,
        "center_x": (img_x_min + img_x_max) / 2.0,
        "center_y": (img_y_min + img_y_max) / 2.0,
        "lon_min": lon_left,
        "lon_max": lon_right,
        "lat_min": lat_bottom,
        "lat_max": lat_top,
    }
    meta_json_path.write_text(json.dumps(meta, indent=2), encoding="utf-8")

    return map_png_path


# 그림 상자 경위도 격자 한 변 마디 수 — 쌍선형 나머지 용화 0.001 · D 0.006 · 큰 상자 0.014 m(589-1).
TM_GRID_NODES = 9


def tm_grid_for_meta(meta: dict, crs_input: str, nodes: int = TM_GRID_NODES) -> dict | None:
    """그림 상자 경위도 고른 마디의 TM 값(589-1 · 화면이 경위도 겹을 TM 자리로 놓음).

    마디 `[j * n + i]` = 경도 `lon_min + i/(n-1)·폭` · 위도 `lat_min + j/(n-1)·높이`(j 는 남 → 북).
    왼위 · 오른아래 마디가 메타 `x_min · y_max` · `x_max · y_min` 과 1 cm 넘게 다르면(메타를 다른
    좌표계로 만든 것) 격자를 주지 않음 — 화면은 옛 같은 네모 식으로 남음.
    """
    keys = ("lon_min", "lon_max", "lat_min", "lat_max", "x_min", "x_max", "y_min", "y_max")
    if any(not isinstance(meta.get(key), (int, float)) for key in keys):
        return None
    transformer = Transformer.from_crs("EPSG:4326", crs_input, always_xy=True)
    lon_span = meta["lon_max"] - meta["lon_min"]
    lat_span = meta["lat_max"] - meta["lat_min"]
    lons = [meta["lon_min"] + i / (nodes - 1) * lon_span for i in range(nodes)]
    lats = [meta["lat_min"] + j / (nodes - 1) * lat_span for j in range(nodes)]
    xs: list[float] = []
    ys: list[float] = []
    for lat in lats:
        x_row, y_row = transformer.transform(lons, [lat] * nodes)
        xs.extend(float(v) for v in x_row)
        ys.extend(float(v) for v in y_row)
    top_left = nodes * (nodes - 1)
    bottom_right = nodes - 1
    off = max(
        abs(xs[top_left] - meta["x_min"]),
        abs(ys[top_left] - meta["y_max"]),
        abs(xs[bottom_right] - meta["x_max"]),
        abs(ys[bottom_right] - meta["y_min"]),
    )
    if not math.isfinite(off) or off > 0.01:
        logger.warning("위성 메타 귀퉁이와 TM 격자가 %.3f m 어긋남 — 격자 안 줌", off)
        return None
    return {"n": nodes, "x": [round(v, 4) for v in xs], "y": [round(v, 4) for v in ys]}