eomsangdonandClaude Opus 5
cb1c5261e2
앞서는 등거리 '점'을 5셀 간격으로 흩뿌려 점 사이가 벌어진 곳에 평탄 삼각형이 남았고 한 단계(+2.5m)만 넣어 그 안에서 계단이 남았다. 사용자 지시(2D에서 사이 등고선을 먼저 만들고 3D화)대로 등거리선을 연속 셀 선으로 뽑아 새 레벨로 등록하고, 그 선들 사이를 다시 이분한다 (SHEET_SURFACE_MIDLINE_ROUNDS=2 → 5m 주곡선이 2.5m·1.25m 가상 등고로). 정점 간격도 5셀에서 2m로 좁혔다. 거리장은 라운드 간 재사용한다. 실측(c1bb453f): 평탄 셀 3.96%→2.48%(최초 9.35%), 경사 2% 미만 셀 27163→17145, 노선 |Δz| 2.56m 유지. 생성 16.4s→36.5s — 라운드 수는 config로 조절 가능. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
582 lines
25 KiB
Python
582 lines
25 KiB
Python
"""도엽등고선 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,
|
||
SHEET_SURFACE_MIDLINE_MAX_M,
|
||
SHEET_SURFACE_MIDLINE_ROUNDS,
|
||
SURFACE_MAX_PREVIEW_VERTICES,
|
||
)
|
||
|
||
logger = logging.getLogger(__name__)
|
||
|
||
# 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값)
|
||
_CONTOUR_FILE = "도엽_등고선.geojson"
|
||
_STREAM_FILE = "도엽_하천중심선.geojson"
|
||
|
||
# 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다.
|
||
SHEET_SOURCE_FILTER = "sheet"
|
||
|
||
|
||
def _load_features_metric(
|
||
processed_dir: Path, epsg: int, filename: str = _CONTOUR_FILE
|
||
) -> list[dict[str, Any]]:
|
||
"""병합 도엽 레이어(WGS84)를 읽어 사업지 CRS(m)로 재투영한다."""
|
||
path = processed_dir / filename
|
||
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 _stream_breakline_vertices(
|
||
spec: Any, burned: np.ndarray, stream_features: list[dict[str, Any]]
|
||
) -> tuple[np.ndarray, np.ndarray] | None:
|
||
"""계곡 기준선(하천중심선)을 따라 등고 교차점 사이를 보간한 정점열을 만든다.
|
||
|
||
구조선 기반 보간(2026-08-30 사용자 확정): 계곡 축을 먼저 세우고, 축이 등고선과
|
||
만나는 점을 Z 앵커로 삼아 앵커 사이를 선 길이 비례로 보간한다. 계곡 바닥이
|
||
등고선 사이에서도 연속으로 내려가는 가상 종단이 되어 골짜기 평탄·역경사가 준다.
|
||
"""
|
||
from shapely import segmentize
|
||
from shapely.geometry import shape
|
||
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import iter_linestrings
|
||
|
||
sample_step_m = 2.0
|
||
vertex_stride = 3 # 2m 샘플 → 6m 간격 정점 (등고 정점 재샘플 밀도와 유사)
|
||
collected_xy: list[np.ndarray] = []
|
||
collected_z: list[np.ndarray] = []
|
||
for feature in stream_features:
|
||
geometry = feature.get("geometry")
|
||
if not geometry:
|
||
continue
|
||
try:
|
||
parsed = shape(geometry)
|
||
except Exception: # noqa: BLE001
|
||
continue
|
||
for line in iter_linestrings(parsed):
|
||
coords = np.asarray(segmentize(line, sample_step_m).coords, dtype=np.float64)
|
||
if len(coords) < 3:
|
||
continue
|
||
xy = coords[:, :2]
|
||
arc = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))])
|
||
row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy())
|
||
inside = (row >= 0) & (col >= 0)
|
||
levels_at = np.full(len(xy), np.nan, dtype=np.float32)
|
||
levels_at[inside] = burned[row[inside], col[inside]]
|
||
hit = np.isfinite(levels_at)
|
||
if hit.sum() < 2:
|
||
continue
|
||
# 등고선 통과 구간(연속 같은 표고 샘플 묶음) 하나 = 앵커 하나.
|
||
anchor_arc: list[float] = []
|
||
anchor_z: list[float] = []
|
||
run_start: int | None = None
|
||
for i in range(len(xy) + 1):
|
||
is_hit = i < len(xy) and hit[i]
|
||
if is_hit and run_start is not None and levels_at[i] != levels_at[run_start]:
|
||
is_hit = False # 표고가 바뀌면 묶음을 끊는다
|
||
if is_hit:
|
||
if run_start is None:
|
||
run_start = i
|
||
elif run_start is not None:
|
||
anchor_arc.append(float(arc[run_start : i if i <= len(xy) else len(xy)].mean()))
|
||
anchor_z.append(float(levels_at[run_start]))
|
||
run_start = i if i < len(xy) and hit[i] else None
|
||
if len(anchor_arc) < 2:
|
||
continue
|
||
# 앵커 사이 구간만 채택 — 첫 앵커 이전·마지막 앵커 이후는 근거가 없다.
|
||
span = (arc >= anchor_arc[0]) & (arc <= anchor_arc[-1]) & inside & ~hit
|
||
take = np.flatnonzero(span)[::vertex_stride]
|
||
if not len(take):
|
||
continue
|
||
collected_xy.append(xy[take])
|
||
collected_z.append(np.interp(arc[take], anchor_arc, anchor_z))
|
||
if not collected_xy:
|
||
return None
|
||
return np.vstack(collected_xy), np.concatenate(collected_z)
|
||
|
||
|
||
def _densify_between_contours(
|
||
spec: Any, burned: np.ndarray, present: list[float], cloud: Any
|
||
) -> Any:
|
||
"""등고선 사이에 2D 가상 등고라인을 반복 이분으로 만들어 정점으로 추가한다.
|
||
|
||
등고선 정점만의 Delaunay TIN은 같은 표고 정점 3개짜리 평탄 삼각형이 굴곡부에
|
||
계단(terrace)을 만든다. 인접 표고 라인 쌍의 등거리선(사이 등고선)을 **연속 셀
|
||
선**으로 뽑아 새 레벨로 등록하고, 그 선들 사이를 다시 이분하는 식으로
|
||
`SHEET_SURFACE_MIDLINE_ROUNDS`회 반복한다 — 5m 주곡선이면 2.5m, 1.25m 가상
|
||
등고가 생겨 TIN이 어디서든 서로 다른 표고를 잇는다(2026-08-30 사용자 지시:
|
||
2D에서 사이 등고선을 먼저 만들고 3D화).
|
||
|
||
등거리 조건은 멀리 있는 다른 표고 쌍에도 우연히 성립하므로(예: 525~530
|
||
골짜기에 555/560 중간점), "그 지점의 최근접 라인이 바로 그 쌍"일 때만 인정한다.
|
||
"""
|
||
from scipy.ndimage import distance_transform_edt
|
||
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ContourCloud
|
||
|
||
if len(present) < 2:
|
||
return cloud
|
||
|
||
# 라운드마다 거리장을 다시 구하되, 직전 라운드 것을 재사용한다 — 새 라인 사이의
|
||
# 이분은 (직전 레벨 거리장, 새 중간선 거리장)만 있으면 되므로 전량 재계산이 필요없다.
|
||
levels = list(present)
|
||
masks = [burned == level for level in levels]
|
||
distances = [
|
||
distance_transform_edt(~mask, sampling=spec.cell_m).astype(np.float32) for mask in masks
|
||
]
|
||
virtual: dict[float, np.ndarray] = {}
|
||
for _ in range(max(SHEET_SURFACE_MIDLINE_ROUNDS, 0)):
|
||
if len(levels) < 2:
|
||
break
|
||
best_distance = np.full(burned.shape, np.inf, dtype=np.float32)
|
||
best_index = np.full(burned.shape, -1, dtype=np.int16)
|
||
for index, distance in enumerate(distances):
|
||
take = distance < best_distance
|
||
best_distance[take] = distance[take]
|
||
best_index[take] = index
|
||
next_levels: list[float] = []
|
||
next_masks: list[np.ndarray] = []
|
||
next_distances: list[np.ndarray] = []
|
||
added_any = False
|
||
for index in range(len(levels)):
|
||
next_levels.append(levels[index])
|
||
next_masks.append(masks[index])
|
||
next_distances.append(distances[index])
|
||
if index == 0:
|
||
continue
|
||
lower, upper = distances[index - 1], distances[index]
|
||
near = (lower < SHEET_SURFACE_MIDLINE_MAX_M) & (upper < SHEET_SURFACE_MIDLINE_MAX_M)
|
||
pair_is_nearest = (best_index == index - 1) | (best_index == index)
|
||
midline = near & pair_is_nearest & (np.abs(lower - upper) <= spec.cell_m)
|
||
midline &= ~masks[index - 1] & ~masks[index]
|
||
if not midline.any():
|
||
continue
|
||
mid_level = (levels[index - 1] + levels[index]) / 2.0
|
||
virtual[mid_level] = midline
|
||
# 새 라인은 마지막(=자기 레벨) 자리 앞에 끼워 정렬을 유지한다.
|
||
next_levels.insert(-1, mid_level)
|
||
next_masks.insert(-1, midline)
|
||
next_distances.insert(
|
||
-1, distance_transform_edt(~midline, sampling=spec.cell_m).astype(np.float32)
|
||
)
|
||
added_any = True
|
||
levels, masks, distances = next_levels, next_masks, next_distances
|
||
if not added_any:
|
||
break
|
||
|
||
if not virtual:
|
||
return cloud
|
||
xs = spec.cell_centers_x()
|
||
ys = spec.cell_centers_y()
|
||
extra_xy: list[np.ndarray] = []
|
||
extra_z: list[np.ndarray] = []
|
||
for level, mask in virtual.items():
|
||
rows, cols = np.nonzero(mask)
|
||
rows, cols = rows[::2], cols[::2] # 1m 셀 → 2m 간격 정점 (TIN 비용 절제)
|
||
extra_xy.append(np.column_stack([xs[cols], ys[rows]]))
|
||
extra_z.append(np.full(len(rows), level))
|
||
added = int(sum(len(a) for a in extra_z))
|
||
logger.info(
|
||
"도엽 서피스: 가상 등고라인 %d단, 정점 %d개 추가 (원 등고 %d단)",
|
||
len(virtual),
|
||
added,
|
||
len(present),
|
||
)
|
||
return ContourCloud(
|
||
xy=np.vstack([cloud.xy, *extra_xy]),
|
||
z=np.concatenate([cloud.z, *extra_z]),
|
||
)
|
||
|
||
|
||
def _cap_flat_summits(surface: np.ndarray, cell_m: float, interval_m: float) -> int:
|
||
"""마지막 등고선 안쪽 평탄 마루를 완만한 능선 돔으로 올린다 (in-place).
|
||
|
||
다음 등고선이 없는 봉우리·능선 마루는 TIN이 마지막 등고 표고로 평탄하게 채운다
|
||
('ㅜ'가 갑자기 'ㅡ'로 변하는 단 — 2026-08-30 사용자 지적). 실제 마루는 그 표고와
|
||
+간격 사이이므로, 마지막 등고선 바깥 사면의 국소 경사를 안쪽으로 연장해 올린다
|
||
(상한 +간격−0.5m — 다음 등고선이 없다는 사실과 모순되지 않게). 바깥 경사를
|
||
구할 수 없으면 +간격/2 돔으로 폴백한다(2026-08-30 사용자 확정 ②a).
|
||
적용 조건: 등표고 평탄 컴포넌트이고 바깥 유효 이웃이 전부 더 낮을 때(=마루).
|
||
골짜기 바닥은 배수 방향을 모르므로 건드리지 않는다. 돋운 컴포넌트 수를 반환.
|
||
"""
|
||
from scipy.ndimage import binary_dilation, distance_transform_edt, label
|
||
|
||
finite = np.isfinite(surface)
|
||
quantized = np.round(surface * 100.0)
|
||
capped = 0
|
||
# 평탄 후보: 상하좌우 이웃 중 같은 표고가 있는 셀만 모아 컴포넌트를 짠다.
|
||
same_right = np.zeros_like(finite)
|
||
same_right[:, :-1] = finite[:, :-1] & finite[:, 1:] & (quantized[:, :-1] == quantized[:, 1:])
|
||
same_down = np.zeros_like(finite)
|
||
same_down[:-1, :] = finite[:-1, :] & finite[1:, :] & (quantized[:-1, :] == quantized[1:, :])
|
||
flat = same_right | same_down
|
||
flat[:, 1:] |= same_right[:, :-1]
|
||
flat[1:, :] |= same_down[:-1, :]
|
||
|
||
# 서로 다른 표고의 평탄면이 맞닿아 있을 수 있으므로 표고값별로 컴포넌트를 짠다.
|
||
# 값별 셀 수 25 미만은 컴포넌트도 25 미만이므로 라벨링 전에 거른다(속도).
|
||
flat_values, flat_counts = np.unique(quantized[flat], return_counts=True)
|
||
for value in flat_values[flat_counts >= 25]:
|
||
value_labels, value_count = label(flat & (quantized == value))
|
||
sizes = np.bincount(value_labels.reshape(-1))
|
||
for component_id in range(1, value_count + 1):
|
||
if sizes[component_id] < 25: # 25㎡ 미만은 시각적으로 단이 아니다
|
||
continue
|
||
rows, cols = np.nonzero(value_labels == component_id)
|
||
if (
|
||
rows.min() == 0
|
||
or cols.min() == 0
|
||
or rows.max() == surface.shape[0] - 1
|
||
or cols.max() == surface.shape[1] - 1
|
||
):
|
||
continue # 격자 가장자리에 닿으면 바깥을 모른다
|
||
# 이후 연산은 컴포넌트 bbox 창(+경사 밴드 여유)에서만 — 전체 격자 반복 회피.
|
||
band_cells = 8
|
||
window = (
|
||
slice(
|
||
max(rows.min() - band_cells, 0),
|
||
min(rows.max() + band_cells + 1, surface.shape[0]),
|
||
),
|
||
slice(
|
||
max(cols.min() - band_cells, 0),
|
||
min(cols.max() + band_cells + 1, surface.shape[1]),
|
||
),
|
||
)
|
||
component = value_labels[window] == component_id
|
||
ring = binary_dilation(component) & ~component & finite[window]
|
||
if not ring.any():
|
||
continue
|
||
level = float(surface[rows[0], cols[0]])
|
||
patch = surface[window]
|
||
if float(patch[ring].max()) >= level - 1e-3:
|
||
continue # 더 높은 이웃이 있으면 마루가 아니다(사면 벤치·골짜기)
|
||
inner = distance_transform_edt(component, sampling=cell_m)
|
||
peak = float(inner.max())
|
||
if peak <= 0.0:
|
||
continue
|
||
# 바깥 사면 경사 추정 — 마지막 등고선 밖 밴드의 (표고 낙차 / 거리) 평균을
|
||
# 안쪽으로 연장한다(2026-08-30 사용자 확정 ②a). 상한 +간격−0.5m:
|
||
# 다음 등고선이 없다는 사실(마루 < L+간격)과 모순되지 않게.
|
||
outer_distance = distance_transform_edt(~component, sampling=cell_m)
|
||
band = (~component) & finite[window] & (outer_distance <= band_cells * cell_m)
|
||
band &= outer_distance > 0
|
||
if band.any():
|
||
drops = level - patch[band].astype(np.float64)
|
||
slope = float(np.mean(drops / outer_distance[band]))
|
||
else:
|
||
slope = 0.0
|
||
if slope > 1e-3:
|
||
rise = np.minimum(slope * inner[component], interval_m - 0.5)
|
||
else: # 경사를 못 구하면 +간격/2 돔 폴백
|
||
rise = (interval_m / 2.0) * (inner[component] / peak)
|
||
patch[component] += rise.astype(surface.dtype)
|
||
capped += 1
|
||
if capped:
|
||
logger.info(
|
||
"도엽 서피스: 평탄 마루 %d곳을 +%.1fm 골격 캡으로 돋움", capped, interval_m / 2.0
|
||
)
|
||
return capped
|
||
|
||
|
||
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_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)
|
||
# 등고 간격(m) — 마루 캡 크기의 근거. 레벨이 하나뿐이면 5m(1:5,000 주곡선) 폴백.
|
||
contour_levels = np.unique(cloud.z)
|
||
interval_m = float(np.diff(contour_levels).min()) if len(contour_levels) > 1 else 5.0
|
||
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ContourCloud
|
||
|
||
burned, levels = rasterize_contours(spec, features, None)
|
||
present = [level for level in sorted(levels) if bool((burned == level).any())]
|
||
# 계단(평탄 삼각형) 해소 — 인접 등고선 사이 중간 보간선을 정점으로 추가.
|
||
cloud = _densify_between_contours(spec, burned, present, cloud)
|
||
# 계곡 구조선 — 하천중심선을 따라 등고 교차점 사이를 보간한 정점 추가.
|
||
stream_features = _load_features_metric(processed_dir, epsg, _STREAM_FILE)
|
||
if stream_features:
|
||
stream_vertices = _stream_breakline_vertices(spec, burned, stream_features)
|
||
if stream_vertices is not None:
|
||
logger.info("도엽 서피스: 계곡 구조선 정점 %d개 추가", len(stream_vertices[1]))
|
||
cloud = ContourCloud(
|
||
xy=np.vstack([cloud.xy, stream_vertices[0]]),
|
||
z=np.concatenate([cloud.z, stream_vertices[1]]),
|
||
)
|
||
surface = interpolate_elevation(spec, cloud) # (R, C), 북→남 행 순서, 외부 NaN
|
||
# 마지막 등고선 안쪽 평탄 마루를 완만한 돔으로 (2026-08-30 사용자 확정).
|
||
_cap_flat_summits(surface, spec.cell_m, interval_m)
|
||
|
||
# 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],
|
||
}
|