eomsangdonandClaude Opus 5.5
66f9c68a1e
- 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>
566 lines
23 KiB
Python
566 lines
23 KiB
Python
"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 608-2 · 616 · 요구 443 · 497).
|
||
|
||
0차곡 = 뚜렷한 유로가 없는 곡두의 집수지형(`노선_배치.md:31` · 교본 `3. 임도망의 배치 절차.md:156`)
|
||
→ 공용 지형 창고 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`(열쇠 = 쓴 도엽 지형 열쇠).
|
||
기본 임계 선 = 도엽마다 미리 뽑아 창고에(`zero_lines_<epsg>_<칸>m.npz` · 도엽 + 여유 200 m 에서 뽑아
|
||
끝점(작은 쪽)이 그 도엽 안인 선만 · 단순화 끝) → 요청 = 도엽 선 이어 붙여 받은 범위로 자름(616-2).
|
||
다른 임계 = 받은 범위 도엽 누적 이어 붙여 칸 고르기 · 선 뽑기(616-1 · 사업 DEM 격자 안 씀).
|
||
다른 임계만 파생 저장본(열쇠 = 받은 범위 · 임계 · 칸 · 쓴 도엽 누적 지문).
|
||
무거운 셈(창고에 없는 지형 · 누적 · 선 · 다른 임계 첫 셈) =
|
||
서버 본체 밖 일꾼 하나(616-3 · `heavy_pool`).
|
||
|
||
계산 자리 = ③ 서버 단독.
|
||
"""
|
||
|
||
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
|
||
|
||
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
|
||
|
||
logger = logging.getLogger(__name__)
|
||
|
||
ZERO_ORDER_KIND = "plan_zero_order"
|
||
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 # 흐름누적 칸
|
||
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))
|
||
|
||
|
||
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
|
||
# 평지 이어 감 — 같은 높이 이웃 표는 한 번 · 아직 안 흐르는 후보 칸만 돌림(방향 · 갱신 차례 그대로)
|
||
eq = [_shift(filled, dr, dc, np.nan) == filled for dr, dc in _D8]
|
||
cand = np.flatnonzero(~done & np.any(eq, axis=0))
|
||
done_f, valid_f, rec_f = done.ravel(), valid.ravel(), receiver.ravel()
|
||
r, c = np.divmod(cand, cols)
|
||
nbr, ok = [], []
|
||
for k, (dr, dc) in enumerate(_D8):
|
||
inb = (r + dr >= 0) & (r + dr < rows) & (c + dc >= 0) & (c + dc < cols)
|
||
ok.append(inb & eq[k].ravel()[cand])
|
||
nbr.append(np.where(inb, cand + dr * cols + dc, 0))
|
||
while True:
|
||
grew = False
|
||
for k, (dr, dc) in enumerate(_D8):
|
||
take = ok[k] & ~done_f[cand] & done_f[nbr[k]] & valid_f[nbr[k]]
|
||
if take.any():
|
||
hit = cand[take]
|
||
rec_f[hit] = hit + dr * cols + dc
|
||
done_f[hit] = True
|
||
grew = True
|
||
if not grew:
|
||
break
|
||
keep = ~done_f[cand]
|
||
cand, nbr, ok = cand[keep], [n[keep] for n in nbr], [o[keep] for o in ok]
|
||
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 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,
|
||
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차곡 선."""
|
||
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"))
|
||
except (OSError, ValueError):
|
||
return None
|
||
|
||
|
||
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"{prefix}_{int(epsg)}_{cell_m:g}m"
|
||
return folder / f"{stem}.npz", folder / f"{stem}.key.json", folder / f".{stem}.lock"
|
||
|
||
|
||
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)
|
||
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 _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"])
|
||
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,
|
||
(iy0 - pad) * cell_m,
|
||
(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, memo)
|
||
if dem is not None and not dem.get("empty"):
|
||
used[no] = dem
|
||
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
|
||
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 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 _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,
|
||
compute: bool = True,
|
||
) -> dict[str, Any] | None:
|
||
"""0차곡 겹 — 받은 범위의 창고 도엽 선(기본 임계) 또는 누적(다른 임계 · 지문 같으면 저장본) ·
|
||
경위도 선 + 가닥 · 길이 합 m. `compute=False` = 셈이 필요하면 None(일꾼으로 넘김)."""
|
||
npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz"
|
||
if file_stamp(npz) is None:
|
||
raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.")
|
||
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
|