feat(B03/B04): LAS 없이 도엽등고선으로 3D 서피스를 만들어 설계한다
- B03: 'LAS 없이 설계' 토글 — LAS 필수카드 비활성화, las_free 플래그로
업로드·완료 검증 면제, 계획노선 CSV를 WF1 분석 입력으로 사용
- B04: 신규 Engine_SheetSurface — 도엽_등고선.geojson을 노선 bbox+300m
직사각형으로 절취, 배수유역 엔진의 등고선 정점구름·Delaunay TIN 보간을
재사용해 dtm_sheet.npz(1m 격자, LAS DTM과 동일 형식) + 프리뷰 glb 생성.
build_surface_sampler('sheet','dtm')로 종·횡단·배수 하류 계산 무수정 동작
- WF1: las_free면 run_sheet_surface_analysis로 분기, sheet/dtm 자동 확정.
VWorld·도엽 확보 블록을 download_geodata()로 추출해 두 경로가 공유
- LAS가 있어도 도엽 서피스를 함께 생성·등록(참고용), B04에
'도엽등고 3D 서피스' 별도 컨테이너로 표시
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -0,0 +1,299 @@
|
||||
"""도엽등고선 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"
|
||||
atomic_npz(
|
||||
model_path,
|
||||
x=x_coords,
|
||||
y=y_coords,
|
||||
z=z_grid,
|
||||
valid_mask=valid_grid,
|
||||
resolution=np.array([SHEET_SURFACE_GRID_M], np.float32),
|
||||
)
|
||||
|
||||
vertices, faces = _preview_mesh(x_coords, y_coords, z_grid, valid_grid)
|
||||
finite_z = z_grid[valid_grid]
|
||||
bounds = np.array(
|
||||
[
|
||||
[x_coords[0], x_coords[-1]],
|
||||
[y_coords[0], y_coords[-1]],
|
||||
[float(finite_z.min()), float(finite_z.max())],
|
||||
]
|
||||
)
|
||||
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],
|
||||
}
|
||||
Reference in New Issue
Block a user