245 lines
9.9 KiB
Python
245 lines
9.9 KiB
Python
from __future__ import annotations
|
|
|
|
import math
|
|
from dataclasses import dataclass
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
|
|
from utils.surface_elevation_sampler import SurfaceElevationSampler
|
|
|
|
|
|
SECTION_SCHEMA_VERSION = 1
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class SectionGenerationOptions:
|
|
station_interval_m: float = 20.0
|
|
cross_half_width_m: float = 15.0
|
|
cross_sample_interval_m: float = 0.5
|
|
long_sample_interval_m: float = 1.0
|
|
include_endpoint: bool = True
|
|
|
|
def validate(self) -> None:
|
|
values = {
|
|
"횡단 측점 간격": self.station_interval_m,
|
|
"횡단 좌우 폭": self.cross_half_width_m,
|
|
"횡단 샘플 간격": self.cross_sample_interval_m,
|
|
"종단 샘플 간격": self.long_sample_interval_m,
|
|
}
|
|
for label, value in values.items():
|
|
if not math.isfinite(value) or value <= 0:
|
|
raise ValueError(f"{label}은 0보다 큰 유한한 값이어야 합니다.")
|
|
if self.cross_sample_interval_m > self.cross_half_width_m * 2:
|
|
raise ValueError("횡단 샘플 간격이 전체 횡단 폭보다 클 수 없습니다.")
|
|
|
|
|
|
def format_station(chainage_m: float) -> str:
|
|
"""WebCAD 도면에서도 재사용할 수 있는 STA.k+mmm.mmm 표기."""
|
|
chainage_m = max(float(chainage_m), 0.0)
|
|
km = int(chainage_m // 1000.0)
|
|
remainder = chainage_m - km * 1000.0
|
|
return f"STA.{km}+{remainder:07.3f}"
|
|
|
|
|
|
def _clean_polyline(polyline: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
|
|
points = np.asarray(polyline, dtype=np.float64)
|
|
if points.ndim != 2 or points.shape[1] < 2 or len(points) < 2:
|
|
raise ValueError("경로는 최소 2개의 (x, y, z) 좌표로 구성되어야 합니다.")
|
|
if points.shape[1] == 2:
|
|
points = np.column_stack([points, np.full(len(points), np.nan)])
|
|
else:
|
|
points = points[:, :3]
|
|
if not np.all(np.isfinite(points[:, :2])):
|
|
raise ValueError("경로 XY 좌표에 NaN 또는 Infinity가 있습니다.")
|
|
|
|
distances = np.hypot(np.diff(points[:, 0]), np.diff(points[:, 1]))
|
|
keep = np.r_[True, distances > 1e-8]
|
|
points = points[keep]
|
|
if len(points) < 2:
|
|
raise ValueError("수평 길이가 있는 경로 구간이 없습니다.")
|
|
distances = np.hypot(np.diff(points[:, 0]), np.diff(points[:, 1]))
|
|
chainage = np.r_[0.0, np.cumsum(distances)]
|
|
return points, chainage
|
|
|
|
|
|
def _chainages(total: float, interval: float, include_endpoint: bool) -> np.ndarray:
|
|
values = np.arange(0.0, total + 1e-9, interval, dtype=np.float64)
|
|
if not len(values) or abs(values[0]) > 1e-9:
|
|
values = np.r_[0.0, values]
|
|
if include_endpoint and total - values[-1] > 1e-6:
|
|
values = np.r_[values, total]
|
|
return values
|
|
|
|
|
|
def _interpolate_xy(points: np.ndarray, chainage: np.ndarray, targets: np.ndarray) -> np.ndarray:
|
|
return np.column_stack([
|
|
np.interp(targets, chainage, points[:, 0]),
|
|
np.interp(targets, chainage, points[:, 1]),
|
|
])
|
|
|
|
|
|
def _tangent_at(points: np.ndarray, chainage: np.ndarray, target: float, probe: float) -> np.ndarray:
|
|
total = float(chainage[-1])
|
|
before = max(0.0, target - probe)
|
|
after = min(total, target + probe)
|
|
if after - before <= 1e-9:
|
|
before = max(0.0, target - 1e-3)
|
|
after = min(total, target + 1e-3)
|
|
pair = _interpolate_xy(points, chainage, np.array([before, after]))
|
|
vector = pair[1] - pair[0]
|
|
length = float(np.hypot(vector[0], vector[1]))
|
|
if length <= 1e-9:
|
|
raise ValueError(f"측점 {target:.3f}m에서 경로 접선 방향을 계산할 수 없습니다.")
|
|
return vector / length
|
|
|
|
|
|
def _float_or_none(value: float) -> float | None:
|
|
return round(float(value), 6) if math.isfinite(float(value)) else None
|
|
|
|
|
|
def generate_sections(
|
|
polyline: np.ndarray | list[list[float]],
|
|
sampler: SurfaceElevationSampler,
|
|
options: SectionGenerationOptions | None = None,
|
|
*,
|
|
source_snapshot: dict[str, Any] | None = None,
|
|
crs: str | None = None,
|
|
) -> dict[str, Any]:
|
|
"""확정 경로로 CAD 인계 가능한 종단·횡단 원시 데이터를 생성한다.
|
|
|
|
CP는 별도 측점을 만들지 않는다. BP(0m)부터 ``station_interval_m`` 간격으로
|
|
진행하고, 설정된 경우 마지막 EP를 추가한다.
|
|
"""
|
|
options = options or SectionGenerationOptions()
|
|
options.validate()
|
|
points, route_chainage = _clean_polyline(np.asarray(polyline, dtype=np.float64))
|
|
total = float(route_chainage[-1])
|
|
|
|
long_chainage = _chainages(total, options.long_sample_interval_m, True)
|
|
long_xy = _interpolate_xy(points, route_chainage, long_chainage)
|
|
long_z, long_valid = sampler.sample_xy(long_xy)
|
|
|
|
station_chainage = _chainages(total, options.station_interval_m, options.include_endpoint)
|
|
station_xy = _interpolate_xy(points, route_chainage, station_chainage)
|
|
tangents = np.vstack([
|
|
_tangent_at(points, route_chainage, float(value), max(options.long_sample_interval_m, 0.5))
|
|
for value in station_chainage
|
|
])
|
|
left_axes = np.column_stack([-tangents[:, 1], tangents[:, 0]])
|
|
|
|
offsets = np.arange(
|
|
-options.cross_half_width_m,
|
|
options.cross_half_width_m + options.cross_sample_interval_m * 0.5,
|
|
options.cross_sample_interval_m,
|
|
dtype=np.float64,
|
|
)
|
|
offsets = offsets[offsets <= options.cross_half_width_m + 1e-9]
|
|
if not np.any(np.isclose(offsets, 0.0, atol=1e-9)):
|
|
offsets = np.sort(np.r_[offsets, 0.0])
|
|
|
|
all_cross_xy = (
|
|
station_xy[:, None, :] + left_axes[:, None, :] * offsets[None, :, None]
|
|
).reshape(-1, 2)
|
|
all_cross_z, all_cross_valid = sampler.sample_xy(all_cross_xy)
|
|
all_cross_z = all_cross_z.reshape(len(station_chainage), len(offsets))
|
|
all_cross_valid = all_cross_valid.reshape(len(station_chainage), len(offsets))
|
|
|
|
stations: list[dict[str, Any]] = []
|
|
cross_sections: list[dict[str, Any]] = []
|
|
for index, value in enumerate(station_chainage):
|
|
station_id = f"station_{int(round(float(value) * 1000)):012d}"
|
|
kind = "bp" if index == 0 else "ep" if abs(float(value) - total) <= 1e-6 else "regular"
|
|
center_index = int(np.argmin(np.abs(offsets)))
|
|
center_z = all_cross_z[index, center_index]
|
|
tangent = tangents[index]
|
|
left = left_axes[index]
|
|
azimuth = (math.degrees(math.atan2(tangent[0], tangent[1])) + 360.0) % 360.0
|
|
frame = {
|
|
"origin": {
|
|
"x": round(float(station_xy[index, 0]), 6),
|
|
"y": round(float(station_xy[index, 1]), 6),
|
|
"z": _float_or_none(center_z),
|
|
},
|
|
"tangent_xy": [round(float(tangent[0]), 9), round(float(tangent[1]), 9)],
|
|
"left_xy": [round(float(left[0]), 9), round(float(left[1]), 9)],
|
|
"up_xyz": [0.0, 0.0, 1.0],
|
|
}
|
|
station = {
|
|
"station_id": station_id,
|
|
"chainage_m": round(float(value), 6),
|
|
"label": format_station(float(value)),
|
|
"kind": kind,
|
|
"center_x": round(float(station_xy[index, 0]), 6),
|
|
"center_y": round(float(station_xy[index, 1]), 6),
|
|
"center_z": _float_or_none(center_z),
|
|
"azimuth_deg": round(azimuth, 6),
|
|
"frame": frame,
|
|
}
|
|
stations.append(station)
|
|
|
|
samples = []
|
|
for offset_index, offset in enumerate(offsets):
|
|
valid = bool(all_cross_valid[index, offset_index])
|
|
xy = all_cross_xy.reshape(len(station_chainage), len(offsets), 2)[index, offset_index]
|
|
samples.append({
|
|
"offset_m": round(float(offset), 6),
|
|
"x": round(float(xy[0]), 6),
|
|
"y": round(float(xy[1]), 6),
|
|
"z": _float_or_none(all_cross_z[index, offset_index]) if valid else None,
|
|
"elevation_m": _float_or_none(all_cross_z[index, offset_index]) if valid else None,
|
|
"valid": valid,
|
|
})
|
|
cross_sections.append({**station, "samples": samples})
|
|
|
|
longitudinal_samples = [
|
|
{
|
|
"chainage_m": round(float(chainage), 6),
|
|
"x": round(float(xy[0]), 6),
|
|
"y": round(float(xy[1]), 6),
|
|
"z": _float_or_none(z) if valid else None,
|
|
"elevation_m": _float_or_none(z) if valid else None,
|
|
"valid": bool(valid),
|
|
}
|
|
for chainage, xy, z, valid in zip(long_chainage, long_xy, long_z, long_valid)
|
|
]
|
|
|
|
finite_z = np.asarray([sample["z"] for sample in longitudinal_samples if sample["z"] is not None])
|
|
datum = math.floor(float(finite_z.min()) / 10.0) * 10.0 if finite_z.size else None
|
|
return {
|
|
"schema_version": SECTION_SCHEMA_VERSION,
|
|
"status": "completed",
|
|
"source": source_snapshot or {},
|
|
"coordinate_reference": {
|
|
"crs": crs,
|
|
"world_axes": {"x": "project_easting", "y": "project_northing", "z": "elevation"},
|
|
"units": {"horizontal": "m", "vertical": "m", "angle": "degree"},
|
|
},
|
|
"cad_exchange": {
|
|
"station_origin": "BP",
|
|
"chainage_direction": "BP_to_EP",
|
|
"cross_offset_sign": {"negative": "right", "positive": "left"},
|
|
"cross_local_axes": {"x": "offset_m", "y": "elevation_m"},
|
|
"recommended_drawing_datum_m": datum,
|
|
},
|
|
"options": {
|
|
"station_interval_m": options.station_interval_m,
|
|
"cross_half_width_m": options.cross_half_width_m,
|
|
"cross_sample_interval_m": options.cross_sample_interval_m,
|
|
"long_sample_interval_m": options.long_sample_interval_m,
|
|
"include_endpoint": options.include_endpoint,
|
|
},
|
|
"longitudinal": {
|
|
"length_m": round(total, 6),
|
|
"samples": longitudinal_samples,
|
|
"stations": stations,
|
|
},
|
|
"cross_sections": cross_sections,
|
|
"summary": {
|
|
"station_count": len(stations),
|
|
"cross_sample_count": int(len(stations) * len(offsets)),
|
|
"invalid_longitudinal_samples": int((~long_valid).sum()),
|
|
"invalid_cross_samples": int((~all_cross_valid).sum()),
|
|
},
|
|
}
|