"""도엽등고선 3D 서피스 — 1:5,000 수치지형도 등고선으로 DTM 격자를 만든다. LAS 없는 설계(2026-08-30 사용자 확정)의 지형 원천이자, LAS가 있어도 참고용으로 같이 만들어 두는 서피스다. 산출 형식은 LAS 파이프라인의 DTM과 완전히 같게 맞춘다 (`dtm_sheet.npz`: x/y/z/valid_mask) — 종·횡단·배수 세부설계가 쓰는 `build_surface_sampler(models_dir, "sheet", "dtm", smooth=False)`가 무수정으로 돈다. 절취 범위는 노선 XY bbox + `SHEET_SURFACE_MARGIN_M`(300m) 직사각형(사용자 확정), 격자는 `SHEET_SURFACE_GRID_M`(1m). 등고선→정점 구름→Delaunay TIN 보간은 배수유역 엔진(`_Watershed_Grid`)의 검증된 경로를 그대로 쓴다. """ import json import logging import time from pathlib import Path from typing import Any import numpy as np from pyproj import Transformer from B04_PreProcess.B04_PreProcess_Engine_ModelContext import ( atomic_npz, clip_and_compact_mesh, grid_faces, grid_vertices, write_glb, ) from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( build_contour_cloud, grid_spec_from_bounds, interpolate_elevation, ) from config.config_system import ( SHEET_SURFACE_GRID_M, SHEET_SURFACE_MARGIN_M, SURFACE_MAX_PREVIEW_VERTICES, ) logger = logging.getLogger(__name__) # 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값) _CONTOUR_FILE = "도엽_등고선.geojson" # 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다. SHEET_SOURCE_FILTER = "sheet" def _load_contour_features_metric(processed_dir: Path, epsg: int) -> list[dict[str, Any]]: """병합 도엽 등고선(WGS84)을 읽어 사업지 CRS(m)로 재투영한다.""" path = processed_dir / _CONTOUR_FILE if not path.exists(): logger.warning("도엽 서피스: 등고선 파일이 없습니다: %s", path) return [] try: with path.open("r", encoding="utf-8") as file: data = json.load(file) except (OSError, json.JSONDecodeError): logger.warning("도엽 서피스: 등고선 GeoJSON을 읽지 못했습니다: %s", path) return [] features = data.get("features") if not isinstance(features, list): return [] transformer = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True) def _map(coords: Any) -> Any: if not isinstance(coords, list): return coords if coords and isinstance(coords[0], (int, float)): x, y = transformer.transform(coords[0], coords[1]) return [x, y, *coords[2:]] return [_map(item) for item in coords] converted: list[dict[str, Any]] = [] for feature in features: geometry = feature.get("geometry") or {} coordinates = _map(geometry.get("coordinates")) if coordinates is None: continue converted.append( { "type": "Feature", "properties": feature.get("properties") or {}, "geometry": {"type": geometry.get("type"), "coordinates": coordinates}, } ) return converted def _preview_mesh( x: np.ndarray, y: np.ndarray, z: np.ndarray, valid: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: """프리뷰용 정점·면 — 정점 수가 상한을 넘으면 격자를 성기게 딴다.""" stride = 1 while (len(x) // stride + 1) * (len(y) // stride + 1) > SURFACE_MAX_PREVIEW_VERTICES: stride += 1 px, py = x[::stride], y[::stride] pz, pv = z[::stride, ::stride], valid[::stride, ::stride] vertices = grid_vertices(px, py, pz.astype(np.float64)) faces = grid_faces(len(py), len(px)) return clip_and_compact_mesh(vertices, faces, pv.reshape(-1)) def build_sheet_surface_model( project_root: Path, processed_dir: Path, models_dir: Path, route_xy: np.ndarray, epsg: int, ) -> dict[str, Any] | None: """도엽등고선으로 DTM npz·프리뷰 glb를 만들고 surface_models 등록용 dict를 돌려준다. `route_xy`: (N, 2) 노선 정점 XY(사업지 CRS, m). 실패하면 None — 호출측은 분석을 계속한다(도엽 미확보 지역 폴백). """ started = time.monotonic() features = _load_contour_features_metric(processed_dir, epsg) if not features: return None x_min = float(np.min(route_xy[:, 0])) - SHEET_SURFACE_MARGIN_M x_max = float(np.max(route_xy[:, 0])) + SHEET_SURFACE_MARGIN_M y_min = float(np.min(route_xy[:, 1])) - SHEET_SURFACE_MARGIN_M y_max = float(np.max(route_xy[:, 1])) + SHEET_SURFACE_MARGIN_M # 저지대 제거(floor) 없이 전부 쓴다 — 종·횡단은 낮은 지반도 필요하다. cloud = build_contour_cloud(features, None, (x_min, y_min, x_max, y_max)) if cloud.is_empty: logger.warning("도엽 서피스: 절취 범위 안에 등고선 정점이 없습니다.") return None spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M) surface = interpolate_elevation(spec, cloud) # (R, C), 북→남 행 순서, 외부 NaN # DtmGridSampler 규약에 맞춰 y 오름차순으로 뒤집어 저장한다. x_coords = spec.cell_centers_x() y_coords = spec.cell_centers_y()[::-1] z_grid = surface[::-1, :].astype(np.float32) valid_grid = np.isfinite(z_grid) if not valid_grid.any(): logger.warning("도엽 서피스: 유효 표고 셀이 없습니다.") return None stem = f"dtm_{SHEET_SOURCE_FILTER}" model_path = models_dir / f"{stem}.npz" preview_path = models_dir / f"{stem}_preview.glb" finite_z = z_grid[valid_grid] # bounds를 npz에 같이 넣는다 — 등고선 API가 이 값을 화면 원점으로 쓴다. 없으면 # LAS structured.npz로 폴백해 메시(glb) 원점과 어긋난다(LAS 없는 설계는 아예 실패). bounds = np.array( [ [x_coords[0], x_coords[-1]], [y_coords[0], y_coords[-1]], [float(finite_z.min()), float(finite_z.max())], ] ) atomic_npz( model_path, x=x_coords, y=y_coords, z=z_grid, valid_mask=valid_grid, bounds=bounds, resolution=np.array([SHEET_SURFACE_GRID_M], np.float32), ) vertices, faces = _preview_mesh(x_coords, y_coords, z_grid, valid_grid) write_glb(preview_path, vertices, faces, bounds) logger.info( "도엽 서피스 생성 완료: %d×%d 격자, 정점 %d개 (%.1fs)", spec.n_rows, spec.n_cols, cloud.xy.shape[0], time.monotonic() - started, ) return { "model_type": "dtm", "source_filter": SHEET_SOURCE_FILTER, "representation": "regular_grid", "model_file_path": str(model_path.relative_to(project_root)).replace("\\", "/"), "resolution_m": SHEET_SURFACE_GRID_M, "generation_params": { "source_filter": SHEET_SOURCE_FILTER, "representation": "regular_grid", "source": "map_sheet_contours", "margin_m": SHEET_SURFACE_MARGIN_M, }, "layers": [ { "layer_name": f"dtm_{SHEET_SOURCE_FILTER}_preview", "geometry_type": "MESH", "file_path": str(preview_path.relative_to(project_root)).replace("\\", "/"), "file_format": "glb", } ], } def build_sheet_surface_from_route( project_root: Path, processed_dir: Path, models_dir: Path ) -> dict[str, Any] | None: """B03 업로드 계획노선 CSV를 찾아 도엽 서피스를 만든다. 없거나 실패하면 None.""" from common_util.common_util_route_geometry import ( find_planned_route_file, read_planned_route_csv, ) route_file = find_planned_route_file(project_root / "B03_FileInput" / "input") if route_file is None: logger.warning("도엽 서피스: 계획 노선 파일이 없습니다.") return None planned = read_planned_route_csv(route_file) if planned is None or len(planned.vertices) < 2: logger.warning("도엽 서피스: 계획 노선 파일을 읽지 못했습니다: %s", route_file.name) return None route_xy = np.array([(v.x, v.y) for v in planned.vertices], dtype=np.float64) return build_sheet_surface_model( project_root, processed_dir, models_dir, route_xy, planned.epsg or 5186 ) def run_sheet_surface_analysis( project_root: Path, route_csv_path: Path, *, on_progress: Any = None, ) -> dict[str, Any]: """LAS 없는 WF1 — 도엽 확보 후 도엽등고선 서피스만으로 분석 결과를 만든다. 반환 형식은 `run_surface_analysis()`와 같다(save_surface_analysis_to_db 호환). """ from common_util.common_util_route_geometry import read_planned_route_csv def _report(percent: int, stage: str, message: str) -> None: if on_progress is not None: on_progress(percent, stage, message) stage_root = project_root / "B04_PreProcess" processed_dir = stage_root / "processed" models_dir = stage_root / "models" processed_dir.mkdir(parents=True, exist_ok=True) models_dir.mkdir(parents=True, exist_ok=True) planned = read_planned_route_csv(route_csv_path) if planned is None or len(planned.vertices) < 2: raise ValueError(f"계획 노선 파일을 읽지 못했습니다: {route_csv_path.name}") epsg = planned.epsg or 5186 route_xy = np.array([(v.x, v.y) for v in planned.vertices], dtype=np.float64) bounds_dict = { "x": [float(route_xy[:, 0].min()), float(route_xy[:, 0].max())], "y": [float(route_xy[:, 1].min()), float(route_xy[:, 1].max())], "z": [ float(min(v.z for v in planned.vertices)), float(max(v.z for v in planned.vertices)), ], } # VWorld 지도·GIS 벡터·도엽 확보 — LAS 경로와 같은 공용 블록 (지연 import로 순환 회피) _report(30, "download_maps", "VWorld 지도 및 수치지형도 도엽 확보 중") from B04_PreProcess.B04_PreProcess_Engine import download_geodata download_geodata( project_root, processed_dir, bounds_dict, route_csv_path.parent, rebuild=False, default_epsg=f"EPSG:{epsg}", report=_report, ) _report(70, "surface_model", "도엽등고선 3D 서피스 생성 중") model = build_sheet_surface_model(project_root, processed_dir, models_dir, route_xy, epsg) if model is None: raise ValueError("도엽등고선으로 지표면을 만들지 못했습니다 — 도엽 확보를 확인하세요.") _report(95, "saving", "결과 저장 중") return { "processed": { "processed_file_path": str( (processed_dir / _CONTOUR_FILE).relative_to(project_root) ).replace("\\", "/"), "converted_file_path": None, "point_count": int(len(route_xy)), "bounds": { "x_min": bounds_dict["x"][0], "x_max": bounds_dict["x"][1], "y_min": bounds_dict["y"][0], "y_max": bounds_dict["y"][1], }, "statistics": { "min_z": bounds_dict["z"][0], "max_z": bounds_dict["z"][1], "mean_z": None, }, }, "ground_summary": {}, "manifest": {"status": "sheet_only"}, "models": [model], }