74 lines
3.5 KiB
Python
74 lines
3.5 KiB
Python
"""G2 — 지표면 TIN과 점 조회."""
|
|
from __future__ import annotations
|
|
|
|
import numpy as np
|
|
from scipy.spatial import Delaunay
|
|
|
|
|
|
class DTM:
|
|
"""Delaunay TIN 지표면. 너무 긴 변의 삼각형은 '자료 없음'으로 본다."""
|
|
|
|
def __init__(self, points: np.ndarray, max_edge: float = 50.0):
|
|
pts = np.asarray(points, dtype=float)
|
|
if pts.ndim != 2 or pts.shape[1] != 3 or len(pts) < 3:
|
|
raise ValueError(
|
|
"DTM에는 지형점이 3개 이상 필요합니다. 종평면도의 표고점·등고선 레이어 또는 GL/DXF 좌표 형식을 확인하세요."
|
|
)
|
|
self.origin = pts[:, :2].mean(axis=0) # 큰 좌표의 수치 오차를 줄이기 위해 원점 이동
|
|
self.xy = pts[:, :2] - self.origin
|
|
self.z = pts[:, 2]
|
|
self.tri = Delaunay(self.xy)
|
|
s = self.tri.simplices
|
|
a, b, c = self.xy[s[:, 0]], self.xy[s[:, 1]], self.xy[s[:, 2]]
|
|
edge = np.maximum.reduce([np.hypot(*(a - b).T), np.hypot(*(b - c).T), np.hypot(*(c - a).T)])
|
|
self.valid = edge <= max_edge
|
|
self.max_edge = max_edge
|
|
|
|
def z_info(self, X: float, Y: float) -> tuple[float | None, float | None]:
|
|
"""(지반고, 그 값을 만든 세 점 중 가장 먼 점까지의 거리). 자료 범위 밖이면 (None, None)."""
|
|
p = np.array([[X, Y]]) - self.origin
|
|
k = int(self.tri.find_simplex(p)[0])
|
|
if k < 0 or not self.valid[k]:
|
|
return None, None
|
|
T = self.tri.transform[k]
|
|
bary2 = T[:2].dot(p[0] - T[2])
|
|
w = np.array([bary2[0], bary2[1], 1 - bary2.sum()])
|
|
s = self.tri.simplices[k]
|
|
far = float(np.hypot(*(self.xy[s] - p[0]).T).max())
|
|
return float(w.dot(self.z[s])), far
|
|
|
|
def z_at(self, X: float, Y: float) -> float | None:
|
|
"""측량 좌표 (X, Y)의 지반고. 자료 범위 밖이면 None."""
|
|
return self.z_info(X, Y)[0]
|
|
|
|
def z_many(self, XY: np.ndarray) -> np.ndarray:
|
|
"""여러 점의 지반고 (자료 없음은 NaN). z_info 와 같은 식을 배열로 계산한다."""
|
|
p = np.asarray(XY, dtype=float).reshape(-1, 2) - self.origin
|
|
k = self.tri.find_simplex(p)
|
|
ok = k >= 0
|
|
ok &= self.valid[np.maximum(k, 0)]
|
|
out = np.full(len(p), np.nan)
|
|
kk = k[ok]
|
|
T = self.tri.transform[kk] # (m, 3, 2)
|
|
b2 = np.einsum("mij,mj->mi", T[:, :2], p[ok] - T[:, 2]) # 첫 두 무게중심 좌표
|
|
s = self.tri.simplices[kk]
|
|
out[ok] = b2[:, 0] * self.z[s[:, 0]] + b2[:, 1] * self.z[s[:, 1]] + (1 - b2.sum(axis=1)) * self.z[s[:, 2]]
|
|
return out
|
|
|
|
def profile(self, start: tuple[float, float], end: tuple[float, float], step: float = 0.5) -> np.ndarray:
|
|
"""start → end 선을 따라 (거리, Z) 배열."""
|
|
(x0, y0), (x1, y1) = start, end
|
|
L = float(np.hypot(x1 - x0, y1 - y0))
|
|
n = max(2, int(L / step) + 1)
|
|
t = np.linspace(0, 1, n)
|
|
XY = np.column_stack([x0 + (x1 - x0) * t, y0 + (y1 - y0) * t])
|
|
return np.column_stack([t * L, self.z_many(XY)])
|
|
|
|
def grid(self, n: int = 200) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
|
"""그림용 격자 (X, Y, Z)."""
|
|
lo = self.xy.min(axis=0) + self.origin
|
|
hi = self.xy.max(axis=0) + self.origin
|
|
gx, gy = np.meshgrid(np.linspace(lo[0], hi[0], n), np.linspace(lo[1], hi[1], n))
|
|
gz = self.z_many(np.column_stack([gx.ravel(), gy.ravel()])).reshape(gx.shape)
|
|
return gx, gy, gz
|