"""보유 도엽 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_PreProcess.B04_PreProcess_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()