123 lines
4.1 KiB
Python
123 lines
4.1 KiB
Python
"""보유 도엽 zip 실측 도곽 vs 도엽번호 알고리즘 대조 검증.
|
|
|
|
Downloads의 (B010)수치지도_{도엽번호}_*.zip에서 도곽선(N3L_A0010000 또는
|
|
N3A_A0010000) SHP를 읽어 EPSG:5187 → WGS84 변환 후, 알고리즘이 예측한
|
|
도곽 범위와 모서리 좌표 오차(m 단위 근사)를 비교한다.
|
|
"""
|
|
|
|
import io
|
|
import re
|
|
import sys
|
|
import zipfile
|
|
from pathlib import Path
|
|
|
|
PROJECT_ROOT = Path(__file__).resolve().parent.parent
|
|
sys.path.insert(0, str(PROJECT_ROOT))
|
|
|
|
import geopandas as gpd
|
|
|
|
from B04_wf1_Surface.B04_wf1_Surface_Engine_MapSheet import (
|
|
latlon_to_sheet5k,
|
|
neighbors_3x3,
|
|
sheet5k_center,
|
|
sheet5k_to_bounds,
|
|
)
|
|
|
|
DOWNLOADS = Path.home() / "Downloads"
|
|
|
|
|
|
def read_sheet_bounds_wgs84(zip_path: Path) -> tuple[float, float, float, float]:
|
|
"""zip 내 도곽선(A0010000) SHP의 전체 bounds를 WGS84로 반환."""
|
|
with zipfile.ZipFile(zip_path) as zf:
|
|
names = zf.namelist()
|
|
shp = [n for n in names if re.search(r"A0010000\.shp$", n, re.IGNORECASE)]
|
|
if not shp:
|
|
raise RuntimeError(f"{zip_path.name}: 도곽선(A0010000) SHP 없음")
|
|
target = shp[0]
|
|
base = target[:-4]
|
|
tmp = {}
|
|
for ext in (".shp", ".shx", ".dbf", ".prj", ".cpg"):
|
|
member = base + ext
|
|
if member in names:
|
|
tmp[ext] = zf.read(member)
|
|
|
|
import tempfile
|
|
|
|
with tempfile.TemporaryDirectory() as td:
|
|
for ext, data in tmp.items():
|
|
(Path(td) / f"sheet{ext}").write_bytes(data)
|
|
gdf = gpd.read_file(Path(td) / "sheet.shp")
|
|
if gdf.crs is None:
|
|
gdf = gdf.set_crs("EPSG:5187")
|
|
gdf = gdf.to_crs("EPSG:4326")
|
|
return tuple(gdf.total_bounds) # (lon_min, lat_min, lon_max, lat_max)
|
|
|
|
|
|
def main():
|
|
zips = sorted(DOWNLOADS.glob("*수치지도_*.zip")) + [
|
|
p for p in (DOWNLOADS / "36906042.zip", DOWNLOADS / "36902002.zip") if p.exists()
|
|
]
|
|
if not zips:
|
|
print("Downloads에 도엽 zip 없음")
|
|
return
|
|
|
|
print(
|
|
f"{'도엽':>10} | {'실측(중심) 위경도':>24} | {'예측번호':>10} | 모서리 최대오차(deg) | 판정"
|
|
)
|
|
all_ok = True
|
|
for zp in zips:
|
|
m = re.search(r"(\d{8})", zp.name)
|
|
if not m:
|
|
continue
|
|
sheet_no = m.group(1)
|
|
try:
|
|
actual = read_sheet_bounds_wgs84(zp)
|
|
except Exception as e:
|
|
print(f"{sheet_no:>10} | 읽기 실패: {e}")
|
|
all_ok = False
|
|
continue
|
|
lon_min, lat_min, lon_max, lat_max = actual
|
|
lat_c = (lat_min + lat_max) / 2
|
|
lon_c = (lon_min + lon_max) / 2
|
|
predicted_no = latlon_to_sheet5k(lat_c, lon_c)
|
|
pb = sheet5k_to_bounds(sheet_no)
|
|
err = max(
|
|
abs(pb[0] - lon_min),
|
|
abs(pb[1] - lat_min),
|
|
abs(pb[2] - lon_max),
|
|
abs(pb[3] - lat_max),
|
|
)
|
|
# 해안 도엽은 바다 쪽 데이터가 잘려 실측 범위가 도곽보다 작을 수 있음
|
|
# → 판정 기준: 번호 일치 + 실측 범위가 예측 도곽 안에 포함(여유 0.0005도)
|
|
tol = 0.0005
|
|
contained = (
|
|
lon_min >= pb[0] - tol
|
|
and lat_min >= pb[1] - tol
|
|
and lon_max <= pb[2] + tol
|
|
and lat_max <= pb[3] + tol
|
|
)
|
|
ok = predicted_no == sheet_no and contained
|
|
all_ok = all_ok and ok
|
|
print(
|
|
f"{sheet_no:>10} | ({lat_c:.5f}, {lon_c:.5f}) | {predicted_no:>10} | "
|
|
f"{err:.6f} | {'OK' if ok else 'MISMATCH'}"
|
|
)
|
|
|
|
print("\n--- 3x3 인접 도엽 산출 테스트 (36906042 중심) ---")
|
|
grid = neighbors_3x3("36906042")
|
|
for i in range(0, 9, 3):
|
|
print(" ", grid[i : i + 3])
|
|
|
|
print("\n--- 구획 경계 이월 테스트 ---")
|
|
for probe in ("36902001", "36902091", "36916100"):
|
|
print(f" {probe} center={sheet5k_center(probe)} -> 3x3:")
|
|
g = neighbors_3x3(probe)
|
|
for i in range(0, 9, 3):
|
|
print(" ", g[i : i + 3])
|
|
|
|
print("\n결과:", "전체 정합" if all_ok else "불일치 존재 — 알고리즘 재검토 필요")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|