eomsangdonandClaude Opus 5
72075563e3
- 지형 파일 개수 제한 해제(정확히 1개 → 1장 이상), 한 카드에 여러 장 담기 (화면 표시 「용화_서편.las 외 1장」, 업로드는 한 장씩 차례로 전송) - 구조화 엔진이 여러 파일을 합친 범위로 한 벌 생성 — WF1 자동 전처리·B04 재분석 모두 프로젝트 지형 파일 전부를 대상으로 실행 - 점이 1억 개를 넘으면 씨닝 — 지면 분류점은 전부 남기고 나머지만 0.5m 칸 최저점으로 축소 (용화 실측: 4,900만점 → 249만점, 지면점 1,140,716개 그대로, 1m 지면격자 표고차 0.0000m) - 머리글만 읽어 5km 넘게 떨어진 파일은 업로드 거부 (임도는 길어도 2~3km) - 화면 조립부 700줄 준수를 위해 terrainCoverage 를 판정 모듈로 이동 Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
339 lines
14 KiB
Python
339 lines
14 KiB
Python
"""B04 LAS/LAZ 고속 구조화 엔진 — 여러 장을 한 벌로 병합한다 (2026-09-06 사용자 확정).
|
|
|
|
드론 라이다는 사업지가 넓으면 도엽별로 여러 장이 온다. 여기서 **합친 범위**로 한 벌을
|
|
만들고, 뒤 단계(지면필터·모델·등고선·배수)는 받는 형식이 그대로라 손대지 않는다.
|
|
|
|
점이 임계를 넘으면 **칸(기본 0.5m)마다 최저점 하나만** 남긴다(씨닝). 설계가 쓰는 격자가
|
|
1m(지면필터 2m·CSF 천 1.5m)라 0.5m 는 설계보다 촘촘해 결과 표고가 사실상 같고, 30GB 두
|
|
장이 메모리 29GB → 1.4GB 로 내려간다. 임계 아래면 원본 점을 그대로 쓴다 — 작은 자료의
|
|
결과는 바뀌지 않는다.
|
|
"""
|
|
|
|
import logging
|
|
import os
|
|
import tempfile
|
|
from collections.abc import Callable, Sequence
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import laspy
|
|
import numpy as np
|
|
|
|
from common_util.common_util_json import replace_with_retry
|
|
from config.config_system import (
|
|
SURFACE_DEFAULT_RGB_VALUE,
|
|
SURFACE_LAS_CHUNK_SIZE,
|
|
SURFACE_MERGE_MAX_GAP_M,
|
|
SURFACE_THIN_CELL_SIZE_M,
|
|
SURFACE_THIN_MAX_CELLS,
|
|
SURFACE_THIN_TRIGGER_POINTS,
|
|
)
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
ProgressCallback = Callable[[int], None]
|
|
PathLike = str | Path
|
|
# 청크에서 점과 함께 옮기는 속성들 — 파일에 없으면 기본값이 남는다.
|
|
_ATTRIBUTES = ("intensity", "rgb", "return_number", "number_of_returns", "classification")
|
|
|
|
|
|
def _as_list(las_path: PathLike | Sequence[PathLike]) -> list[Path]:
|
|
if isinstance(las_path, (str, Path)):
|
|
return [Path(las_path)]
|
|
return [Path(item) for item in las_path]
|
|
|
|
|
|
def point_cloud_extent(path: PathLike) -> tuple[int, tuple[float, float, float, float]]:
|
|
"""머리글만 읽어 점 수와 XY 범위를 돌려준다 — 파일 크기와 무관하게 즉시 끝난다."""
|
|
with laspy.open(Path(path)) as las_file:
|
|
header = las_file.header
|
|
return int(header.point_count), (
|
|
float(header.mins[0]),
|
|
float(header.mins[1]),
|
|
float(header.maxs[0]),
|
|
float(header.maxs[1]),
|
|
)
|
|
|
|
|
|
def merge_gap_error(
|
|
paths: Sequence[PathLike], gap_m: float = SURFACE_MERGE_MAX_GAP_M
|
|
) -> str | None:
|
|
"""서로 멀리 떨어진 지형 파일이 섞였는지 — 문제면 안내 문구, 없으면 None.
|
|
|
|
다른 사업지 파일이나 좌표계가 다른 파일이 섞이면 합친 범위가 통째로 어긋나 격자가
|
|
터진다. 도엽으로 나뉜 자료는 경계가 맞닿으므로 여유를 두고 **어느 파일과도 만나지
|
|
않는 파일**만 걸러 낸다.
|
|
"""
|
|
sources = _as_list(paths)
|
|
if len(sources) < 2:
|
|
return None
|
|
boxes = [(path, point_cloud_extent(path)[1]) for path in sources]
|
|
for index, (path, box) in enumerate(boxes):
|
|
near = any(
|
|
box[0] - gap_m <= other[2]
|
|
and other[0] - gap_m <= box[2]
|
|
and box[1] - gap_m <= other[3]
|
|
and other[1] - gap_m <= box[3]
|
|
for other_index, (_, other) in enumerate(boxes)
|
|
if other_index != index
|
|
)
|
|
if not near:
|
|
return (
|
|
f"지형 파일 「{path.name}」의 좌표가 다른 파일과 "
|
|
f"{int(gap_m):,}m 넘게 떨어져 있습니다."
|
|
" 같은 사업지의 파일인지, 좌표계가 같은지 확인해 주십시오."
|
|
)
|
|
return None
|
|
|
|
|
|
class _Merged:
|
|
"""합친 점을 담는 그릇 — 원본 유지형과 씨닝형이 같은 모양으로 낸다."""
|
|
|
|
def __init__(self, capacity: int) -> None:
|
|
self.xyz = np.empty((capacity, 3), dtype=np.float64)
|
|
self.intensity = np.zeros(capacity, dtype=np.uint16)
|
|
self.rgb = np.full((capacity, 3), SURFACE_DEFAULT_RGB_VALUE, dtype=np.uint8)
|
|
self.return_number = np.ones(capacity, dtype=np.uint8)
|
|
self.number_of_returns = np.ones(capacity, dtype=np.uint8)
|
|
self.classification = np.zeros(capacity, dtype=np.uint8)
|
|
self.size = 0
|
|
|
|
def arrays(self) -> dict[str, np.ndarray]:
|
|
end = self.size
|
|
return {
|
|
"xyz": self.xyz[:end],
|
|
"intensity": self.intensity[:end],
|
|
"rgb": self.rgb[:end],
|
|
"return_number": self.return_number[:end],
|
|
"number_of_returns": self.number_of_returns[:end],
|
|
"classification": self.classification[:end],
|
|
}
|
|
|
|
def append(self, columns: dict[str, np.ndarray]) -> None:
|
|
count = len(columns["x"])
|
|
section = slice(self.size, self.size + count)
|
|
self.xyz[section, 0] = columns["x"]
|
|
self.xyz[section, 1] = columns["y"]
|
|
self.xyz[section, 2] = columns["z"]
|
|
for key in _ATTRIBUTES:
|
|
if key in columns:
|
|
getattr(self, key)[section] = columns[key]
|
|
self.size += count
|
|
|
|
|
|
def _chunk_columns(chunk: Any, dimensions: set[str]) -> dict[str, np.ndarray]:
|
|
"""청크에서 쓸 값만 꺼낸다. 파일에 없는 항목은 키를 빼서 기본값이 남게 한다."""
|
|
columns: dict[str, np.ndarray] = {
|
|
"x": np.asarray(chunk.x, dtype=np.float64),
|
|
"y": np.asarray(chunk.y, dtype=np.float64),
|
|
"z": np.asarray(chunk.z, dtype=np.float64),
|
|
}
|
|
if "intensity" in dimensions:
|
|
columns["intensity"] = np.asarray(chunk.intensity, dtype=np.uint16)
|
|
if {"red", "green", "blue"}.issubset(dimensions):
|
|
colors = np.stack(
|
|
[
|
|
np.asarray(chunk.red, dtype=np.float64),
|
|
np.asarray(chunk.green, dtype=np.float64),
|
|
np.asarray(chunk.blue, dtype=np.float64),
|
|
],
|
|
axis=1,
|
|
)
|
|
if colors.size and float(colors.max()) > 255.0:
|
|
colors /= 256.0
|
|
columns["rgb"] = colors.clip(0, 255).astype(np.uint8)
|
|
if {"return_number", "number_of_returns"}.issubset(dimensions):
|
|
columns["return_number"] = np.asarray(chunk.return_number, dtype=np.uint8)
|
|
columns["number_of_returns"] = np.asarray(chunk.number_of_returns, dtype=np.uint8)
|
|
if "classification" in dimensions:
|
|
columns["classification"] = np.asarray(chunk.classification, dtype=np.uint8)
|
|
return columns
|
|
|
|
|
|
class _ThinGrid:
|
|
"""씨닝형 — **지면 분류점은 전부** 남기고, 나머지는 칸마다 최저점 하나만 남긴다.
|
|
|
|
설계 지표면을 만드는 것은 지면점이다(업체가 분류해 준 ASPRS class 2). 그 점을 하나도
|
|
버리지 않으므로 **지면 결과는 씨닝 전과 완전히 같다**(2026-09-06 용화 실측: 1m 지면
|
|
격자 141,969칸 전부 표고 차이 0). 지면점은 원본의 2~3%뿐이라 남겨도 가볍다.
|
|
|
|
나머지(수목·구조물·잡음)는 칸마다 최저점만 남긴다 — 분류가 없는 자료에서 CSF·PMF가
|
|
지면을 찾을 밑그림으로 충분하다(CSF 천 간격 1.5m > 칸 0.5m).
|
|
|
|
한 청크 안에서 같은 칸이 여러 번 나오면 뒤에 쓴 값이 이겨 최저점이 아니게 된다.
|
|
그래서 청크를 (칸, 표고)로 정렬해 **칸마다 첫 점**만 골라 낸 뒤 격자와 견준다.
|
|
"""
|
|
|
|
#: ASPRS 지면 분류 코드.
|
|
GROUND_CLASS = 2
|
|
|
|
def __init__(self, bounds: np.ndarray, cell_size: float) -> None:
|
|
self.cell_size = cell_size
|
|
self.x_min = float(bounds[0, 0])
|
|
self.y_min = float(bounds[1, 0])
|
|
self.width = int(np.ceil((float(bounds[0, 1]) - self.x_min) / cell_size)) + 1
|
|
self.height = int(np.ceil((float(bounds[1, 1]) - self.y_min) / cell_size)) + 1
|
|
cells = self.width * self.height
|
|
if cells > SURFACE_THIN_MAX_CELLS:
|
|
raise ValueError(
|
|
"지형 자료의 합친 범위가 너무 넓습니다."
|
|
" 같은 사업지의 파일인지, 좌표계가 같은지 확인해 주십시오."
|
|
)
|
|
self.best_z = np.full(cells, np.inf, dtype=np.float64)
|
|
self.x = np.zeros(cells, dtype=np.float64)
|
|
self.y = np.zeros(cells, dtype=np.float64)
|
|
self.intensity = np.zeros(cells, dtype=np.uint16)
|
|
self.rgb = np.full((cells, 3), SURFACE_DEFAULT_RGB_VALUE, dtype=np.uint8)
|
|
self.return_number = np.ones(cells, dtype=np.uint8)
|
|
self.number_of_returns = np.ones(cells, dtype=np.uint8)
|
|
self.classification = np.zeros(cells, dtype=np.uint8)
|
|
# 그대로 남길 지면점 — 청크마다 모아 두었다가 마지막에 잇는다.
|
|
self.ground: list[dict[str, np.ndarray]] = []
|
|
|
|
def add(self, columns: dict[str, np.ndarray]) -> None:
|
|
classification = columns.get("classification")
|
|
if classification is not None:
|
|
is_ground = classification == self.GROUND_CLASS
|
|
if is_ground.any():
|
|
self.ground.append({key: value[is_ground] for key, value in columns.items()})
|
|
keep = ~is_ground
|
|
columns = {key: value[keep] for key, value in columns.items()}
|
|
x, y, z = columns["x"], columns["y"], columns["z"]
|
|
if not len(x):
|
|
return
|
|
grid_x = np.clip(((x - self.x_min) / self.cell_size).astype(np.int64), 0, self.width - 1)
|
|
grid_y = np.clip(((y - self.y_min) / self.cell_size).astype(np.int64), 0, self.height - 1)
|
|
cell = grid_y * self.width + grid_x
|
|
order = np.lexsort((z, cell))
|
|
sorted_cell = cell[order]
|
|
first = np.ones(len(order), dtype=bool)
|
|
first[1:] = sorted_cell[1:] != sorted_cell[:-1]
|
|
candidate = order[first]
|
|
candidate_cell = cell[candidate]
|
|
better = z[candidate] < self.best_z[candidate_cell]
|
|
chosen = candidate[better]
|
|
target = candidate_cell[better]
|
|
self.best_z[target] = z[chosen]
|
|
self.x[target] = x[chosen]
|
|
self.y[target] = y[chosen]
|
|
for key in _ATTRIBUTES:
|
|
if key in columns:
|
|
getattr(self, key)[target] = columns[key][chosen]
|
|
|
|
def collect(self) -> _Merged:
|
|
occupied = np.flatnonzero(np.isfinite(self.best_z))
|
|
ground_count = sum(len(item["x"]) for item in self.ground)
|
|
merged = _Merged(len(occupied) + ground_count)
|
|
merged.xyz[: len(occupied), 0] = self.x[occupied]
|
|
merged.xyz[: len(occupied), 1] = self.y[occupied]
|
|
merged.xyz[: len(occupied), 2] = self.best_z[occupied]
|
|
for key in _ATTRIBUTES:
|
|
getattr(merged, key)[: len(occupied)] = getattr(self, key)[occupied]
|
|
merged.size = len(occupied)
|
|
for item in self.ground:
|
|
merged.append(item)
|
|
return merged
|
|
|
|
|
|
def _headers(sources: list[Path]) -> tuple[int, np.ndarray, bool]:
|
|
"""전체 점 수·합친 범위(3x2)·색 보유 여부를 머리글만 읽어 구한다."""
|
|
total = 0
|
|
has_rgb = False
|
|
mins = np.full(3, np.inf, dtype=np.float64)
|
|
maxs = np.full(3, -np.inf, dtype=np.float64)
|
|
for source in sources:
|
|
with laspy.open(source) as las_file:
|
|
header = las_file.header
|
|
total += int(header.point_count)
|
|
mins = np.minimum(mins, np.asarray(header.mins, dtype=np.float64))
|
|
maxs = np.maximum(maxs, np.asarray(header.maxs, dtype=np.float64))
|
|
dimensions = set(header.point_format.dimension_names)
|
|
has_rgb = has_rgb or {"red", "green", "blue"}.issubset(dimensions)
|
|
if not np.isfinite(mins).all():
|
|
mins = np.zeros(3, dtype=np.float64)
|
|
maxs = np.zeros(3, dtype=np.float64)
|
|
return total, np.column_stack((mins, maxs)), has_rgb
|
|
|
|
|
|
def _merge_sources(
|
|
sources: list[Path],
|
|
total_points: int,
|
|
bounds: np.ndarray,
|
|
thin: bool,
|
|
progress_callback: ProgressCallback | None,
|
|
) -> _Merged:
|
|
grid = _ThinGrid(bounds, SURFACE_THIN_CELL_SIZE_M) if thin else None
|
|
merged = _Merged(total_points) if grid is None else None
|
|
done = 0
|
|
for source in sources:
|
|
with laspy.open(source) as las_file:
|
|
dimensions = set(las_file.header.point_format.dimension_names)
|
|
for chunk in las_file.chunk_iterator(SURFACE_LAS_CHUNK_SIZE):
|
|
columns = _chunk_columns(chunk, dimensions)
|
|
if grid is not None:
|
|
grid.add(columns)
|
|
else:
|
|
merged.append(columns)
|
|
done += len(columns["x"])
|
|
if progress_callback:
|
|
progress_callback(int(done / total_points * 100) if total_points else 100)
|
|
return grid.collect() if grid is not None else merged
|
|
|
|
|
|
def structurize_las(
|
|
las_path: PathLike | Sequence[PathLike],
|
|
output_dir: str | Path,
|
|
progress_callback: ProgressCallback | None = None,
|
|
) -> Path:
|
|
"""지형 파일 한 장 또는 여러 장을 청크로 읽어 B04 structured.npz로 원자적 저장한다."""
|
|
sources = _as_list(las_path)
|
|
if not sources:
|
|
raise ValueError("구조화할 지형 파일이 없습니다.")
|
|
target_dir = Path(output_dir)
|
|
target_dir.mkdir(parents=True, exist_ok=True)
|
|
target = target_dir / "structured.npz"
|
|
|
|
total_points, bounds, has_rgb = _headers(sources)
|
|
thin = total_points > SURFACE_THIN_TRIGGER_POINTS
|
|
merged = _merge_sources(sources, total_points, bounds, thin, progress_callback)
|
|
logger.info(
|
|
"B04 구조화: 파일 %d장 원본 %d점 → 저장 %d점 (씨닝 %s)",
|
|
len(sources),
|
|
total_points,
|
|
merged.size,
|
|
f"{SURFACE_THIN_CELL_SIZE_M}m 칸" if thin else "없음",
|
|
)
|
|
|
|
temporary_path: Path | None = None
|
|
try:
|
|
with tempfile.NamedTemporaryFile(
|
|
mode="wb",
|
|
dir=target_dir,
|
|
prefix=".structured.",
|
|
suffix=".npz.tmp",
|
|
delete=False,
|
|
) as temporary:
|
|
temporary_path = Path(temporary.name)
|
|
np.savez_compressed(
|
|
temporary,
|
|
**merged.arrays(),
|
|
bounds=bounds,
|
|
total_points=np.array([merged.size], dtype=np.int64),
|
|
source_point_count=np.array([total_points], dtype=np.int64),
|
|
source_file_count=np.array([len(sources)], dtype=np.int64),
|
|
thinned=np.array([int(thin)], dtype=np.int8),
|
|
has_rgb=np.array([int(has_rgb)], dtype=np.int8),
|
|
)
|
|
temporary.flush()
|
|
os.fsync(temporary.fileno())
|
|
# 이 교체가 윈도우 공유 위반으로 죽으면 WF1 분석 전체가 실패한다 — 재시도한다.
|
|
replace_with_retry(temporary_path, target)
|
|
temporary_path = None
|
|
finally:
|
|
if temporary_path is not None:
|
|
temporary_path.unlink(missing_ok=True)
|
|
|
|
if progress_callback:
|
|
progress_callback(100)
|
|
return target
|