feat(B04): 도엽 서피스를 구조선(계곡·능선) 기반 보간으로 개선한다
5m 등고선만 보고 보간하던 것을 구조선 기준으로 보강(2026-08-30 사용자 확정): - 계곡: 도엽_하천중심선을 기준선으로, 등고선 교차점을 Z 앵커 삼아 선 길이 비례로 보간한 정점열을 TIN에 추가 — 계곡 바닥이 등고선 사이에서도 연속으로 내려간다 (실데이터 605정점, 480~599m) - 능선: 마루 캡을 고정 +2.5m에서 마지막 등고선 바깥 사면의 국소 경사 연장으로 교체(상한 +간격−0.5m, 경사 불명 시 돔 폴백) — 마루가 주변 사면 기울기를 따라 이어진다 실측: 평탄 셀 5.14%→3.96%(최초 9.35%), 노선 |Δz| 불변(2.57m), 생성 16.4s. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -42,16 +42,19 @@ logger = logging.getLogger(__name__)
|
||||
|
||||
# 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값)
|
||||
_CONTOUR_FILE = "도엽_등고선.geojson"
|
||||
_STREAM_FILE = "도엽_하천중심선.geojson"
|
||||
|
||||
# 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다.
|
||||
SHEET_SOURCE_FILTER = "sheet"
|
||||
|
||||
|
||||
def _load_contour_features_metric(processed_dir: Path, epsg: int) -> list[dict[str, Any]]:
|
||||
"""병합 도엽 등고선(WGS84)을 읽어 사업지 CRS(m)로 재투영한다."""
|
||||
path = processed_dir / _CONTOUR_FILE
|
||||
def _load_features_metric(
|
||||
processed_dir: Path, epsg: int, filename: str = _CONTOUR_FILE
|
||||
) -> list[dict[str, Any]]:
|
||||
"""병합 도엽 레이어(WGS84)를 읽어 사업지 CRS(m)로 재투영한다."""
|
||||
path = processed_dir / filename
|
||||
if not path.exists():
|
||||
logger.warning("도엽 서피스: 등고선 파일이 없습니다: %s", path)
|
||||
logger.warning("도엽 서피스: 도엽 레이어 파일이 없습니다: %s", path)
|
||||
return []
|
||||
try:
|
||||
with path.open("r", encoding="utf-8") as file:
|
||||
@@ -102,7 +105,77 @@ def _preview_mesh(
|
||||
return clip_and_compact_mesh(vertices, faces, pv.reshape(-1))
|
||||
|
||||
|
||||
def _densify_between_contours(spec: Any, features: list[dict[str, Any]], cloud: Any) -> Any:
|
||||
def _stream_breakline_vertices(
|
||||
spec: Any, burned: np.ndarray, stream_features: list[dict[str, Any]]
|
||||
) -> tuple[np.ndarray, np.ndarray] | None:
|
||||
"""계곡 기준선(하천중심선)을 따라 등고 교차점 사이를 보간한 정점열을 만든다.
|
||||
|
||||
구조선 기반 보간(2026-08-30 사용자 확정): 계곡 축을 먼저 세우고, 축이 등고선과
|
||||
만나는 점을 Z 앵커로 삼아 앵커 사이를 선 길이 비례로 보간한다. 계곡 바닥이
|
||||
등고선 사이에서도 연속으로 내려가는 가상 종단이 되어 골짜기 평탄·역경사가 준다.
|
||||
"""
|
||||
from shapely import segmentize
|
||||
from shapely.geometry import shape
|
||||
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import iter_linestrings
|
||||
|
||||
sample_step_m = 2.0
|
||||
vertex_stride = 3 # 2m 샘플 → 6m 간격 정점 (등고 정점 재샘플 밀도와 유사)
|
||||
collected_xy: list[np.ndarray] = []
|
||||
collected_z: list[np.ndarray] = []
|
||||
for feature in stream_features:
|
||||
geometry = feature.get("geometry")
|
||||
if not geometry:
|
||||
continue
|
||||
try:
|
||||
parsed = shape(geometry)
|
||||
except Exception: # noqa: BLE001
|
||||
continue
|
||||
for line in iter_linestrings(parsed):
|
||||
coords = np.asarray(segmentize(line, sample_step_m).coords, dtype=np.float64)
|
||||
if len(coords) < 3:
|
||||
continue
|
||||
xy = coords[:, :2]
|
||||
arc = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))])
|
||||
row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy())
|
||||
inside = (row >= 0) & (col >= 0)
|
||||
levels_at = np.full(len(xy), np.nan, dtype=np.float32)
|
||||
levels_at[inside] = burned[row[inside], col[inside]]
|
||||
hit = np.isfinite(levels_at)
|
||||
if hit.sum() < 2:
|
||||
continue
|
||||
# 등고선 통과 구간(연속 같은 표고 샘플 묶음) 하나 = 앵커 하나.
|
||||
anchor_arc: list[float] = []
|
||||
anchor_z: list[float] = []
|
||||
run_start: int | None = None
|
||||
for i in range(len(xy) + 1):
|
||||
is_hit = i < len(xy) and hit[i]
|
||||
if is_hit and run_start is not None and levels_at[i] != levels_at[run_start]:
|
||||
is_hit = False # 표고가 바뀌면 묶음을 끊는다
|
||||
if is_hit:
|
||||
if run_start is None:
|
||||
run_start = i
|
||||
elif run_start is not None:
|
||||
anchor_arc.append(float(arc[run_start : i if i <= len(xy) else len(xy)].mean()))
|
||||
anchor_z.append(float(levels_at[run_start]))
|
||||
run_start = i if i < len(xy) and hit[i] else None
|
||||
if len(anchor_arc) < 2:
|
||||
continue
|
||||
# 앵커 사이 구간만 채택 — 첫 앵커 이전·마지막 앵커 이후는 근거가 없다.
|
||||
span = (arc >= anchor_arc[0]) & (arc <= anchor_arc[-1]) & inside & ~hit
|
||||
take = np.flatnonzero(span)[::vertex_stride]
|
||||
if not len(take):
|
||||
continue
|
||||
collected_xy.append(xy[take])
|
||||
collected_z.append(np.interp(arc[take], anchor_arc, anchor_z))
|
||||
if not collected_xy:
|
||||
return None
|
||||
return np.vstack(collected_xy), np.concatenate(collected_z)
|
||||
|
||||
|
||||
def _densify_between_contours(
|
||||
spec: Any, burned: np.ndarray, present: list[float], cloud: Any
|
||||
) -> Any:
|
||||
"""인접 표고 등고선 쌍의 등거리 중간선을 중간 표고 정점으로 추가한다.
|
||||
|
||||
등고선 정점만으로 Delaunay TIN을 치면 같은 표고 정점 3개짜리 평탄 삼각형이
|
||||
@@ -113,11 +186,8 @@ def _densify_between_contours(spec: Any, features: list[dict[str, Any]], cloud:
|
||||
"""
|
||||
from scipy.ndimage import distance_transform_edt
|
||||
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ContourCloud
|
||||
|
||||
burned, levels = rasterize_contours(spec, features, None)
|
||||
present = [level for level in sorted(levels) if bool((burned == level).any())]
|
||||
if len(present) < 2:
|
||||
return cloud
|
||||
|
||||
@@ -169,8 +239,9 @@ def _cap_flat_summits(surface: np.ndarray, cell_m: float, interval_m: float) ->
|
||||
|
||||
다음 등고선이 없는 봉우리·능선 마루는 TIN이 마지막 등고 표고로 평탄하게 채운다
|
||||
('ㅜ'가 갑자기 'ㅡ'로 변하는 단 — 2026-08-30 사용자 지적). 실제 마루는 그 표고와
|
||||
+간격 사이이므로, 경계에서 가장 먼 골격선이 +간격/2가 되도록 거리 비례로 올린다
|
||||
(기대값 캡, 오차 상한 ±간격/2 = 원천 데이터 불확실성 이내).
|
||||
+간격 사이이므로, 마지막 등고선 바깥 사면의 국소 경사를 안쪽으로 연장해 올린다
|
||||
(상한 +간격−0.5m — 다음 등고선이 없다는 사실과 모순되지 않게). 바깥 경사를
|
||||
구할 수 없으면 +간격/2 돔으로 폴백한다(2026-08-30 사용자 확정 ②a).
|
||||
적용 조건: 등표고 평탄 컴포넌트이고 바깥 유효 이웃이 전부 더 낮을 때(=마루).
|
||||
골짜기 바닥은 배수 방향을 모르므로 건드리지 않는다. 돋운 컴포넌트 수를 반환.
|
||||
"""
|
||||
@@ -205,24 +276,46 @@ def _cap_flat_summits(surface: np.ndarray, cell_m: float, interval_m: float) ->
|
||||
or cols.max() == surface.shape[1] - 1
|
||||
):
|
||||
continue # 격자 가장자리에 닿으면 바깥을 모른다
|
||||
# 이후 연산은 컴포넌트 bbox(+1) 창에서만 — 전체 격자 반복을 피한다.
|
||||
# 이후 연산은 컴포넌트 bbox 창(+경사 밴드 여유)에서만 — 전체 격자 반복 회피.
|
||||
band_cells = 8
|
||||
window = (
|
||||
slice(rows.min() - 1, rows.max() + 2),
|
||||
slice(cols.min() - 1, cols.max() + 2),
|
||||
slice(
|
||||
max(rows.min() - band_cells, 0),
|
||||
min(rows.max() + band_cells + 1, surface.shape[0]),
|
||||
),
|
||||
slice(
|
||||
max(cols.min() - band_cells, 0),
|
||||
min(cols.max() + band_cells + 1, surface.shape[1]),
|
||||
),
|
||||
)
|
||||
component = value_labels[window] == component_id
|
||||
ring = binary_dilation(component) & ~component & finite[window]
|
||||
if not ring.any():
|
||||
continue
|
||||
level = float(surface[rows[0], cols[0]])
|
||||
if float(surface[window][ring].max()) >= level - 1e-3:
|
||||
patch = surface[window]
|
||||
if float(patch[ring].max()) >= level - 1e-3:
|
||||
continue # 더 높은 이웃이 있으면 마루가 아니다(사면 벤치·골짜기)
|
||||
inner = distance_transform_edt(component, sampling=cell_m)
|
||||
peak = float(inner.max())
|
||||
if peak <= 0.0:
|
||||
continue
|
||||
patch = surface[window]
|
||||
patch[component] += (interval_m / 2.0) * (inner[component] / peak).astype(surface.dtype)
|
||||
# 바깥 사면 경사 추정 — 마지막 등고선 밖 밴드의 (표고 낙차 / 거리) 평균을
|
||||
# 안쪽으로 연장한다(2026-08-30 사용자 확정 ②a). 상한 +간격−0.5m:
|
||||
# 다음 등고선이 없다는 사실(마루 < L+간격)과 모순되지 않게.
|
||||
outer_distance = distance_transform_edt(~component, sampling=cell_m)
|
||||
band = (~component) & finite[window] & (outer_distance <= band_cells * cell_m)
|
||||
band &= outer_distance > 0
|
||||
if band.any():
|
||||
drops = level - patch[band].astype(np.float64)
|
||||
slope = float(np.mean(drops / outer_distance[band]))
|
||||
else:
|
||||
slope = 0.0
|
||||
if slope > 1e-3:
|
||||
rise = np.minimum(slope * inner[component], interval_m - 0.5)
|
||||
else: # 경사를 못 구하면 +간격/2 돔 폴백
|
||||
rise = (interval_m / 2.0) * (inner[component] / peak)
|
||||
patch[component] += rise.astype(surface.dtype)
|
||||
capped += 1
|
||||
if capped:
|
||||
logger.info(
|
||||
@@ -244,7 +337,7 @@ def build_sheet_surface_model(
|
||||
분석을 계속한다(도엽 미확보 지역 폴백).
|
||||
"""
|
||||
started = time.monotonic()
|
||||
features = _load_contour_features_metric(processed_dir, epsg)
|
||||
features = _load_features_metric(processed_dir, epsg)
|
||||
if not features:
|
||||
return None
|
||||
|
||||
@@ -263,8 +356,24 @@ def build_sheet_surface_model(
|
||||
# 등고 간격(m) — 마루 캡 크기의 근거. 레벨이 하나뿐이면 5m(1:5,000 주곡선) 폴백.
|
||||
contour_levels = np.unique(cloud.z)
|
||||
interval_m = float(np.diff(contour_levels).min()) if len(contour_levels) > 1 else 5.0
|
||||
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ContourCloud
|
||||
|
||||
burned, levels = rasterize_contours(spec, features, None)
|
||||
present = [level for level in sorted(levels) if bool((burned == level).any())]
|
||||
# 계단(평탄 삼각형) 해소 — 인접 등고선 사이 중간 보간선을 정점으로 추가.
|
||||
cloud = _densify_between_contours(spec, features, cloud)
|
||||
cloud = _densify_between_contours(spec, burned, present, cloud)
|
||||
# 계곡 구조선 — 하천중심선을 따라 등고 교차점 사이를 보간한 정점 추가.
|
||||
stream_features = _load_features_metric(processed_dir, epsg, _STREAM_FILE)
|
||||
if stream_features:
|
||||
stream_vertices = _stream_breakline_vertices(spec, burned, stream_features)
|
||||
if stream_vertices is not None:
|
||||
logger.info("도엽 서피스: 계곡 구조선 정점 %d개 추가", len(stream_vertices[1]))
|
||||
cloud = ContourCloud(
|
||||
xy=np.vstack([cloud.xy, stream_vertices[0]]),
|
||||
z=np.concatenate([cloud.z, stream_vertices[1]]),
|
||||
)
|
||||
surface = interpolate_elevation(spec, cloud) # (R, C), 북→남 행 순서, 외부 NaN
|
||||
# 마지막 등고선 안쪽 평탄 마루를 완만한 돔으로 (2026-08-30 사용자 확정).
|
||||
_cap_flat_summits(surface, spec.cell_m, interval_m)
|
||||
|
||||
Reference in New Issue
Block a user