Aislo/0_old/utils/download_gis_vectors.py

import sys
import urllib.request
import urllib.parse
from pathlib import Path
import json

# 루트 경로 추가
root_dir = Path(__file__).resolve().parent.parent
sys.path.append(str(root_dir))

try:
    import config
    VWORLD_API_KEY = config.VWORLD_API_KEY
except ImportError:
    VWORLD_API_KEY = "3DBD7306-7DBD-38BB-B292-267C5ED7AC6B"


def request_vworld_wfs(typename: str, filename: str, bounds_wgs84: dict, output_dir: Path):
    """브이월드 WFS API를 사용하여 지정된 Bounds 영역 내의 공간 벡터 데이터를 GeoJSON으로 다운로드합니다."""
    # 4326(위경도) bbox 파라미터 구성 (BBOX=minx,miny,maxx,maxy)
    bbox_str = f"{bounds_wgs84['min_lon']},{bounds_wgs84['min_lat']},{bounds_wgs84['max_lon']},{bounds_wgs84['max_lat']}"
    
    # WFS 파라미터 조립
    params = {
        "key": VWORLD_API_KEY,
        "domain": "localhost",
        "service": "WFS",
        "version": "1.1.0",
        "request": "GetFeature",
        "typename": typename,
        "bbox": bbox_str,
        "output": "application/json", # GeoJSON 포맷 요청
        "srsName": "EPSG:4326"
    }
    
    url = "http://api.vworld.kr/req/wfs?" + urllib.parse.urlencode(params)
    headers = {"User-Agent": "Mozilla/5.0"}
    
    req = urllib.request.Request(url, headers=headers)
    try:
        with urllib.request.urlopen(req, timeout=15) as res:
            content = res.read()
            
        # 결과가 JSON 데이터인지 검증
        data = json.loads(content.decode("utf-8"))
        
        # 파일 쓰기
        output_path = output_dir / filename
        output_path.write_text(json.dumps(data, indent=2, ensure_ascii=False), encoding="utf-8")
        
        features_cnt = len(data.get("features", []))
        print(f"  -> {filename} 다운로드 완료 (피처 수: {features_cnt}개)")
    except Exception as e:
        print(f"  -> {filename} 다운로드 실패: {e}")


