eomsangdon/Aislo
1
0
Fork 0
Code Issues Pull Requests 2 Actions Packages Projects Releases Wiki Activity

fix(B04): 등고선 보간을 등간격으로 — 얇은 래스터·앵커 제거·최소 다듬기

Browse Source 
사용자 지적(1m 보간선이 등간격 아님)의 원인 세 가지를 잡았다.
- 등고 라인 래스터가 all_touched=True라 2px 두께였다. 그 폭만큼 정확히 5m
  배수 표고인 평탄 띠가 생겨(격자의 7.6%) 사이 1m 선 간격이 찌그러졌다.
  서피스 전용 얇은 래스터화(all_touched=False)로 바꾸고, 봉우리 폐합 링이
  사라지지 않도록 길이 필터도 뺐다.
- 계곡 앵커를 1m로 양자화해 제약에 섞은 탓에 제약 표고가 29단→121단이 되어
  계곡 주변만 1m 간격이 됐다. 앵커를 뺐다 — V자 등고선이 계곡 하강을 이미
  담고 있어 거리 보간만으로 충분하다.
- 등고선 고정 완화를 수렴시키면 harmonic 해가 되어 마루가 눌린다. z=r은
  biharmonic이라서다 — ANUDEM/Topo to Raster가 라플라스가 아니라 thin plate
  spline을 쓰는 이유(Hutchinson 1988/89, 조사 결과). 다듬기를 omega=1·5회로
  줄였다. 합성 원뿔 검증: 거리보간 오차 0.121m·간격 CV 2.18 → 5회 0.118m·1.43,
  20회 이상 CV 50↑ 악화. 원뿔 회귀 테스트를 남겼다.

