eomsangdonandClaude Opus 5.5
9c2b5e151b
- 도엽이 위성 상자 밖 0.8 폭까지 나가 9×9 격자 늘려 씀 자리에서 0.33 m 어긋남(화면 실측) - 격자 범위 = 상자 ± 한 폭(3 × 3 폭) · 상자 안 칸 크기 그대로(8 등분) · 격자 자기 경위도 범위 칸 더함 - 지킴 시험 = 상자 밖 0.9 폭까지 pyproj 0.02 m 안 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
251 lines
11 KiB
Python
251 lines
11 KiB
Python
"""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.8 폭까지 · 3 × 3 폭) · 한 변 25 마디
|
|
# = 상자 안 8 등분 · 쌍선형 나머지 용화 0.001 · D 0.006 · 큰 상자 0.014 m(589-1).
|
|
TM_GRID_NODES = 25
|
|
TM_GRID_MARGIN = 1.0
|
|
|
|
|
|
def tm_grid_for_meta(meta: dict, crs_input: str, nodes: int = TM_GRID_NODES) -> dict | None:
|
|
"""그림 상자 둘레 경위도 고른 마디의 TM 값(589-1 · 화면이 경위도 겹을 TM 자리로 놓음).
|
|
|
|
격자 범위 = 상자 경위도 범위를 둘레로 `TM_GRID_MARGIN` 폭씩 넓힌 것 · 마디 `[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)
|
|
(x0, x1), (y1, y0) = transformer.transform(
|
|
[meta["lon_min"], meta["lon_max"]], [meta["lat_max"], meta["lat_min"]]
|
|
)
|
|
off = max(
|
|
abs(x0 - meta["x_min"]),
|
|
abs(y1 - meta["y_max"]),
|
|
abs(x1 - meta["x_max"]),
|
|
abs(y0 - meta["y_min"]),
|
|
)
|
|
if not math.isfinite(off) or off > 0.01:
|
|
logger.warning("위성 메타 귀퉁이와 TM 격자가 %.3f m 어긋남 — 격자 안 줌", off)
|
|
return None
|
|
lon_pad = (meta["lon_max"] - meta["lon_min"]) * TM_GRID_MARGIN
|
|
lat_pad = (meta["lat_max"] - meta["lat_min"]) * TM_GRID_MARGIN
|
|
lon_min, lon_max = meta["lon_min"] - lon_pad, meta["lon_max"] + lon_pad
|
|
lat_min, lat_max = meta["lat_min"] - lat_pad, meta["lat_max"] + lat_pad
|
|
lons = [lon_min + i / (nodes - 1) * (lon_max - lon_min) for i in range(nodes)]
|
|
xs: list[float] = []
|
|
ys: list[float] = []
|
|
for j in range(nodes):
|
|
lat = lat_min + j / (nodes - 1) * (lat_max - lat_min)
|
|
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)
|
|
return {
|
|
"n": nodes,
|
|
"lon_min": lon_min,
|
|
"lon_max": lon_max,
|
|
"lat_min": lat_min,
|
|
"lat_max": lat_max,
|
|
"x": [round(v, 4) for v in xs],
|
|
"y": [round(v, 4) for v in ys],
|
|
}
|