From 0dae8e301632ada7625de0f3d88f17f7dfe37884 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Sun, 30 Aug 2026 15:10:04 +0900 Subject: [PATCH] =?UTF-8?q?feat(B04):=20=EB=8F=84=EC=97=BD=20=EC=84=9C?= =?UTF-8?q?=ED=94=BC=EC=8A=A4=EB=A5=BC=20=EA=B5=AC=EC=A1=B0=EC=84=A0(?= =?UTF-8?q?=EA=B3=84=EA=B3=A1=C2=B7=EB=8A=A5=EC=84=A0)=20=EA=B8=B0?= =?UTF-8?q?=EB=B0=98=20=EB=B3=B4=EA=B0=84=EC=9C=BC=EB=A1=9C=20=EA=B0=9C?= =?UTF-8?q?=EC=84=A0=ED=95=9C=EB=8B=A4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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) --- .../B04_PreProcess_Engine_SheetSurface.py | 145 +++++++++++++++--- 1 file changed, 127 insertions(+), 18 deletions(-) diff --git a/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py b/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py index 14bdbd95..a06831ab 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py +++ b/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py @@ -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)