def download_all_gis_vectors(prj_path: Path, bounds_meter: dict, output_dir: Path):
    """지형의 로컬 미터단위 bounds 정보를 위경도로 변환 후 5대 국가 GIS 데이터를 다운로드합니다."""
    from pyproj import Transformer
    from utils.vworld_downloader import get_epsg_from_prj
    
    output_dir.mkdir(parents=True, exist_ok=True)
    
    # 1. vworld_white_meta.json이 생성되어 있다면 직접 위경도 경계 획득하여 pyproj 축왜곡 방지
    meta_path = prj_path.parent / "vworld_white_meta.json"
    if not meta_path.exists():
        meta_path = prj_path.parent / "vworld_satellite_meta.json"
        
    if meta_path.exists():
        try:
            with open(meta_path, "r", encoding="utf-8") as f:
                meta_data = json.load(f)
            # 주변 데이터를 충분히 확보하기 위해 수집 반경(마진 패딩)을 0.002도에서 0.010도(약 1km 반경)로 확장
            bounds_wgs84 = {
                "min_lon": meta_data["lon_min"] - 0.010,
                "max_lon": meta_data["lon_max"] + 0.010,
                "min_lat": meta_data["lat_min"] - 0.010,
                "max_lat": meta_data["lat_max"] + 0.010
            }
            print(f"WFS 영역 요청 WGS84 Bounds (메타파일 참조 - 마진확장): {bounds_wgs84}")
        except Exception as ex:
            print(f"메타파일 좌표 파싱 오류: {ex}. pyproj 폴백 사용.")
            # 폴백용 pyproj 변환
            src_epsg = "EPSG:5186"
            if prj_path.exists():
                src_epsg = get_epsg_from_prj(prj_path.read_text(encoding="utf-8", errors="ignore"))
            transformer = Transformer.from_crs(src_epsg, "EPSG:4326", always_xy=True)
            x_min, x_max = bounds_meter["x"][0], bounds_meter["x"][1]
            y_min, y_max = bounds_meter["y"][0], bounds_meter["y"][1]
            lon_min, lat_min = transformer.transform(x_min, y_min)
            lon_max, lat_max = transformer.transform(x_max, y_max)
            bounds_wgs84 = {
                "min_lon": min(lon_min, lon_max) - 0.002,
                "max_lon": max(lon_min, lon_max) + 0.002,
                "min_lat": min(lat_min, lat_max) - 0.002,
                "max_lat": max(lat_min, lat_max) + 0.002
            }
    else:
        # 폴백용 pyproj 변환
        src_epsg = "EPSG:5186"
        if prj_path.exists():
            src_epsg = get_epsg_from_prj(prj_path.read_text(encoding="utf-8", errors="ignore"))
        transformer = Transformer.from_crs(src_epsg, "EPSG:4326", always_xy=True)
        x_min, x_max = bounds_meter["x"][0], bounds_meter["x"][1]
        y_min, y_max = bounds_meter["y"][0], bounds_meter["y"][1]
        lon_min, lat_min = transformer.transform(x_min, y_min)
        lon_max, lat_max = transformer.transform(x_max, y_max)
        bounds_wgs84 = {
            "min_lon": min(lon_min, lon_max) - 0.010,
            "max_lon": max(lon_min, lon_max) + 0.010,
            "min_lat": min(lat_min, lat_max) - 0.010,
            "max_lat": max(lat_min, lat_max) + 0.010
        }
    
    # 5대 레이어 다운로드
    request_vworld_wfs("lt_c_landinfobasemap", "연속지적도_bounds.geojson", bounds_wgs84, output_dir)
    request_vworld_wfs("lt_c_uq111", "용도지역도_bounds.geojson", bounds_wgs84, output_dir)
    request_vworld_wfs("lt_c_adsigg", "행정구역_시군구_bounds.geojson", bounds_wgs84, output_dir)
    request_vworld_wfs("lt_c_adri", "행정구역_읍면동_bounds.geojson", bounds_wgs84, output_dir)
    request_vworld_wfs("lt_c_gimshydro", "수계망_물줄기_bounds.geojson", bounds_wgs84, output_dir)
    request_vworld_wfs("lt_c_kfdrssigugrade", "산사태위험등급_bounds.geojson", bounds_wgs84, output_dir)

    # 6. gpkg 등고선(national_contours.gpkg) 데이터에서 기준 영역 크롭하여 등고선_bounds.geojson 생성
    gpkg_path = root_dir / "source" / "grobal_contours" / "national_contours.gpkg"
    if gpkg_path.exists():
        print(f"  -> 등고선 추출 시작 ({gpkg_path.name})")
        try:
            import geopandas as gpd
            from shapely.geometry import box
            
            # WGS84 영역 박스를 GPKG의 EPSG:5179 좌표계로 변환하여 로드하기 위해 
            # transformer 재설정 (4326 -> 5179)
            from pyproj import Transformer
            t_to_5179 = Transformer.from_crs("EPSG:4326", "EPSG:5179", always_xy=True)
            min_x_5179, min_y_5179 = t_to_5179.transform(bounds_wgs84["min_lon"], bounds_wgs84["min_lat"])
            max_x_5179, max_y_5179 = t_to_5179.transform(bounds_wgs84["max_lon"], bounds_wgs84["max_lat"])
            
            bbox = box(min_x_5179, min_y_5179, max_x_5179, max_y_5179)
            
            # GeoPackage에서 bbox 교차 영역만 스트리밍 로드
            gdf = gpd.read_file(gpkg_path, bbox=bbox)
            if not gdf.empty:
                # WGS84(4326)로 변환
                gdf_wgs84 = gdf.to_crs("EPSG:4326")
                geojson_out = output_dir / "등고선_bounds.geojson"
                gdf_wgs84.to_file(geojson_out, driver="GeoJSON")
                print(f"  -> 등고선 추출 및 등고선_bounds.geojson 저장 완료 (피처 수: {len(gdf_wgs84)}개)")
            else:
                print("  -> 해당 영역의 등고선 데이터가 GPKG에 존재하지 않습니다.")
        except Exception as e:
            print(f"  -> 등고선 GPKG 크롭 연산 실패: {e}")
    else:
        print(f"  -> {gpkg_path} 파일이 존재하지 않아 등고선 추출을 생략합니다.")


if __name__ == "__main__":
    # 수동 검증용 (강원도 정선군 신동읍 일대)
    dummy_bounds = {
        "x": [183433.605, 183805.800],
        "y": [489168.796, 489570.168]
    }
    prj_file = root_dir / "samples" / "step1_scan" / "result.prj"
    download_all_gis_vectors(prj_file, dummy_bounds, root_dir / "instance" / "sample")