Files
Aislo/scratch/test_mapsheet_verification.py
T
eomsangdonandClaude Fable 5 f7528a4aa4 refactor(B04): B04_wf1_Surface -> B04_PreProcess 전면 개명
- 폴더·내부 파일 51개 접두사 개명 (git mv, 이력 보존)
- 저장소 전체 참조 치환 67파일: import 경로, 라우트 슬러그(b04-preprocess),
  라우트 키(B04_PREPROCESS), storage 경로 상수, locale, SQL 주석
- 로직 변경 없음 (기계적 치환). typecheck·백엔드 import 검증 통과

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-08-08 10:01:36 +09:00

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_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()