diff --git a/B04_wf1_Surface/B04_wf1_Surface_Engine_MapSheet.py b/B04_wf1_Surface/B04_wf1_Surface_Engine_MapSheet.py new file mode 100644 index 00000000..e1b0c142 --- /dev/null +++ b/B04_wf1_Surface/B04_wf1_Surface_Engine_MapSheet.py @@ -0,0 +1,90 @@ +# B04_wf1_Surface_Engine_MapSheet.py +# 1:5,000 수치지형도 도엽번호 산출 엔진. +# 국토지리정보원 도엽번호 체계: +# - 경위도 1°x1° 구획: 위도 2자리 + 경도 끝 1자리 (예: 369 = 위도 36~37°, 경도 129~130°) +# - 1:50,000 = 1° 구획을 4x4(15'x15') 분할, NW부터 행 우선 01~16 (예: 36906) +# - 1:5,000 = 1:50,000 도엽을 10x10(1.5'x1.5') 분할, NW부터 행 우선 001~100 +# (예: 36906042 = 36906 도엽의 042 구획) + +from __future__ import annotations + +# 1:5,000 도엽 한 변 크기 (도 단위): 15' / 10 = 1.5' +SHEET5K_SIZE_DEG = 0.025 +# 1:50,000 도엽 한 변 크기 (도 단위): 15' +SHEET50K_SIZE_DEG = 0.25 + + +def _lon_digit_to_degree(digit: int) -> int: + """경도 끝자리 → 경도 정수부 복원. 한반도 범위(124~132°E)에서 유일.""" + if digit >= 4: + return 120 + digit # 4~9 → 124~129 + return 130 + digit # 0~2 → 130~132 + + +def latlon_to_sheet5k(lat: float, lon: float) -> str: + """WGS84 위경도 → 1:5,000 도엽번호(8자리 문자열).""" + lat_deg = int(lat) + lon_deg = int(lon) + cell = f"{lat_deg:02d}{lon_deg % 10}" + + # 1° 구획 내 오프셋 (북서 원점 기준: 북쪽에서 남쪽으로 행 증가) + lat_frac = (lat_deg + 1) - lat # 구획 상단으로부터 남쪽 거리 + lon_frac = lon - lon_deg + + row50 = min(int(lat_frac / SHEET50K_SIZE_DEG), 3) + col50 = min(int(lon_frac / SHEET50K_SIZE_DEG), 3) + idx50 = row50 * 4 + col50 + 1 # 01~16 + + # 1:50,000 도엽 내 오프셋 + lat_in_50 = lat_frac - row50 * SHEET50K_SIZE_DEG + lon_in_50 = lon_frac - col50 * SHEET50K_SIZE_DEG + row5 = min(int(lat_in_50 / SHEET5K_SIZE_DEG), 9) + col5 = min(int(lon_in_50 / SHEET5K_SIZE_DEG), 9) + idx5 = row5 * 10 + col5 + 1 # 001~100 + + return f"{cell}{idx50:02d}{idx5:03d}" + + +def sheet5k_to_bounds(sheet_no: str) -> tuple[float, float, float, float]: + """1:5,000 도엽번호 → WGS84 도곽 범위 (lon_min, lat_min, lon_max, lat_max).""" + s = str(sheet_no).strip() + if len(s) != 8 or not s.isdigit(): + raise ValueError(f"1:5,000 도엽번호는 8자리 숫자여야 함: {sheet_no!r}") + + lat_deg = int(s[0:2]) + lon_deg = _lon_digit_to_degree(int(s[2])) + idx50 = int(s[3:5]) + idx5 = int(s[5:8]) + if not (1 <= idx50 <= 16): + raise ValueError(f"1:50,000 구획번호(01~16) 범위 밖: {sheet_no!r}") + if not (1 <= idx5 <= 100): + raise ValueError(f"1:5,000 구획번호(001~100) 범위 밖: {sheet_no!r}") + + row50, col50 = divmod(idx50 - 1, 4) + row5, col5 = divmod(idx5 - 1, 10) + + lat_max = (lat_deg + 1) - row50 * SHEET50K_SIZE_DEG - row5 * SHEET5K_SIZE_DEG + lon_min = lon_deg + col50 * SHEET50K_SIZE_DEG + col5 * SHEET5K_SIZE_DEG + return (lon_min, lat_max - SHEET5K_SIZE_DEG, lon_min + SHEET5K_SIZE_DEG, lat_max) + + +def sheet5k_center(sheet_no: str) -> tuple[float, float]: + """도엽 중심 위경도 (lat, lon).""" + lon_min, lat_min, lon_max, lat_max = sheet5k_to_bounds(sheet_no) + return ((lat_min + lat_max) / 2.0, (lon_min + lon_max) / 2.0) + + +def neighbors_3x3(sheet_no: str) -> list[str]: + """중심 도엽 + 인접 8매 = 3x3 도엽번호 목록 (NW부터 행 우선 9매). + + 구획(1:50,000, 1°) 경계를 넘어가는 인접 도엽은 중심 좌표에 + 도엽 크기만큼 오프셋을 더해 재산출하므로 번호 이월이 자동 처리된다. + """ + lat_c, lon_c = sheet5k_center(sheet_no) + result = [] + for dr in (1, 0, -1): # 북 → 남 + for dc in (-1, 0, 1): # 서 → 동 + result.append( + latlon_to_sheet5k(lat_c + dr * SHEET5K_SIZE_DEG, lon_c + dc * SHEET5K_SIZE_DEG) + ) + return result diff --git a/B04_wf1_Surface/B04_wf1_Surface_Engine_SheetStore.py b/B04_wf1_Surface/B04_wf1_Surface_Engine_SheetStore.py new file mode 100644 index 00000000..46618b3d --- /dev/null +++ b/B04_wf1_Surface/B04_wf1_Surface_Engine_SheetStore.py @@ -0,0 +1,181 @@ +# B04_wf1_Surface_Engine_SheetStore.py +# 1:5,000 수치지형도 도엽 zip 전역 영구저장소 관리. +# +# 조사 결과(2026-07-26): +# - 브이월드 오픈 API에는 도엽 색인·수치지형도 레이어가 없어 도엽번호 산출은 +# 자체 알고리즘(B04_wf1_Surface_Engine_MapSheet) 사용 — 보유 도엽 10매 실측 정합 확인. +# - 브이월드 데이터센터 수치지형도 v2.0(1:5000, dsId=30205)은 시군구 단위 zip이며 +# 다운로드(GET /dtmk/downloadResourceFile.do?ds_id=30205&fileNo=..)는 OIDC SSO +# 로그인 세션 필수(비로그인 시 200/0바이트) → 서버 자동 다운로드 불가. +# - 따라서 취득은 사용자 수동 다운로드 + 업로드(인제스트) 폴백으로 확정하고, +# 확보된 zip은 resources/map_sheets/{도엽번호}.zip 으로 영구 보관한다. + +from __future__ import annotations + +import datetime +import json +import re +import shutil +import tempfile +import zipfile +from pathlib import Path + +from config.config_system import MAP_SHEETS_DIR, MAP_SHEETS_INDEX_PATH + +from .B04_wf1_Surface_Engine_MapSheet import latlon_to_sheet5k, sheet5k_to_bounds + +# 도곽선 레이어 코드 (수치지형도 v2.0 도엽본) +_SHEET_FRAME_CODE = "A0010000" +_SHEET_NO_RE = re.compile(r"(\d{8})") + +# 도엽본 좌표계 후보: 동부원점(5187)이 기본, 오라벨 zip은 중부원점(5186) 선언 사례 실측됨 +_CRS_CANDIDATES = ("EPSG:5187", "EPSG:5186", "EPSG:5185", "EPSG:5188") + +VWORLD_DATASET_URL = "https://www.vworld.kr/dtmk/dtmk_ntads_s002.do?dsId=30205" + + +def _load_index() -> dict: + p = Path(MAP_SHEETS_INDEX_PATH) + if p.exists(): + return json.loads(p.read_text(encoding="utf-8")) + return {"sheets": {}} + + +def _save_index(index: dict) -> None: + Path(MAP_SHEETS_DIR).mkdir(parents=True, exist_ok=True) + Path(MAP_SHEETS_INDEX_PATH).write_text( + json.dumps(index, ensure_ascii=False, indent=2), encoding="utf-8" + ) + + +def get_sheet_path(sheet_no: str) -> Path | None: + """영구저장소에 보관된 도엽 zip 경로. 없으면 None.""" + entry = _load_index()["sheets"].get(str(sheet_no)) + if not entry: + return None + path = Path(MAP_SHEETS_DIR) / entry["file"] + return path if path.exists() else None + + +def missing_sheets(sheet_nos: list[str]) -> list[str]: + """요청 도엽 중 영구저장소에 없는 번호 목록.""" + index = _load_index()["sheets"] + return [ + s + for s in sheet_nos + if str(s) not in index or not (Path(MAP_SHEETS_DIR) / index[str(s)]["file"]).exists() + ] + + +def acquisition_guidance(missing: list[str]) -> dict: + """미보유 도엽 취득 안내 (수동 다운로드 폴백).""" + return { + "missing": missing, + "vworld_dataset_url": VWORLD_DATASET_URL, + "note": ( + "브이월드 데이터센터(수치지형도 v2.0 1:5,000)는 로그인 후 시군구 단위 zip, " + "국토정보플랫폼은 도엽 단위 zip((B010)수치지도_도엽번호_*.zip)을 제공한다. " + "다운로드한 zip을 업로드하면 도곽 검증 후 영구저장소에 등록된다." + ), + } + + +def _read_frame_bounds_wgs84(zip_path: Path) -> tuple[tuple[float, float, float, float], str]: + """zip 내 도곽선(A0010000) SHP bounds를 WGS84로 반환. (bounds, 사용된 CRS). + + 선언된 PRJ가 오라벨(중부/동부원점 뒤바뀜)인 실측 사례가 있어, 파일명 도엽번호 + 기준 예측 도곽과 대조해 정합한 CRS를 자동 판별한다. + """ + import geopandas as gpd + + with zipfile.ZipFile(zip_path) as zf: + names = zf.namelist() + shp = [n for n in names if n.upper().endswith(f"{_SHEET_FRAME_CODE}.SHP")] + if not shp: + raise ValueError(f"{zip_path.name}: 도곽선({_SHEET_FRAME_CODE}) SHP 없음") + base = shp[0][:-4] + members = {} + for ext in (".shp", ".shx", ".dbf", ".prj", ".cpg"): + for n in names: + if n.lower() == (base + ext).lower(): + members[ext] = zf.read(n) + + m = _SHEET_NO_RE.search(zip_path.name) + expected = sheet5k_to_bounds(m.group(1)) if m else None + + with tempfile.TemporaryDirectory() as td: + for ext, data in members.items(): + (Path(td) / f"frame{ext}").write_bytes(data) + shp_path = Path(td) / "frame.shp" + + for crs in _CRS_CANDIDATES: + gdf = gpd.read_file(shp_path) + gdf = gdf.set_crs(crs, allow_override=True) + b = tuple(gdf.to_crs("EPSG:4326").total_bounds) + if expected is None: + return b, crs + tol = 0.002 + if ( + b[0] >= expected[0] - tol + and b[1] >= expected[1] - tol + and b[2] <= expected[2] + tol + and b[3] <= expected[3] + tol + ): + return b, crs + raise ValueError(f"{zip_path.name}: 어떤 좌표계 후보로도 도엽번호 도곽과 정합하지 않음") + + +def ingest_zip(src_path: str | Path, sheet_no: str | None = None) -> dict: + """도엽 zip을 검증 후 영구저장소에 등록한다. + + - 도엽번호: 인자 우선, 없으면 파일명에서 8자리 추출. + - 도곽선 실측 bounds가 도엽번호 예측 도곽 안에 들어야 등록(해안 도엽은 부분 도곽 허용). + - PRJ 오라벨은 CRS 후보 대조로 자동 교정하고 crs 필드에 기록. + """ + src = Path(src_path) + if not src.exists(): + raise FileNotFoundError(str(src)) + + if sheet_no is None: + m = _SHEET_NO_RE.search(src.name) + if not m: + raise ValueError(f"파일명에서 도엽번호(8자리)를 찾을 수 없음: {src.name}") + sheet_no = m.group(1) + sheet_no = str(sheet_no) + + bounds, crs = _read_frame_bounds_wgs84(src) + center_no = latlon_to_sheet5k((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2) + if center_no != sheet_no: + raise ValueError(f"도곽 실측 중심의 도엽번호({center_no})가 요청 번호({sheet_no})와 불일치") + + Path(MAP_SHEETS_DIR).mkdir(parents=True, exist_ok=True) + dest_name = f"{sheet_no}.zip" + dest = Path(MAP_SHEETS_DIR) / dest_name + if src.resolve() != dest.resolve(): + shutil.copy2(src, dest) + + index = _load_index() + index["sheets"][sheet_no] = { + "file": dest_name, + "original_name": src.name, + "crs": crs, + "bounds_wgs84": [round(v, 8) for v in bounds], + "acquired": datetime.date.today().isoformat(), + "source": "vworld/ngii 수동 다운로드", + } + _save_index(index) + return index["sheets"][sheet_no] + + +def ensure_sheets(sheet_nos: list[str]) -> dict: + """도엽 목록의 확보 상태 요약: 보유 경로 + 미보유 취득 안내.""" + available = {} + for s in sheet_nos: + p = get_sheet_path(s) + if p: + available[str(s)] = str(p) + missing = [str(s) for s in sheet_nos if str(s) not in available] + result = {"available": available, "missing": missing} + if missing: + result["guidance"] = acquisition_guidance(missing) + return result diff --git a/config/config_system.py b/config/config_system.py index ce2ca9a7..0418cba3 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -430,6 +430,10 @@ FOREST_ROAD_PROFILE_ALIGNMENT = { # ───────────────────────────────────────────────────────────────────────── STORAGE_BASE_DIR = os.path.join(os.path.dirname(os.path.dirname(__file__)), "storage") +# 1:5,000 수치지형도 도엽 전역 영구저장소 (국가 데이터, 프로젝트 간 공유 — 재다운로드 금지) +MAP_SHEETS_DIR = os.path.join(os.path.dirname(os.path.dirname(__file__)), "resources", "map_sheets") +MAP_SHEETS_INDEX_PATH = os.path.join(MAP_SHEETS_DIR, "map_sheets_index.json") + # 시스템 리소스 로그 (루트 log 폴더에 단일 파일, 1개월 보관) LOG_BASE_DIR = os.path.join(os.path.dirname(os.path.dirname(__file__)), "log") RESOURCE_LOG_PATH = os.path.join(LOG_BASE_DIR, "system_resources.log") diff --git a/scratch/test_mapsheet_verification.py b/scratch/test_mapsheet_verification.py new file mode 100644 index 00000000..be8f1978 --- /dev/null +++ b/scratch/test_mapsheet_verification.py @@ -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() diff --git a/scratch/test_sheetstore_ingest.py b/scratch/test_sheetstore_ingest.py new file mode 100644 index 00000000..100dfd6b --- /dev/null +++ b/scratch/test_sheetstore_ingest.py @@ -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()