auto: 2026-07-26 18:05 (ESD_LAPTOP)

2026-07-26 18:05:36 +09:00
parent 0a1913eb6c
commit 97f5c824a9
5 changed files with 444 additions and 0 deletions
@@ -0,0 +1,122 @@
+"""보유 도엽 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()
@@ -0,0 +1,47 @@
+"""도엽 영구저장소(SheetStore) 인제스트 실테스트.
+
+Downloads의 보유 도엽 zip 전량을 resources/map_sheets/에 등록하고,
+CRS 오라벨 zip(36902002.zip) 자동 교정과 ensure_sheets 3x3 흐름을 확인한다.
+"""
+
+import sys
+from pathlib import Path
+
+PROJECT_ROOT = Path(__file__).resolve().parent.parent
+sys.path.insert(0, str(PROJECT_ROOT))
+
+from B04_wf1_Surface.B04_wf1_Surface_Engine_MapSheet import neighbors_3x3
+from B04_wf1_Surface.B04_wf1_Surface_Engine_SheetStore import (
+    ensure_sheets,
+    get_sheet_path,
+    ingest_zip,
+)
+
+DOWNLOADS = Path.home() / "Downloads"
+
+
+def main():
+    zips = sorted(DOWNLOADS.glob("*수치지도_*.zip")) + [
+        p for p in (DOWNLOADS / "36906042.zip", DOWNLOADS / "36902002.zip") if p.exists()
+    ]
+    for zp in zips:
+        try:
+            entry = ingest_zip(zp)
+            print(f"OK   {zp.name} -> {entry['file']} crs={entry['crs']}")
+        except Exception as e:
+            print(f"FAIL {zp.name}: {e}")
+
+    print("\n--- ensure_sheets: 36906042 중심 3x3 ---")
+    result = ensure_sheets(neighbors_3x3("36906042"))
+    print("보유:", sorted(result["available"].keys()))
+    print("미보유:", result["missing"])
+    if "guidance" in result:
+        print("안내 URL:", result["guidance"]["vworld_dataset_url"])
+
+    print("\n--- get_sheet_path ---")
+    print("36906042:", get_sheet_path("36906042"))
+    print("99999999:", get_sheet_path("99999999"))
+
+
+if __name__ == "__main__":
+    main()