"""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