실측(c1bb453f): 평탄 셀 0.00%, 5m 구간 표고 분포 8.4~10.3%(균등),
생성 6.3s, 노선 |Δz| 2.64m.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
eomsangdon
2026-08-30 15:51:08 +09:00
co-authored by Claude Opus 5
parent 31e3eb4ecf
commit 15ce46d54d
2 changed files with 112 additions and 97 deletions
Download Patch File Download Diff File Expand all files Collapse all files
B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py
+108 -91
View File
@@ -9,8 +9,9 @@ LAS 없는 설계(2026-08-30 사용자 확정)의 지형 원천이자, LAS가
격자는 `SHEET_SURFACE_GRID_M`(1m).
순서는 **2D 먼저, 메시는 맨 마지막**이다(2026-08-30 사용자 지시):
① 등고 라인을 격자에 굽고 ② 계곡 구조선 앵커를 제약으로 더한 뒤
③ 등고선 사이 거리 비례 보간으로 표고 격자를 만들고 ④ 마루를 캡한다.
① 등고 라인을 격자에 굽고 ② 등고선 사이 거리 비례 보간으로 표고 격자를 만든 뒤
③ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장하고 ④ 라인을 고정한 채
완화(라플라스)해 등고 간격을 고르게 한다.
격자에서 뽑는 1m 등고선이 곧 2D 보간선이며, 메시(glb)는 그 격자의 표현일 뿐이다.
"""
@@ -34,6 +35,7 @@ from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import grid_spec_from_b
from config.config_system import (
SHEET_SURFACE_GRID_M,
SHEET_SURFACE_MARGIN_M,
SHEET_SURFACE_RELAX_ITERATIONS,
SURFACE_MAX_PREVIEW_VERTICES,
)
@@ -41,7 +43,6 @@ 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"
@@ -104,72 +105,107 @@ def _preview_mesh(
return clip_and_compact_mesh(vertices, faces, pv.reshape(-1))
def _stream_breakline_vertices(
spec: Any, burned: np.ndarray, stream_features: list[dict[str, Any]]
) -> tuple[np.ndarray, np.ndarray] | None:
"""계곡 기준선(하천중심선)을 따라 등고 교차점 사이를 보간한 정점열을 만든다.
def _rasterize_contour_levels(spec: Any, features: list[dict[str, Any]]) -> np.ndarray:
"""등고 라인을 격자에 굽는다 — 셀 = 그 위를 지나는 라인의 표고, 그 외 NaN.
구조선 기반 보간(2026-08-30 사용자 확정): 계곡 축을 먼저 세우고, 축이 등고선과
만나는 점을 Z 앵커로 삼아 앵커 사이를 선 길이 비례로 보간한다. 계곡 바닥이
등고선 사이에서도 연속으로 내려가는 가상 종단이 되어 골짜기 평탄·역경사가 준다.
배수유역의 `rasterize_contours()`와 두 가지가 다르다(둘 다 서피스 품질 때문이다):
· `all_touched=False` — 스치는 셀까지 칠하면 라인이 2px 두께가 되고, 그 폭만큼
정확히 등고 표고인 **평탄 띠**가 생겨 사이 1m 등고선 간격이 찌그러진다
(2026-08-30 사용자 지적: 보간선이 등간격이 아님).
· 길이 필터 없음 — 봉우리 폐합 링 같은 짧은 등고선을 버리면 그 일대가 통째로
평평해진다. 배수유역은 노이즈를 버려야 하지만 지형면은 다 있어야 한다.
같은 셀을 두 표고가 지나면 낮은 쪽을 남긴다(배수유역과 같은 규칙).
"""
from shapely import segmentize
from rasterio.features import rasterize
from shapely.geometry import shape
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import iter_linestrings
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
ELEVATION_KEYS,
grid_transform,
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:
by_level: dict[float, list[Any]] = {}
for feature in features:
geometry = feature.get("geometry")
if not geometry:
continue
properties = feature.get("properties") or {}
elevation = next(
(float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None),
None,
)
if elevation is None:
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)
by_level.setdefault(elevation, []).append(line)
burned = np.full((spec.n_rows, spec.n_cols), np.nan, dtype=np.float32)
transform = grid_transform(spec)
for elevation in sorted(by_level, reverse=True):
stamp = rasterize(
[(line, 1) for line in by_level[elevation]],
out_shape=(spec.n_rows, spec.n_cols),
transform=transform,
fill=0,
dtype="uint8",
all_touched=False,
).astype(bool)
burned[stamp] = elevation
logger.info(
"도엽 서피스: 등고 라인 %d단을 격자에 굽어 %d셀",
len(by_level),
int(np.isfinite(burned).sum()),
)
return burned
def _relax_surface(surface: np.ndarray, fixed: np.ndarray, iterations: int) -> None:
"""등고 라인을 고정한 채 이웃 평균으로 표고를 **살짝만** 다듬는다 (in-place).
거리 보간은 두 라인의 굽은 정도가 다르면 능선·계곡 중심축에 각진 자국이 남는다.
이웃 평균 몇 회로 그 자국만 지운다.
**끝까지 수렴시키면 안 된다.** 수렴한 해는 라플라스 방정식의 해(harmonic)인데,
원뿔·능선 형상 z=r은 harmonic이 아니라 biharmonic이라 harmonic 해로 끌고 가면
마루가 눌리고 평탄 셀이 늘어난다. 등고선 DEM의 표준인 ANUDEM/Topo to Raster가
라플라스가 아니라 **thin plate spline(biharmonic)** 을 쓰는 이유가 이것이다
(Hutchinson 1988/89). 거리 보간은 원뿔을 정확히 재현하므로 그것을 바탕으로 두고
다듬기만 한다.
합성 원뿔 검증(2026-08-30): 거리보간만 오차 0.121m·간격 CV 2.18 →
5회 다듬기 0.118m·1.43(개선) → 20회 이상 CV 50↑·평탄 셀 증가(악화).
"""
if iterations <= 0:
return
free = ~fixed & np.isfinite(surface)
if not free.any():
return
# 체커보드(red-black) 순서 — 한 색을 갱신할 때 이웃(다른 색)은 그대로라 제자리
# 갱신이 안전하다. 과완화는 쓰지 않는다(omega=1): 빨리 수렴시킬수록 harmonic 해에
# 가까워져 능선이 눌린다. 여기서 원하는 건 수렴이 아니라 자국 제거다.
omega = np.float32(1.0)
rows, cols = np.indices(surface.shape)
red = free & (((rows + cols) & 1) == 0)
black = free & ~red
padded = np.zeros((surface.shape[0] + 2, surface.shape[1] + 2), dtype=np.float32)
for _ in range(iterations):
for colour in (red, black):
padded[1:-1, 1:-1] = surface
padded[0, 1:-1] = surface[0]
padded[-1, 1:-1] = surface[-1]
padded[1:-1, 0] = surface[:, 0]
padded[1:-1, -1] = surface[:, -1]
neighbours = (
padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:]
) * np.float32(0.25)
surface[colour] += omega * (neighbours[colour] - surface[colour])
logger.info("도엽 서피스: 등고선 고정 완화 %d회 (자유 셀 %d개)", iterations, int(free.sum()))
def _interpolate_between_contours(burned: np.ndarray, cell_m: float) -> np.ndarray:
@@ -229,27 +265,9 @@ def _interpolate_between_contours(burned: np.ndarray, cell_m: float) -> np.ndarr
return surface
def _burn_stream_anchors(
spec: Any, burned: np.ndarray, stream_features: list[dict[str, Any]]
) -> int:
"""계곡 구조선 앵커 보간값을 등고 격자에 제약으로 굽는다. 구운 셀 수 반환."""
vertices = _stream_breakline_vertices(spec, burned, stream_features)
if vertices is None:
return 0
xy, z = vertices
row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy())
inside = (row >= 0) & (col >= 0)
row, col, z = row[inside], col[inside], z[inside]
free = ~np.isfinite(burned[row, col]) # 등고 라인 위는 덮지 않는다
# 1m로 양자화 — 거리장을 표고별로 한 번씩 도는 방식이라 레벨 수가 곧 비용이다.
# 원천이 5m 주곡선이므로 0.5m 이내 반올림은 정확도에 영향이 없다.
burned[row[free], col[free]] = np.round(z[free]).astype(np.float32)
return int(free.sum())
def _resolve_enclosed_interiors(
burned: np.ndarray, present: list[float], surface: np.ndarray, cell_m: float, interval_m: float
) -> int:
) -> np.ndarray:
"""폐합 등고선 안쪽(봉우리·웅덩이)을 바깥 사면 경사로 연장한다 (in-place).
거리 보간은 "가장 가까운 서로 다른 두 라인 사이"를 채우므로, 마지막 등고선
@@ -260,7 +278,7 @@ def _resolve_enclosed_interiors(
폐합 라인 내부에 다른 표고 제약이 하나도 없으면 그 안이 마루(바깥이 낮을 때)
또는 웅덩이(바깥이 높을 때)다. 바깥 사면의 국소 경사를 안쪽으로 연장하되
±(간격−0.5m)로 제한한다 — 다음 등고선이 없다는 사실과 모순되지 않는 범위다.
처리한 영역 수를 반환한다.
처리한 영역의 마스크를 반환한다 — 이어지는 완화에서 함께 고정해야 한다.
"""
from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label
@@ -268,6 +286,7 @@ def _resolve_enclosed_interiors(
band_cells = 8
limit = max(interval_m - 0.5, 0.5)
resolved = 0
handled = np.zeros(burned.shape, dtype=bool)
for level in present:
mask = burned == level
interior = binary_fill_holes(mask) & ~mask
@@ -301,10 +320,11 @@ def _resolve_enclosed_interiors(
peak = float(inner.max())
offset = (limit / 2.0) * (inner[component] / peak) if peak > 0 else 0.0
surface[component] = level + direction * offset
handled |= component
resolved += 1
if resolved:
logger.info("도엽 서피스: 폐합 등고선 내부 %d곳을 사면 경사로 연장", resolved)
return resolved
return handled
def build_sheet_surface_model(
@@ -331,28 +351,25 @@ def build_sheet_surface_model(
spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M)
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours
# ① 2D — 등고 라인을 격자에 굽는다(라인 셀 = 표고, 그 외 NaN).
burned, levels = rasterize_contours(spec, features, None)
present = sorted({level for level in levels if bool((burned == level).any())})
burned = _rasterize_contour_levels(spec, features)
present = sorted(np.unique(burned[np.isfinite(burned)]).tolist())
if len(present) < 2:
logger.warning("도엽 서피스: 절취 범위 안에 등고선이 부족합니다.")
return None
# 등고 간격(m) — 마루 캡 크기의 근거. 레벨이 하나뿐이면 5m(1:5,000 주곡선) 폴백.
# 등고 간격(m) — 마루 연장 상한의 근거. 레벨이 하나뿐이면 5m(주곡선) 폴백.
interval_m = float(np.diff(np.array(present)).min()) if len(present) > 1 else 5.0
# ② 2D — 계곡 구조선(하천중심선) 앵커 보간값을 같은 격자에 제약으로 굽는다.
stream_features = _load_features_metric(processed_dir, epsg, _STREAM_FILE)
if stream_features:
burned_stream = _burn_stream_anchors(spec, burned, stream_features)
logger.info("도엽 서피스: 계곡 구조선 제약 %d셀", burned_stream)
# ③ 2D — 등고선 사이 거리 비례 보간(메시 없음). 여기서 나온 격자에서 1m 등고선을
# ② 2D — 등고선 사이 거리 비례 보간(메시 없음). 여기서 나온 격자에서 1m 등고선을
# 뽑으므로 화면 등고선이 곧 2D 보간선이다(2026-08-30 사용자 지시).
# 제약은 **등고선만** 쓴다 — 계곡 구조선 앵커를 1m로 양자화해 섞었더니 제약
# 표고가 29단→121단이 되어 계곡 주변만 1m 간격이 되고 보간선이 등간격을
# 잃었다(2026-08-30 사용자 지적). V자 등고선이 계곡 하강을 이미 담고 있다.
surface = _interpolate_between_contours(burned, spec.cell_m)
# ④ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장 (2026-08-30 사용자 확정).
_resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m)
# ③ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장 (2026-08-30 사용자 확정).
summits = _resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m)
# ④ 등고 라인(과 마루)을 고정한 채 완화 — 중심축 뭉침을 풀어 간격을 고르게 한다.
_relax_surface(surface, np.isfinite(burned) | summits, SHEET_SURFACE_RELAX_ITERATIONS)
# DtmGridSampler 규약에 맞춰 y 오름차순으로 뒤집어 저장한다.
x_coords = spec.cell_centers_x()