feat(D03 겹): 노선분석 지도 단층 · 지질구조선 · 0차곡 겹 — 묶음 셋째 · 처음 끔 (PLAN 584-1 · 요구 443)
- 서버: 선 겹 = 브이월드 WFS lt_l_gimsfault · lt_l_gimslinea 받은 범위로 잘라 가닥 · 길이 · 덮으면 다시 안 받음 · 실패 = 502 · 새 선 겹은 LINE_KINDS 한 줄 - 서버: 0차곡 = 사업 DEM 2 m 솎음 D8 흐름누적(웅덩이 메움 · 평지 넘김 · numpy) · 위 집수 0.1~1 ha 칸 세선화 선 · 10 m 미만 토막 거름 · 파생 저장본(격자 지문 · 임계) · 임계 = 요청 min_m2 · max_m2 - 화면: 단추 셋(bolt · lineDashed · arrowsJoin 부품 상자에 더함) · 색 --draw-map-fault · lineament · zero-order · 범례 = 가닥 · km · 첫 셈 늦으면 알림 · 바탕 캔버스 - 선 긋기 · 그림 겹 놓기 · 그림 읽기를 LayerDefs 로 옮김(700줄 한도) - test_553_4: 겹 목록에 faults · lineaments · zeroOrder 더함(새 겹 · 기대 목록 갱신) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -0,0 +1,183 @@
|
||||
"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 요구 443).
|
||||
|
||||
0차곡 = 뚜렷한 유로가 없는 곡두의 집수지형(`노선_배치.md:31` · 교본 `3. 임도망의 배치 절차.md:156`)
|
||||
→ 사업 DEM(평활 안 한 `models_plan/dtm_sheet_<방법>.npz`)을 2 m 로 솎아 D8 흐름누적 ·
|
||||
위 집수 넓이가 `min_m2` 이상 · `max_m2` 미만인 칸 = 0차곡(그보다 크면 1차 계류) →
|
||||
세선화 · 선 뽑기(B05 `_mask_to_polylines`) · 경위도 선 + 가닥 · 길이 합.
|
||||
임계값 = 지식DB 수치 없음 → 정함(브레인) 0.1 ha · 1 ha · 바꿈 = 요청 `min_m2` · `max_m2`.
|
||||
파생 저장본(열쇠 = 격자 지문 · 임계 · 칸).
|
||||
|
||||
계산 자리 = ③ 서버 단독.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
from common_util.common_util_derived_store import (
|
||||
derived_key,
|
||||
file_stamp,
|
||||
read_derived,
|
||||
write_derived,
|
||||
)
|
||||
from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD
|
||||
|
||||
ZERO_ORDER_KIND = "plan_zero_order"
|
||||
ZERO_ORDER_VERSION = 1
|
||||
ZERO_ORDER_MIN_M2 = 1_000.0 # 집수 시작(0.1 ha)
|
||||
ZERO_ORDER_MAX_M2 = 10_000.0 # 이 넓이부터 1차 계류(1 ha)
|
||||
ZERO_ORDER_CELL_M = 2.0 # 흐름누적 칸(1 m 격자는 솎음)
|
||||
ZERO_ORDER_MIN_LEN_M = 10.0 # 세선화 갈림목 토막 거름
|
||||
SIMPLIFY_M = 1.0
|
||||
|
||||
_D8 = ((-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1))
|
||||
|
||||
|
||||
def _shift(grid: np.ndarray, dr: int, dc: int, fill: Any) -> np.ndarray:
|
||||
"""칸마다 (dr, dc) 이웃 값 — 격자 밖 = `fill`."""
|
||||
rows, cols = grid.shape
|
||||
pad = np.pad(grid, 1, constant_values=fill)
|
||||
return pad[1 + dr : 1 + dr + rows, 1 + dc : 1 + dc + cols]
|
||||
|
||||
|
||||
def d8_accumulation(z: np.ndarray, valid: np.ndarray) -> np.ndarray:
|
||||
"""D8 흐름누적(칸 수 · 자기 칸 포함 · 빈 칸 0).
|
||||
|
||||
웅덩이 메움(형태 재구성 · 가장자리 · 빈 칸 곁 = 출구) → 가장 가파른 아래 이웃 하나로 흐름 →
|
||||
평지(등고선 칸 · 메운 웅덩이)는 이미 흐르는 같은 높이 이웃으로 이어 감 → 위에서 아래로 더함.
|
||||
"""
|
||||
from skimage.morphology import reconstruction
|
||||
|
||||
rows, cols = z.shape
|
||||
height = np.where(valid, np.asarray(z, float), np.nan)
|
||||
floor = np.nanmin(height) if valid.any() else 0.0
|
||||
edge = valid & ~np.all([_shift(valid, dr, dc, False) for dr, dc in _D8], axis=0)
|
||||
ground = np.where(valid, height, floor)
|
||||
seed = np.where(edge | ~valid, ground, ground.max())
|
||||
filled = np.where(valid, reconstruction(seed, ground, method="erosion"), np.nan)
|
||||
index = np.arange(rows * cols).reshape(rows, cols)
|
||||
best = np.zeros(z.shape)
|
||||
receiver = np.full(z.shape, -1, np.int64)
|
||||
with np.errstate(invalid="ignore"):
|
||||
for dr, dc in _D8:
|
||||
drop = (filled - _shift(filled, dr, dc, np.nan)) / math.hypot(dr, dc)
|
||||
better = drop > best
|
||||
best = np.where(better, drop, best)
|
||||
receiver = np.where(better, index + dr * cols + dc, receiver)
|
||||
done = (receiver >= 0) | edge | ~valid
|
||||
while True:
|
||||
grew = False
|
||||
for dr, dc in _D8:
|
||||
take = ~done & _shift(done & valid, dr, dc, False)
|
||||
take &= _shift(filled, dr, dc, np.nan) == filled
|
||||
if take.any():
|
||||
receiver[take] = index[take] + dr * cols + dc
|
||||
done |= take
|
||||
grew = True
|
||||
if not grew:
|
||||
break
|
||||
down = receiver.ravel().tolist()
|
||||
acc = valid.astype(float).ravel().tolist()
|
||||
upstream = np.bincount(receiver[receiver >= 0], minlength=rows * cols).tolist()
|
||||
stack = [i for i, n in enumerate(upstream) if n == 0]
|
||||
while stack:
|
||||
i = stack.pop()
|
||||
j = down[i]
|
||||
if j >= 0:
|
||||
acc[j] += acc[i]
|
||||
upstream[j] -= 1
|
||||
if upstream[j] == 0:
|
||||
stack.append(j)
|
||||
return np.asarray(acc).reshape(rows, cols)
|
||||
|
||||
|
||||
def zero_order_mask(
|
||||
acc: np.ndarray, valid: np.ndarray, cell_m: float, min_m2: float, max_m2: float
|
||||
) -> np.ndarray:
|
||||
"""위 집수 넓이(칸 수 × 칸 넓이)가 [min_m2, max_m2) 인 칸."""
|
||||
area = acc * cell_m * cell_m
|
||||
return valid & (area >= min_m2) & (area < max_m2)
|
||||
|
||||
|
||||
def zero_order_lines(
|
||||
z: np.ndarray,
|
||||
valid: np.ndarray,
|
||||
x: np.ndarray,
|
||||
y: np.ndarray,
|
||||
cell_m: float,
|
||||
min_m2: float = ZERO_ORDER_MIN_M2,
|
||||
max_m2: float = ZERO_ORDER_MAX_M2,
|
||||
) -> list[list[tuple[float, float]]]:
|
||||
"""격자(행 = 남 → 북) → 0차곡 선(작업 좌표계 m 점 목록 · `ZERO_ORDER_MIN_LEN_M` 미만 거름)."""
|
||||
from B05_Profile.B05_Profile_Engine_Skeleton import _mask_to_polylines
|
||||
|
||||
mask = zero_order_mask(d8_accumulation(z, valid), valid, cell_m, min_m2, max_m2)
|
||||
lines = [[(p[0], p[1]) for p in line["polyline"]] for line in _mask_to_polylines(mask, x, y, z)]
|
||||
return [
|
||||
line
|
||||
for line in lines
|
||||
if sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= ZERO_ORDER_MIN_LEN_M
|
||||
]
|
||||
|
||||
|
||||
def zero_order_layer(
|
||||
project_root: Path,
|
||||
models_dir: Path,
|
||||
epsg: int,
|
||||
min_m2: float = ZERO_ORDER_MIN_M2,
|
||||
max_m2: float = ZERO_ORDER_MAX_M2,
|
||||
) -> dict[str, Any]:
|
||||
"""0차곡 겹 — 지문(격자 파일 · 임계 · 칸) 같으면 저장본 · 경위도 선 + 가닥 · 길이 합 m."""
|
||||
from shapely.geometry import LineString, mapping
|
||||
from shapely.ops import transform
|
||||
|
||||
from D03_PlanAnalysis.D03_PlanAnalysis_Layers import _transforms
|
||||
|
||||
npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz"
|
||||
stamp = file_stamp(npz)
|
||||
if stamp is None:
|
||||
raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.")
|
||||
key = derived_key(
|
||||
ZERO_ORDER_KIND,
|
||||
ZERO_ORDER_VERSION,
|
||||
npz.name,
|
||||
stamp,
|
||||
int(epsg),
|
||||
float(min_m2),
|
||||
float(max_m2),
|
||||
ZERO_ORDER_CELL_M,
|
||||
ZERO_ORDER_MIN_LEN_M,
|
||||
SIMPLIFY_M,
|
||||
)
|
||||
saved = read_derived(project_root, ZERO_ORDER_KIND, key)
|
||||
if saved is not None:
|
||||
return saved["value"]
|
||||
with np.load(npz) as dem:
|
||||
x, y = np.asarray(dem["x"], float), np.asarray(dem["y"], float)
|
||||
cell = float(dem["resolution"][0]) if "resolution" in dem.files else float(x[1] - x[0])
|
||||
step = max(1, round(ZERO_ORDER_CELL_M / cell))
|
||||
z = np.asarray(dem["z"], float)[::step, ::step]
|
||||
valid = np.asarray(dem["valid_mask"], bool)[::step, ::step]
|
||||
lines = zero_order_lines(z, valid, x[::step], y[::step], cell * step, min_m2, max_m2)
|
||||
_, to_ll = _transforms(epsg)
|
||||
shapes = [LineString(line).simplify(SIMPLIFY_M) for line in lines]
|
||||
value = {
|
||||
"min_m2": float(min_m2),
|
||||
"max_m2": float(max_m2),
|
||||
"cell_m": cell * step,
|
||||
"count": len(shapes),
|
||||
"length_m": round(sum(s.length for s in shapes), 1),
|
||||
"collection": {
|
||||
"type": "FeatureCollection",
|
||||
"features": [
|
||||
{"type": "Feature", "geometry": mapping(transform(to_ll, s)), "properties": {}}
|
||||
for s in shapes
|
||||
],
|
||||
},
|
||||
}
|
||||
write_derived(project_root, ZERO_ORDER_KIND, key, value)
|
||||
return value
|
||||
Reference in New Issue
Block a user