HYUNJUNGLEE
b9342f6726
S-CANVAS (Saman Corp.) — DXF + DEM + AI 기반 3D 조감도 생성 엔진. ~24k LOC Python (scanvas_maker.py 7072 LOC GUI + 구조물 파서/빌더 다수). 이 커밋은 7-iter cleanup이 적용된 상태로 import: - F821 8 + B023 6: 비동기 lambda + except/loop 변수 캡처 NameError (Py3.13에서 reproduce 확인된 진짜 버그) - RUF012 4 + RUF013 1: ClassVar / implicit Optional 명시화 - F811/B905/B904/F401/F841/W293/F541/UP/SIM/RUF/PLR 700+ cleanup/modernization 신규 파일: - ruff.toml: target=py313, Korean unicode/저자 스타일/도메인 복잡도 무력화 - requirements-py313.txt: pyproj>=3.7, scipy>=1.14, numpy>=2.0.2 (Py3.13 wheel) - .gitignore: gcp-key.json, 캐시, 백업, 생성 이미지 제외 검증: ruff 0 errors, py_compile 0 errors, import 33/33 OK on Py3.13.13. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
868 lines
40 KiB
Python
868 lines
40 KiB
Python
"""실제 지형 DEM으로 DXF TIN 외곽을 확장하는 유틸리티.
|
||
|
||
동기:
|
||
사용자가 제공한 DXF 평면도의 등고선 범위를 벗어나면 3D 장면에서 지형이
|
||
절벽처럼 끊긴다. 조감도에서 이 구역을 AI 프롬프트로 메우면 사실 기반이
|
||
아니라 AI 상상이라 결과가 가변적이다. 이 모듈은 DXF 범위 바깥에 대해
|
||
실제 공개 DEM(AWS Open Terrain Tiles, terrarium PNG 포맷)을 받아 외곽 링
|
||
메시를 만들고, 기존 TIN과 Z-연속성을 맞춰 이어 붙인다.
|
||
|
||
기본 소스:
|
||
AWS Open Terrain Tiles (terrarium): API 키 없음, 글로벌, ~30m 수직 정확도.
|
||
URL: https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png
|
||
디코딩: elev_m = (R*256 + G + B/256) - 32768
|
||
|
||
옵션 소스:
|
||
- 로컬 GeoTIFF: `cache/dem/local.tif` 가 존재하고 rasterio가 설치되어
|
||
있으면 우선 사용. 한국이면 NGII 5m DEM을 여기 놓아두면 정확도 상승.
|
||
|
||
좌표계:
|
||
- TIN/메시 프레임: 투영 CRS(예: EPSG:5187) - origin (zero-basing)
|
||
- DEM fetch: WGS84 (EPSG:4326)
|
||
- 본 모듈은 투영 CRS와 origin만 받으면 내부에서 pyproj로 변환 처리.
|
||
|
||
Seam 처리:
|
||
1. 수직 datum 자동 보정: inner 경계 근처 TIN Z와 DEM Z의 평균 차이로
|
||
global offset을 잡아 DEM 전체 Z에서 차감.
|
||
2. Feathering: inner 경계로부터 feather_m 이내 도넛 점은 가장 가까운
|
||
TIN 정점의 Z와 보정된 DEM Z를 선형 블렌드.
|
||
"""
|
||
from __future__ import annotations
|
||
|
||
import hashlib
|
||
import io
|
||
import math
|
||
from dataclasses import dataclass
|
||
from pathlib import Path
|
||
from collections.abc import Callable
|
||
|
||
import numpy as np
|
||
import pyproj
|
||
import pyvista as pv
|
||
import requests
|
||
from PIL import Image
|
||
from scipy.spatial import Delaunay, cKDTree
|
||
|
||
|
||
# --- AWS Open Terrain Tiles (terrarium) ---------------------------------
|
||
_AWS_TERRARIUM_URL = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png"
|
||
_HEADERS = {
|
||
"User-Agent": "Mozilla/5.0 (S-CANVAS DEM extender)",
|
||
"Accept": "image/png,image/*;q=0.8",
|
||
}
|
||
|
||
|
||
@dataclass
|
||
class DemExtendResult:
|
||
mesh: pv.PolyData
|
||
n_points: int
|
||
n_faces: int
|
||
buffer_m: float
|
||
source: str
|
||
vertical_offset_m: float = 0.0
|
||
zoom: int = 0
|
||
grid_step_m: float = 0.0
|
||
feather_m: float = 0.0
|
||
info: str = ""
|
||
|
||
|
||
def _latlon_to_tile(lat: float, lon: float, zoom: int) -> tuple[int, int]:
|
||
lat_rad = math.radians(lat)
|
||
n = 2 ** zoom
|
||
x = int((lon + 180.0) / 360.0 * n)
|
||
y = int((1.0 - math.log(math.tan(lat_rad) + 1.0 / math.cos(lat_rad)) / math.pi) / 2.0 * n)
|
||
return x, y
|
||
|
||
|
||
def _tile_to_latlon(x: int, y: int, zoom: int) -> tuple[float, float]:
|
||
n = 2 ** zoom
|
||
lon = x / n * 360.0 - 180.0
|
||
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
|
||
return math.degrees(lat_rad), lon
|
||
|
||
|
||
def _terrarium_decode(rgb: np.ndarray) -> np.ndarray:
|
||
"""(H,W,3) uint8 RGB → (H,W) float elevation in meters."""
|
||
r = rgb[..., 0].astype(np.float64)
|
||
g = rgb[..., 1].astype(np.float64)
|
||
b = rgb[..., 2].astype(np.float64)
|
||
return (r * 256.0 + g + b / 256.0) - 32768.0
|
||
|
||
|
||
def _bbox_cache_key(min_lat: float, min_lon: float, max_lat: float, max_lon: float, zoom: int) -> str:
|
||
raw = f"{min_lat:.6f}_{min_lon:.6f}_{max_lat:.6f}_{max_lon:.6f}_z{zoom}"
|
||
return hashlib.sha1(raw.encode("utf-8")).hexdigest()[:16]
|
||
|
||
|
||
def fetch_terrarium_grid(min_lat: float, min_lon: float, max_lat: float, max_lon: float,
|
||
zoom: int = 13,
|
||
cache_dir: str | Path = "cache/dem",
|
||
log_fn: Callable[[str], None] = print,
|
||
timeout: float = 15.0) -> tuple[np.ndarray, tuple[float, float, float, float]]:
|
||
"""AWS terrarium 타일을 BBOX 범위로 받아 합쳐 (elev_grid, bounds) 반환.
|
||
|
||
Returns:
|
||
elev_grid: (H,W) float64, 위→아래로 lat 감소, 좌→우로 lon 증가.
|
||
bounds: (grid_lat_max, grid_lon_min, grid_lat_min, grid_lon_max) — 타일
|
||
경계에 정렬된 실제 격자 범위.
|
||
"""
|
||
cache_dir = Path(cache_dir)
|
||
cache_dir.mkdir(parents=True, exist_ok=True)
|
||
cache_key = _bbox_cache_key(min_lat, min_lon, max_lat, max_lon, zoom)
|
||
cache_path = cache_dir / f"terrarium_{cache_key}.png"
|
||
bounds_path = cache_dir / f"terrarium_{cache_key}.bounds.txt"
|
||
|
||
x_min, y_min = _latlon_to_tile(max_lat, min_lon, zoom)
|
||
x_max, y_max = _latlon_to_tile(min_lat, max_lon, zoom)
|
||
cols = x_max - x_min + 1
|
||
rows = y_max - y_min + 1
|
||
|
||
grid_lat_max, grid_lon_min = _tile_to_latlon(x_min, y_min, zoom)
|
||
grid_lat_min, grid_lon_max = _tile_to_latlon(x_max + 1, y_max + 1, zoom)
|
||
bounds = (grid_lat_max, grid_lon_min, grid_lat_min, grid_lon_max)
|
||
|
||
if cache_path.exists() and bounds_path.exists():
|
||
try:
|
||
img = Image.open(cache_path).convert("RGB")
|
||
arr = np.asarray(img, dtype=np.uint8)
|
||
elev = _terrarium_decode(arr)
|
||
log_fn(f" [DEM] 캐시 사용: {cache_path.name} ({arr.shape[1]}x{arr.shape[0]}px, {cols}x{rows} tiles z{zoom})")
|
||
return elev, bounds
|
||
except Exception as e:
|
||
log_fn(f" [DEM] 캐시 로드 실패 ({e}), 재다운로드합니다.")
|
||
|
||
tile_size = 256
|
||
merged = Image.new("RGB", (cols * tile_size, rows * tile_size), (128, 0, 0))
|
||
ok = 0
|
||
fail = 0
|
||
for ty in range(y_min, y_max + 1):
|
||
for tx in range(x_min, x_max + 1):
|
||
url = _AWS_TERRARIUM_URL.format(z=zoom, x=tx, y=ty)
|
||
try:
|
||
resp = requests.get(url, headers=_HEADERS, timeout=timeout)
|
||
if resp.status_code == 200 and len(resp.content) > 200:
|
||
tile_img = Image.open(io.BytesIO(resp.content)).convert("RGB")
|
||
merged.paste(tile_img, ((tx - x_min) * tile_size, (ty - y_min) * tile_size))
|
||
ok += 1
|
||
else:
|
||
fail += 1
|
||
except requests.exceptions.RequestException:
|
||
fail += 1
|
||
|
||
log_fn(f" [DEM] 타일: 성공 {ok}장, 실패 {fail}장 (z{zoom}, {cols}x{rows})")
|
||
if ok == 0:
|
||
raise RuntimeError("AWS terrarium 타일을 한 장도 받지 못했습니다. 네트워크를 확인하세요.")
|
||
|
||
merged.save(cache_path)
|
||
bounds_path.write_text(
|
||
f"{grid_lat_max} {grid_lon_min} {grid_lat_min} {grid_lon_max} z{zoom}\n",
|
||
encoding="utf-8",
|
||
)
|
||
|
||
arr = np.asarray(merged, dtype=np.uint8)
|
||
elev = _terrarium_decode(arr)
|
||
return elev, bounds
|
||
|
||
|
||
def _sample_grid_bilinear(elev: np.ndarray,
|
||
bounds: tuple[float, float, float, float],
|
||
lats: np.ndarray, lons: np.ndarray) -> np.ndarray:
|
||
"""(H,W) 격자에서 (lats, lons) 점들의 Z를 bilinear 샘플링."""
|
||
grid_lat_max, grid_lon_min, grid_lat_min, grid_lon_max = bounds
|
||
h, w = elev.shape
|
||
fx = (lons - grid_lon_min) / (grid_lon_max - grid_lon_min) * (w - 1)
|
||
fy = (grid_lat_max - lats) / (grid_lat_max - grid_lat_min) * (h - 1)
|
||
fx = np.clip(fx, 0.0, w - 1 - 1e-9)
|
||
fy = np.clip(fy, 0.0, h - 1 - 1e-9)
|
||
x0 = np.floor(fx).astype(np.int64); x1 = x0 + 1
|
||
y0 = np.floor(fy).astype(np.int64); y1 = y0 + 1
|
||
x1 = np.clip(x1, 0, w - 1)
|
||
y1 = np.clip(y1, 0, h - 1)
|
||
tx = fx - x0
|
||
ty = fy - y0
|
||
v00 = elev[y0, x0]; v01 = elev[y0, x1]
|
||
v10 = elev[y1, x0]; v11 = elev[y1, x1]
|
||
v0 = v00 * (1 - tx) + v01 * tx
|
||
v1 = v10 * (1 - tx) + v11 * tx
|
||
return v0 * (1 - ty) + v1 * ty
|
||
|
||
|
||
def _try_local_geotiff(cache_dir: Path) -> tuple[np.ndarray, tuple[float, float, float, float], str] | None:
|
||
"""cache/dem/local.tif 가 있으면 읽어 (elev, (lat_max,lon_min,lat_min,lon_max), label) 반환.
|
||
|
||
rasterio가 없거나 파일 없으면 None.
|
||
"""
|
||
tif = cache_dir / "local.tif"
|
||
if not tif.exists():
|
||
return None
|
||
try:
|
||
import rasterio # type: ignore
|
||
from rasterio.warp import transform_bounds # type: ignore
|
||
except ImportError:
|
||
return None
|
||
try:
|
||
with rasterio.open(tif) as ds:
|
||
arr = ds.read(1).astype(np.float64)
|
||
# nodata 처리
|
||
nodata = ds.nodata
|
||
if nodata is not None:
|
||
arr = np.where(arr == nodata, np.nan, arr)
|
||
src_bounds = ds.bounds
|
||
src_crs = ds.crs
|
||
lon_min, lat_min, lon_max, lat_max = transform_bounds(
|
||
src_crs, "EPSG:4326",
|
||
src_bounds.left, src_bounds.bottom, src_bounds.right, src_bounds.top,
|
||
densify_pts=21,
|
||
)
|
||
# rasterio는 top-to-bottom: 첫 행이 lat_max
|
||
bounds = (lat_max, lon_min, lat_min, lon_max)
|
||
return arr, bounds, "local_geotiff"
|
||
except Exception:
|
||
return None
|
||
|
||
|
||
def _generate_ring_points(inner_xyxy: tuple[float, float, float, float],
|
||
outer_xyxy: tuple[float, float, float, float],
|
||
grid_step: float) -> np.ndarray:
|
||
"""도넛(ring) 영역을 규칙 격자로 채운 XY 좌표 반환. inner bbox 내부는 제외.
|
||
|
||
Legacy(사각 bbox 경계): hull을 못 얻는 경우의 폴백.
|
||
"""
|
||
ix0, iy0, ix1, iy1 = inner_xyxy
|
||
ox0, oy0, ox1, oy1 = outer_xyxy
|
||
xs = np.arange(ox0, ox1 + grid_step * 0.5, grid_step)
|
||
ys = np.arange(oy0, oy1 + grid_step * 0.5, grid_step)
|
||
gx, gy = np.meshgrid(xs, ys)
|
||
pts = np.column_stack([gx.ravel(), gy.ravel()])
|
||
inside = (pts[:, 0] > ix0) & (pts[:, 0] < ix1) & (pts[:, 1] > iy0) & (pts[:, 1] < iy1)
|
||
ring_pts = pts[~inside]
|
||
# inner 경계 자체에 정확히 놓이는 점들도 추가 (seam 정합용)
|
||
ex = np.arange(ix0, ix1 + grid_step * 0.5, grid_step)
|
||
edge_top = np.column_stack([ex, np.full_like(ex, iy1)])
|
||
edge_bot = np.column_stack([ex, np.full_like(ex, iy0)])
|
||
ey = np.arange(iy0, iy1 + grid_step * 0.5, grid_step)
|
||
edge_left = np.column_stack([np.full_like(ey, ix0), ey])
|
||
edge_right = np.column_stack([np.full_like(ey, ix1), ey])
|
||
return np.vstack([ring_pts, edge_top, edge_bot, edge_left, edge_right])
|
||
|
||
|
||
def _compute_tin_hull_projected(tin_xyz_zerobased: np.ndarray,
|
||
origin: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
|
||
"""TIN 점들의 XY 컨벡스 헐을 투영 CRS 좌표로 반환.
|
||
|
||
Returns:
|
||
hull_xy_abs: (N,2) 반시계 순 헐 정점 XY (projected CRS, 절대좌표)
|
||
hull_z_abs: (N,) 해당 정점의 TIN Z (절대 높이)
|
||
"""
|
||
from scipy.spatial import ConvexHull
|
||
xy_abs = tin_xyz_zerobased[:, :2] + origin[:2]
|
||
hull = ConvexHull(xy_abs)
|
||
hull_xy_abs = xy_abs[hull.vertices]
|
||
hull_z_abs = tin_xyz_zerobased[hull.vertices, 2] + origin[2]
|
||
return hull_xy_abs, hull_z_abs
|
||
|
||
|
||
def _densify_polygon(poly_xy: np.ndarray, max_seg_len: float) -> np.ndarray:
|
||
"""다각형 엣지를 max_seg_len 이하로 선형 보간해 세분화."""
|
||
n = len(poly_xy)
|
||
out = []
|
||
for i in range(n):
|
||
p0 = poly_xy[i]
|
||
p1 = poly_xy[(i + 1) % n]
|
||
seg = p1 - p0
|
||
L = float(np.hypot(seg[0], seg[1]))
|
||
k = max(1, int(np.ceil(L / max(max_seg_len, 1e-6))))
|
||
out.extend(p0 + seg * t for t in np.linspace(0.0, 1.0, k, endpoint=False))
|
||
return np.asarray(out, dtype=np.float64)
|
||
|
||
|
||
def _generate_ring_points_hull(outer_xyxy: tuple[float, float, float, float],
|
||
inner_hull_xy: np.ndarray,
|
||
grid_step: float,
|
||
precomputed_boundary: np.ndarray | None = None,
|
||
) -> tuple[np.ndarray, int]:
|
||
"""hull 바깥쪽 격자점 + hull 엣지 세분화 점들을 결합해 링 XY 반환.
|
||
|
||
Args:
|
||
precomputed_boundary: (N,2) TIN mesh에 이미 존재하는 bbox 변 정점 XY.
|
||
제공되면 densify를 건너뛰고 이 점들을 그대로 inner 경계로 사용 →
|
||
TIN 정점과 DEM 링 정점이 **공유 XY**가 되어 T-vertex/fin 원천 제거.
|
||
|
||
Returns:
|
||
pts: (M, 2) 전체 링 점
|
||
n_hull_boundary: 끝 n개가 hull 경계 점. Z 오버라이드/feather용.
|
||
"""
|
||
ox0, oy0, ox1, oy1 = outer_xyxy
|
||
xs = np.arange(ox0, ox1 + grid_step * 0.5, grid_step)
|
||
ys = np.arange(oy0, oy1 + grid_step * 0.5, grid_step)
|
||
gx, gy = np.meshgrid(xs, ys)
|
||
pts = np.column_stack([gx.ravel(), gy.ravel()])
|
||
# inner_hull = TIN bbox 4 모서리. bbox 내부(margin 포함)의 격자점 제거.
|
||
# 여유 margin = grid_step*0.5 — bbox 선에 너무 가까운 격자점도 제거해 링 삼각형이
|
||
# bbox 내부로 밀고 들어가는 현상(error1.png 톱니 fin) 차단.
|
||
ix0 = float(np.min(inner_hull_xy[:, 0])); ix1 = float(np.max(inner_hull_xy[:, 0]))
|
||
iy0 = float(np.min(inner_hull_xy[:, 1])); iy1 = float(np.max(inner_hull_xy[:, 1]))
|
||
guard = grid_step * 0.5
|
||
near_inner = (
|
||
(pts[:, 0] > ix0 - guard) & (pts[:, 0] < ix1 + guard)
|
||
& (pts[:, 1] > iy0 - guard) & (pts[:, 1] < iy1 + guard)
|
||
)
|
||
outside_grid = pts[~near_inner]
|
||
|
||
if precomputed_boundary is not None and len(precomputed_boundary) >= 3:
|
||
# **경계 densify 추가** — TIN bbox 공유정점 사이 간격이 외곽 grid_step 보다
|
||
# 크면 링 Delaunay 가 bbox 를 가로질러 큰 삼각형을 만들어 TIN 영역을 덮음.
|
||
# 따라서 공유정점을 유지하되, 인접 공유정점 사이 gap > grid_step 구간은
|
||
# 외곽 격자와 같은 밀도로 세분화해 삽입.
|
||
raw = np.asarray(precomputed_boundary, dtype=np.float64)
|
||
b_tol = max(ix1 - ix0, iy1 - iy0) * 1e-4 + 1e-3
|
||
on_bot = np.abs(raw[:, 1] - iy0) < b_tol
|
||
on_top = np.abs(raw[:, 1] - iy1) < b_tol
|
||
on_lft = np.abs(raw[:, 0] - ix0) < b_tol
|
||
on_rgt = np.abs(raw[:, 0] - ix1) < b_tol
|
||
parts: list[np.ndarray] = []
|
||
for mask, sort_col, fixed_col, fixed_val in [
|
||
(on_bot, 0, 1, iy0), (on_top, 0, 1, iy1),
|
||
(on_lft, 1, 0, ix0), (on_rgt, 1, 0, ix1),
|
||
]:
|
||
if mask.sum() == 0:
|
||
continue
|
||
side = raw[mask]
|
||
# 변 sort + 중복 제거
|
||
side = side[np.argsort(side[:, sort_col])]
|
||
_, uq = np.unique(np.round(side[:, sort_col], 3), return_index=True)
|
||
side = side[np.sort(uq)]
|
||
parts.append(side)
|
||
if len(side) < 2:
|
||
continue
|
||
main = side[:, sort_col]
|
||
gaps = np.diff(main)
|
||
for k, g in enumerate(gaps):
|
||
if g > grid_step:
|
||
n_add = int(np.ceil(g / grid_step)) - 1
|
||
if n_add < 1:
|
||
continue
|
||
mids = np.linspace(main[k], main[k + 1], n_add + 2)[1:-1]
|
||
add = np.zeros((len(mids), 2))
|
||
add[:, sort_col] = mids
|
||
add[:, fixed_col] = fixed_val
|
||
parts.append(add)
|
||
boundary_pts = np.unique(np.round(np.vstack(parts), 3), axis=0) if parts else raw
|
||
else:
|
||
boundary_pts = _densify_polygon(inner_hull_xy, grid_step)
|
||
|
||
# 최종 중복 제거 (outside_grid 점 중 boundary 와 매우 가까운 것 제거 → sliver 방지)
|
||
if len(boundary_pts) > 0 and len(outside_grid) > 0:
|
||
from scipy.spatial import cKDTree as _cKDT
|
||
bt = _cKDT(boundary_pts)
|
||
d, _ = bt.query(outside_grid, k=1)
|
||
outside_grid = outside_grid[d > grid_step * 0.25]
|
||
|
||
n_hull_boundary = len(boundary_pts)
|
||
all_pts = np.vstack([outside_grid, boundary_pts])
|
||
return all_pts, n_hull_boundary
|
||
|
||
|
||
def build_extended_terrain_ring(projected_bounds: tuple[float, float, float, float],
|
||
origin: np.ndarray,
|
||
src_crs: str,
|
||
buffer_m: float = 1000.0,
|
||
tin_xyz_zerobased: np.ndarray | None = None,
|
||
grid_step_m: float | None = None,
|
||
feather_m: float = 80.0,
|
||
zoom: int = 13,
|
||
cache_dir: str | Path = "cache/dem",
|
||
source: str = "auto",
|
||
use_hull_boundary: bool = True,
|
||
datum_offset_override: float | None = None,
|
||
elev_grid_override: np.ndarray | None = None,
|
||
grid_bounds_override: tuple[float, float, float, float] | None = None,
|
||
log_fn: Callable[[str], None] = print) -> DemExtendResult:
|
||
"""DXF bbox 외곽을 실제 DEM으로 확장한 메시 생성.
|
||
|
||
Args:
|
||
projected_bounds: (xmin, ymin, xmax, ymax) 투영 CRS 원본 좌표.
|
||
origin: (ox, oy, oz) — TIN 생성 시 zero-basing에 사용한 offset.
|
||
src_crs: 투영 CRS 문자열 (예: "EPSG:5187").
|
||
buffer_m: 외곽으로 확장할 거리(미터).
|
||
tin_xyz_zerobased: (N,3) zero-based TIN 정점. 제공되면 수직 datum 보정
|
||
+ feathering + (use_hull_boundary=True면) 컨벡스 헐 경계 스냅에 사용.
|
||
None이면 DEM Z 원본을 그대로 사용하고 bbox 경계 사용.
|
||
grid_step_m: 외곽 링 격자 간격. None이면 buffer/25 기준 자동 (최소 10m).
|
||
feather_m: inner 경계에서 이 거리 이내는 TIN Z로 블렌드.
|
||
zoom: AWS terrarium zoom level (13 ≈ 19m/px at mid-lat).
|
||
cache_dir: 캐시 디렉토리.
|
||
source: "auto" | "aws_terrain" | "local_geotiff".
|
||
use_hull_boundary: True면 링 inner 경계를 TIN 컨벡스 헐로 스냅. 헐 정점은
|
||
TIN Z를 그대로 써 공유 정점으로 삼아 gap/T-vertex 제거. False면 사각 bbox
|
||
경계(legacy).
|
||
datum_offset_override: None이 아니면 이 값을 DEM Z에서 차감(vertical datum 보정).
|
||
TIN/내부채움과 **같은 offset**을 써야 bbox seam에 Z 단차가 없음. 호출자가
|
||
create_tin_from_dxf의 `_dem_datum_offset`을 그대로 넘기는 것을 권장.
|
||
elev_grid_override, grid_bounds_override: None이 아니면 DEM 타일을 **재사용**
|
||
(TIN 생성 시 이미 받은 격자 재활용 → datum 일관성 + 속도).
|
||
log_fn: 로그 callback.
|
||
|
||
Returns:
|
||
DemExtendResult: 메시 + 메타데이터.
|
||
"""
|
||
cache_dir = Path(cache_dir)
|
||
cache_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
xmin, ymin, xmax, ymax = projected_bounds
|
||
inner = (xmin, ymin, xmax, ymax)
|
||
outer = (xmin - buffer_m, ymin - buffer_m, xmax + buffer_m, ymax + buffer_m)
|
||
|
||
if grid_step_m is None:
|
||
grid_step_m = max(10.0, buffer_m / 25.0)
|
||
|
||
# --- 좌표계 진단 로그 -------------------------------------------------
|
||
# TIN/DXF는 src_crs(예: EPSG:5187), DEM 타일은 WGS84로 받아와 Z만 샘플링.
|
||
# 결과 mesh의 XY는 src_crs 평면 그대로이므로 TIN과 동일 CRS. 혹시 변환
|
||
# 오차가 의심되면 왕복 변환 거리로 확인 가능.
|
||
log_fn(f" [DEM] 좌표계: DXF={src_crs} → WGS84로 DEM 샘플링 → 결과 mesh는 {src_crs} 평면")
|
||
try:
|
||
_mid_x = 0.5 * (xmin + xmax); _mid_y = 0.5 * (ymin + ymax)
|
||
_to_wgs_t = pyproj.Transformer.from_crs(src_crs, "EPSG:4326", always_xy=True)
|
||
_to_proj_t = pyproj.Transformer.from_crs("EPSG:4326", src_crs, always_xy=True)
|
||
_lon, _lat = _to_wgs_t.transform(_mid_x, _mid_y)
|
||
_rx, _ry = _to_proj_t.transform(_lon, _lat)
|
||
_rt_err = float(((_rx - _mid_x) ** 2 + (_ry - _mid_y) ** 2) ** 0.5)
|
||
log_fn(f" [DEM] TIN bbox: X=[{xmin:.1f}, {xmax:.1f}] Y=[{ymin:.1f}, {ymax:.1f}] "
|
||
f"(중심 WGS84 Lon={_lon:.5f} Lat={_lat:.5f})")
|
||
log_fn(f" [DEM] 좌표 왕복 변환 오차: {_rt_err:.4f}m (중심점 기준, 0에 가까워야 정상)")
|
||
except Exception as _ce:
|
||
log_fn(f" [DEM] 좌표계 진단 실패: {_ce}")
|
||
|
||
# 1) 링 점 생성 — **TIN bbox를 inner polygon으로 사용** (convex hull 대신).
|
||
# 이전 hull 모드는 TIN pts의 bbox ≠ convex hull 이라 두 메시가 hull 외곽과
|
||
# bbox 사이 얇은 영역을 중복 덮거나 T-vertex/fin을 만들었다. bbox를 inner로
|
||
# 쓰면 DEM 링 inner 경계가 TIN bbox 선과 정확히 맞닿아 중복/fin 제거.
|
||
hull_xy_abs = None
|
||
n_hull_boundary = 0
|
||
boundary_mode = "bbox"
|
||
if use_hull_boundary and tin_xyz_zerobased is not None and len(tin_xyz_zerobased) >= 3:
|
||
try:
|
||
tin_abs = tin_xyz_zerobased + origin # (N,3) 절대좌표
|
||
tbx0 = float(tin_abs[:, 0].min()); tbx1 = float(tin_abs[:, 0].max())
|
||
tby0 = float(tin_abs[:, 1].min()); tby1 = float(tin_abs[:, 1].max())
|
||
hull_xy_abs = np.array([
|
||
[tbx0, tby0], [tbx1, tby0], [tbx1, tby1], [tbx0, tby1],
|
||
], dtype=np.float64)
|
||
# **TIN의 실제 bbox 변 정점**만 선별해 DEM 링 inner 경계로 사용 →
|
||
# 두 메시가 동일 XY의 공유 정점을 가지므로 T-vertex/fin 원천 제거.
|
||
bbox_tol_tin = max(tbx1 - tbx0, tby1 - tby0) * 1e-4 + 1e-3
|
||
on_bbox_mask = (
|
||
(np.abs(tin_abs[:, 0] - tbx0) < bbox_tol_tin)
|
||
| (np.abs(tin_abs[:, 0] - tbx1) < bbox_tol_tin)
|
||
| (np.abs(tin_abs[:, 1] - tby0) < bbox_tol_tin)
|
||
| (np.abs(tin_abs[:, 1] - tby1) < bbox_tol_tin)
|
||
)
|
||
tin_bbox_vertices_xy = tin_abs[on_bbox_mask, :2]
|
||
if len(tin_bbox_vertices_xy) < 4:
|
||
# bbox 변에 정점이 너무 적으면 폴백(densify 기본 로직)
|
||
tin_bbox_vertices_xy = None
|
||
|
||
ring_xy, n_hull_boundary = _generate_ring_points_hull(
|
||
outer, hull_xy_abs, grid_step_m,
|
||
precomputed_boundary=tin_bbox_vertices_xy,
|
||
)
|
||
n_shared = len(tin_bbox_vertices_xy) if tin_bbox_vertices_xy is not None else 0
|
||
boundary_mode = (f"TIN bbox 공유정점 {n_shared}개" if n_shared > 0
|
||
else "TIN bbox 4모서리 densify 폴백")
|
||
except Exception as e:
|
||
log_fn(f" [DEM] bbox 경계 구성 실패 ({e}) — 단순 bbox 링 폴백")
|
||
ring_xy = _generate_ring_points(inner, outer, grid_step_m)
|
||
hull_xy_abs = None; n_hull_boundary = 0
|
||
else:
|
||
ring_xy = _generate_ring_points(inner, outer, grid_step_m)
|
||
log_fn(f" [DEM] 링 격자점: {len(ring_xy)}개 (step={grid_step_m:.1f}m, "
|
||
f"buffer={buffer_m:.0f}m, 경계={boundary_mode})")
|
||
|
||
# 2) WGS84로 변환
|
||
to_wgs = pyproj.Transformer.from_crs(src_crs, "EPSG:4326", always_xy=True)
|
||
lons, lats = to_wgs.transform(ring_xy[:, 0], ring_xy[:, 1])
|
||
lons = np.asarray(lons); lats = np.asarray(lats)
|
||
|
||
# 3) DEM 소스 결정 및 fetch
|
||
resolved_source = source
|
||
elev_grid = None
|
||
grid_bounds = None
|
||
|
||
# 3-0) 호출자 override — TIN densify에서 이미 받은 격자 재사용 (datum 일관성)
|
||
if elev_grid_override is not None and grid_bounds_override is not None:
|
||
elev_grid = elev_grid_override
|
||
grid_bounds = grid_bounds_override
|
||
resolved_source = "reused_from_caller"
|
||
log_fn(f" [DEM] 호출자 override 격자 재사용 ({elev_grid.shape[1]}x{elev_grid.shape[0]}) "
|
||
f"— TIN 생성 때 받은 것과 **동일 datum 보장**")
|
||
|
||
if elev_grid is None and source in ("auto", "local_geotiff"):
|
||
local = _try_local_geotiff(cache_dir)
|
||
if local is not None:
|
||
elev_grid, grid_bounds, _ = local
|
||
resolved_source = "local_geotiff"
|
||
log_fn(f" [DEM] 로컬 GeoTIFF 사용: cache/dem/local.tif ({elev_grid.shape[1]}x{elev_grid.shape[0]})")
|
||
|
||
if elev_grid is None:
|
||
ox0, oy0, ox1, oy1 = outer
|
||
margin = max(grid_step_m * 2, 50.0)
|
||
corners_x = np.array([ox0 - margin, ox1 + margin, ox0 - margin, ox1 + margin])
|
||
corners_y = np.array([oy0 - margin, oy0 - margin, oy1 + margin, oy1 + margin])
|
||
cx_lon, cx_lat = to_wgs.transform(corners_x, corners_y)
|
||
min_lat, max_lat = float(np.min(cx_lat)), float(np.max(cx_lat))
|
||
min_lon, max_lon = float(np.min(cx_lon)), float(np.max(cx_lon))
|
||
elev_grid, grid_bounds = fetch_terrarium_grid(
|
||
min_lat, min_lon, max_lat, max_lon,
|
||
zoom=zoom, cache_dir=cache_dir, log_fn=log_fn,
|
||
)
|
||
resolved_source = "aws_terrain"
|
||
|
||
# 4) 링 점 고도 샘플링
|
||
z_dem_raw = _sample_grid_bilinear(elev_grid, grid_bounds, lats, lons)
|
||
if np.any(np.isnan(z_dem_raw)):
|
||
med = float(np.nanmedian(z_dem_raw))
|
||
z_dem_raw = np.where(np.isnan(z_dem_raw), med, z_dem_raw)
|
||
|
||
# 4a) 전역 outlier 클립 + **NaN 즉시 채움**.
|
||
# terrarium 디코딩 실패·타일 엣지에서 극단값/NaN 이 링 삼각형을 수십 m
|
||
# 솟게 만들었던 현상 방지. IQR 기반 + 절대범위(−500m~9000m) 가드.
|
||
finite = np.isfinite(z_dem_raw)
|
||
if finite.any():
|
||
vals = z_dem_raw[finite]
|
||
q1, q3 = np.percentile(vals, [5, 95])
|
||
iqr = max(q3 - q1, 1.0)
|
||
lo = max(-500.0, q1 - iqr * 3.0)
|
||
hi = min(9000.0, q3 + iqr * 3.0)
|
||
med_all = float(np.median(vals))
|
||
z_dem_raw = np.where(finite, z_dem_raw, med_all)
|
||
z_dem_raw = np.clip(z_dem_raw, lo, hi)
|
||
else:
|
||
z_dem_raw = np.zeros_like(z_dem_raw)
|
||
log_fn(" [DEM] 경고: 샘플 전체가 NaN — 0으로 채움 (DEM 타일 fetch 점검 필요)")
|
||
|
||
# 4b) 국소 스파이크 필터 — 각 점의 Z가 반경 R 내 이웃 median 대비 3×MAD 이상
|
||
# 벗어나면 median으로 치환. 경계 값이 튀는 원인(인접 타일 디코딩 차이)을 직접 제거.
|
||
try:
|
||
radius = max(grid_step_m * 3.5, 60.0)
|
||
tree_local = cKDTree(ring_xy)
|
||
neigh_idxs = tree_local.query_ball_point(ring_xy, r=radius)
|
||
cleaned = z_dem_raw.copy()
|
||
spike_count = 0
|
||
for i, nbrs in enumerate(neigh_idxs):
|
||
if len(nbrs) < 5:
|
||
continue
|
||
neigh_z = z_dem_raw[nbrs]
|
||
med_ = float(np.median(neigh_z))
|
||
mad_ = float(np.median(np.abs(neigh_z - med_))) + 1e-6
|
||
if abs(z_dem_raw[i] - med_) > 4.0 * mad_:
|
||
cleaned[i] = med_
|
||
spike_count += 1
|
||
if spike_count:
|
||
log_fn(f" [DEM] 스파이크 {spike_count}개 제거 (radius={radius:.0f}m, 4×MAD)")
|
||
z_dem_raw = cleaned
|
||
except Exception as _se:
|
||
log_fn(f" [DEM] 스파이크 필터 경고: {_se}")
|
||
|
||
# 5) 수직 datum 보정 + feathering
|
||
# **override 우선** — TIN 생성 시 계산한 offset을 그대로 쓰면 bbox seam에 Z 단차 0.
|
||
vertical_offset = 0.0
|
||
z_final = z_dem_raw.copy()
|
||
if datum_offset_override is not None:
|
||
vertical_offset = float(datum_offset_override)
|
||
z_final = z_dem_raw - vertical_offset
|
||
log_fn(f" [DEM] 수직 datum **override** offset={vertical_offset:+.2f}m "
|
||
f"(TIN 생성 시 계산한 값 재사용 → bbox seam Z 단차 0)")
|
||
if tin_xyz_zerobased is not None and len(tin_xyz_zerobased) > 0:
|
||
tin_abs = tin_xyz_zerobased + origin
|
||
tree = cKDTree(tin_abs[:, :2])
|
||
ring_abs = np.column_stack([ring_xy[:, 0], ring_xy[:, 1]])
|
||
dists, idxs = tree.query(ring_abs, k=1)
|
||
near_mask = dists < max(feather_m * 1.5, 50.0)
|
||
if np.any(near_mask) and datum_offset_override is None:
|
||
tin_z_near = tin_abs[idxs[near_mask], 2]
|
||
dem_z_near = z_dem_raw[near_mask]
|
||
diffs = dem_z_near - tin_z_near
|
||
vertical_offset = float(np.median(diffs))
|
||
log_fn(f" [DEM] 수직 datum 자동 보정: offset={vertical_offset:+.2f}m "
|
||
f"(비교 점 {int(near_mask.sum())}개)")
|
||
z_final = z_dem_raw - vertical_offset
|
||
if feather_m > 0 and np.any(near_mask):
|
||
# smoothstep 블렌드 (선형보다 부드러운 C1 전이 — 경계 Z 단차 제거).
|
||
# t=0(경계)에서 TIN Z, t=1(feather_m 바깥)에서 DEM Z. 사이 값은
|
||
# 3t²−2t³ 곡선으로 양 끝에서 미분 0 → 절벽 같은 Z 점프 제거.
|
||
raw_t = np.clip(dists / feather_m, 0.0, 1.0)
|
||
t = raw_t * raw_t * (3.0 - 2.0 * raw_t) # smoothstep
|
||
tin_z_all = tin_abs[idxs, 2]
|
||
blended = (1.0 - t) * tin_z_all + t * z_final
|
||
z_final = np.where(dists < feather_m, blended, z_final)
|
||
|
||
# 5b) **Outer smooth blend** — outer 경계 근처를 이웃 평균으로 부드럽게 평활.
|
||
# 이전 "하위 10% 분위수" ramp-down 은 점마다 강제 하강량이 달라 톱니 fin 을
|
||
# 만들었음(error1.png). 이번에는:
|
||
# - 이웃 **평균** Z 를 타겟으로(분위수 NO)
|
||
# - smoothstep 가중치 (경계 1, ramp 밖 0)
|
||
# - 이웃 반경 = grid_step*6 로 넓게 → fin 유발하는 국소 편차 제거
|
||
# - 최대 하강/상승 폭 제한 (|Δ| ≤ grid_step*0.8) → 튀는 값 차단
|
||
try:
|
||
ox0, oy0, ox1, oy1 = outer
|
||
dist_to_outer = np.minimum.reduce([
|
||
ring_xy[:, 0] - ox0,
|
||
ox1 - ring_xy[:, 0],
|
||
ring_xy[:, 1] - oy0,
|
||
oy1 - ring_xy[:, 1],
|
||
])
|
||
dist_to_outer = np.maximum(dist_to_outer, 0.0)
|
||
ramp_width = max(float(grid_step_m) * 4.0, float(buffer_m) * 0.20)
|
||
in_ramp = dist_to_outer < ramp_width
|
||
n_ramp = int(in_ramp.sum())
|
||
if n_ramp > 0:
|
||
tree_outer = cKDTree(ring_xy)
|
||
r_neigh = max(float(grid_step_m) * 6.0, 120.0)
|
||
idxs_ramp = np.where(in_ramp)[0]
|
||
nbrs_ramp = tree_outer.query_ball_point(ring_xy[idxs_ramp], r=r_neigh)
|
||
z_target = z_final[idxs_ramp].copy()
|
||
for k, nb in enumerate(nbrs_ramp):
|
||
if len(nb) < 5:
|
||
continue
|
||
z_target[k] = float(np.mean(z_final[nb]))
|
||
# smoothstep weight: outer(t=0) → w=1, ramp 밖(t=1) → w=0
|
||
t = np.clip(dist_to_outer[idxs_ramp] / ramp_width, 0.0, 1.0)
|
||
u = 1.0 - t
|
||
w = u * u * (3.0 - 2.0 * u)
|
||
z_before = z_final[idxs_ramp].copy()
|
||
blended = (1.0 - w) * z_before + w * z_target
|
||
# 튀는 값 가드 — 이웃 평균과의 편차를 grid_step*0.8 이내로 제한
|
||
cap = float(grid_step_m) * 0.8
|
||
blended = np.clip(blended, z_before - cap, z_before + cap)
|
||
z_final[idxs_ramp] = blended
|
||
log_fn(
|
||
f" [DEM] outer smooth blend {n_ramp}개 점 "
|
||
f"(ramp폭={ramp_width:.0f}m, r={r_neigh:.0f}m, |Δ|≤{cap:.1f}m) "
|
||
f"— 이웃 평균 블렌드, fin/톱니 제거"
|
||
)
|
||
except Exception as _rd:
|
||
log_fn(f" [DEM] outer smooth blend 경고: {_rd}")
|
||
|
||
# 5c) bbox 경계 densified 점 Z 설정.
|
||
# (1) **TIN 공유 정점** (= 실제 TIN vertex 그대로) 은 TIN Z 그대로.
|
||
# (2) **경계 보간 densify** 신규 점(= TIN vertex 없는 bbox 변 위 점)은 양옆
|
||
# TIN vertex 사이 선형 보간 Z (거리 가중 k=2) 로 설정 → TIN 경계선이
|
||
# 그대로 연장되어 seam 에서 Z 점프 0.
|
||
if n_hull_boundary > 0 and hull_xy_abs is not None and tin_xyz_zerobased is not None:
|
||
n_total = len(ring_xy)
|
||
boundary_start = n_total - n_hull_boundary
|
||
boundary_xy = ring_xy[boundary_start:]
|
||
tin_abs_full = tin_xyz_zerobased + origin
|
||
tin_tree_full = cKDTree(tin_abs_full[:, :2])
|
||
d_near, i_near = tin_tree_full.query(boundary_xy, k=2)
|
||
# k=1 이면 거의 일치하는 TIN vertex → 그대로 사용
|
||
# k=2 가중평균으로 densify 중간 점 Z 를 선형 보간
|
||
exact = d_near[:, 0] < 1e-3
|
||
z_b = np.empty(len(boundary_xy), dtype=np.float64)
|
||
z_b[exact] = tin_abs_full[i_near[exact, 0], 2]
|
||
if (~exact).any():
|
||
w0 = 1.0 / np.maximum(d_near[~exact, 0], 1e-6) ** 2
|
||
w1 = 1.0 / np.maximum(d_near[~exact, 1], 1e-6) ** 2
|
||
z0 = tin_abs_full[i_near[~exact, 0], 2]
|
||
z1 = tin_abs_full[i_near[~exact, 1], 2]
|
||
z_b[~exact] = (w0 * z0 + w1 * z1) / (w0 + w1)
|
||
z_final[boundary_start:] = z_b
|
||
|
||
# 5d) **최종 NaN/무한 가드** — 어떤 경로로 들어온 NaN도 여기서 확정적으로 제거.
|
||
# 인접 이웃의 median 으로 채움. 이웃 없으면 전체 median.
|
||
bad = ~np.isfinite(z_final)
|
||
if bad.any():
|
||
good_idx = np.where(~bad)[0]
|
||
if len(good_idx) > 0:
|
||
tree_nan = cKDTree(ring_xy[good_idx])
|
||
_, j = tree_nan.query(ring_xy[bad], k=1)
|
||
z_final[bad] = z_final[good_idx[j]]
|
||
else:
|
||
z_final[bad] = 0.0
|
||
log_fn(f" [DEM] NaN/Inf {int(bad.sum())}개 이웃 값으로 채움")
|
||
|
||
# 5e) **링 Z Laplacian 1pass** — 톱니(fin) 마지막 평활. TIN 공유 경계 정점은
|
||
# 고정(Laplacian 에서 제외) → TIN–링 seam Z 는 변경되지 않음.
|
||
try:
|
||
tree_ring = cKDTree(ring_xy)
|
||
r_sm = max(float(grid_step_m) * 2.2, 40.0)
|
||
fixed_mask = np.zeros(len(ring_xy), dtype=bool)
|
||
if n_hull_boundary > 0:
|
||
fixed_mask[len(ring_xy) - n_hull_boundary:] = True
|
||
movable = np.where(~fixed_mask)[0]
|
||
nbrs = tree_ring.query_ball_point(ring_xy[movable], r=r_sm)
|
||
z_new = z_final.copy()
|
||
for k, nb in enumerate(nbrs):
|
||
if len(nb) < 4:
|
||
continue
|
||
vid = movable[k]
|
||
nb_arr = np.asarray(nb, dtype=np.int64)
|
||
nb_arr = nb_arr[nb_arr != vid]
|
||
if len(nb_arr) == 0:
|
||
continue
|
||
z_new[vid] = 0.5 * z_final[vid] + 0.5 * float(np.mean(z_final[nb_arr]))
|
||
z_final = z_new
|
||
except Exception as _ls:
|
||
log_fn(f" [DEM] 링 Laplacian 평활 경고: {_ls}")
|
||
|
||
# 6) Zero-basing
|
||
ring_xyz = np.column_stack([
|
||
ring_xy[:, 0] - origin[0],
|
||
ring_xy[:, 1] - origin[1],
|
||
z_final - origin[2],
|
||
])
|
||
|
||
# 7) Delaunay (XY 기준) → 도넛 triangulation
|
||
try:
|
||
tri = Delaunay(ring_xyz[:, :2])
|
||
except Exception as e:
|
||
raise RuntimeError(f"링 메시 Delaunay 실패: {e}") from e
|
||
|
||
# inner(TIN bbox) 내부 삼각형 제거 — **3중 가드**로 침범 원천 차단.
|
||
# (a) centroid 가 inner bbox 내부
|
||
# (b) **세 정점 중 하나라도 inner bbox strict 내부** (≥ strict_tol 안쪽)
|
||
# (c) **엣지 중점** 셋 중 하나라도 inner bbox strict 내부
|
||
# 이전에는 (a)+(b) 만으로 bbox 선을 따라 "세 정점 전부 경계 위"인 납작한
|
||
# 삼각형이 남아 엣지 중점이 bbox 내부로 진입하며 TIN 위를 덮었음(error.png
|
||
# 침범 증상). (c) 추가로 완전 차단.
|
||
cx0 = inner[0] - origin[0]; cy0 = inner[1] - origin[1]
|
||
cx1 = inner[2] - origin[0]; cy1 = inner[3] - origin[1]
|
||
bbox_w = max(cx1 - cx0, cy1 - cy0)
|
||
strict_tol = bbox_w * 1e-4 + 1e-3 # 수치 오차 가드 (≈mm 단위)
|
||
|
||
v0 = ring_xyz[tri.simplices[:, 0], :2]
|
||
v1 = ring_xyz[tri.simplices[:, 1], :2]
|
||
v2 = ring_xyz[tri.simplices[:, 2], :2]
|
||
centroids = (v0 + v1 + v2) / 3.0
|
||
mid01 = (v0 + v1) * 0.5
|
||
mid12 = (v1 + v2) * 0.5
|
||
mid20 = (v2 + v0) * 0.5
|
||
|
||
def _strict_in_bbox(p: np.ndarray) -> np.ndarray:
|
||
return (
|
||
(p[:, 0] > cx0 + strict_tol) & (p[:, 0] < cx1 - strict_tol)
|
||
& (p[:, 1] > cy0 + strict_tol) & (p[:, 1] < cy1 - strict_tol)
|
||
)
|
||
|
||
# (a) centroid
|
||
inside_inner_centroid = _strict_in_bbox(centroids - np.array([[0, 0]])) # strict only
|
||
# centroid 는 strict 내부까지 너무 보수적이면 seam 삼각형도 지울 수 있으므로
|
||
# 너그럽게: bbox 완전 내부만.
|
||
inside_inner_centroid = (
|
||
(centroids[:, 0] > cx0 + strict_tol) & (centroids[:, 0] < cx1 - strict_tol)
|
||
& (centroids[:, 1] > cy0 + strict_tol) & (centroids[:, 1] < cy1 - strict_tol)
|
||
)
|
||
|
||
# (b) 세 정점 중 하나라도 strict 내부
|
||
vx = ring_xyz[:, 0]; vy = ring_xyz[:, 1]
|
||
strict_inner_vert = (
|
||
(vx > cx0 + strict_tol) & (vx < cx1 - strict_tol)
|
||
& (vy > cy0 + strict_tol) & (vy < cy1 - strict_tol)
|
||
)
|
||
has_strict_vertex = (
|
||
strict_inner_vert[tri.simplices[:, 0]]
|
||
| strict_inner_vert[tri.simplices[:, 1]]
|
||
| strict_inner_vert[tri.simplices[:, 2]]
|
||
)
|
||
|
||
# (c) 엣지 중점 중 하나라도 strict 내부
|
||
edge_mid_inside = (
|
||
_strict_in_bbox(mid01) | _strict_in_bbox(mid12) | _strict_in_bbox(mid20)
|
||
)
|
||
|
||
inside_inner = inside_inner_centroid | has_strict_vertex | edge_mid_inside
|
||
n_cent = int(inside_inner_centroid.sum())
|
||
n_strict = int(has_strict_vertex.sum())
|
||
n_emid = int(edge_mid_inside.sum())
|
||
tris_keep = tri.simplices[~inside_inner]
|
||
log_fn(f" [DEM] inner 제거 삼각형: centroid {n_cent}개 + strict-vertex {n_strict}개 "
|
||
f"+ edge-midpoint {n_emid}개 (tol={strict_tol:.3f}m) — TIN 침범 3중 차단")
|
||
|
||
# 링 삼각형 벽 컷 — **2단계 보호 + 추가 완화**.
|
||
# 이전: slope_ratio>2.5 & Z스팬>10m 도 실제 72° 산사면/급한 골짜기를 자르면서
|
||
# 접합점에 삼각 구멍(= error.png의 두 번째 이슈)을 냈다.
|
||
# 이번:
|
||
# (1) **TIN 공유 경계 정점을 가진 삼각형은 컷 제외** — bbox seam 삼각형은
|
||
# DEM 급경사 한 번에 날아가면 보이는 **접합점 구멍**이 된다. 보호.
|
||
# (2) 기준 상향: slope_ratio > 4.0 (≈76°) AND Z스팬 > 30m — 자연 급사면보다
|
||
# 가파르고, 30m 이상 Z 점프가 있어야만 진짜 "수직 벽" 후보로 간주.
|
||
# (3) 추가 가드: 에지 길이 e_max가 너무 짧지 않아야 함(<5m는 제외: 극미세
|
||
# sliver는 별도 수치 오차, 벽 아님).
|
||
if len(tris_keep) > 0:
|
||
zs = ring_xyz[tris_keep][:, :, 2]
|
||
z_span_tri = zs.max(axis=1) - zs.min(axis=1)
|
||
pxy = ring_xyz[tris_keep][:, :, :2]
|
||
e01 = np.linalg.norm(pxy[:, 0] - pxy[:, 1], axis=1)
|
||
e12 = np.linalg.norm(pxy[:, 1] - pxy[:, 2], axis=1)
|
||
e20 = np.linalg.norm(pxy[:, 2] - pxy[:, 0], axis=1)
|
||
e_max_tri = np.maximum(np.maximum(e01, e12), e20)
|
||
e_max_safe = np.maximum(e_max_tri, 1e-6)
|
||
slope_ratio_tri = z_span_tri / e_max_safe
|
||
|
||
# (1) seam 보호 마스크 — TIN bbox 변 공유 정점(= boundary_pts, 마지막
|
||
# n_hull_boundary 개)을 가진 삼각형은 "접합부"라 컷 제외.
|
||
protect_vertex = np.zeros(len(ring_xyz), dtype=bool)
|
||
if n_hull_boundary > 0:
|
||
protect_vertex[len(ring_xyz) - n_hull_boundary:] = True
|
||
pv0 = protect_vertex[tris_keep[:, 0]]
|
||
pv1 = protect_vertex[tris_keep[:, 1]]
|
||
pv2 = protect_vertex[tris_keep[:, 2]]
|
||
seam_tri = pv0 | pv1 | pv2
|
||
|
||
# (2+3) 완화된 기준
|
||
wall_mask = (slope_ratio_tri > 4.0) & (z_span_tri > 30.0) & (e_max_tri > 5.0) & (~seam_tri)
|
||
if wall_mask.any():
|
||
tris_keep = tris_keep[~wall_mask]
|
||
log_fn(f" [DEM] 링 벽 삼각형 {int(wall_mask.sum())}개 제거 "
|
||
f"(slope_ratio>4.0(≈76°) & Z스팬>30m, seam {int(seam_tri.sum())}개 보호) "
|
||
f"— 자연 급사면 보존, 접합점 구멍 방지")
|
||
else:
|
||
log_fn(f" [DEM] 링 벽 컷: 대상 없음 (seam {int(seam_tri.sum())}개 보호) "
|
||
f"— 자연 경사면 유지")
|
||
|
||
if len(tris_keep) == 0:
|
||
raise RuntimeError("링 삼각형이 하나도 남지 않았습니다 (buffer_m/grid_step_m 확인).")
|
||
|
||
faces = np.column_stack([np.full(len(tris_keep), 3), tris_keep])
|
||
mesh = pv.PolyData(ring_xyz, faces)
|
||
mesh["Elevation"] = ring_xyz[:, 2]
|
||
|
||
info = (f"source={resolved_source}, buffer={buffer_m:.0f}m, step={grid_step_m:.1f}m, "
|
||
f"feather={feather_m:.0f}m, offset={vertical_offset:+.2f}m, "
|
||
f"boundary={boundary_mode}, pts={len(ring_xyz)}, tris={len(tris_keep)}")
|
||
log_fn(f" [DEM] 링 메시 완료: {info}")
|
||
|
||
return DemExtendResult(
|
||
mesh=mesh,
|
||
n_points=len(ring_xyz),
|
||
n_faces=len(tris_keep),
|
||
buffer_m=buffer_m,
|
||
source=resolved_source,
|
||
vertical_offset_m=vertical_offset,
|
||
zoom=zoom,
|
||
grid_step_m=grid_step_m,
|
||
feather_m=feather_m,
|
||
info=info,
|
||
)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
# 스모크 테스트: 사연댐 대략 좌표(WGS84) 주변
|
||
# 투영 좌표로 변환해서 테스트
|
||
to_proj = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:5187", always_xy=True)
|
||
cx, cy = to_proj.transform(129.10, 35.54) # 대략 사연댐 부근
|
||
half = 1000.0
|
||
bounds = (cx - half, cy - half, cx + half, cy + half)
|
||
origin = np.array([cx - half, cy - half, 100.0])
|
||
result = build_extended_terrain_ring(
|
||
projected_bounds=bounds,
|
||
origin=origin,
|
||
src_crs="EPSG:5187",
|
||
buffer_m=500.0,
|
||
grid_step_m=50.0,
|
||
feather_m=0.0,
|
||
log_fn=print,
|
||
)
|
||
print(result.info)
|