feat(D 예정노선): 0차곡 흐름누적 = 도엽별 공용 창고 (PLAN 608-2 · 요구 486)
- 도엽 + 여유 200 m(창고에 있는 도엽만 · 없는 이웃 = 출구)로 셈 · storage/_shared/terrain/<도엽>/zero_acc_<epsg>_2m.npz - 열쇠 = 쓴 도엽 지형 열쇠 · 이웃 지형이 새로 생기면 다시 셈 - 0차곡 겹 = 도엽 누적 이어 붙임 · 임계 바꿈은 칸 고르기 · 선 뽑기만 · 창고 격자 아닌 옛 DEM 은 통째 셈 - 새 시험 608 = 조각 대 통째 0차곡 칸 차 0 · 저장 · 이웃 들어옴 다시 · 임계 바꿈 누적 그대로 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -1,37 +1,46 @@
|
||||
"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 요구 443).
|
||||
"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 요구 443 · 486).
|
||||
|
||||
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`.
|
||||
파생 저장본(열쇠 = 격자 지문 · 임계 · 칸).
|
||||
파생 저장본(열쇠 = 격자 지문 · 임계 · 칸 · 쓴 도엽 누적 지문).
|
||||
흐름누적 = 공용 지형 창고 도엽마다(도엽 + 여유 200 m · 창고에 있는 이웃 도엽만) 한 번 셈 →
|
||||
`storage/_shared/terrain/<도엽>/zero_acc_<epsg>_<칸>m.npz`(열쇠 = 쓴 도엽 지형 열쇠) ·
|
||||
임계 바꿈 = 저장 누적에서 칸 고르기 · 선 뽑기만 · 사업 DEM 이 창고 격자가 아니면(옛 받기) 통째 셈.
|
||||
|
||||
계산 자리 = ③ 서버 단독.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import json
|
||||
import math
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
from common_util import common_util_terrain_store as terrain_store
|
||||
from common_util.common_util_atomic import atomic_write_npz
|
||||
from common_util.common_util_derived_store import (
|
||||
derived_key,
|
||||
file_stamp,
|
||||
read_derived,
|
||||
write_derived,
|
||||
)
|
||||
from common_util.common_util_json import atomic_write_json
|
||||
from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD
|
||||
|
||||
ZERO_ORDER_KIND = "plan_zero_order"
|
||||
ZERO_ORDER_VERSION = 1
|
||||
ZERO_ORDER_VERSION = 2
|
||||
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 # 세선화 갈림목 토막 거름
|
||||
ZERO_ORDER_MARGIN_M = 200.0 # 도엽 누적 여유(통째 셈과 0차곡 칸 차 0 · 분석 486 3-1)
|
||||
ACC_CAP = 65535 # 저장 누적 상한(칸 수 · 2 m 면 26 ha · 임계 1 ha 보다 넉넉)
|
||||
SIMPLIFY_M = 1.0
|
||||
|
||||
_D8 = ((-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1))
|
||||
@@ -103,6 +112,28 @@ def zero_order_mask(
|
||||
return valid & (area >= min_m2) & (area < max_m2)
|
||||
|
||||
|
||||
def mask_lines(
|
||||
acc: np.ndarray,
|
||||
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(acc, 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_lines(
|
||||
z: np.ndarray,
|
||||
valid: np.ndarray,
|
||||
@@ -112,16 +143,108 @@ def zero_order_lines(
|
||||
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
|
||||
"""격자(행 = 남 → 북) 통째 흐름누적 → 0차곡 선."""
|
||||
return mask_lines(d8_accumulation(z, valid), z, valid, x, y, cell_m, min_m2, max_m2)
|
||||
|
||||
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 _read_json(path: Path) -> Any:
|
||||
try:
|
||||
return json.loads(path.read_text(encoding="utf-8"))
|
||||
except (OSError, ValueError):
|
||||
return None
|
||||
|
||||
|
||||
def _acc_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]:
|
||||
folder = terrain_store.TERRAIN_ROOT / str(sheet_no)
|
||||
stem = f"zero_acc_{int(epsg)}_{cell_m:g}m"
|
||||
return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock"
|
||||
|
||||
|
||||
def _dem_key(sheet_no: str, epsg: int, cell_m: float) -> dict[str, Any] | None:
|
||||
"""창고 지형 열쇠(셈 시간 뺌) — 안 셌거나 원천이 바뀌었으면 None."""
|
||||
stored = terrain_store._stored_ok(sheet_no, epsg, cell_m)
|
||||
return None if stored is None else {k: v for k, v in stored.items() if k != "seconds"}
|
||||
|
||||
|
||||
def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL_M) -> Path | None:
|
||||
"""도엽 하나의 흐름누적(창고) — 도엽 + 여유 200 m 지형(창고에 있는 도엽만 · 없는 이웃 =
|
||||
출구)으로 셈 · 도엽 칸만 저장(uint16 · `ACC_CAP` 로 자름) · 열쇠 = 쓴 도엽 지형 열쇠 ·
|
||||
지형이 없으면 None."""
|
||||
own = _dem_key(sheet_no, epsg, cell_m)
|
||||
dem_npz = terrain_store._names(sheet_no, epsg, cell_m)[0]
|
||||
if own is None or own.get("empty") or not dem_npz.is_file():
|
||||
return None
|
||||
with np.load(dem_npz) as data:
|
||||
ix0, iy0 = int(data["ix0"]), int(data["iy0"])
|
||||
rows, cols = data["z"].shape
|
||||
pad = math.ceil(ZERO_ORDER_MARGIN_M / cell_m)
|
||||
box = (
|
||||
(ix0 - pad) * cell_m,
|
||||
(iy0 - pad) * cell_m,
|
||||
(ix0 + cols + pad) * cell_m,
|
||||
(iy0 + rows + pad) * cell_m,
|
||||
)
|
||||
used = {}
|
||||
for no in terrain_store.sheets_for_box(box, epsg):
|
||||
dem = _dem_key(no, epsg, cell_m)
|
||||
if dem is not None and not dem.get("empty"):
|
||||
used[no] = dem
|
||||
key = {"version": ZERO_ORDER_VERSION, "margin_m": ZERO_ORDER_MARGIN_M, "dem": used}
|
||||
npz, key_path, lock_path = _acc_names(sheet_no, epsg, cell_m)
|
||||
if npz.is_file() and _read_json(key_path) == key:
|
||||
return npz
|
||||
with terrain_store._Lock(lock_path):
|
||||
if npz.is_file() and _read_json(key_path) == key: # 기다리는 동안 남이 셌을 수 있음
|
||||
return npz
|
||||
grid = terrain_store.load_terrain(epsg, sheets=sorted(used), box=box, cell_m=cell_m)
|
||||
acc = d8_accumulation(grid["z"], grid["valid_mask"])[pad : pad + rows, pad : pad + cols]
|
||||
key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게
|
||||
atomic_write_npz(
|
||||
npz,
|
||||
acc=np.minimum(acc, ACC_CAP).astype(np.uint16),
|
||||
ix0=np.int64(ix0),
|
||||
iy0=np.int64(iy0),
|
||||
cell_m=np.float64(cell_m),
|
||||
)
|
||||
atomic_write_json(key_path, key)
|
||||
return npz
|
||||
|
||||
|
||||
def store_accumulation_paths(
|
||||
x: np.ndarray, y: np.ndarray, epsg: int, cell_m: float
|
||||
) -> list[Path] | None:
|
||||
"""사업 DEM 격자(창고 격자 위)를 덮는 도엽 누적 — 창고 격자 아님 · 지형 안 센 도엽 = None."""
|
||||
fx, fy = x[0] / cell_m - 0.5, y[0] / cell_m - 0.5
|
||||
if cell_m != ZERO_ORDER_CELL_M or abs(fx - round(fx)) > 1e-6 or abs(fy - round(fy)) > 1e-6:
|
||||
return None
|
||||
half = cell_m / 2
|
||||
box = (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half)
|
||||
paths: list[Path] = []
|
||||
for no in terrain_store.sheets_for_box(box, epsg):
|
||||
dem = _dem_key(no, epsg, cell_m)
|
||||
if dem is None:
|
||||
return None
|
||||
path = None if dem.get("empty") else sheet_accumulation(no, epsg, cell_m)
|
||||
if path is not None:
|
||||
paths.append(path)
|
||||
return paths
|
||||
|
||||
|
||||
def paste_accumulation(
|
||||
paths: list[Path], x: np.ndarray, y: np.ndarray, cell_m: float
|
||||
) -> np.ndarray:
|
||||
"""도엽 누적을 사업 격자(y 오름차순)에 이어 붙임 — 없는 칸 0."""
|
||||
ix0, iy0 = round(x[0] / cell_m - 0.5), round(y[0] / cell_m - 0.5)
|
||||
acc = np.zeros((len(y), len(x)))
|
||||
for path in paths:
|
||||
with np.load(path) as data:
|
||||
px, py, part = int(data["ix0"]), int(data["iy0"]), data["acc"]
|
||||
r0, c0 = max(py, iy0), max(px, ix0)
|
||||
r1, c1 = min(py + part.shape[0], iy0 + len(y)), min(px + part.shape[1], ix0 + len(x))
|
||||
if r0 < r1 and c0 < c1:
|
||||
src = part[r0 - py : r1 - py, c0 - px : c1 - px]
|
||||
np.copyto(acc[r0 - iy0 : r1 - iy0, c0 - ix0 : c1 - ix0], src, where=src > 0)
|
||||
return acc
|
||||
|
||||
|
||||
def zero_order_layer(
|
||||
@@ -131,7 +254,8 @@ def zero_order_layer(
|
||||
min_m2: float = ZERO_ORDER_MIN_M2,
|
||||
max_m2: float = ZERO_ORDER_MAX_M2,
|
||||
) -> dict[str, Any]:
|
||||
"""0차곡 겹 — 지문(격자 파일 · 임계 · 칸) 같으면 저장본 · 경위도 선 + 가닥 · 길이 합 m."""
|
||||
"""0차곡 겹 — 지문(격자 파일 · 임계 · 칸 · 도엽 누적) 같으면 저장본 · 경위도 선 + 가닥 ·
|
||||
길이 합 m."""
|
||||
from shapely.geometry import LineString, mapping
|
||||
from shapely.ops import transform
|
||||
|
||||
@@ -141,6 +265,12 @@ def zero_order_layer(
|
||||
stamp = file_stamp(npz)
|
||||
if stamp is None:
|
||||
raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.")
|
||||
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))
|
||||
x, y, cell_m = x[::step], y[::step], cell * step
|
||||
paths = store_accumulation_paths(x, y, epsg, cell_m)
|
||||
key = derived_key(
|
||||
ZERO_ORDER_KIND,
|
||||
ZERO_ORDER_VERSION,
|
||||
@@ -152,23 +282,25 @@ def zero_order_layer(
|
||||
ZERO_ORDER_CELL_M,
|
||||
ZERO_ORDER_MIN_LEN_M,
|
||||
SIMPLIFY_M,
|
||||
None if paths is None else [f"{p.parent.name}:{file_stamp(p)}" for p in paths],
|
||||
)
|
||||
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)
|
||||
if paths is None:
|
||||
lines = zero_order_lines(z, valid, x, y, cell_m, min_m2, max_m2)
|
||||
else:
|
||||
acc = paste_accumulation(paths, x, y, cell_m)
|
||||
lines = mask_lines(acc, z, valid, x, y, cell_m, 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,
|
||||
"cell_m": cell_m,
|
||||
"count": len(shapes),
|
||||
"length_m": round(sum(s.length for s in shapes), 1),
|
||||
"collection": {
|
||||
|
||||
@@ -0,0 +1,104 @@
|
||||
"""D 0차곡 흐름누적 = 도엽별 공용 창고(PLAN 608-2) — 도엽 + 여유 200 m 조각 = 통째 셈 ·
|
||||
처음만 셈 · 이웃 도엽 지형이 새로 생기면 다시 · 임계 바꿈은 누적 그대로 · 창고 격자 아닌 DEM = 통째.
|
||||
|
||||
./venv/Scripts/python.exe -m pytest resources/tester/test_608_plan_zero_order_store.py -q
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from common_util import common_util_terrain_store as store
|
||||
from common_util.common_util_atomic import atomic_write_npz
|
||||
from common_util.common_util_json import atomic_write_json
|
||||
from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zo
|
||||
|
||||
EPSG = 5187
|
||||
WEST, EAST = "36901074", "36901075"
|
||||
CELL = 20.0 # 시험은 굵은 칸(여유 200 m = 10 칸)
|
||||
|
||||
|
||||
def _height(x, y):
|
||||
"""기운 면 + 골 · 마루(파장 240 m) — 웅덩이 없음 · 작은 집수."""
|
||||
return 0.08 * x + 6.0 * np.sin(y / 38.0) + 2.0 * np.sin(x / 55.0)
|
||||
|
||||
|
||||
def _put_sheet(no: str) -> None:
|
||||
ix0, iy0, own = store._owned(no, EPSG, CELL)
|
||||
rows, cols = own.shape
|
||||
gx, gy = np.meshgrid((np.arange(cols) + ix0 + 0.5) * CELL, (np.arange(rows) + iy0 + 0.5) * CELL)
|
||||
z = _height(gx - 300000.0, gy).astype(np.float32)
|
||||
z[~own] = np.nan
|
||||
npz, key_path, _ = store._names(no, EPSG, CELL)
|
||||
atomic_write_npz(npz, z=z, ix0=np.int64(ix0), iy0=np.int64(iy0), cell_m=np.float64(CELL))
|
||||
atomic_write_json(key_path, {**store._key(EPSG, CELL, {}), "source": "gpkg", "seconds": 1})
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def shared(tmp_path, monkeypatch):
|
||||
monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path / "terrain")
|
||||
calls = []
|
||||
real = zo.d8_accumulation
|
||||
monkeypatch.setattr(zo, "d8_accumulation", lambda *a: calls.append(1) or real(*a))
|
||||
return calls
|
||||
|
||||
|
||||
def _acc(no: str):
|
||||
with np.load(zo.sheet_accumulation(no, EPSG, CELL)) as data:
|
||||
return int(data["ix0"]), int(data["iy0"]), np.asarray(data["acc"], float)
|
||||
|
||||
|
||||
def test_tiles_equal_whole_and_saved(shared) -> None:
|
||||
calls = shared
|
||||
_put_sheet(WEST)
|
||||
assert zo.sheet_accumulation(EAST, EPSG, CELL) is None # 지형 안 센 도엽
|
||||
lone = _acc(WEST)[2]
|
||||
assert len(calls) == 1
|
||||
_acc(WEST)
|
||||
assert len(calls) == 1 # 두 번째 = 저장본
|
||||
_put_sheet(EAST) # 이웃 지형이 생김 = 다시 셈
|
||||
pieces = [_acc(WEST), _acc(EAST)]
|
||||
assert len(calls) == 3
|
||||
|
||||
whole = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL)
|
||||
ix0, iy0 = round(whole["x"][0] / CELL - 0.5), round(whole["y"][0] / CELL - 0.5)
|
||||
acc = zo.d8_accumulation(whole["z"], whole["valid_mask"])
|
||||
tiled = np.zeros_like(acc)
|
||||
for px, py, part in pieces:
|
||||
dst = tiled[py - iy0 : py - iy0 + part.shape[0], px - ix0 : px - ix0 + part.shape[1]]
|
||||
np.copyto(dst, part, where=part > 0)
|
||||
valid = whole["valid_mask"]
|
||||
cells = CELL * CELL
|
||||
for low, high in ((2, 8), (1, 10)): # 0차곡 칸 = 집수 칸 수 [low, high)
|
||||
a = zo.zero_order_mask(acc, valid, CELL, low * cells, high * cells)
|
||||
b = zo.zero_order_mask(tiled, valid, CELL, low * cells, high * cells)
|
||||
assert a.any() and np.array_equal(a, b)
|
||||
# 이웃 없이 셌던 서쪽 도엽은 동쪽 경계에서 달랐음(이웃이 들어와 고쳐짐)
|
||||
assert not np.array_equal(lone, pieces[0][2])
|
||||
|
||||
|
||||
def test_layer_uses_store_and_threshold_reuses(shared, tmp_path, monkeypatch) -> None:
|
||||
"""사업 DEM 이 창고 격자 위면 도엽 누적을 이어 붙임 · 임계 바꿈 = 누적 다시 안 셈."""
|
||||
calls = shared
|
||||
monkeypatch.setattr(zo, "ZERO_ORDER_CELL_M", CELL)
|
||||
_put_sheet(WEST)
|
||||
_put_sheet(EAST)
|
||||
grid = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL)
|
||||
x, y, z = grid["x"][5:-5], grid["y"][5:-5], grid["z"][5:-5, 5:-5]
|
||||
models = tmp_path / "models"
|
||||
models.mkdir()
|
||||
np.savez(
|
||||
models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz",
|
||||
x=x,
|
||||
y=y,
|
||||
z=np.nan_to_num(z, nan=0.0),
|
||||
valid_mask=np.isfinite(z),
|
||||
resolution=np.array([CELL], np.float32),
|
||||
)
|
||||
cells = CELL * CELL
|
||||
first = zo.zero_order_layer(tmp_path, models, EPSG, 2 * cells, 8 * cells)
|
||||
assert first["count"] >= 1 and first["cell_m"] == CELL
|
||||
assert len(calls) == 2 # 도엽 둘 한 번씩
|
||||
zo.zero_order_layer(tmp_path, models, EPSG, 1 * cells, 10 * cells)
|
||||
assert len(calls) == 2 # 임계 바꿈 = 저장 누적에서 칸 고르기만
|
||||
# 창고 격자에서 반 칸 어긋난 DEM = 통째 셈
|
||||
assert zo.store_accumulation_paths(x + CELL / 2, y, EPSG, CELL) is None
|
||||
Reference in New Issue
Block a user