"""B04 지표면 분석 엔진 오케스트레이터. 원본 LAS를 구조화하고 지면 필터를 실행한 뒤 지표면 5종 모델을 빌드하는 동기 계산 파이프라인. 라우터에서 asyncio.to_thread로 호출한다. """ import json import logging import time from collections.abc import Callable from pathlib import Path from typing import Any import numpy as np from B04_PreProcess.B04_PreProcess_Engine_Ground import ( build_ground_masks, detect_extra_filters, run_ground_filter, summarize_masks, ) from B04_PreProcess.B04_PreProcess_Engine_Pipeline import build_all_terrain_models from B04_PreProcess.B04_PreProcess_Engine_Structurize import structurize_las from common_util.common_util_atomic import atomic_write_npz from common_util.common_util_json import atomic_write_json from config.config_system import SURFACE_GROUND_RATIO_WARN, build_surface_model_config logger = logging.getLogger(__name__) # 진행 콜백 시그니처: (진행률 0~100, 현재 단계 키, 메시지) ProgressCallback = Callable[[int, str, str], None] GROUND_POINT_SAMPLE_LIMIT = 500_000 GROUND_POINT_CACHE_VERSION = 2 def _source_identity(las_path: Path) -> dict[str, Any]: """입력 LAS의 정체성(이름·크기·수정시각)으로 캐시 세대를 식별한다 (PLAN B-2).""" stat = las_path.stat() return { "filename": las_path.name, "size_bytes": int(stat.st_size), "mtime": float(stat.st_mtime), } def _relative_to_project(project_root: Path, path: Path) -> str: """프로젝트 루트 기준 posix 상대 경로 문자열.""" return path.relative_to(project_root).as_posix() def cache_ground_points( structured_path: Path, filter_key: str, mask: np.ndarray | None = None, ) -> Path: """필터링된 지면 포인트 미리보기 캐시를 생성하고 경로를 반환한다. mask 미전달 시 영구 저장된 mask_{filter}.npy를 우선 재사용한다 (PLAN A-2). """ with np.load(structured_path) as structured: xyz = np.asarray(structured["xyz"], dtype=np.float32) mask_path = structured_path.parent / f"mask_{filter_key}.npy" if mask is None and mask_path.is_file(): stored_mask = np.load(mask_path) if len(stored_mask) == len(xyz): mask = np.asarray(stored_mask, dtype=bool) if mask is None: mask = build_ground_masks(structured, [filter_key])[filter_key] np.save(mask_path, np.asarray(mask, dtype=bool)) ground_indexes = np.flatnonzero(mask) ground_points = xyz[ground_indexes] ground_point_count = int(len(ground_points)) if ground_point_count: data_bounds = np.column_stack((ground_points.min(axis=0), ground_points.max(axis=0))) else: data_bounds = np.zeros((3, 2), dtype=np.float64) if ground_point_count > GROUND_POINT_SAMPLE_LIMIT: rng = np.random.default_rng(20260717) sample_indexes = rng.choice( ground_point_count, GROUND_POINT_SAMPLE_LIMIT, replace=False ) ground_indexes = ground_indexes[sample_indexes] points = ground_points[sample_indexes] else: points = ground_points arrays = { "xyz": points, "bounds": np.asarray(structured["bounds"], dtype=np.float64), "data_bounds": np.asarray(data_bounds, dtype=np.float64), "cache_version": np.asarray(GROUND_POINT_CACHE_VERSION, dtype=np.int16), "point_count": np.asarray(ground_point_count, dtype=np.int64), "sampled_count": np.asarray(len(points), dtype=np.int64), } if "rgb" in structured: rgb = np.asarray(structured["rgb"]) if rgb.ndim > 0 and len(rgb) == len(xyz): arrays["rgb"] = rgb[ground_indexes] cache_path = structured_path.parent / f"ground_points_{filter_key}.npz" atomic_write_npz(cache_path, **arrays) return cache_path def run_surface_analysis( project_root: Path, las_path: Path, *, source_filters: list[str], methods: list[str], force: bool = False, on_progress: ProgressCallback | None = None, ) -> dict[str, Any]: """구조화→필터→모델 빌드를 수행하고 산출 메타데이터를 반환한다. 반환 dict: - processed: {processed_file_path, converted_file_path, point_count, bounds, statistics} - ground_summary: 필터별 지면 포인트 요약 - manifest: 지표면 모델 파이프라인 manifest - models: [{model_type, model_file_path, resolution_m, generation_params, layers}] """ def _report(percent: int, stage: str, message: str) -> None: if on_progress is not None: on_progress(percent, stage, message) total_started = time.monotonic() 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) # 0. 입력 세대 검증: LAS가 바뀌었으면 모든 캐시를 재계산한다 (PLAN B-2) identity_path = processed_dir / "source_identity.json" current_identity = _source_identity(las_path) stored_identity: dict[str, Any] | None = None if identity_path.is_file(): try: stored_identity = json.loads(identity_path.read_text(encoding="utf-8")) except (OSError, ValueError): stored_identity = None rebuild = force or stored_identity != current_identity # 1. LAS 구조화 (structured.npz) — 캐시가 유효하면 스킵 (PLAN B-1) structured_path = processed_dir / "structured.npz" if rebuild or not structured_path.is_file(): _report(10, "structurize", "LAS 구조화 중") step_started = time.monotonic() structured_path = structurize_las(las_path, processed_dir) atomic_write_json(identity_path, current_identity) logger.info( "B04 LAS 구조화 완료: %s (%.1fs)", las_path.name, time.monotonic() - step_started ) else: _report(10, "structurize", "구조화 캐시 재사용") logger.info("B04 구조화 캐시 재사용: %s", structured_path.name) with np.load(structured_path) as structured: xyz = structured["xyz"] bounds = structured["bounds"] total_points = int(len(xyz)) stats = { "min_z": float(bounds[2, 0]), "max_z": float(bounds[2, 1]), "mean_z": float(np.mean(xyz[:, 2])) if total_points else None, } bounds_dict = { "x_min": float(bounds[0, 0]), "x_max": float(bounds[0, 1]), "y_min": float(bounds[1, 0]), "y_max": float(bounds[1, 1]), } data = {"xyz": xyz, "bounds": bounds, "classification": structured["classification"]} # 2. 지면 필터 실행 — mask_{filter}.npy 영구 캐시 우선 재사용 (PLAN A-1) _report(40, "ground_filter", "지면 필터 적용 중") # LAS가 지면분류를 싣고 왔으면 필터 하나로 더 쓴다 — 미분류 LAS면 조용히 빠진다. extra_filters = detect_extra_filters(data, source_filters) if extra_filters: source_filters = list(source_filters) + extra_filters logger.info("B04 LAS 자체 지면분류 감지 — 필터 추가: %s", ", ".join(extra_filters)) masks: dict[str, np.ndarray] = {} for filter_key in source_filters: mask_path = processed_dir / f"mask_{filter_key}.npy" mask: np.ndarray | None = None if not rebuild and mask_path.is_file(): stored_mask = np.load(mask_path) if len(stored_mask) == total_points: mask = np.asarray(stored_mask, dtype=bool) logger.info("B04 지면 필터 캐시 재사용: %s", mask_path.name) if mask is None: step_started = time.monotonic() mask = np.asarray(run_ground_filter(filter_key, data), dtype=bool) np.save(mask_path, mask) logger.info( "B04 지면 필터 계산 완료: %s (%.1fs)", filter_key, time.monotonic() - step_started, ) masks[filter_key] = mask ground_summary = summarize_masks(data, masks) # 지면점 비율은 필터가 조용히 실패해도 유일하게 드러나는 신호다 — 반드시 남긴다. for filter_key, entry in ground_summary.items(): ratio = float(entry["ground_ratio"]) log = logger.warning if ratio < SURFACE_GROUND_RATIO_WARN else logger.info log( "B04 지면점 %s: %d / %d (%.2f%%)", filter_key, entry["ground_point_count"], entry["total_point_count"], ratio * 100, ) for filter_key, mask in masks.items(): cache_path = processed_dir / f"ground_points_{filter_key}.npz" if not rebuild and cache_path.is_file(): with np.load(cache_path) as cached: if ( "cache_version" in cached and int(cached["cache_version"]) == GROUND_POINT_CACHE_VERSION ): continue cache_ground_points(structured_path, filter_key, mask) # 3. 지표면 5종 모델 빌드 _report(70, "surface_model", "지표면 모델 생성 중") config = build_surface_model_config() config["source_filters"] = list(source_filters) config["precompute"] = list(methods) step_started = time.monotonic() def _model_progress(percent: int, detail: str) -> None: _report(70 + int(percent * 0.2), "surface_model", detail) manifest = build_all_terrain_models( data, masks, models_dir, config, force=rebuild, progress=_model_progress ) logger.info( "B04 지표면 모델 빌드 완료: status=%s (%.1fs)", manifest.get("status"), time.monotonic() - step_started, ) # 3-2·3-3. VWorld 지도·국가 GIS 벡터·수치지형도 도엽 (공용 블록 — 도엽 서피스도 사용) _report(90, "download_maps", "VWorld 지도 및 GIS 벡터 데이터 다운로드 중") las_bounds_dict = { "x": [float(bounds[0, 0]), float(bounds[0, 1])], "y": [float(bounds[1, 0]), float(bounds[1, 1])], "z": [float(bounds[2, 0]), float(bounds[2, 1])], } download_geodata( project_root, processed_dir, las_bounds_dict, las_path.parent, rebuild, report=_report ) # 3-4. 도엽등고선 3D 서피스 — LAS가 있어도 참고용으로 같이 만들어 영구저장한다 # (2026-08-30 사용자 확정). 실패해도 분석은 계속한다. _report(94, "surface_model", "도엽등고선 3D 서피스 생성 중") sheet_models: list[dict[str, Any]] = [] try: from B04_PreProcess.B04_PreProcess_Engine_SheetSurface import ( build_sheet_surface_from_route, ) sheet_models = build_sheet_surface_from_route(project_root, processed_dir, models_dir) except Exception as exc: logger.warning("도엽등고선 서피스 생성 실패: %s", exc) _report(95, "saving", "결과 저장 중") return _collect_analysis_result( project_root, models_dir, structured_path, bounds_dict, stats, total_points, ground_summary, manifest, sheet_models, total_started, ) def download_geodata( project_root: Path, processed_dir: Path, las_bounds_dict: dict[str, list[float]], prj_search_dir: Path, rebuild: bool, *, default_epsg: str = "EPSG:5186", report: Any = None, ) -> None: """VWorld 지도·국가 GIS 벡터·수치지형도 도엽 확보 (공용 블록). LAS 분석(run_surface_analysis)과 LAS 없는 도엽 서피스 분석이 같이 쓴다. 실패해도 예외를 밖으로 던지지 않는다 — 분석 본체를 막지 않는다. """ try: # 입력 파일과 같은 폴더의 .prj를 우선 사용 (PLAN E-1: 전체 재귀 glob 제거) prj_candidates = sorted(prj_search_dir.glob("*.prj")) or sorted( project_root.glob("B03_FileInput/**/*.prj") ) prj_path = prj_candidates[0] if prj_candidates else project_root / "result.prj" from B04_PreProcess.B04_PreProcess_Engine_Extent import ( download_extent, map_meta_covers, satellite_extent, ) from B04_PreProcess.B04_PreProcess_Engine_GisVector import download_all_gis_vectors from B04_PreProcess.B04_PreProcess_Engine_VWorld import ( download_vworld_satellite_map, get_epsg_from_prj, ) project_epsg = default_epsg if prj_path.exists(): project_epsg = get_epsg_from_prj(prj_path.read_text(encoding="utf-8", errors="ignore")) # 국가 GIS 벡터는 라이다∪계획노선 범위로 받는다. bounds_dict_for_download = download_extent(project_root, las_bounds_dict, project_epsg) # 배경 지도(위성·하이브리드·백지도)는 계획노선이 걸치는 기준 도엽의 도곽 범위로 받는다 # — 수치지형도 도엽과 같은 눈금. 주변 도엽은 받지 않는다(2026-08-01 사용자 지시). map_bounds_for_download = satellite_extent(project_root, las_bounds_dict, project_epsg) # VWorld 지도 및 GIS 데이터 저장 위치는 B04_PreProcess/processed에 보관. layers = [ {"layer": "Satellite", "ext": "jpeg"}, {"layer": "Hybrid", "ext": "png"}, {"layer": "white", "ext": "png"}, ] for item in layers: meta_path = processed_dir / f"vworld_{item['layer'].lower()}_meta.json" # 파일이 있어도 계획노선·여유 셀이 바뀌어 범위를 못 덮으면 다시 받는다. # (B03 업로드가 부르는 경로는 rebuild=False라, 존재 여부만 보면 영영 갱신되지 않는다.) if not rebuild and map_meta_covers(meta_path, map_bounds_for_download): continue try: step_started = time.monotonic() download_vworld_satellite_map( prj_path, map_bounds_for_download, processed_dir, layer_name=item["layer"], ext=item["ext"], ) logger.info( "B04 VWorld %s 지도 다운로드 완료 (%.1fs)", item["layer"], time.monotonic() - step_started, ) except Exception as exc: logger.warning("B04 VWorld %s 지도 다운로드 실패: %s", item["layer"], exc) # 등고선은 국가 gpkg 추가 시점이 전처리보다 늦을 수 있으므로 개별로 존재 여부를 확인한다. contour_geojson = processed_dir / "등고선_bounds.geojson" if ( rebuild or not any(processed_dir.glob("*_bounds.geojson")) or not contour_geojson.exists() ): try: step_started = time.monotonic() download_all_gis_vectors(prj_path, bounds_dict_for_download, processed_dir) logger.info( "B04 국가 GIS 벡터 다운로드 완료 (%.1fs)", time.monotonic() - step_started ) except Exception as exc: logger.warning("B04 국가 GIS 벡터 다운로드 실패: %s", exc) # 3-3. 1:5,000 수치지형도 도엽 확보 → 프로젝트 영구저장소 # 기준은 계획노선 시점·종점 (같은 도엽이면 9매, 이웃 도엽에 걸치면 12매). # (실패해도 분석은 계속 — 폴백은 수동 다운로드 + 인제스트) if report is not None: report(92, "download_maps", "수치지형도 도엽 확보 중") try: from B04_PreProcess.B04_PreProcess_Engine_Extent import sheet_reference_points_wgs84 from B04_PreProcess.B04_PreProcess_Engine_MapSheet import neighbors_for_points from B04_PreProcess.B04_PreProcess_Engine_SheetStore import ( ensure_sheets, get_project_map_sheets_dir, prune_sheets, ) step_started = time.monotonic() # 도엽 기준은 계획노선 시점·종점 — 노선이 두 도엽에 걸치면 양쪽 주변까지 확보한다. # 계획노선이 없으면 라이다 범위 중심으로 되돌아간다(2026-08-01 사용자 지시). reference_points = sheet_reference_points_wgs84( project_root, las_bounds_dict, project_epsg ) sheet_grid = neighbors_for_points(reference_points) sheet_store = get_project_map_sheets_dir(project_root) sheet_result = ensure_sheets(sheet_store, sheet_grid) if sheet_result["failed"]: logger.warning("B04 수치지형도 도엽 일부 확보 실패: %s", sheet_result["failed"]) logger.info( "B04 수치지형도 도엽 확보 완료: %d매 (%.1fs)", len(sheet_result["available"]), time.monotonic() - step_started, ) # 선정에서 빠진 zip은 쓰이지 않으므로 정리한다(다른 지역 잔재 포함). prune_sheets(sheet_store, sheet_grid) # 확보된 도엽을 레이어별 병합 GeoJSON으로 산출 (기존 산출물 있으면 스킵) if sheet_result["available"] and ( rebuild or not any(processed_dir.glob("도엽_*.geojson")) ): from B04_PreProcess.B04_PreProcess_Engine_SheetParser import ( parse_and_merge_sheets, ) step_started = time.monotonic() merged = parse_and_merge_sheets( sheet_store, list(sheet_result["available"]), processed_dir ) logger.info( "B04 도엽 레이어 병합 완료: %s (%.1fs)", {k: v["count"] for k, v in merged.items()}, time.monotonic() - step_started, ) except Exception as exc: logger.warning("B04 수치지형도 도엽 확보 실패: %s", exc) except Exception as exc: logger.warning("B04 지도·GIS 다운로드 단계 실패: %s", exc) def _collect_analysis_result( project_root: Path, models_dir: Path, structured_path: Path, bounds_dict: dict[str, float], stats: dict[str, Any], total_points: int, ground_summary: dict[str, Any], manifest: dict[str, Any], sheet_models: list[dict[str, Any]], total_started: float, ) -> dict[str, Any]: """manifest에서 모델 목록을 추려 분석 결과 dict를 조립한다.""" processed = { "processed_file_path": _relative_to_project(project_root, structured_path), "converted_file_path": None, "point_count": total_points, "bounds": bounds_dict, "statistics": stats, } # manifest에서 저장된 모델별 정보 추출 (필터별 대표 모델) models: list[dict[str, Any]] = [] for filter_key, filter_entry in manifest.get("source_filters", {}).items(): for method, meta in filter_entry.get("methods", {}).items(): if meta.get("status") != "completed": continue model_file = meta.get("model_file") model_path = (models_dir / model_file) if model_file else None layers: list[dict[str, Any]] = [] if meta.get("preview_file"): layers.append( { "layer_name": f"{method}_{filter_key}_preview", "geometry_type": "MESH" if method != "meshfree" else "POINTCLOUD", "file_path": _relative_to_project( project_root, models_dir / meta["preview_file"] ), "file_format": "glb" if method != "meshfree" else "ply", } ) models.append( { "model_type": method, "source_filter": filter_key, "representation": meta.get("representation"), "model_file_path": _relative_to_project(project_root, model_path) if model_path else None, "resolution_m": meta.get("grid_resolution_meters"), "generation_params": { "source_filter": filter_key, "representation": meta.get("representation"), "footprint_area_m2": meta.get("footprint_area_m2"), }, "layers": layers, } ) models.extend(sheet_models) logger.info( "B04 WF1 분석 완료: 모델 %d개, 총 %.1fs", len(models), time.monotonic() - total_started ) return { "processed": processed, "ground_summary": ground_summary, "manifest": manifest, "models": models, }