perf(D03): 0차곡 = 창고 격자 · 도엽 선 미리 뽑아 이어 붙임 · 무거운 셈 일꾼 (PLAN 616-1 · 616-2 · 616-3)

- 0차곡 셈 = 사업 DEM 대신 창고 격자(받은 범위 도엽 누적) · 옛 격자 통째 갈래 걷음
- 기본 임계 선 = 도엽마다 창고 `zero_lines_<epsg>_2m.npz`(여유 200 m 에서 뽑아 끝점 도엽만 · 단순화 끝) · 요청 = 이어 붙여 받은 범위로 자름(파생 저장본 안 씀) · 다른 임계 = 누적 이어 붙여 셈 · 저장본
- 창고 누적 열쇠 판 그대로(ACC_VERSION 2 · 전국 누적 다시 셈 없음) · 파생 판 3
- compute=False = 셈 필요하면 None · 서버 본체 밖 일꾼 하나(`heavy_pool` · 띄울 때 빈 __main__ · 서버 끝나면 따라 끝남)
- 미리 계산 단계 zero_lines 덧(zero 뒤 · 기본 단계에 듦)
- 시험 test_616 새로 · test_608 반 칸 어긋난 DEM 통째 셈 확인 걷음(갈래 없어짐) · test_584_1 사업 DEM 통째 0차곡 겹 시험 걷음(창고 길로 바뀜 · 616 시험이 대신)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
2026-10-31 19:18:43 +09:00
co-authored by Claude Opus 5.5
parent 2648c92b70
commit 66f9c68a1e
5 changed files with 551 additions and 138 deletions
+341 -103
View File
@@ -1,14 +1,18 @@
"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 요구 443 · 486).
"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 616 · 요구 443 · 497).
0차곡 = 뚜렷한 유로가 없는 곡두의 집수지형(`노선_배치.md:31` · 교본 `3. 임도망의 배치 절차.md:156`)
→ 사업 DEM(평활 안 한 `models_plan/dtm_sheet_<방법>.npz`)을 2 m 로 솎아 D8 흐름누적 ·
→ 공용 지형 창고 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 이 창고 격자가 아니면(옛 받기) 통째 셈.
흐름누적 = 창고 도엽마다(도엽 + 여유 200 m · 창고에 있는 이웃 도엽만) 한 번 셈 →
`storage/_shared/terrain/<도엽>/zero_acc_<epsg>_<칸>m.npz`(열쇠 = 쓴 도엽 지형 열쇠).
기본 임계 선 = 도엽마다 미리 뽑아 창고에(`zero_lines_<epsg>_<칸>m.npz` · 도엽 + 여유 200 m 에서 뽑아
끝점(작은 쪽)이 그 도엽 안인 선만 · 단순화 끝) → 요청 = 도엽 선 이어 붙여 받은 범위로 자름(616-2).
다른 임계 = 받은 범위 도엽 누적 이어 붙여 칸 고르기 · 선 뽑기(616-1 · 사업 DEM 격자 안 씀).
다른 임계만 파생 저장본(열쇠 = 받은 범위 · 임계 · 칸 · 쓴 도엽 누적 지문).
무거운 셈(창고에 없는 지형 · 누적 · 선 · 다른 임계 첫 셈) =
서버 본체 밖 일꾼 하나(616-3 · `heavy_pool`).
계산 자리 = ③ 서버 단독.
"""
@@ -16,7 +20,14 @@
from __future__ import annotations
import json
import logging
import math
import os
import sys
import threading
import time
import types
from concurrent.futures import ProcessPoolExecutor
from pathlib import Path
from typing import Any
@@ -33,13 +44,17 @@ from common_util.common_util_derived_store import (
from common_util.common_util_json import atomic_write_json
from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD
logger = logging.getLogger(__name__)
ZERO_ORDER_KIND = "plan_zero_order"
ZERO_ORDER_VERSION = 2
ZERO_ORDER_VERSION = 3 # 창고 격자 · 도엽 선(616)
ACC_VERSION = 2 # 창고 누적 열쇠 판 — 올리면 전국 누적을 다시 셈
LINES_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_CELL_M = 2.0 # 흐름누적 칸
ZERO_ORDER_MIN_LEN_M = 10.0 # 세선화 갈림목 토막 거름
ZERO_ORDER_MARGIN_M = 200.0 # 도엽 누적 여유(통째 셈과 0차곡 칸 차 0 · 분석 486 3-1)
ZERO_ORDER_MARGIN_M = 200.0 # 도엽 누적 · 선 여유(통째 셈과 0차곡 칸 차 0 · 분석 486 3-1)
ACC_CAP = 65535 # 저장 누적 상한(칸 수 · 2 m 면 26 ha · 임계 1 ha 보다 넉넉)
SIMPLIFY_M = 1.0
@@ -159,6 +174,13 @@ def zero_order_lines(
return mask_lines(d8_accumulation(z, valid), z, valid, x, y, cell_m, min_m2, max_m2)
class NeedsCompute(Exception):
"""`compute=False` 인데 창고에 없거나 낡아 셈이 필요함 — 부르는 쪽이 일꾼으로 넘김."""
Memo = dict[tuple[str, int, float], Any]
def _read_json(path: Path) -> Any:
try:
return json.loads(path.read_text(encoding="utf-8"))
@@ -166,29 +188,62 @@ def _read_json(path: Path) -> Any:
return None
def _acc_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]:
def _store_names(prefix: str, 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"
stem = f"{prefix}_{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."""
def _acc_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]:
return _store_names("zero_acc", sheet_no, epsg, cell_m)
def _lines_names(sheet_no: str, epsg: int, cell_m: float) -> tuple[Path, Path, Path]:
return _store_names("zero_lines", sheet_no, epsg, cell_m)
def _dem_key(sheet_no: str, epsg: int, cell_m: float, memo: Memo | None = None) -> dict | None:
"""창고 지형 열쇠(셈 시간 뺌) — 안 셌거나 원천이 바뀌었으면 None · `memo` = 한 요청 안."""
slot = (str(sheet_no), int(epsg), float(cell_m))
if memo is not None and slot in memo:
return memo[slot]
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"}
key = None if stored is None else {k: v for k, v in stored.items() if k != "seconds"}
if memo is not None:
memo[slot] = key
return key
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:
def _has_terrain(
sheet_no: str, epsg: int, cell_m: float, compute: bool = True, memo: Memo | None = None
) -> bool:
"""창고 지형이 있음(없으면 셈) — 빈 도엽 · 받기 실패 = False."""
dem = _dem_key(sheet_no, epsg, cell_m, memo)
if dem is None:
if not compute:
raise NeedsCompute(sheet_no)
try:
terrain_store.ensure_sheet_terrain(sheet_no, epsg, cell_m)
except RuntimeError as error:
logger.warning("0차곡 지형 없음 %s: %s", sheet_no, error)
return False
if memo is not None:
memo.pop((str(sheet_no), int(epsg), float(cell_m)), None)
dem = _dem_key(sheet_no, epsg, cell_m, memo)
return dem is not None and not dem.get("empty")
def _sheet_frame(path: Path, cell_m: float) -> tuple[int, int, int, int, tuple[float, ...]]:
"""창고 도엽 파일(지형 · 누적) → (ix0, iy0, 행, 열, 여유 200 m 상자) — 큰 칸은 머리만 읽음."""
import zipfile
with np.load(path) as data:
ix0, iy0 = int(data["ix0"]), int(data["iy0"])
rows, cols = data["z"].shape
name = "z.npy" if "z" in data.files else "acc.npy"
with zipfile.ZipFile(path) as archive, archive.open(name) as head:
version = np.lib.format.read_magic(head)
read = np.lib.format.read_array_header_1_0
rows, cols = (read if version == (1, 0) else np.lib.format.read_array_header_2_0)(head)[0]
pad = math.ceil(ZERO_ORDER_MARGIN_M / cell_m)
box = (
(ix0 - pad) * cell_m,
@@ -196,15 +251,36 @@ def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL
(ix0 + cols + pad) * cell_m,
(iy0 + rows + pad) * cell_m,
)
return ix0, iy0, rows, cols, box
def sheet_accumulation(
sheet_no: str,
epsg: int,
cell_m: float = ZERO_ORDER_CELL_M,
compute: bool = True,
memo: Memo | None = None,
) -> Path | None:
"""도엽 하나의 흐름누적(창고) — 도엽 + 여유 200 m 지형(창고에 있는 도엽만 · 없는 이웃 =
출구)으로 셈 · 도엽 칸만 저장(uint16 · `ACC_CAP` 로 자름) · 열쇠 = 쓴 도엽 지형 열쇠 ·
지형이 없으면 None."""
own = _dem_key(sheet_no, epsg, cell_m, memo)
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
ix0, iy0, rows, cols, box = _sheet_frame(dem_npz, cell_m)
pad = math.ceil(ZERO_ORDER_MARGIN_M / cell_m)
used = {}
for no in terrain_store.sheets_for_box(box, epsg):
dem = _dem_key(no, epsg, cell_m)
dem = _dem_key(no, epsg, cell_m, memo)
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}
key = {"version": ACC_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
if not compute:
raise NeedsCompute(sheet_no)
with terrain_store._Lock(lock_path):
if npz.is_file() and _read_json(key_path) == key: # 기다리는 동안 남이 셌을 수 있음
return npz
@@ -222,30 +298,10 @@ def sheet_accumulation(sheet_no: str, epsg: int, cell_m: float = ZERO_ORDER_CELL
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."""
"""도엽 누적을 창고 격자(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:
@@ -259,69 +315,251 @@ def paste_accumulation(
return acc
def _box_lines(
acc_paths: list[Path],
epsg: int,
box: tuple[float, ...],
cell_m: float,
min_m2: float,
max_m2: float,
) -> list[Any]:
"""누적 이어 붙여 상자 안 창고 격자에서 선 뽑기 → 단순화한 shapely 선(작업 좌표계 m)."""
import shapely
sheets = [p.parent.name for p in acc_paths]
grid = terrain_store.load_terrain(epsg, sheets=sheets, box=box, cell_m=cell_m)
x, y = grid["x"], grid["y"]
if not len(x) or not len(y):
return []
acc = paste_accumulation(acc_paths, x, y, cell_m)
lines = mask_lines(acc, grid["z"], grid["valid_mask"], x, y, cell_m, min_m2, max_m2)
if not lines:
return []
index = np.repeat(np.arange(len(lines)), [len(line) for line in lines])
geoms = shapely.linestrings(np.concatenate(lines), indices=index)
return list(shapely.simplify(geoms, SIMPLIFY_M))
def _owned_lines(lines: list[Any], sheet_no: str, epsg: int) -> list[Any]:
"""끝점 둘 중 작은 쪽(x, y 차례)이 그 도엽(경위도) 안인 선 — 이웃과 겹쳐 뽑힌 선은 한 곳만."""
import shapely
from pyproj import Transformer
from B04_PreProcess.B04_PreProcess_Engine_MapSheet import sheet5k_to_bounds
if not lines:
return []
first = shapely.get_coordinates(shapely.get_point(lines, 0))
last = shapely.get_coordinates(shapely.get_point(lines, -1))
swap = (last[:, 0] < first[:, 0]) | ((last[:, 0] == first[:, 0]) & (last[:, 1] < first[:, 1]))
ends = np.where(swap[:, None], last, first)
lon, lat = Transformer.from_crs(
terrain_store._crs(epsg), "EPSG:4326", always_xy=True
).transform(ends[:, 0], ends[:, 1])
lon0, lat0, lon1, lat1 = sheet5k_to_bounds(sheet_no)
keep = (lon >= lon0) & (lon < lon1) & (lat >= lat0) & (lat < lat1)
return [g for g, k in zip(lines, keep) if k]
def sheet_lines(
sheet_no: str,
epsg: int,
cell_m: float = ZERO_ORDER_CELL_M,
compute: bool = True,
memo: Memo | None = None,
) -> Path | None:
"""도엽 하나의 기본 임계 0차곡 선(창고 · 616-2) — 도엽 + 여유 200 m 에서 뽑아(창고에 있는 이웃
누적만) 끝점이 그 도엽 안인 선만 · 단순화 끝 · 열쇠 = 쓴 누적 지문 · 지형 없으면 None."""
import shapely
own = sheet_accumulation(sheet_no, epsg, cell_m, compute, memo)
if own is None:
return None
box = _sheet_frame(own, cell_m)[4]
used = {}
for no in terrain_store.sheets_for_box(box, epsg):
path = own if no == sheet_no else _acc_names(no, epsg, cell_m)[0]
if (stamp := file_stamp(path)) is not None:
used[no] = stamp
key = {
"version": LINES_VERSION,
"min_m2": ZERO_ORDER_MIN_M2,
"max_m2": ZERO_ORDER_MAX_M2,
"min_len_m": ZERO_ORDER_MIN_LEN_M,
"simplify_m": SIMPLIFY_M,
"margin_m": ZERO_ORDER_MARGIN_M,
"acc": used,
}
npz, key_path, lock_path = _lines_names(sheet_no, epsg, cell_m)
if npz.is_file() and _read_json(key_path) == key:
return npz
if not compute:
raise NeedsCompute(sheet_no)
with terrain_store._Lock(lock_path):
if npz.is_file() and _read_json(key_path) == key: # 기다리는 동안 남이 셌을 수 있음
return npz
paths = [_acc_names(no, epsg, cell_m)[0] for no in sorted(used)]
lines = _box_lines(paths, epsg, box, cell_m, ZERO_ORDER_MIN_M2, ZERO_ORDER_MAX_M2)
mine = _owned_lines(lines, sheet_no, epsg)
coords, index = shapely.get_coordinates(mine, return_index=True)
key_path.unlink(missing_ok=True) # 쓰는 도중 끊기면 옛 열쇠가 새 파일을 덮지 않게
atomic_write_npz(npz, coords=coords, index=index.astype(np.int32))
atomic_write_json(key_path, key)
return npz
def _read_lines(path: Path) -> list[Any]:
import shapely
with np.load(path) as data:
coords, index = data["coords"], data["index"]
return list(shapely.linestrings(coords, indices=index)) if len(coords) else []
def _project_box(npz: Path) -> tuple[float, float, float, float]:
"""사업 DEM 범위(칸 모서리) — 창고 격자에서 이 상자 안 칸을 씀."""
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])
half = cell / 2
return (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half)
def _collection(shapes: list[Any], epsg: int) -> dict[str, Any]:
"""작업 좌표계 선 → 경위도 GeoJSON(좌표 변환 한 번)."""
import shapely
from pyproj import Transformer
features: list[dict[str, Any]] = []
if shapes:
coords, index = shapely.get_coordinates(shapes, return_index=True)
lon, lat = Transformer.from_crs(
terrain_store._crs(epsg), "EPSG:4326", always_xy=True
).transform(coords[:, 0], coords[:, 1])
points = np.column_stack([lon, lat])
for part in np.split(points, np.flatnonzero(np.diff(index)) + 1):
geometry = {"type": "LineString", "coordinates": part.tolist()}
features.append({"type": "Feature", "geometry": geometry, "properties": {}})
return {"type": "FeatureCollection", "features": features}
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
compute: bool = True,
) -> dict[str, Any] | None:
"""0차곡 겹 — 받은 범위의 창고 도엽 선(기본 임계) 또는 누적(다른 임계 · 지문 같으면 저장본) ·
경위도 선 + 가닥 · 길이 합 m. `compute=False` = 셈이 필요하면 None(일꾼으로 넘김)."""
npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz"
stamp = file_stamp(npz)
if stamp is None:
if file_stamp(npz) 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,
npz.name,
stamp,
int(epsg),
float(min_m2),
float(max_m2),
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:
z = np.asarray(dem["z"], float)[::step, ::step]
valid = np.asarray(dem["valid_mask"], bool)[::step, ::step]
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_m,
"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
],
},
}
box = _project_box(npz)
cell_m = ZERO_ORDER_CELL_M
default = (float(min_m2), float(max_m2)) == (ZERO_ORDER_MIN_M2, ZERO_ORDER_MAX_M2)
memo: Memo = {}
try:
paths: dict[str, Path] = {}
for no in terrain_store.sheets_for_box(box, epsg):
if _has_terrain(no, epsg, cell_m, compute, memo):
if (got := sheet_accumulation(no, epsg, cell_m, compute, memo)) is not None:
paths[no] = got
if default: # 누적을 모두 갖춘 뒤 선(이웃 누적이 선 열쇠에 듦)
for no in list(paths):
paths[no] = sheet_lines(
no, epsg, cell_m, compute, memo
) # 창고 도엽 선 = 저장본 · 이어 붙여 자르기만(파생 저장본 안 씀)
return _layer_value(_clipped_lines(list(paths.values()), box), epsg, min_m2, max_m2)
key = derived_key(
ZERO_ORDER_KIND,
ZERO_ORDER_VERSION,
[round(v, 3) for v in box],
int(epsg),
float(min_m2),
float(max_m2),
cell_m,
ZERO_ORDER_MIN_LEN_M,
SIMPLIFY_M,
[f"{p.name}:{no}:{file_stamp(p)}" for no, p in sorted(paths.items())],
)
saved = read_derived(project_root, ZERO_ORDER_KIND, key)
if saved is not None:
return saved["value"]
if not compute:
raise NeedsCompute("임계")
except NeedsCompute:
return None
shapes = _box_lines([paths[no] for no in sorted(paths)], epsg, box, cell_m, min_m2, max_m2)
value = _layer_value(shapes, epsg, min_m2, max_m2)
write_derived(project_root, ZERO_ORDER_KIND, key, value)
return value
def _clipped_lines(paths: list[Path], box: tuple[float, ...]) -> list[Any]:
"""도엽 선 이어 붙여 받은 범위로 자름 — 잘려 짧아진 토막(`ZERO_ORDER_MIN_LEN_M` 미만) 거름."""
import shapely
lines = [g for p in paths for g in _read_lines(p)]
if not lines:
return []
parts = shapely.get_parts(shapely.clip_by_rect(lines, *box))
keep = (shapely.get_type_id(parts) == 1) & (shapely.length(parts) >= ZERO_ORDER_MIN_LEN_M)
return list(parts[keep])
def _layer_value(shapes: list[Any], epsg: int, min_m2: float, max_m2: float) -> dict[str, Any]:
import shapely
return {
"min_m2": float(min_m2),
"max_m2": float(max_m2),
"cell_m": ZERO_ORDER_CELL_M,
"count": len(shapes),
"length_m": round(float(shapely.length(shapes).sum()), 1) if shapes else 0.0,
"collection": _collection(shapes, epsg),
}
# ── 무거운 셈 일꾼(616-3) ─────────────────────────────────────────────
_POOL: ProcessPoolExecutor | None = None
def _watch_parent(pid: int) -> None:
"""일꾼 시작 — 서버 프로세스가 끝나면(다시 켬 · 강제 종료) 따라 끝남(고아 일꾼 막음)."""
import psutil
parent = psutil.Process(pid)
def loop() -> None:
while parent.is_running():
time.sleep(2.0)
os._exit(0)
threading.Thread(target=loop, daemon=True).start()
def _warm() -> None:
"""일꾼 띄울 때 무거운 모듈을 미리 읽음(첫 셈이 읽기를 안 기다리게)."""
import shapely # noqa: F401
import skimage.morphology # noqa: F401
from B05_Profile import B05_Profile_Engine_Skeleton # noqa: F401
def heavy_pool() -> ProcessPoolExecutor:
"""서버 본체 밖 일꾼 하나 — 처음 부를 때 띄워 둠 · 죽었으면 새로(이벤트 루프 실에서만 부름)."""
global _POOL
if _POOL is None or getattr(_POOL, "_broken", False):
_POOL = ProcessPoolExecutor(1, initializer=_watch_parent, initargs=(os.getpid(),))
# ponytail: spawn 일꾼이 main.py(라우터 전부 · 약 10 s)를 다시 읽지 않게 띄우는 순간만
# 빈 __main__ 으로 바꿈 — 일꾼은 이 모듈만 읽음 · 파이썬 판이 바뀌어 안 먹히면 느려질 뿐
main = sys.modules["__main__"]
sys.modules["__main__"] = types.ModuleType("__main__")
try:
_POOL.submit(_warm)
finally:
sys.modules["__main__"] = main
return _POOL
@@ -182,31 +182,6 @@ def test_zero_order_lines_head_only() -> None:
assert sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= zo.ZERO_ORDER_MIN_LEN_M
def test_zero_order_layer_saved_and_threshold_key(tmp_path, monkeypatch) -> None:
n = 60
models = tmp_path / "models"
models.mkdir()
np.savez(
models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz",
x=230000.0 + np.arange(n) * 2.0,
y=450000.0 + np.arange(n) * 2.0,
z=_valley(n),
valid_mask=np.ones((n, n), bool),
resolution=np.array([2.0, 2.0]),
)
calls = []
real = zo.zero_order_lines
monkeypatch.setattr(zo, "zero_order_lines", lambda *a, **k: calls.append(1) or real(*a, **k))
first = zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 2000.0)
assert first["count"] >= 1 and first["length_m"] > 0 and first["cell_m"] == 2.0
assert first["collection"]["features"][0]["geometry"]["type"] == "LineString"
again = zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 2000.0)
assert (again["count"], again["length_m"]) == (first["count"], first["length_m"])
assert len(calls) == 1 # 두 번째 = 저장본
zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 3000.0) # 임계 바뀜 = 다시 셈
assert len(calls) == 2
def test_zero_order_no_dem(tmp_path) -> None:
with pytest.raises(FileNotFoundError):
zo.zero_order_layer(tmp_path, tmp_path, EPSG)
@@ -1,5 +1,5 @@
"""D 0차곡 흐름누적 = 도엽별 공용 창고(PLAN 608-2) — 도엽 + 여유 200 m 조각 = 통째 셈 ·
처음만 셈 · 이웃 도엽 지형이 새로 생기면 다시 · 임계 바꿈은 누적 그대로 · 창고 격자 아닌 DEM = 통째.
처음만 셈 · 이웃 도엽 지형이 새로 생기면 다시 · 임계 바꿈은 누적 그대로(사업 DEM 격자 무관 = 616).
./venv/Scripts/python.exe -m pytest resources/tester/test_608_plan_zero_order_store.py -q
"""
@@ -77,7 +77,7 @@ def test_tiles_equal_whole_and_saved(shared) -> None:
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)
@@ -100,5 +100,3 @@ def test_layer_uses_store_and_threshold_reuses(shared, tmp_path, monkeypatch) ->
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
@@ -0,0 +1,189 @@
"""PLAN 616 — 0차곡 첫 계산 빠르게: 창고 격자에서 셈(616-1) · 도엽 선 미리 뽑아 이어 붙임(616-2) ·
무거운 셈 = 서버 밖 일꾼(616-3) · 급경사 값 그림 바이너리(616-4).
./venv/Scripts/python.exe -m pytest resources/tester/test_616_zero_order_fast.py -q
"""
from __future__ import annotations
import base64
import os
import sys
from pathlib import Path
import numpy as np
import pytest
import shapely
ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))
from common_util import common_util_terrain_store as store # noqa: E402
from common_util.common_util_atomic import atomic_write_npz # noqa: E402
from common_util.common_util_json import atomic_write_json # noqa: E402
from D03_PlanAnalysis import D03_PlanAnalysis_Layers as layers # noqa: E402
from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zo # noqa: E402
from tools import precompute_terrain as tool # noqa: E402
EPSG = 5187
WEST, EAST = "36901074", "36901075"
CELL = 10.0 # 시험은 굵은 칸(여유 200 m = 20 칸 · 기본 임계 = 10 ~ 100 칸)
def _height(x, y):
"""기운 면 + 골 · 마루 — 웅덩이 없음 · 작은 집수."""
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):
"""임시 창고 · 칸 10 m · d8 부름 셈."""
monkeypatch.setattr(store, "TERRAIN_ROOT", tmp_path / "terrain")
monkeypatch.setattr(zo, "ZERO_ORDER_CELL_M", CELL)
calls = []
real = zo.d8_accumulation
monkeypatch.setattr(zo, "d8_accumulation", lambda *a: calls.append(1) or real(*a))
return calls
def _models(tmp_path: Path, shift: float = 0.0) -> Path:
"""두 도엽 가운데를 덮는 사업 DEM(칸 = 시험 칸 · 값은 안 씀 · 범위만)."""
grid = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL)
x, y = grid["x"][30:-30] + shift, grid["y"][30:-30]
models = tmp_path / f"models{shift:g}"
models.mkdir()
np.savez(
models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz",
x=x,
y=y,
z=np.zeros((len(y), len(x))),
valid_mask=np.ones((len(y), len(x)), bool),
resolution=np.array([CELL], np.float32),
)
return models
def test_sheet_lines_join_equals_box(shared, tmp_path) -> None:
"""도엽 선(여유 200 m 에서 뽑아 끝점 도엽만) 이어 붙임 = 두 도엽 통째 상자에서 뽑은 선."""
_put_sheet(WEST)
_put_sheet(EAST)
paths = [zo.sheet_lines(no, EPSG, CELL) for no in (WEST, EAST)]
joined = [g for p in paths for g in zo._read_lines(p)]
grid = store.load_terrain(EPSG, sheets=[WEST, EAST], cell_m=CELL)
half = CELL / 2
box = (grid["x"][0] - half, grid["y"][0] - half, grid["x"][-1] + half, grid["y"][-1] + half)
accs = [zo._acc_names(no, EPSG, CELL)[0] for no in (WEST, EAST)]
whole = zo._box_lines(accs, EPSG, box, CELL, zo.ZERO_ORDER_MIN_M2, zo.ZERO_ORDER_MAX_M2)
assert len(whole) >= 20
# 여유 200 m 보다 긴 선(시험 지형 = 동서로 긴 골)만 창 끝에서 잘려 다를 수 있음 — 1 % 안
assert abs(len(joined) - len(whole)) <= 0.01 * len(whole)
total = sum(g.length for g in whole)
assert abs(sum(g.length for g in joined) - total) <= 0.01 * total
keys = {tuple(map(tuple, shapely.get_coordinates(g).round(3))) for g in joined}
assert len(keys) == len(joined) # 겹쳐 뽑힌 선은 한 도엽에만
def test_layer_default_reads_lines_light(shared, tmp_path) -> None:
"""기본 임계 = 창고 도엽 선 이어 붙여 자름 — 창고에 없으면 compute=False 가 None(일꾼 몫) ·
있으면 d8 · 선 뽑기 없이 바로."""
calls = shared
_put_sheet(WEST)
_put_sheet(EAST)
models = _models(tmp_path)
assert zo.zero_order_layer(tmp_path, models, EPSG, compute=False) is None
assert calls == []
first = zo.zero_order_layer(tmp_path, models, EPSG)
assert first["count"] >= 10 and first["cell_m"] == CELL
assert len(calls) == 2 # 도엽 둘 누적 한 번씩
light = zo.zero_order_layer(tmp_path, models, EPSG, compute=False)
assert (light["count"], light["length_m"]) == (first["count"], first["length_m"])
assert len(calls) == 2
assert light["collection"]["features"][0]["geometry"]["type"] == "LineString"
# 받은 범위 밖 선은 잘림 — 모든 점이 범위 안
with np.load(models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz") as dem:
x0, x1 = dem["x"][0] - CELL / 2, dem["x"][-1] + CELL / 2
from pyproj import Transformer
to_m = Transformer.from_crs("EPSG:4326", f"EPSG:{EPSG}", always_xy=True)
pts = np.array(
[p for f in light["collection"]["features"] for p in f["geometry"]["coordinates"]]
)
mx, _ = to_m.transform(pts[:, 0], pts[:, 1])
assert mx.min() >= x0 - 0.01 and mx.max() <= x1 + 0.01
def test_layer_custom_threshold_saved_and_grid_free(shared, tmp_path) -> None:
"""다른 임계 = 창고 누적에서 셈(사업 DEM 격자 안 씀 · 반 칸 어긋나도 창고) · 둘째 = 저장본."""
calls = shared
_put_sheet(WEST)
_put_sheet(EAST)
models = _models(tmp_path)
cells = CELL * CELL
assert zo.zero_order_layer(tmp_path, models, EPSG, 5 * cells, 40 * cells, False) is None
first = zo.zero_order_layer(tmp_path, models, EPSG, 5 * cells, 40 * cells)
assert first["count"] >= 1 and len(calls) == 2
again = zo.zero_order_layer(tmp_path, models, EPSG, 5 * cells, 40 * cells, False)
assert again == first # 저장본
shifted = zo.zero_order_layer(
tmp_path, _models(tmp_path, CELL / 2), EPSG, 5 * cells, 40 * cells
)
assert shifted["count"] >= 1 and len(calls) == 2 # 옛 격자 갈래 없음 = 누적 다시 안 셈
def test_lines_recomputed_when_neighbor_acc_appears(shared) -> None:
_put_sheet(WEST)
_put_sheet(EAST)
path = zo.sheet_lines(WEST, EPSG, CELL) # 동쪽 누적 없음 = 빼고 셈
stamp = path.stat().st_mtime_ns
assert zo.sheet_lines(WEST, EPSG, CELL, compute=False) == path
zo.sheet_accumulation(EAST, EPSG, CELL)
with pytest.raises(zo.NeedsCompute):
zo.sheet_lines(WEST, EPSG, CELL, compute=False)
assert zo.sheet_lines(WEST, EPSG, CELL).stat().st_mtime_ns != stamp
def test_acc_key_version_unchanged() -> None:
"""누적 열쇠 판을 올리면 전국 미리 계산 누적을 통째 다시 셈 — 616 은 안 올림."""
assert zo.ACC_VERSION == 2
def test_precompute_zero_lines_step(shared, monkeypatch) -> None:
_put_sheet(WEST)
_put_sheet(EAST)
monkeypatch.setattr(tool, "CELL_M", CELL)
assert "zero_lines" in tool.STEPS and tool.STEPS.index("zero_lines") == 2
_, _, info = tool.run_step("zero_lines", WEST)
assert info["ok"] and info["fresh"] and not info["empty"]
_, _, again = tool.run_step("zero_lines", WEST)
assert again["ok"] and not again["fresh"]
def test_raster_bytes_roundtrip() -> None:
png = b"\x89PNG\r\n\x1a\nfake"
value = {"image": "data:image/png;base64," + base64.b64encode(png).decode(), "counts": [1, 2]}
got, meta = layers.raster_bytes(value)
assert got == png and meta == {"counts": [1, 2]}
def test_heavy_pool_runs_outside() -> None:
"""일꾼 = 다른 프로세스 · 띄울 때 __main__ 을 되돌려 놓음."""
main = sys.modules["__main__"]
pool = zo.heavy_pool()
try:
assert sys.modules["__main__"] is main
assert pool.submit(os.getpid).result(timeout=120) != os.getpid()
assert zo.heavy_pool() is pool
finally:
pool.shutdown()
zo._POOL = None
+19 -6
View File
@@ -2,7 +2,8 @@
"""전국 미리 계산 일괄 도구(PLAN 612-1 · 요구 493) — 공용 지형 창고에 도엽마다 미리 셈.
실행: ./venv/Scripts/python.exe tools/precompute_terrain.py --sheets 36901 --workers 4
--steps dem,zero,e07 단계(차례 고정 = 지형 → 0차곡 누적 → E07 · 뒤에 tiles)
--steps dem,zero,zero_lines,e07 단계(차례 고정 = 지형 → 0차곡 누적 → 0차곡 선 → E07 ·
뒤에 tiles · zero_lines = 두 PC zero 단계가 다 끝난 뒤 · 이웃 누적 없으면 그 이웃 빼고 셈)
--sheets 36901,36712089 도엽 번호 또는 앞자리(1:50,000 묶음) · 주면 높이 차 거름 안 씀
--box lon0,lat0,lon1,lat1 도엽 가운데가 이 경위도 상자 안
--min-relief 150 도엽 안 높이 차(m) 이상만(산지) · --limit N 앞 N 장만
@@ -14,8 +15,8 @@
- 도엽 목록 = 전국 등고선 파일 꼭짓점이 닿는 도엽(0.025°) — `_precompute_sheets.json`(전국 파일
지문이 같으면 다시 안 셈) · 차례 = 등고선 길이 큰 것(산지)부터 · 같으면 높이 차.
- 좌표계 = 도엽 가운데 경도의 TM 띠(`tm_epsg_for_lon` · 5185~5187) · 격자 2 m.
- 셈 = 공용 함수만 부름(창고 지형 · D 0차곡 누적 · E07 도엽) — 저장 열쇠가 같으면 함수가 건너뜀 ·
쓰기 = 함수 쪽 원자 쓰기(끊겨도 반쯤 쓴 파일 없음).
- 셈 = 공용 함수만 부름(창고 지형 · D 0차곡 누적 · 선 · E07 도엽) — 저장 열쇠가 같으면 함수가
건너뜀 · 쓰기 = 함수 쪽 원자 쓰기(끊겨도 반쯤 쓴 파일 없음).
- 이어 하기 = PC 별 `_precompute_progress_<컴퓨터 이름>.json`(단계 · 도엽마다 성공 · 시간 ·
시놀로지가 두 PC 사이 옮김 → 한 파일은 한 PC 만 씀) · 옛 `_precompute_progress.json` 의 성공은
읽어서 이어 받음(안 씀) · 성공한 도엽은 다시 안 부름 · 실패는 다음에 다시 · 전국 파일이 바뀌면
@@ -59,7 +60,7 @@ from common_util import common_util_terrain_store as store # noqa: E402
from common_util.common_util_crs import tm_epsg_for_lon # noqa: E402
from common_util.common_util_json import atomic_write_json # noqa: E402
STEPS = ("dem", "zero", "e07")
STEPS = ("dem", "zero", "zero_lines", "e07")
CELL_M = 2.0
PROGRESS_NAME = "_precompute_progress.json" # 옛 한 벌 기록(읽기만)
PC_NAME = platform.node() or "pc"
@@ -70,7 +71,12 @@ INDEX_CHUNK = 50_000 # 가닥 묶음(약 0.5 GB 읽기)
INDEX_STRIDE = 8 # 꼭짓점 솎음(도엽 2 km 에 비해 넉넉)
SAVE_EVERY_S = 10.0
# 단계별 창고 파일 앞머리(크기 셈)
STEP_FILES = {"dem": ("dem_",), "zero": ("zero_acc_",), "e07": ("e07_", "e07si_")}
STEP_FILES = {
"dem": ("dem_",),
"zero": ("zero_acc_",),
"zero_lines": ("zero_lines_",),
"e07": ("e07_", "e07si_"),
}
_ENVELOPE = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}
@@ -322,6 +328,13 @@ def run_step(step: str, job: str, part: Part | None = None) -> tuple[str, str, d
before = _mtime(target)
empty = zero.sheet_accumulation(sheet_no, epsg, CELL_M) is None
fresh = _mtime(target) != before
elif step == "zero_lines":
from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zero
target = zero._lines_names(sheet_no, epsg, CELL_M)[0]
before = _mtime(target)
empty = zero.sheet_lines(sheet_no, epsg, CELL_M) is None
fresh = _mtime(target) != before
elif step == "e07":
from E07_StabArea import E07_StabArea_Router as e07_router
from E07_StabArea import E07_StabArea_Sheet as e07
@@ -558,7 +571,7 @@ def summary(sheets: list[str], progress: dict[str, Any]) -> dict[str, Any]:
def main(argv: list[str] | None = None) -> int:
parser = argparse.ArgumentParser(description="공용 지형 창고 전국 미리 계산(PLAN 612)")
parser.add_argument("--steps", default="dem,zero,e07")
parser.add_argument("--steps", default="dem,zero,zero_lines,e07")
parser.add_argument("--sheets", default="")
parser.add_argument("--box", default="")
parser.add_argument("--min-relief", type=float, default=150.0)