"""D 예정노선 지형 받기 — 사업 전체 한 벌(도엽 · 위성 · DEM) (요구 345 · 389 · PLAN 530-1 · 561-2). 지형은 D 프로젝트 storage 안 `B04_PreProcess/` 배치 그대로(지형 코드 주인 = B04 · 고칠 곳 0): - 범위 = 전체 시점 · 종점(+ 그은 점) 상자 + 여유 · 받은 범위가 있으면 그 범위와 합침(줄지 않음) - 도엽 = 범위 상자에 닿는 도엽 모두 → `ensure_sheets`(없는 것만 · prune 없음) - 병합 = 지난번 병합 도엽 ∪ 이번 도엽(옛 도엽 유지 · 열쇠 같으면 건너뜀) - 위성 = 범위 도엽 도곽 ∪ 옛 위성 범위(덮으면 건너뜀) - DEM = 한 벌 `models_plan/`(연도는 D03 나눔 뒤에 정해짐) · 받은 범위 = `plan_terrain.json` - 작업 좌표계 = 지점 epsg · `B03_FileInput/input/prj/plan.prj` 한 장(B04 · B05 좌표계 창구가 읽음) - 받기 전 범위 보기 = `plan_extent`(내려받기 없음) · 받기는 사용자가 누를 때만 계산 자리 = ③ 서버 단독. DB 행(`surface_models`)은 라우터가 씀. """ from __future__ import annotations import json import logging import math from pathlib import Path from typing import Any, Callable import numpy as np from B04_PreProcess.B04_PreProcess_Engine_MapSheet import ( SHEET5K_SIZE_DEG, sheet5k_to_bounds, sheets_for_points, ) from B04_PreProcess.B04_PreProcess_Engine_SheetParser import PARSER_VERSION, parse_and_merge_sheets from B04_PreProcess.B04_PreProcess_Engine_SheetStore import ( ensure_sheets, get_project_map_sheets_dir, get_sheet_path, ) from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import build_sheet_surface_model from B04_PreProcess.B04_PreProcess_Engine_VWorld import download_vworld_satellite_map from common_util.common_util_json import atomic_write_json from config.config_system_terrain import ( SHEET_SURFACE_DEFAULT_METHOD, SHEET_SURFACE_GRID_M, SHEET_SURFACE_MARGIN_M, ) logger = logging.getLogger(__name__) # 탐색 복도 반폭(D 만 · B 기본 100 m) — 실측(530-2): 100 · 200 m 는 못 지킨 곳 남음 · 300 m 0~1 PLAN_CORRIDOR_M = 300.0 # DEM 이 복도 밖으로 더 덮는 폭(복도 끝 칸 기울기 · 횡단 표본 ±20 m) DEM_EXTRA_M = 50.0 # 여유 처음 값 = 복도 + 더 덮는 폭(옛 연도 구간 DEM 과 같은 폭) · 바닥 = B 절취 여유 DEFAULT_MARGIN_M = PLAN_CORRIDOR_M + DEM_EXTRA_M MAX_MARGIN_M = 5000.0 STAGE = Path("B04_PreProcess") PLAN_MODELS = "models_plan" TERRAIN_META = "plan_terrain.json" PRJ_NAME = "plan.prj" Progress = Callable[[int, str], None] Box = tuple[float, float, float, float] # x0, y0, x1, y1(작업 좌표계 m) XY = tuple[float, float] def models_dir_of(project_root: Path, year: int | None = None) -> Path: """사업 전체 DEM 자리 — `year` 는 옛 연도 분석 부르는 꼴 맞춤(안 씀).""" return project_root / STAGE / PLAN_MODELS def clamp_margin(margin_m: float | None) -> float: """여유 — 빈 값 = 처음 값 · B 절취 여유 ~ 5 km.""" if margin_m is None or not math.isfinite(margin_m): return DEFAULT_MARGIN_M return min(MAX_MARGIN_M, max(SHEET_SURFACE_MARGIN_M, float(margin_m))) def received_box(project_root: Path) -> Box | None: """받은 DEM 범위 — 받은 적 없거나 읽지 못하면 None.""" meta = _read_json(models_dir_of(project_root) / TERRAIN_META) or {} try: x0, y0, x1, y1 = (float(v) for v in meta["box"]) except (KeyError, TypeError, ValueError): return None return x0, y0, x1, y1 def wanted_box(points: list[XY], margin_m: float, have: Box | None) -> Box: """받을 범위 = 점 상자 + 여유 ∪ 받은 범위.""" xs, ys = [p[0] for p in points], [p[1] for p in points] box = (min(xs) - margin_m, min(ys) - margin_m, max(xs) + margin_m, max(ys) + margin_m) if have is None: return box return min(box[0], have[0]), min(box[1], have[1]), max(box[2], have[2]), max(box[3], have[3]) def covers(have: Box | None, box: Box) -> bool: return ( have is not None and all(have[i] <= box[i] for i in (0, 1)) and all(have[i] >= box[i] for i in (2, 3)) ) def box_ring(box: Box, epsg: int) -> list[list[float]]: """상자 네 귀 경위도([경도, 위도] · 왼아래부터 반시계).""" from common_util.common_util_crs import tm_to_lonlat x0, y0, x1, y1 = box return [list(tm_to_lonlat(x, y, epsg)) for x, y in ((x0, y0), (x1, y0), (x1, y1), (x0, y1))] def sheets_for_box(box: Box, epsg: int) -> list[str]: """상자에 닿는 도엽 — 경위도 외접 상자를 도엽 반 칸 간격으로 짚음.""" ring = box_ring(box, epsg) lons, lats = [p[0] for p in ring], [p[1] for p in ring] step = SHEET5K_SIZE_DEG / 2 def ticks(lo: float, hi: float) -> list[float]: return [lo + step * i for i in range(int((hi - lo) / step) + 1)] + [hi] grid = [ (lat, lon) for lat in ticks(min(lats), max(lats)) for lon in ticks(min(lons), max(lons)) ] return sorted(set(sheets_for_points(grid))) def plan_extent( project_root: Path, points: list[XY], epsg: int, margin_m: float | None = None ) -> dict[str, Any]: """받기 전 범위 보기 — 내려받기 없음. `points` = 전체 시점 · 종점(+ 그은 점 · 작업 좌표계 m).""" margin = clamp_margin(margin_m) have = received_box(project_root) box = wanted_box(points, margin, have) store = get_project_map_sheets_dir(project_root) sheets = [] for no in sheets_for_box(box, epsg): lon0, lat0, lon1, lat1 = sheet5k_to_bounds(no) sheets.append( { "no": no, "ring": [[lon0, lat0], [lon1, lat0], [lon1, lat1], [lon0, lat1]], "stored": get_sheet_path(store, no) is not None, } ) return { "epsg": int(epsg), "margin_m": margin, "box_m": list(box), "box": box_ring(box, epsg), "have_m": list(have) if have else None, "have": box_ring(have, epsg) if have else None, "covered": covers(have, box), "sheets": sheets, # DEM 칸 수 어림 — 받기 시간 · 메모리 가늠 "cells": int((box[2] - box[0]) * (box[3] - box[1]) / SHEET_SURFACE_GRID_M**2), } def write_plan_prj(project_root: Path, epsg: int) -> Path: """작업 좌표계 .prj — 내용이 같으면 그대로.""" from pyproj import CRS path = project_root / "B03_FileInput" / "input" / "prj" / PRJ_NAME text = CRS.from_epsg(int(epsg)).to_wkt("WKT1_ESRI") if not path.is_file() or path.read_text(encoding="utf-8", errors="ignore") != text: path.parent.mkdir(parents=True, exist_ok=True) path.write_text(text, encoding="utf-8") return path def _read_json(path: Path) -> dict[str, Any] | None: try: return json.loads(path.read_text(encoding="utf-8")) except (OSError, ValueError): return None def _merge_sheets(processed: Path, store: Path, available: list[str]) -> list[str]: """병합 = 지난번 병합 도엽(zip 남은 것) ∪ 이번 도엽 · 열쇠 같으면 건너뜀.""" key_path = processed / "sheet_merge.key.json" old = (_read_json(key_path) or {}).get("sheets") or [] sheets = sorted({*available, *(s for s in old if get_sheet_path(store, s) is not None)}) key = {"sheets": sheets, "parser": PARSER_VERSION} if sheets and (_read_json(key_path) != key or not any(processed.glob("도엽_*.geojson"))): parse_and_merge_sheets(store, sheets, processed) atomic_write_json(key_path, key) return sheets def _satellite(processed: Path, prj: Path, points: list[tuple[float, float]], epsg: int) -> str: """위성 = 기준 도엽 도곽 ∪ 옛 범위 · 이미 덮으면 건너뜀. 돌려줌 = skip · done · failed.""" from pyproj import Transformer to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True) xs: list[float] = [] ys: list[float] = [] for sheet in sheets_for_points(points): lon0, lat0, lon1, lat1 = sheet5k_to_bounds(sheet) for lon, lat in ((lon0, lat0), (lon0, lat1), (lon1, lat0), (lon1, lat1)): x, y = to_metric.transform(lon, lat) xs.append(x) ys.append(y) need = {"x": [min(xs), max(xs)], "y": [min(ys), max(ys)]} meta = _read_json(processed / "vworld_satellite_meta.json") try: old = [float(meta[k]) for k in ("x_min", "x_max", "y_min", "y_max")] if meta else None except (KeyError, TypeError, ValueError): old = None if old and old[0] <= need["x"][0] and old[1] >= need["x"][1]: if old[2] <= need["y"][0] and old[3] >= need["y"][1]: return "skip" if old: need = { "x": [min(old[0], need["x"][0]), max(old[1], need["x"][1])], "y": [min(old[2], need["y"][0]), max(old[3], need["y"][1])], } try: done = download_vworld_satellite_map(prj, need, processed, layer_name="Satellite") except Exception as exc: # noqa: BLE001 — 위성은 배경 · 지형 받기를 막지 않음 logger.warning("D 위성 지도 받기 실패: %s", exc) return "failed" return "done" if done else "failed" def ensure_plan_terrain( project_root: Path, points: list[XY], epsg: int, margin_m: float | None = None, on_progress: Progress | None = None, ) -> dict[str, Any]: """사업 전체 지형 받기 — `points` = 전체 시점 · 종점(+ 그은 점 · 작업 좌표계 m). 돌려줌 = 도엽(받은 · 실패 · 병합) · 위성 상태 · 받은 범위 · 모델 등록용 dict 목록. """ def report(percent: int, message: str) -> None: if on_progress is not None: on_progress(percent, message) if not points: raise ValueError("시점 · 종점이 없습니다.") processed = project_root / STAGE / "processed" processed.mkdir(parents=True, exist_ok=True) prj = write_plan_prj(project_root, epsg) margin = clamp_margin(margin_m) box = wanted_box(points, margin, received_box(project_root)) report(5, "수치지형도 도엽 고르는 중") wanted = sheets_for_box(box, epsg) store = get_project_map_sheets_dir(project_root) before = {s for s in wanted if get_sheet_path(store, s) is not None} report(10, f"수치지형도 도엽 {len(wanted) - len(before)}매 받는 중") got = ensure_sheets(store, wanted) if not got["available"]: raise ValueError(f"수치지형도 도엽을 하나도 받지 못했습니다: {got['failed']}") report(45, "도엽 레이어 병합 중") merged = _merge_sheets(processed, store, list(got["available"])) report(60, "위성 지도 받는 중") satellite = _satellite(processed, prj, [(p[1], p[0]) for p in box_ring(box, epsg)], epsg) report(65, "사업 지형면 만드는 중") models_dir = models_dir_of(project_root) models_dir.mkdir(parents=True, exist_ok=True) # B 절취 = 넘긴 꼭짓점 상자 + `SHEET_SURFACE_MARGIN_M` → 그만큼 안쪽으로 넘김(여유 ≥ 그 값) inner = SHEET_SURFACE_MARGIN_M outline = np.array( [[box[0] + inner, box[1] + inner], [box[2] - inner, box[3] - inner]], dtype=np.float64 ) models = build_sheet_surface_model( project_root, processed, models_dir, outline, f"EPSG:{int(epsg)}", [SHEET_SURFACE_DEFAULT_METHOD], ) if not models: raise ValueError( "사업 지형면을 만들지 못했습니다 — 범위 안 도엽 등고선 높이 확인" f"(병합 도엽 {len(merged)}매)" ) atomic_write_json(models_dir / TERRAIN_META, {"box": list(box), "margin_m": margin}) report(100, "지형 받기 끝") return { "sheets": { "wanted": len(wanted), "downloaded": sorted(set(got["available"]) - before), "failed": got["failed"], "merged": merged, }, "satellite": satellite, "box_m": list(box), "models": models, }