Aislo/B04_wf1_Surface/B04_wf1_Surface_Engine_VWorld.py

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

from __future__ import annotations

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

from PIL import Image
from pyproj import Transformer

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

    VWORLD_API_KEY = getattr(
        config_system, "VWORLD_API_KEY", "3DBD7306-7DBD-38BB-B292-267C5ED7AC6B"
    )
except ImportError:
    VWORLD_API_KEY = "3DBD7306-7DBD-38BB-B292-267C5ED7AC6B"


def get_epsg_from_prj(prj_content: str) -> str:
    """PRJ 파일 텍스트 내용을 해석하여 적합한 EPSG 좌표계 코드를 추정합니다.

    기본 한반도 주요 좌표계 매핑:
      - Korea Central Belt (중부원점 20만, 60만) -> EPSG:5186 (Korea 2000 Central)
      - Korea East Belt (동부원점) -> EPSG:5187
      - Korea West Belt (서부원점) -> EPSG:5185
      - WGS84 UTM 52N -> EPSG:32652
    """
    prj_upper = prj_content.upper()
    if "EAST" in prj_upper:
        return "EPSG:5187"
    elif "WEST" in prj_upper:
        return "EPSG:5185"
    elif "CENTRAL" in prj_upper:
        return "EPSG:5186"
    elif "UTM" in prj_upper and "52N" in prj_upper:
        return "EPSG:32652"

    # 디폴트는 한국 중부 2000 좌표계 적용
    return "EPSG:5186"


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 18 초고해상도 적용)
    zoom = 18

    # 영역을 포괄하는 좌상단 타일, 우하단 타일 인덱스 산출
    x1, y1 = latlon_to_tile(lat_max, lon_min, zoom)
    x2, y2 = latlon_to_tile(lat_min, lon_max, zoom)

    # 타일 경계 마진 패딩
    x_start = min(x1, x2) - 1
    x_end = max(x1, x2) + 1
    y_start = min(y1, y2) - 1
    y_end = max(y1, y2) + 1

    tile_w = x_end - x_start + 1
    tile_h = y_end - y_start + 1

    # 너무 많은 타일을 다운로드하여 IP 차단되는 것을 방지
    if tile_w > 15:
        tile_w = 15
        x_end = x_start + 14
    if tile_h > 15:
        tile_h = 15
        y_end = y_start + 14

    # 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