eomsangdonandClaude Opus 5.5
42b63663fa
- 열쇠 상자를 소수 셋째 자리로 줄여 받은 범위를 못 덮는다고 봄 → 실제 값 그대로 저장 - 시험: 소수 끝자리까지 같은 상자 두 번 = 한 번만 받음 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
260 lines
11 KiB
Python
260 lines
11 KiB
Python
"""D 노선분석 지도 겹 — 도엽 겹 · 지적도 · 급경사 그림 (요구 395-2 · PLAN 569-1).
|
|
|
|
- 도엽 겹(하천중심선 · 도로중심선 …) = 지형 받기 병합본 · 없으면 지난 병합 도엽에서 그 겹만 풀어 냄
|
|
- 지적도 = 브이월드 WFS LX맵(`lt_c_landinfobasemap`) · B04 `request_vworld_wfs` 를 1 km 조각으로 ·
|
|
한 번 1,000 필지 한도 → 꽉 차면 넷으로 다시 쪼갬 · `pnu` 로 겹친 필지 거름 ·
|
|
저장 = `processed/연속지적도_bounds.geojson`(B04 `/geojson?layer=지적도` 자리) + 곁 열쇠
|
|
- 급경사 = 평활 안 한 사업 DEM(`models_plan/dtm_sheet_<방법>.npz`) · E07 `slope_rad` ·
|
|
15° · 35° · 40° 단계 PNG · 파생 저장본(열쇠 = 격자 지문)
|
|
|
|
계산 자리 = ③ 서버 단독.
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import base64
|
|
import io
|
|
import json
|
|
import logging
|
|
import math
|
|
import tempfile
|
|
from concurrent.futures import ThreadPoolExecutor
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
|
|
from B04_PreProcess.B04_PreProcess_Engine_GisVector import request_vworld_wfs
|
|
from B04_PreProcess.B04_PreProcess_Engine_SheetParser import SHEET_LAYERS, parse_and_merge_sheets
|
|
from B04_PreProcess.B04_PreProcess_Engine_SheetStore import get_project_map_sheets_dir
|
|
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__)
|
|
|
|
Box = tuple[float, float, float, float] # x0, y0, x1, y1(작업 좌표계 m)
|
|
|
|
CADASTRAL_TYPENAME = "lt_c_landinfobasemap"
|
|
CADASTRAL_FILE = "연속지적도_bounds.geojson"
|
|
CADASTRAL_KEY = "연속지적도_bounds.key.json"
|
|
TILE_M = 1000.0
|
|
WFS_LIMIT = 1000
|
|
MAX_SPLIT = 3 # 1 km → 125 m 까지 쪼갬
|
|
WFS_WORKERS = 4
|
|
|
|
SLOPE_KIND = "plan_slope"
|
|
SLOPE_VERSION = 1
|
|
# (아래 °, 위 °, 색 RGBA) — 임도망_편성 완 · 중 · 급경사지 · 타당성평가 35° · 노선_배치 40° 연속
|
|
SLOPE_BANDS = (
|
|
(15.0, 35.0, (250, 204, 21, 110)),
|
|
(35.0, 40.0, (249, 115, 22, 140)),
|
|
(40.0, 90.0, (220, 38, 38, 160)),
|
|
)
|
|
MAX_SIDE_PX = 4096 # ponytail: 한 변 넘으면 솎음 · 넓은 사업은 조각 그림으로
|
|
|
|
|
|
def _read_json(path: Path) -> Any:
|
|
try:
|
|
return json.loads(path.read_text(encoding="utf-8"))
|
|
except (OSError, ValueError):
|
|
return None
|
|
|
|
|
|
# ── 도엽 겹 ─────────────────────────────────────────────
|
|
|
|
|
|
def ensure_sheet_layer(project_root: Path, key: str) -> Path:
|
|
"""`processed/도엽_{key}.geojson` — 없으면 지난 병합 도엽에서 그 겹만 풀어 냄(0 건 = 빈 것)."""
|
|
if key not in SHEET_LAYERS:
|
|
raise KeyError(key)
|
|
processed = project_root / "B04_PreProcess" / "processed"
|
|
path = processed / f"도엽_{key}.geojson"
|
|
if path.is_file():
|
|
return path
|
|
sheets = (_read_json(processed / "sheet_merge.key.json") or {}).get("sheets") or []
|
|
if not sheets:
|
|
raise FileNotFoundError("병합한 도엽이 없습니다 — 지형 받기 뒤 보입니다.")
|
|
parse_and_merge_sheets(get_project_map_sheets_dir(project_root), sheets, processed, [key])
|
|
if not path.is_file():
|
|
atomic_write_json(path, {"type": "FeatureCollection", "features": []})
|
|
return path
|
|
|
|
|
|
# ── 지적도 ──────────────────────────────────────────────
|
|
|
|
|
|
def tiles(box: Box, size: float = TILE_M) -> list[Box]:
|
|
"""상자를 `size` 칸으로 쪼갬(끝 칸은 남은 만큼)."""
|
|
x0, y0, x1, y1 = box
|
|
nx, ny = max(1, math.ceil((x1 - x0) / size)), max(1, math.ceil((y1 - y0) / size))
|
|
return [
|
|
(x0 + i * size, y0 + j * size, min(x1, x0 + (i + 1) * size), min(y1, y0 + (j + 1) * size))
|
|
for j in range(ny)
|
|
for i in range(nx)
|
|
]
|
|
|
|
|
|
def quarter(box: Box) -> list[Box]:
|
|
x0, y0, x1, y1 = box
|
|
xm, ym = (x0 + x1) / 2, (y0 + y1) / 2
|
|
return [(x0, y0, xm, ym), (xm, y0, x1, ym), (x0, ym, xm, y1), (xm, ym, x1, y1)]
|
|
|
|
|
|
def _lonlat_bounds(box: Box, epsg: int) -> dict[str, float]:
|
|
from common_util.common_util_crs import tm_to_lonlat
|
|
|
|
x0, y0, x1, y1 = box
|
|
ring = [tm_to_lonlat(x, y, epsg) for x, y in ((x0, y0), (x1, y0), (x1, y1), (x0, y1))]
|
|
lons, lats = [p[0] for p in ring], [p[1] for p in ring]
|
|
return {"min_lon": min(lons), "min_lat": min(lats), "max_lon": max(lons), "max_lat": max(lats)}
|
|
|
|
|
|
def _wfs(box: Box, epsg: int) -> list[dict[str, Any]] | None:
|
|
"""조각 하나 — B04 WFS 받기(실패 = 파일 없음 · 키 값은 남기지 않음). None = 실패."""
|
|
with tempfile.TemporaryDirectory() as folder:
|
|
request_vworld_wfs(
|
|
CADASTRAL_TYPENAME, "tile.geojson", _lonlat_bounds(box, epsg), Path(folder)
|
|
)
|
|
data = _read_json(Path(folder) / "tile.geojson")
|
|
if not isinstance(data, dict) or not isinstance(data.get("features"), list):
|
|
return None
|
|
return data["features"]
|
|
|
|
|
|
def fetch_tile(box: Box, epsg: int, depth: int = 0) -> tuple[list[dict[str, Any]], int]:
|
|
"""(필지, 실패 조각 수) — 1,000 꽉 차면 넷으로 다시 쪼갬."""
|
|
features = _wfs(box, epsg)
|
|
if features is None:
|
|
return [], 1
|
|
if len(features) < WFS_LIMIT or depth >= MAX_SPLIT:
|
|
return features, 0
|
|
out: list[dict[str, Any]] = []
|
|
failed = 0
|
|
for part in quarter(box):
|
|
got, bad = fetch_tile(part, epsg, depth + 1)
|
|
out.extend(got)
|
|
failed += bad
|
|
return out, failed
|
|
|
|
|
|
def dedup_parcels(features: list[dict[str, Any]]) -> list[dict[str, Any]]:
|
|
"""조각 경계에 걸친 필지 — `pnu`(없으면 지물 id) 첫 것만."""
|
|
seen: set[str] = set()
|
|
out = []
|
|
for feature in features:
|
|
props = feature.get("properties") or {}
|
|
key = str(props.get("pnu") or feature.get("id") or json.dumps(feature.get("geometry")))
|
|
if key not in seen:
|
|
seen.add(key)
|
|
out.append(feature)
|
|
return out
|
|
|
|
|
|
def _covers(have: Any, box: Box) -> bool:
|
|
try:
|
|
return have[0] <= box[0] and have[1] <= box[1] and have[2] >= box[2] and have[3] >= box[3]
|
|
except (TypeError, IndexError):
|
|
return False
|
|
|
|
|
|
def ensure_cadastral(processed: Path, box: Box, epsg: int) -> dict[str, Any]:
|
|
"""지적도 — 열쇠 상자가 덮고 실패 조각 없으면 그대로 · 아니면 상자 통째로 다시 받음."""
|
|
path, key_path = processed / CADASTRAL_FILE, processed / CADASTRAL_KEY
|
|
key = _read_json(key_path) or {}
|
|
if path.is_file() and not key.get("failed") and _covers(key.get("box"), box):
|
|
return {"path": path, "fetched": False, "count": key.get("count"), "failed": 0}
|
|
parts = tiles(box)
|
|
with ThreadPoolExecutor(max_workers=WFS_WORKERS) as pool:
|
|
results = list(pool.map(lambda part: fetch_tile(part, epsg), parts))
|
|
failed = sum(bad for _, bad in results)
|
|
if failed == len(parts):
|
|
raise RuntimeError(f"지적도를 받지 못했습니다 — 브이월드 WFS 조각 {failed}개 실패")
|
|
features = dedup_parcels([f for got, _ in results for f in got])
|
|
processed.mkdir(parents=True, exist_ok=True)
|
|
atomic_write_json(path, {"type": "FeatureCollection", "features": features})
|
|
atomic_write_json(
|
|
key_path,
|
|
{
|
|
"box": [float(v) for v in box],
|
|
"epsg": int(epsg),
|
|
"count": len(features),
|
|
"failed": failed,
|
|
},
|
|
)
|
|
logger.info("D 지적도 받기: 조각 %d · 실패 %d · 필지 %d", len(parts), failed, len(features))
|
|
return {"path": path, "fetched": True, "count": len(features), "failed": failed}
|
|
|
|
|
|
# ── 급경사 ──────────────────────────────────────────────
|
|
|
|
|
|
def slope_bands_png(z: np.ndarray, valid: np.ndarray, cell_m: float) -> dict[str, Any]:
|
|
"""격자(행 = 남 → 북) → 단계 PNG(위 = 북 · 15° 미만 · 빈 칸 투명) + 단계 칸 비율."""
|
|
from PIL import Image
|
|
|
|
from E07_StabArea.E07_StabArea_Engine import slope_rad
|
|
|
|
stride = max(1, math.ceil(max(z.shape) / MAX_SIDE_PX))
|
|
grid = np.where(valid, np.asarray(z, float), np.nan)[::stride, ::stride]
|
|
deg = np.degrees(slope_rad(grid, cell_m * stride))
|
|
known = np.isfinite(deg)
|
|
index = np.zeros(deg.shape, np.uint8)
|
|
for i, (low, high, _) in enumerate(SLOPE_BANDS, start=1):
|
|
index[known & (deg >= low) & (deg < high)] = i
|
|
total = max(1, int(known.sum()))
|
|
image = Image.fromarray(np.flipud(index), mode="P")
|
|
palette = [0, 0, 0] + [c for _, _, rgba in SLOPE_BANDS for c in rgba[:3]]
|
|
image.putpalette(palette + [0] * (768 - len(palette)))
|
|
alpha = bytes([0] + [rgba[3] for _, _, rgba in SLOPE_BANDS])
|
|
buffer = io.BytesIO()
|
|
image.save(buffer, format="PNG", optimize=True, transparency=alpha)
|
|
return {
|
|
"width": image.width,
|
|
"height": image.height,
|
|
"image": "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode("ascii"),
|
|
"bands": [
|
|
{
|
|
"min_deg": low,
|
|
"max_deg": high,
|
|
"color": "#{:02x}{:02x}{:02x}".format(*rgba[:3]),
|
|
"ratio": round(int((index == i).sum()) / total, 4),
|
|
}
|
|
for i, (low, high, rgba) in enumerate(SLOPE_BANDS, start=1)
|
|
],
|
|
}
|
|
|
|
|
|
def slope_layer(project_root: Path, models_dir: Path, epsg: int) -> dict[str, Any]:
|
|
"""급경사 겹 — 지문(격자 파일 · 단계 · 판) 같으면 저장본 그대로."""
|
|
from common_util.common_util_crs import tm_to_lonlat
|
|
|
|
npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz"
|
|
stamp = file_stamp(npz)
|
|
if stamp is None:
|
|
raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.")
|
|
key = derived_key(
|
|
SLOPE_KIND, SLOPE_VERSION, npz.name, stamp, int(epsg), SLOPE_BANDS, MAX_SIDE_PX
|
|
)
|
|
saved = read_derived(project_root, SLOPE_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])
|
|
value = slope_bands_png(dem["z"], np.asarray(dem["valid_mask"], bool), cell)
|
|
half = cell / 2
|
|
box = (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half)
|
|
value["box_m"] = [float(v) for v in box]
|
|
value["box"] = [
|
|
list(tm_to_lonlat(px, py, epsg))
|
|
for px, py in ((box[0], box[1]), (box[2], box[1]), (box[2], box[3]), (box[0], box[3]))
|
|
]
|
|
write_derived(project_root, SLOPE_KIND, key, value)
|
|
return value
|