eomsangdonandClaude Opus 5
0a1ebcfbb5
① 3D 서피스 겹쳐보기 어긋남 — `write_glb()`는 넘겨받은 상자의 중심을 원점으로 삼는데 모델마다 자기 상자를 넘겨 도엽 서피스와 라이다 지표면이 다른 원점에 섰다 (2f940d8a 실측: 수평 7.85m·높이 1.28m). 도엽 서피스가 라이다 포인트 상자 (`structured.npz`, 화면 `setReferenceBounds`와 같은 값)를 화면 원점으로 쓰게 했다. 표고 격자 상자(`bounds`)는 절취 범위 그대로 두고 `scene_bounds`를 npz에 따로 남긴다 — 등고선 API도 이 값을 먼저 써서 등고선이 메시 위에 얹힌다. 라이다 없는 사업지는 기준이 자기뿐이라 종전대로 자기 상자를 쓴다. ② 세부유역 저장본 좌표계 — `04_detailed_basins.geojson`에 변환에 쓴 좌표계를 `crs_input`으로 남기고, B07 유역도가 그 값으로 되돌린다. 기록이 없는 옛 저장본은 노선 CSV의 EPSG 라벨로 쓰였으므로 그 라벨로 되돌린다(경고 로그 + 재확정 안내). 검증: tmp/tests 42 passed (신규 test_scene_origin_and_basin_crs.py 5건). 패치 코드로 도엽 서피스를 다시 만들어 GLB 정점 상자를 대조 — sheet [-426.58, -68.95, -386.41]~[411.42, 71.52, 379.59], csf [-172.12, -44.06, -199.69]~[184.11, 3.22, 165.72] 로 같은 원점. 수정 전 sheet는 [-419.5, -70.24, -383.0]~[418.5, 70.24, 383.0] (자기 중심)였다. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
585 lines
25 KiB
Python
585 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).
|
||
|
||
순서는 **2D 먼저, 메시는 맨 마지막**이다(2026-08-30 사용자 지시):
|
||
① 등고 라인을 격자에 굽고 ② 등고선 사이 거리 비례 보간으로 표고 격자를 만든 뒤
|
||
③ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장하고 ④ 라인을 고정한 채
|
||
완화(라플라스)해 등고 간격을 고르게 한다.
|
||
격자에서 뽑는 1m 등고선이 곧 2D 보간선이며, 메시(glb)는 그 격자의 표현일 뿐이다.
|
||
"""
|
||
|
||
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_SheetMethods import (
|
||
SHEET_METHOD_BUILDERS,
|
||
SHEET_METHOD_LABELS,
|
||
)
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import grid_spec_from_bounds
|
||
from config.config_system import (
|
||
SHEET_SURFACE_GRID_M,
|
||
SHEET_SURFACE_MARGIN_M,
|
||
SHEET_SURFACE_METHODS,
|
||
SURFACE_MAX_PREVIEW_VERTICES,
|
||
SURFACE_SMOOTHING_DTM_PREVIEW_RESOLUTION_M,
|
||
SURFACE_SMOOTHING_DTM_SIGMA_M,
|
||
SURFACE_SMOOTHING_DTM_SPLINE_SMOOTH,
|
||
)
|
||
|
||
logger = logging.getLogger(__name__)
|
||
|
||
# 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값)
|
||
_CONTOUR_FILE = "도엽_등고선.geojson"
|
||
|
||
# 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다.
|
||
SHEET_SOURCE_FILTER = "sheet"
|
||
|
||
|
||
def _load_features_metric(
|
||
processed_dir: Path, crs: str, 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", crs, 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 _scene_bounds(project_root: Path, bounds: np.ndarray) -> np.ndarray:
|
||
"""프리뷰 메시를 놓을 **화면 원점 기준 상자**.
|
||
|
||
`write_glb()`는 넘겨받은 상자의 중심을 원점으로 삼아 정점을 옮긴다. 모델마다 자기
|
||
범위를 주면 도엽 서피스와 라이다 지표면이 서로 다른 원점에 서서, 겹쳐 보기·계획선이
|
||
어긋난다(2026-09-01 실측: 수평 7.85m·높이 1.28m). 그래서 라이다 포인트 상자가 있으면
|
||
그 상자를 같이 쓴다 — 화면이 기준으로 삼는 상자(`setReferenceBounds`)와 같은 값이다.
|
||
라이다가 없는 사업지(도엽만)는 기준이 이 서피스뿐이라 자기 상자를 그대로 쓴다.
|
||
"""
|
||
structured = project_root / "B04_PreProcess" / "processed" / "structured.npz"
|
||
if not structured.is_file():
|
||
return bounds
|
||
try:
|
||
with np.load(structured) as data:
|
||
if "bounds" not in data:
|
||
return bounds
|
||
return np.asarray(data["bounds"], dtype=float)
|
||
except (OSError, ValueError) as exc:
|
||
logger.warning("도엽 서피스: 라이다 상자를 읽지 못해 자기 범위로 놓습니다 — %s", exc)
|
||
return bounds
|
||
|
||
|
||
def _preview_mesh(
|
||
x: np.ndarray, y: np.ndarray, z: np.ndarray, valid: np.ndarray
|
||
) -> tuple[np.ndarray, np.ndarray]:
|
||
"""프리뷰용 정점·면 — 정점 수가 상한을 넘으면 격자를 성기게 딴다.
|
||
|
||
격자를 그대로 잇는다. 한때 NURBS 곡면을 걸었으나(2026-08-30) DTM 스무딩이
|
||
들어오면서 곡면 적합이 이중으로 걸려 되돌렸다 — 스무딩은 표고 정본(npz)에서
|
||
한 번만 한다.
|
||
"""
|
||
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 _write_smoothed(
|
||
models_dir: Path,
|
||
stem: str,
|
||
x: np.ndarray,
|
||
y: np.ndarray,
|
||
z: np.ndarray,
|
||
valid: np.ndarray,
|
||
bounds: np.ndarray,
|
||
scene: np.ndarray,
|
||
) -> None:
|
||
"""`{stem}_smooth.npz`·`_smooth_preview.glb`를 만든다 — LAS DTM 스무딩과 같은 절차.
|
||
|
||
`B04_PreProcess_Engine_Smooth.smooth_dtm()`은 `TerrainContext`(라이다 발자국)를
|
||
받으므로 그대로 못 쓴다. 그래서 계수는 **config 값을 그대로** 두고 같은 두 단계만
|
||
옮긴다(2026-08-30 사용자 지시 — 계수 변경 금지):
|
||
|
||
① 무효 영역이 번지지 않는 정규화 가우시안 (`smoothing_dtm_sigma_meters`)
|
||
② C² 바이큐빅 B-spline 재평가 (`kx=ky=3`, `s=smoothing_dtm_spline_smooth`)를
|
||
`smoothing_dtm_preview_resolution_meters` 격자에서
|
||
|
||
화면 스무딩 토글과 확정 스냅샷이 이 파일을 찾으므로 이름 규칙을 지켜야 한다.
|
||
"""
|
||
from scipy.interpolate import RectBivariateSpline
|
||
|
||
from B04_PreProcess.B04_PreProcess_Engine_Smooth import _masked_gaussian_filter
|
||
|
||
cell_m = float(x[1] - x[0]) if len(x) > 1 else SHEET_SURFACE_GRID_M
|
||
sigma_pixels = SURFACE_SMOOTHING_DTM_SIGMA_M / cell_m if cell_m > 0 else 0.0
|
||
# 결측이 하나라도 있으면 스플라인 결과가 통째로 NaN이 된다(TIN·TIN 곡면은 볼록껍질
|
||
# 밖이 결측이다). 최근접 표고로 메워 적합하고 아래에서 원래 마스크로 되돌린다.
|
||
filled = z.astype(np.float64)
|
||
if not valid.all():
|
||
from scipy.ndimage import distance_transform_edt
|
||
|
||
_, (near_row, near_col) = distance_transform_edt(~valid, return_indices=True)
|
||
filled = filled[near_row, near_col]
|
||
z_pre = _masked_gaussian_filter(filled, valid, sigma_pixels)
|
||
try:
|
||
spline = RectBivariateSpline(y, x, z_pre, kx=3, ky=3, s=SURFACE_SMOOTHING_DTM_SPLINE_SMOOTH)
|
||
except Exception as exc: # noqa: BLE001 — 스무딩 실패가 원본 산출을 막으면 안 된다
|
||
logger.warning("도엽 서피스(%s): 스무딩 스플라인 실패(%s) — 건너뜁니다.", stem, exc)
|
||
return
|
||
|
||
# 재평가 격자는 config의 프리뷰 해상도를 쓰되 원본보다 성기게 잡지 않는다 —
|
||
# 이 npz는 화면용이자 **스무딩 확정 시 종·횡단이 샘플링하는 표고 정본**이라,
|
||
# 정점 상한(화면 사정)으로 해상도를 깎으면 설계 정밀도가 같이 깎인다.
|
||
# 메시 정점 수는 _preview_mesh가 알아서 성기게 딴다.
|
||
step = max(SURFACE_SMOOTHING_DTM_PREVIEW_RESOLUTION_M, cell_m)
|
||
sx = np.arange(x[0], x[-1] + step * 0.5, step)
|
||
sy = np.arange(y[0], y[-1] + step * 0.5, step)
|
||
sz = np.asarray(spline(sy, sx), dtype=np.float32)
|
||
# 원본 유효 마스크를 최근접으로 옮겨 무효 영역을 그대로 지킨다.
|
||
col = np.clip(np.searchsorted(x, sx) - 1, 0, len(x) - 1)
|
||
row = np.clip(np.searchsorted(y, sy) - 1, 0, len(y) - 1)
|
||
svalid = valid[np.ix_(row, col)]
|
||
sz[~svalid] = np.nan
|
||
|
||
atomic_npz(
|
||
models_dir / f"{stem}_smooth.npz",
|
||
x=sx,
|
||
y=sy,
|
||
z=sz,
|
||
valid_mask=svalid,
|
||
bounds=bounds,
|
||
scene_bounds=scene,
|
||
resolution=np.array([step], np.float32),
|
||
)
|
||
vertices, faces = _preview_mesh(sx, sy, np.nan_to_num(sz, nan=float(bounds[2, 0])), svalid)
|
||
write_glb(models_dir / f"{stem}_smooth_preview.glb", vertices, faces, scene)
|
||
|
||
|
||
def _rasterize_contour_levels(spec: Any, features: list[dict[str, Any]]) -> np.ndarray:
|
||
"""등고 라인을 격자에 굽는다 — 셀 = 그 위를 지나는 라인의 표고, 그 외 NaN.
|
||
|
||
배수유역의 `rasterize_contours()`와 두 가지가 다르다(둘 다 서피스 품질 때문이다):
|
||
· `all_touched=False` — 스치는 셀까지 칠하면 라인이 2px 두께가 되고, 그 폭만큼
|
||
정확히 등고 표고인 **평탄 띠**가 생겨 사이 1m 등고선 간격이 찌그러진다
|
||
(2026-08-30 사용자 지적: 보간선이 등간격이 아님).
|
||
· 길이 필터 없음 — 봉우리 폐합 링 같은 짧은 등고선을 버리면 그 일대가 통째로
|
||
평평해진다. 배수유역은 노이즈를 버려야 하지만 지형면은 다 있어야 한다.
|
||
|
||
같은 셀을 두 표고가 지나면 낮은 쪽을 남긴다(배수유역과 같은 규칙).
|
||
"""
|
||
from rasterio.features import rasterize
|
||
from shapely.geometry import shape
|
||
|
||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
|
||
ELEVATION_KEYS,
|
||
grid_transform,
|
||
iter_linestrings,
|
||
)
|
||
|
||
by_level: dict[float, list[Any]] = {}
|
||
for feature in features:
|
||
geometry = feature.get("geometry")
|
||
if not geometry:
|
||
continue
|
||
properties = feature.get("properties") or {}
|
||
elevation = next(
|
||
(float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None),
|
||
None,
|
||
)
|
||
if elevation is None:
|
||
continue
|
||
try:
|
||
parsed = shape(geometry)
|
||
except Exception: # noqa: BLE001
|
||
continue
|
||
for line in iter_linestrings(parsed):
|
||
by_level.setdefault(elevation, []).append(line)
|
||
|
||
burned = np.full((spec.n_rows, spec.n_cols), np.nan, dtype=np.float32)
|
||
transform = grid_transform(spec)
|
||
for elevation in sorted(by_level, reverse=True):
|
||
stamp = rasterize(
|
||
[(line, 1) for line in by_level[elevation]],
|
||
out_shape=(spec.n_rows, spec.n_cols),
|
||
transform=transform,
|
||
fill=0,
|
||
dtype="uint8",
|
||
all_touched=False,
|
||
).astype(bool)
|
||
burned[stamp] = elevation
|
||
logger.info(
|
||
"도엽 서피스: 등고 라인 %d단을 격자에 굽어 %d셀",
|
||
len(by_level),
|
||
int(np.isfinite(burned).sum()),
|
||
)
|
||
return burned
|
||
|
||
|
||
def _resolve_enclosed_interiors(
|
||
burned: np.ndarray, present: list[float], surface: np.ndarray, cell_m: float, interval_m: float
|
||
) -> np.ndarray:
|
||
"""폐합 등고선 안쪽(봉우리·웅덩이)을 바깥 사면 경사로 연장한다 (in-place).
|
||
|
||
거리 보간은 "가장 가까운 서로 다른 두 라인 사이"를 채우므로, 마지막 등고선
|
||
안쪽에는 더 높은 라인이 없어 아래쪽 라인 쪽으로 끌려 **분화구처럼 파인다**.
|
||
등고선이 'ㅜ'에서 점점 짧아지다 사라지는 마루가 바로 이 자리다(2026-08-30
|
||
사용자 지적).
|
||
|
||
폐합 라인 내부에 다른 표고 제약이 하나도 없으면 그 안이 마루(바깥이 낮을 때)
|
||
또는 웅덩이(바깥이 높을 때)다. 바깥 사면의 국소 경사를 안쪽으로 연장하되
|
||
±(간격−0.5m)로 제한한다 — 다음 등고선이 없다는 사실과 모순되지 않는 범위다.
|
||
처리한 영역의 마스크를 반환한다 — 이어지는 완화에서 함께 고정해야 한다.
|
||
"""
|
||
from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label
|
||
|
||
constrained = np.isfinite(burned)
|
||
band_cells = 8
|
||
limit = max(interval_m - 0.5, 0.5)
|
||
resolved = 0
|
||
handled = np.zeros(burned.shape, dtype=bool)
|
||
for level in present:
|
||
mask = burned == level
|
||
interior = binary_fill_holes(mask) & ~mask
|
||
if not interior.any():
|
||
continue
|
||
components, count = label(interior)
|
||
for component_id in range(1, count + 1):
|
||
component = components == component_id
|
||
if (constrained & component).any():
|
||
continue # 안에 다른 제약이 있으면 마루가 아니다(보통의 감싸는 링)
|
||
ring = binary_dilation(binary_fill_holes(mask) | mask) & ~component & ~mask
|
||
ring &= np.isfinite(surface)
|
||
if not ring.any():
|
||
continue
|
||
outside_mean = float(surface[ring].mean())
|
||
direction = 1.0 if outside_mean < level else -1.0
|
||
inner = distance_transform_edt(component, sampling=cell_m)
|
||
# 바깥 사면 경사 — 라인 밖 band_cells 이내 유효 셀의 (낙차 / 거리) 평균.
|
||
outer_distance = distance_transform_edt(~(mask | component), sampling=cell_m)
|
||
band = (outer_distance > 0) & (outer_distance <= band_cells * cell_m)
|
||
band &= np.isfinite(surface) & ~component & ~mask
|
||
if band.any():
|
||
slope = float(
|
||
np.mean(np.abs(level - surface[band].astype(np.float64)) / outer_distance[band])
|
||
)
|
||
else:
|
||
slope = 0.0
|
||
if slope > 1e-3:
|
||
offset = np.minimum(slope * inner[component], limit)
|
||
else:
|
||
peak = float(inner.max())
|
||
offset = (limit / 2.0) * (inner[component] / peak) if peak > 0 else 0.0
|
||
surface[component] = level + direction * offset
|
||
handled |= component
|
||
resolved += 1
|
||
if resolved:
|
||
logger.info("도엽 서피스: 폐합 등고선 내부 %d곳을 사면 경사로 연장", resolved)
|
||
return handled
|
||
|
||
|
||
def _write_method_model(
|
||
project_root: Path,
|
||
models_dir: Path,
|
||
spec: Any,
|
||
surface: np.ndarray,
|
||
method_key: str,
|
||
) -> dict[str, Any] | None:
|
||
"""방식 하나의 표고 격자를 npz·프리뷰 glb로 저장하고 등록용 dict를 만든다."""
|
||
# 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("도엽 서피스(%s): 유효 표고 셀이 없습니다.", method_key)
|
||
return None
|
||
|
||
source_filter = f"{SHEET_SOURCE_FILTER}_{method_key}"
|
||
stem = f"dtm_{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())],
|
||
]
|
||
)
|
||
# 표고 격자 상자(`bounds`)는 절취 범위 그대로 두고, 화면 원점만 라이다와 맞춘다.
|
||
scene = _scene_bounds(project_root, bounds)
|
||
atomic_npz(
|
||
model_path,
|
||
x=x_coords,
|
||
y=y_coords,
|
||
z=z_grid,
|
||
valid_mask=valid_grid,
|
||
bounds=bounds,
|
||
scene_bounds=scene,
|
||
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, scene)
|
||
_write_smoothed(models_dir, stem, x_coords, y_coords, z_grid, valid_grid, bounds, scene)
|
||
return {
|
||
"model_type": "dtm",
|
||
"source_filter": 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": source_filter,
|
||
"representation": "regular_grid",
|
||
"source": "map_sheet_contours",
|
||
"interpolation": method_key,
|
||
"interpolation_label": SHEET_METHOD_LABELS.get(method_key, method_key),
|
||
"margin_m": SHEET_SURFACE_MARGIN_M,
|
||
},
|
||
"layers": [
|
||
{
|
||
"layer_name": f"{stem}_preview",
|
||
"geometry_type": "MESH",
|
||
"file_path": str(preview_path.relative_to(project_root)).replace("\\", "/"),
|
||
"file_format": "glb",
|
||
}
|
||
],
|
||
}
|
||
|
||
|
||
def build_sheet_surface_model(
|
||
project_root: Path,
|
||
processed_dir: Path,
|
||
models_dir: Path,
|
||
route_xy: np.ndarray,
|
||
crs: str,
|
||
methods: list[str] | None = None,
|
||
) -> list[dict[str, Any]]:
|
||
"""도엽등고선으로 방식별 DTM npz·프리뷰 glb를 만들고 등록용 dict 목록을 돌려준다.
|
||
|
||
방식을 하나로 고르지 않고 전부 만들어 두는 이유: 문헌상 지형에 따라 우열이 갈려
|
||
화면에서 바꿔 보며 정해야 한다(2026-08-30 사용자 지시). 실패하면 빈 목록 —
|
||
호출측은 분석을 계속한다(도엽 미확보 지역 폴백).
|
||
|
||
`route_xy`: (N, 2) 노선 정점 XY(사업지 CRS, m), `crs`: 그 좌표계(pyproj 입력 문자열).
|
||
"""
|
||
started = time.monotonic()
|
||
features = _load_features_metric(processed_dir, crs)
|
||
if not features:
|
||
return []
|
||
|
||
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
|
||
spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M)
|
||
|
||
# ① 2D — 등고 라인을 격자에 굽는다(라인 셀 = 표고, 그 외 NaN).
|
||
burned = _rasterize_contour_levels(spec, features)
|
||
present = sorted(np.unique(burned[np.isfinite(burned)]).tolist())
|
||
if len(present) < 2:
|
||
logger.warning("도엽 서피스: 절취 범위 안에 등고선이 부족합니다.")
|
||
return []
|
||
# 등고 간격(m) — 마루 연장 상한의 근거. 레벨이 하나뿐이면 5m(주곡선) 폴백.
|
||
interval_m = float(np.diff(np.array(present)).min()) if len(present) > 1 else 5.0
|
||
|
||
selected = methods or list(SHEET_SURFACE_METHODS)
|
||
models: list[dict[str, Any]] = []
|
||
for method_key in selected:
|
||
builder = SHEET_METHOD_BUILDERS.get(method_key)
|
||
if builder is None:
|
||
logger.warning("도엽 서피스: 알 수 없는 보간 방식 %s — 건너뜁니다.", method_key)
|
||
continue
|
||
step_started = time.monotonic()
|
||
try:
|
||
# ② 2D 보간 — 여기서 나온 격자에서 1m 등고선을 뽑으므로 화면 등고선이 곧
|
||
# 2D 보간선이다. 메시(glb)는 그 격자의 표현일 뿐이다(사용자 지시).
|
||
surface = builder(spec, burned, features, spec.cell_m)
|
||
# ③ 폐합 등고선 안쪽(마루·웅덩이)은 방식과 무관하게 같은 규칙으로 채운다.
|
||
_resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m)
|
||
except Exception as exc: # noqa: BLE001 — 한 방식이 죽어도 나머지는 만든다
|
||
logger.warning("도엽 서피스(%s) 생성 실패: %s", method_key, exc)
|
||
continue
|
||
model = _write_method_model(project_root, models_dir, spec, surface, method_key)
|
||
if model is not None:
|
||
models.append(model)
|
||
logger.info("도엽 서피스(%s) 완료 (%.1fs)", method_key, time.monotonic() - step_started)
|
||
|
||
logger.info(
|
||
"도엽 서피스 생성 완료: %d×%d 격자, 등고 %d단, 방식 %d개 (%.1fs)",
|
||
spec.n_rows,
|
||
spec.n_cols,
|
||
len(present),
|
||
len(models),
|
||
time.monotonic() - started,
|
||
)
|
||
return models
|
||
|
||
|
||
def build_sheet_surface_from_route(
|
||
project_root: Path,
|
||
processed_dir: Path,
|
||
models_dir: Path,
|
||
methods: list[str] | None = None,
|
||
) -> list[dict[str, Any]]:
|
||
"""B03 업로드 계획노선을 찾아 지정한 방식의 도엽 서피스를 만든다. 없으면 빈 목록.
|
||
|
||
`methods`를 주지 않으면 `SHEET_SURFACE_METHODS` 전체를 만든다 — 관리자가 화면에서
|
||
한 방식을 요청할 때 그 목록만 넘긴다.
|
||
"""
|
||
from B04_PreProcess.B04_PreProcess_Engine_Extent import project_epsg_from_prj
|
||
from common_util.common_util_route_geometry import load_design_route
|
||
|
||
# 노선 CSV의 `crs_epsg` 열은 표시용 라벨이라 LAS(.prj) 좌표계와 다를 수 있다
|
||
# (2026-09-01 실측: 라벨 5179, 실제 5176 — 도엽 서피스만 딴 자리에 만들어졌다).
|
||
# `load_design_route()`가 .prj 좌표계로 재투영해 주므로 라이다 지표면과 한 자리에 선다.
|
||
planned = load_design_route(project_root)
|
||
if planned is None or len(planned.vertices) < 2:
|
||
logger.warning("도엽 서피스: 계획 노선 파일을 읽지 못했습니다.")
|
||
return []
|
||
crs = planned.crs_input or project_epsg_from_prj(project_root)
|
||
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, crs, methods
|
||
)
|
||
|
||
|
||
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
|
||
|
||
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(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 서피스 생성 중")
|
||
models = build_sheet_surface_model(
|
||
project_root, processed_dir, models_dir, route_xy, f"EPSG:{epsg}"
|
||
)
|
||
if not models:
|
||
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": models,
|
||
}
|