eomsangdonandClaude Opus 5.5
edc20df849
- 서버: 15 · 35 · 40° 단계 PNG 길을 걷고 기울기 % 값 그림 한 장(칸 = % 내림 0~250 · 빈 칸 255) + 1 % 칸 수 · 파생 저장본 판 2(열쇠 = 격자 지문) - 화면: 고른 범위대로 값 그림을 칠함(서버 왕복 없음) · 범례 비율 = 칸 수로 바로 셈 · 범례 급경사 줄 ⚙ → 모달(아래 % · 위 % · 색 · ° 환산 · 범위 추가 · 지움 · 기본값 · 빈 값 · 겹침 막음) - 기본 = 지식DB 15 · 35 · 40° 의 % 26.8 · 70.0 · 83.9 · 기억 = 화면 취향 aislo:pref:d03.slope-bands(계정) - 시험: 569_1 급경사 시험을 값 그림 · 칸 수로 고침, 581_1 은 bands 대신 counts(옛 단계 응답이 없어짐) · 새 603 시험(칸 수 합 = 유효 칸 · 기본 비율 = 옛 15 · 35 · 40° 비율) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
617 lines
27 KiB
Python
617 lines
27 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` ·
|
|
기울기 % 값 그림 한 장(1 % 단계 · 회색) + 1 % 칸 수 · 범위 색칠 · 비율 = 화면(PLAN 603) ·
|
|
파생 저장본(열쇠 = 격자 지문)
|
|
- 음영기복 = 같은 DEM · B07 `hillshade_png`(PLAN 581-1)
|
|
- 보호 · 규제 구역 = 브이월드 WFS 여덟 종(겹 이름 = `_분석_440-C` 표 1) ·
|
|
받은 범위로 잘라 종류 · 이름 · 넓이만 · 저장 = `processed/보호구역_bounds.geojson` + 곁 열쇠 ·
|
|
종류 하나 실패 = 그 종류만 빠짐
|
|
- 산림입지도 = 브이월드 WFS `lt_c_fsdifrsts`(원자료 산림청 산림입지토양도 · PLAN 581-2) ·
|
|
받은 범위로 잘라 토양형 코드 · 이름 · 넓이만 ·
|
|
저장 = `processed/산림입지도_bounds.geojson` + 곁 열쇠
|
|
- 선 겹 = 브이월드 WFS 단층 `lt_l_gimsfault` · 지질구조선 `lt_l_gimslinea`(PLAN 584-1) ·
|
|
받은 범위로 잘라 가닥 · 길이만 · 저장 = `processed/<이름>_bounds.geojson` + 곁 열쇠
|
|
- 0차곡 = `D03_PlanAnalysis_ZeroOrder`
|
|
- 계류 100 m · 도로 300 m 띠 = 도엽 하천중심선(세류 뺌) · 도로중심선(미분류 · 소로 ·
|
|
면리간도로 뺌) 합친 띠 · 받은 범위로 자름 · 파생 저장본(열쇠 = 두 파일 지문)
|
|
|
|
계산 자리 = ③ 서버 단독.
|
|
"""
|
|
|
|
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"
|
|
HILLSHADE_KIND = "plan_hillshade"
|
|
HILLSHADE_VERSION = 2 # 빈 칸 투명
|
|
RASTER_KINDS = ("slope", "hillshade")
|
|
SLOPE_VERSION = 2 # 값 그림 + 1 % 칸 수(PLAN 603)
|
|
SLOPE_MAX_PCT = 250 # 칸 값 0~250 = 기울기 % 내림(250 = 250 % 이상)
|
|
SLOPE_EMPTY = 255 # 빈 칸
|
|
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, typename: str) -> list[dict[str, Any]] | None:
|
|
"""조각 하나 — B04 WFS 받기(실패 = 파일 없음 · 키 값은 남기지 않음). None = 실패."""
|
|
with tempfile.TemporaryDirectory() as folder:
|
|
request_vworld_wfs(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, typename: str = CADASTRAL_TYPENAME
|
|
) -> tuple[list[dict[str, Any]], int]:
|
|
"""(지물, 실패 조각 수) — 1,000 꽉 차면 넷으로 다시 쪼갬."""
|
|
features = _wfs(box, epsg, typename)
|
|
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, typename)
|
|
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]:
|
|
"""지적도 — 열쇠 상자가 덮고 실패 조각 없으면 그대로 · 덮는데 실패 조각이 남았으면 그 조각만
|
|
다시(PLAN 595-1) · 못 덮으면 상자 통째로 다시 받음."""
|
|
path, key_path = processed / CADASTRAL_FILE, processed / CADASTRAL_KEY
|
|
key = _read_json(key_path) or {}
|
|
covered = path.is_file() and _covers(key.get("box"), box)
|
|
if covered and not key.get("failed"):
|
|
return {"path": path, "fetched": False, "count": key.get("count"), "failed": 0}
|
|
retry = key.get("failed_tiles") if covered else None # 옛 열쇠(조각 기록 없음) = 통째로
|
|
old = ((_read_json(path) or {}).get("features") or []) if retry else []
|
|
frame = key["box"] if retry else [float(v) for v in box]
|
|
parts = [tuple(t) for t in retry] if retry else tiles(box)
|
|
with ThreadPoolExecutor(max_workers=WFS_WORKERS) as pool:
|
|
results = list(pool.map(lambda part: fetch_tile(part, epsg), parts))
|
|
bad = [list(part) for part, (_, n) in zip(parts, results) if n]
|
|
if len(bad) == len(parts) and not old:
|
|
raise RuntimeError(f"지적도를 받지 못했습니다 — 브이월드 WFS 조각 {len(bad)}개 실패")
|
|
features = dedup_parcels([*old, *(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": frame,
|
|
"epsg": int(epsg),
|
|
"count": len(features),
|
|
"failed": len(bad),
|
|
"failed_tiles": bad,
|
|
},
|
|
)
|
|
logger.info("D 지적도 받기: 조각 %d · 실패 %d · 필지 %d", len(parts), len(bad), len(features))
|
|
return {"path": path, "fetched": True, "count": len(features), "failed": len(bad)}
|
|
|
|
|
|
# ── 급경사 ──────────────────────────────────────────────
|
|
|
|
|
|
def slope_values_png(z: np.ndarray, valid: np.ndarray, cell_m: float) -> dict[str, Any]:
|
|
"""격자(행 = 남 → 북) → 기울기 % 값 그림(회색 · 위 = 북 · 칸 값 = % 내림 0~250 · 빈 칸 255)
|
|
+ 1 % 칸 수(`counts[v]` = 값 v 칸 수 · 합 = 유효 칸)."""
|
|
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]
|
|
pct = np.tan(slope_rad(grid, cell_m * stride)) * 100.0
|
|
known = np.isfinite(pct)
|
|
value = np.full(pct.shape, SLOPE_EMPTY, np.uint8)
|
|
value[known] = np.minimum(np.floor(pct[known]), SLOPE_MAX_PCT).astype(np.uint8)
|
|
image = Image.fromarray(np.flipud(value), mode="L")
|
|
buffer = io.BytesIO()
|
|
image.save(buffer, format="PNG", optimize=True)
|
|
counts = np.bincount(value[known], minlength=SLOPE_MAX_PCT + 1)
|
|
return {
|
|
"width": image.width,
|
|
"height": image.height,
|
|
"image": "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode("ascii"),
|
|
"counts": [int(n) for n in counts],
|
|
}
|
|
|
|
|
|
def hillshade_layer_png(z: np.ndarray, valid: np.ndarray, cell_m: float) -> dict[str, Any]:
|
|
"""음영기복 PNG(B07 그대로 · 빈 칸 흰색 255 → 투명 · 값 칸은 90~250)."""
|
|
from PIL import Image
|
|
|
|
from B07_DesignDetail.B07_DesignDetail_Engine_Cad_Lidar import hillshade_png
|
|
|
|
image, width, height = hillshade_png(z, valid, cell_m)
|
|
gray = Image.open(io.BytesIO(base64.b64decode(image.split(",", 1)[1]))).convert("L")
|
|
gray.putalpha(gray.point(lambda v: 0 if v == 255 else 255))
|
|
buffer = io.BytesIO()
|
|
gray.save(buffer, format="PNG", optimize=True)
|
|
data = "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode("ascii")
|
|
return {"width": width, "height": height, "image": data, "bands": []}
|
|
|
|
|
|
def slope_layer(
|
|
project_root: Path, models_dir: Path, epsg: int, kind: str = "slope"
|
|
) -> dict[str, Any]:
|
|
"""급경사(`slope`) · 음영기복(`hillshade`) — 지문(격자 파일 · 판) 같으면 저장본."""
|
|
npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz"
|
|
stamp = file_stamp(npz)
|
|
if stamp is None:
|
|
raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.")
|
|
store_kind, version, build = (
|
|
(SLOPE_KIND, SLOPE_VERSION, slope_values_png)
|
|
if kind == "slope"
|
|
else (HILLSHADE_KIND, HILLSHADE_VERSION, hillshade_layer_png)
|
|
)
|
|
key = derived_key(store_kind, version, npz.name, stamp, int(epsg), MAX_SIDE_PX)
|
|
saved = read_derived(project_root, store_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 = build(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"] = lonlat_corners(box, epsg)
|
|
write_derived(project_root, store_kind, key, value)
|
|
return value
|
|
|
|
|
|
def lonlat_corners(box: Any, epsg: int) -> list[list[float]]:
|
|
"""TM 상자 → 경위도 네 귀퉁이(왼쪽 아래부터 돎) — GL 바탕 그림 자리(586)."""
|
|
from common_util.common_util_crs import tm_to_lonlat
|
|
|
|
x0, y0, x1, y1 = (float(v) for v in box)
|
|
return [list(tm_to_lonlat(px, py, epsg)) for px, py in ((x0, y0), (x1, y0), (x1, y1), (x0, y1))]
|
|
|
|
|
|
# ── 보호 · 규제 구역 ─────────────────────────────────────
|
|
|
|
# (열쇠, 이름, 브이월드 WFS 겹 이름) — `_분석_440-C` 표 1 그대로
|
|
ZONE_KINDS: tuple[tuple[str, str, tuple[str, ...]], ...] = (
|
|
("forest", "산림보호구역", ("lt_c_uf151",)),
|
|
("baekdu", "백두대간보호지역", ("lt_c_uf901",)),
|
|
("water", "상수원보호구역", ("lt_c_um710",)),
|
|
("park", "자연공원", ("lt_c_wgisnpgug", "lt_c_wgisnpdo", "lt_c_wgisnpgun")),
|
|
("heritage", "국가유산 지정 · 보호구역", ("lt_c_uo301",)),
|
|
("wetland", "습지보호지역", ("lt_c_um901",)),
|
|
("wildlife", "야생생물보호구역", ("lt_c_um221",)),
|
|
("greenbelt", "개발제한구역", ("lt_c_ud801",)),
|
|
)
|
|
ZONE_FILE = "보호구역_bounds.geojson"
|
|
ZONE_KEY = "보호구역_bounds.key.json"
|
|
ZONE_VERSION = 1
|
|
ZONE_NAME_PROPS = ("alias", "park_name", "uname")
|
|
|
|
|
|
def _transforms(epsg: int) -> tuple[Any, Any]:
|
|
"""(경위도 → m, m → 경위도) shapely 변환 함수."""
|
|
from pyproj import Transformer
|
|
|
|
to_m = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True).transform
|
|
to_ll = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True).transform
|
|
return to_m, to_ll
|
|
|
|
|
|
def _clip_metric(feature: dict[str, Any], frame: Any, to_m: Any) -> Any:
|
|
"""지물 면(경위도) → 작업 좌표계 m 로 바꿔 범위로 자름 · 못 읽거나 빈 것 = None."""
|
|
from shapely.geometry import shape
|
|
from shapely.ops import transform
|
|
|
|
try:
|
|
metric = transform(to_m, shape(feature["geometry"]))
|
|
except (KeyError, TypeError, ValueError, AttributeError):
|
|
return None
|
|
if not metric.is_valid:
|
|
metric = metric.buffer(0)
|
|
part = metric.intersection(frame)
|
|
return None if part.is_empty or part.area <= 0 else part
|
|
|
|
|
|
def clip_zones(
|
|
found: dict[str, list[dict[str, Any]]], box: Box, epsg: int
|
|
) -> tuple[list[dict[str, Any]], dict[str, dict[str, Any]]]:
|
|
"""종류별 지물 → 받은 범위로 자른 지물(종류 · 이름) + 종류별 {이름, 개수, 넓이 m²}."""
|
|
from shapely.geometry import box as rect
|
|
from shapely.geometry import mapping
|
|
from shapely.ops import transform
|
|
|
|
to_m, to_ll = _transforms(epsg)
|
|
frame = rect(*box)
|
|
features: list[dict[str, Any]] = []
|
|
summary: dict[str, dict[str, Any]] = {}
|
|
for kind, label, _ in ZONE_KINDS:
|
|
count, area = 0, 0.0
|
|
for feature in dedup_parcels(found.get(kind) or []):
|
|
part = _clip_metric(feature, frame, to_m)
|
|
if part is None:
|
|
continue
|
|
props = feature.get("properties") or {}
|
|
name = next((str(props[k]) for k in ZONE_NAME_PROPS if props.get(k)), label)
|
|
count += 1
|
|
area += part.area
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"geometry": mapping(transform(to_ll, part)),
|
|
"properties": {"kind": kind, "name": name, "area_m2": round(part.area, 1)},
|
|
}
|
|
)
|
|
summary[kind] = {"label": label, "count": count, "area_m2": round(area, 1)}
|
|
return features, summary
|
|
|
|
|
|
def ensure_zones(processed: Path, box: Box, epsg: int) -> dict[str, Any]:
|
|
"""보호 · 규제 구역 — 열쇠 상자가 덮고 실패 종류 없으면 그대로 · 아니면 통째로 다시 받음."""
|
|
path, key_path = processed / ZONE_FILE, processed / ZONE_KEY
|
|
key = _read_json(key_path) or {}
|
|
if (
|
|
path.is_file()
|
|
and key.get("version") == ZONE_VERSION
|
|
and not key.get("failed")
|
|
and _covers(key.get("box"), box)
|
|
):
|
|
return {"path": path, "fetched": False, "failed": []}
|
|
jobs = [(kind, name) for kind, _, names in ZONE_KINDS for name in names]
|
|
# ponytail: 겹마다 범위 통째 한 번 · 1,000 꽉 차면 넷 — 구역은 필지보다 훨씬 적음
|
|
with ThreadPoolExecutor(max_workers=WFS_WORKERS) as pool:
|
|
results = list(pool.map(lambda job: fetch_tile(box, epsg, 0, job[1]), jobs))
|
|
failed = sorted({kind for (kind, _), (_, bad) in zip(jobs, results) if bad})
|
|
if len(failed) == len(ZONE_KINDS):
|
|
raise RuntimeError("보호 · 규제 구역을 받지 못했습니다 — 브이월드 WFS 실패")
|
|
found: dict[str, list[dict[str, Any]]] = {}
|
|
for (kind, _), (got, _) in zip(jobs, results):
|
|
found.setdefault(kind, []).extend(got)
|
|
features, summary = clip_zones(found, box, epsg)
|
|
for kind in failed:
|
|
summary[kind]["failed"] = True
|
|
processed.mkdir(parents=True, exist_ok=True)
|
|
atomic_write_json(path, {"type": "FeatureCollection", "features": features, "summary": summary})
|
|
atomic_write_json(
|
|
key_path,
|
|
{
|
|
"version": ZONE_VERSION,
|
|
"box": [float(v) for v in box],
|
|
"epsg": int(epsg),
|
|
"failed": failed,
|
|
},
|
|
)
|
|
logger.info("D 보호구역 받기: 지물 %d · 실패 종류 %s", len(features), failed)
|
|
return {"path": path, "fetched": True, "failed": failed}
|
|
|
|
|
|
# ── 산림입지도 ──────────────────────────────────────────
|
|
|
|
SOIL_TYPENAME = "lt_c_fsdifrsts"
|
|
SOIL_FILE = "산림입지도_bounds.geojson"
|
|
SOIL_KEY = "산림입지도_bounds.key.json"
|
|
SOIL_VERSION = 1
|
|
|
|
|
|
def clip_soils(
|
|
found: list[dict[str, Any]], box: Box, epsg: int
|
|
) -> tuple[list[dict[str, Any]], list[dict[str, Any]]]:
|
|
"""산림입지도 지물 → 받은 범위로 자른 지물(토양형 코드 · 이름 · 넓이) + 토양형별 {코드, 이름,
|
|
개수, 넓이 m²}(코드 차례). 이름 없으면 토양형 글자 그대로."""
|
|
from shapely.geometry import box as rect
|
|
from shapely.geometry import mapping
|
|
from shapely.ops import transform
|
|
|
|
to_m, to_ll = _transforms(epsg)
|
|
frame = rect(*box)
|
|
features: list[dict[str, Any]] = []
|
|
summary: dict[str, dict[str, Any]] = {}
|
|
for feature in dedup_parcels(found):
|
|
part = _clip_metric(feature, frame, to_m)
|
|
if part is None:
|
|
continue
|
|
props = feature.get("properties") or {}
|
|
code = str(props.get("toyanghyun") or props.get("sltp_cd") or "?")
|
|
name = str(props.get("name") or code)
|
|
row = summary.setdefault(
|
|
code,
|
|
{"code": code, "sort": str(props.get("sltp_cd") or code), "name": name, "count": 0},
|
|
)
|
|
row["count"] += 1
|
|
row["area_m2"] = row.get("area_m2", 0.0) + part.area
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"geometry": mapping(transform(to_ll, part)),
|
|
"properties": {"code": code, "name": name, "area_m2": round(part.area, 1)},
|
|
}
|
|
)
|
|
rows = sorted(summary.values(), key=lambda row: (row["sort"], row["code"]))
|
|
for row in rows:
|
|
del row["sort"]
|
|
row["area_m2"] = round(row["area_m2"], 1)
|
|
return features, rows
|
|
|
|
|
|
def ensure_soils(processed: Path, box: Box, epsg: int) -> dict[str, Any]:
|
|
"""산림입지도 — 열쇠 상자가 덮으면 그대로 · 아니면 범위 통째로 다시 받음 · 실패 = 오류."""
|
|
path, key_path = processed / SOIL_FILE, processed / SOIL_KEY
|
|
key = _read_json(key_path) or {}
|
|
if path.is_file() and key.get("version") == SOIL_VERSION and _covers(key.get("box"), box):
|
|
return {"path": path, "fetched": False}
|
|
got, failed = fetch_tile(box, epsg, 0, SOIL_TYPENAME)
|
|
if failed:
|
|
raise RuntimeError("산림입지도를 받지 못했습니다 — 브이월드 WFS 실패")
|
|
features, summary = clip_soils(got, box, epsg)
|
|
processed.mkdir(parents=True, exist_ok=True)
|
|
atomic_write_json(path, {"type": "FeatureCollection", "features": features, "summary": summary})
|
|
atomic_write_json(
|
|
key_path, {"version": SOIL_VERSION, "box": [float(v) for v in box], "epsg": int(epsg)}
|
|
)
|
|
logger.info("D 산림입지도 받기: 지물 %d · 토양형 %d", len(features), len(summary))
|
|
return {"path": path, "fetched": True}
|
|
|
|
|
|
# ── 선 겹(단층 · 지질구조선) ─────────────────────────────
|
|
|
|
# 열쇠 → (이름, 브이월드 WFS 겹 이름) — 새 선 겹은 여기 한 줄(PLAN 584-1)
|
|
LINE_KINDS: dict[str, tuple[str, str]] = {
|
|
"faults": ("단층", "lt_l_gimsfault"),
|
|
"lineaments": ("지질구조선", "lt_l_gimslinea"),
|
|
}
|
|
LINE_VERSION = 1
|
|
|
|
|
|
def clip_lines(
|
|
found: list[dict[str, Any]], box: Box, epsg: int
|
|
) -> tuple[list[dict[str, Any]], dict[str, Any]]:
|
|
"""선 지물 → 받은 범위로 자른 선(길이 m · 원 속성 `legend`) + {가닥, 길이 m}."""
|
|
from shapely.geometry import box as rect
|
|
from shapely.geometry import mapping, shape
|
|
from shapely.ops import transform
|
|
|
|
to_m, to_ll = _transforms(epsg)
|
|
frame = rect(*box)
|
|
features: list[dict[str, Any]] = []
|
|
total = 0.0
|
|
for feature in dedup_parcels(found):
|
|
try:
|
|
part = transform(to_m, shape(feature["geometry"])).intersection(frame)
|
|
except (KeyError, TypeError, ValueError, AttributeError):
|
|
continue
|
|
if part.is_empty or part.length <= 0:
|
|
continue
|
|
total += part.length
|
|
props = feature.get("properties") or {}
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"geometry": mapping(transform(to_ll, part)),
|
|
"properties": {"legend": props.get("legend"), "length_m": round(part.length, 1)},
|
|
}
|
|
)
|
|
return features, {"count": len(features), "length_m": round(total, 1)}
|
|
|
|
|
|
def ensure_lines(processed: Path, box: Box, epsg: int, kind: str) -> dict[str, Any]:
|
|
"""선 겹 — 열쇠 상자가 덮으면 그대로 · 아니면 범위 통째로 다시 받음 · 실패 = 오류."""
|
|
label, typename = LINE_KINDS[kind]
|
|
path, key_path = processed / f"{label}_bounds.geojson", processed / f"{label}_bounds.key.json"
|
|
key = _read_json(key_path) or {}
|
|
if path.is_file() and key.get("version") == LINE_VERSION and _covers(key.get("box"), box):
|
|
return {"path": path, "fetched": False}
|
|
got, failed = fetch_tile(box, epsg, 0, typename)
|
|
if failed:
|
|
raise RuntimeError(f"{label}을 받지 못했습니다 — 브이월드 WFS 실패")
|
|
features, summary = clip_lines(got, box, epsg)
|
|
processed.mkdir(parents=True, exist_ok=True)
|
|
atomic_write_json(path, {"type": "FeatureCollection", "features": features, "summary": summary})
|
|
atomic_write_json(
|
|
key_path, {"version": LINE_VERSION, "box": [float(v) for v in box], "epsg": int(epsg)}
|
|
)
|
|
logger.info("D %s 받기: 선 %d · %.0f m", label, summary["count"], summary["length_m"])
|
|
return {"path": path, "fetched": True}
|
|
|
|
|
|
# ── 계류 · 도로 띠 ───────────────────────────────────────
|
|
|
|
BAND_KIND = "plan_band"
|
|
BAND_VERSION = 1
|
|
# 타당성평가 「산지계류 100 m 이내」 · 「도로법 도로 300 m 이내」
|
|
STREAM_BAND_M = 100.0
|
|
ROAD_BAND_M = 300.0
|
|
# ponytail: 원문 = 1/25,000 주요 수계 · 도로법 도로 — 1/5,000 도엽에서 세류 · 소로급을 뺀
|
|
# 어림(협의 값)
|
|
SMALL_STREAMS = frozenset({"세류"})
|
|
LOCAL_ROADS = frozenset({"미분류", "소로", "면리간도로"})
|
|
BAND_SEGMENTS = 8 # 둥근 끝 사분원 꼭짓점
|
|
|
|
|
|
def band_union(lines: list[Any], dist_m: float, box: Box) -> Any:
|
|
"""선(작업 좌표계 m) → 거리 `dist_m` 띠 하나(겹친 곳 합침 · 범위로 자름) · 선 없으면 None."""
|
|
from shapely import union_all
|
|
from shapely.geometry import box as rect
|
|
|
|
if not lines:
|
|
return None
|
|
band = union_all(lines).buffer(dist_m, quad_segs=BAND_SEGMENTS).intersection(rect(*box))
|
|
return None if band.is_empty else band
|
|
|
|
|
|
def _band(
|
|
path: Path, epsg: int, box: Box, dist_m: float, prop: str, skip: frozenset[str]
|
|
) -> dict[str, Any]:
|
|
from shapely.geometry import mapping, shape
|
|
from shapely.ops import transform
|
|
|
|
to_m, to_ll = _transforms(epsg)
|
|
features = (_read_json(path) or {}).get("features") or []
|
|
lines = [
|
|
transform(to_m, shape(f["geometry"]))
|
|
for f in features
|
|
if f.get("geometry") and (f.get("properties") or {}).get(prop) not in skip
|
|
]
|
|
band = band_union(lines, dist_m, box)
|
|
frame = (box[2] - box[0]) * (box[3] - box[1])
|
|
area = band.area if band is not None else 0.0
|
|
shown = [] if band is None else [band.simplify(0.5)]
|
|
return {
|
|
"dist_m": dist_m,
|
|
"lines": len(lines),
|
|
"area_m2": round(area, 1),
|
|
"ratio": round(area / frame, 4) if frame > 0 else 0.0,
|
|
"collection": {
|
|
"type": "FeatureCollection",
|
|
"features": [
|
|
{"type": "Feature", "geometry": mapping(transform(to_ll, g)), "properties": {}}
|
|
for g in shown
|
|
],
|
|
},
|
|
}
|
|
|
|
|
|
def band_layer(project_root: Path, epsg: int, box: Box) -> dict[str, Any]:
|
|
"""계류 100 m · 도로 300 m 띠 — 두 도엽 겹 지문 · 범위 같으면 저장본 그대로."""
|
|
stream = ensure_sheet_layer(project_root, "하천중심선")
|
|
road = ensure_sheet_layer(project_root, "도로중심선")
|
|
key = derived_key(
|
|
BAND_KIND,
|
|
BAND_VERSION,
|
|
file_stamp(stream),
|
|
file_stamp(road),
|
|
int(epsg),
|
|
[round(v, 1) for v in box],
|
|
STREAM_BAND_M,
|
|
ROAD_BAND_M,
|
|
sorted(SMALL_STREAMS),
|
|
sorted(LOCAL_ROADS),
|
|
)
|
|
saved = read_derived(project_root, BAND_KIND, key)
|
|
if saved is not None:
|
|
return saved["value"]
|
|
value = {
|
|
"stream": _band(stream, epsg, box, STREAM_BAND_M, "구분", SMALL_STREAMS),
|
|
"road": _band(road, epsg, box, ROAD_BAND_M, "도로구분", LOCAL_ROADS),
|
|
}
|
|
write_derived(project_root, BAND_KIND, key, value)
|
|
return value
|