"""D 노선분석 지도 암상 지질도 겹(PLAN 760-1 · 요구 597-3 · 598-3 · 근거 `_분석_597-3`). - 출처 = 한국지질자원연구원 지오빅데이터 오픈API(수치지질도 1:50,000 · CC BY-NC) · 키 = `.env` `GEOBIGDATA_API_KEY`(쿼리 `key` · 없으면 `MissingKeyError` · 값은 어떤 로그 · 응답 · 오류 글에도 안 냄) - 그림 = WMS `L_50K_Geology_Map` 조각 대리(`common_util_wms_tiles` 「geology」) - 도형 = WFS `geoOpen:l_50k_geology_litho_latest`(GeoJSON · 경위도로 받아 작업 좌표계로 옮김) · 한글 칸이 「?」로 깨짐 → 기호(`symnum`)마다 그 면 안 점 하나 GetFeatureInfo 로 시대 · 지층 · 대표암상 · 기호를 채움(대문자 전체 인자 · 1.1.1) - 받은 범위로 자르고 기호별 넓이 · 바위 갈래(`ROCK_GROUPS` · 색 = 화면 `--draw-map-geology-*`) · 저장 = `<받은 겹 자리>/암상지질도_bounds.json`(상자가 덮으면 다시 안 받음) - 노선 따라 구간 = `route_segments`(1 m 걸음 표본 · 참고 표시 — 토사 · 암 판정 아님) 계산 자리 = ③ 서버 단독(브라우저는 KIGAM 을 부르지 않음). """ from __future__ import annotations import json import logging import os import re import time import urllib.parse import urllib.request from pathlib import Path from typing import Any import numpy as np from common_util.common_util_json import atomic_write_json logger = logging.getLogger(__name__) Box = tuple[float, float, float, float] KEY_ENV = "GEOBIGDATA_API_KEY" WMS_URL = "https://data.kigam.re.kr/openapi/wms" WFS_URL = "https://data.kigam.re.kr/openapi/wfs" MAP_LAYER = "L_50K_Geology_Map" LITHO_TYPE = "geoOpen:l_50k_geology_litho_latest" FILE = "암상지질도_bounds.json" VERSION = 1 MAX_FEATURES = 1000 TIMEOUT_S = 30 STEP_M = 1.0 SOURCE = "한국지질자원연구원 지오빅데이터 · 수치지질도 1:50,000 (CC BY-NC)" # 바위 갈래(차례 = 화면 `--draw-map-geology-1..6`) — 대표암상의 첫 「…암」 낱말로 고름(없으면 글 통째) ROCK_GROUPS: list[tuple[str, tuple[str, ...]]] = [ ("심성암(화강암류)", ("화강", "섬록", "반려", "섬장", "몬조", "페그마")), ("화산암 · 맥암", ("안산", "유문", "현무", "응회", "조면", "반암", "맥암", "화산", "석영맥")), ("퇴적암", ("사암", "셰일", "세일", "이암", "역암", "석회", "미사", "이질", "처트", "백운")), ( "변성암", ("편마", "편암", "규암", "천매", "대리", "혼펠", "각섬암", "점판", "미그마", "결정질"), ), ("미고결 퇴적층", ("충적", "자갈", "모래", "점토", "퇴적층", "붕적", "선상지")), ] OTHER = len(ROCK_GROUPS) # 그 밖 META = 3 # 변성암 — 글 어디든 낱말이 있으면 먼저(「화강편마암」 · 「결정질 석회암」) ROCK_GROUP_NAMES = [name for name, _ in ROCK_GROUPS] + ["그 밖"] SYMBOL = re.compile(r"^[A-Z][A-Za-z0-9-]*$") class MissingKeyError(RuntimeError): """`.env` 에 키 없음.""" def rock_group(litho: str) -> int: """대표암상 → 바위 갈래 번호(0 ~ `OTHER`).""" words = litho.split() head = next((w for w in words if w.endswith(("암", "층")) and not w.endswith("암질")), litho) if any(key in litho for key in ROCK_GROUPS[META][1]): return META for text in (head, litho): for index, (_, keys) in enumerate(ROCK_GROUPS): if any(key in text for key in keys): return index return OTHER def _get(url: str, params: dict[str, Any]) -> bytes: """KIGAM 요청 — 키 없음 = `MissingKeyError` · 실패 = `RuntimeError`(오류 글에 주소 · 키 없음).""" key = os.environ.get(KEY_ENV, "").strip() if not key: raise MissingKeyError(f"암상 지질도 키 없음 — .env {KEY_ENV}") query = urllib.parse.urlencode({"key": key, **params}) request = urllib.request.Request(f"{url}?{query}", headers={"User-Agent": "Mozilla/5.0"}) try: with urllib.request.urlopen(request, timeout=TIMEOUT_S) as response: return response.read() except OSError as error: raise RuntimeError(f"암상 지질도를 받지 못했습니다 — {type(error).__name__}") from None def _json(body: bytes) -> Any: """JSON 몸(같은 이름 칸이 겹쳐도 다 남김 — 깨진 한글 칸 「??」) · 아니면 `RuntimeError`.""" try: return json.loads(body.decode("utf-8"), object_pairs_hook=list) except ValueError: raise RuntimeError("암상 지질도를 받지 못했습니다 — JSON 이 아닌 응답") from None def _pairs_get(pairs: Any, name: str) -> Any: return next((v for k, v in pairs if k == name), None) if isinstance(pairs, list) else None def parse_litho(body: bytes) -> list[dict[str, Any]]: """WFS GeoJSON → [{symnum, symbol, geometry(경위도 GeoJSON)}] — 기호 = 깨진 칸 중 영문 기호 값.""" features = _pairs_get(_json(body), "features") or [] out = [] for feature in features: props = _pairs_get(feature, "properties") or [] geometry = _pairs_get(feature, "geometry") if not geometry: continue symnum = str(_pairs_get(props, "symnum") or "") symbol = next( ( str(v) for k, v in props if k not in ("symnum", "mapname", "mapidx") and SYMBOL.match(str(v)) ), "", ) out.append({"symnum": symnum or symbol, "symbol": symbol, "geometry": _plain(geometry)}) return out def _plain(value: Any) -> Any: """`object_pairs_hook=list` 로 읽은 값 → 보통 dict · list.""" if ( isinstance(value, list) and value and all(isinstance(p, tuple) and len(p) == 2 for p in value) ): return {k: _plain(v) for k, v in value} if isinstance(value, list): return [_plain(v) for v in value] return value def fetch_litho(bounds: dict[str, float]) -> list[dict[str, Any]]: """받은 범위(경위도 상자) 안 암상 면.""" body = _get( WFS_URL, { "service": "WFS", "version": "1.1.0", "request": "GetFeature", "typeName": LITHO_TYPE, "maxFeatures": MAX_FEATURES, "outputFormat": "application/json", "srsName": "EPSG:4326", "bbox": f"{bounds['min_lon']},{bounds['min_lat']},{bounds['max_lon']},{bounds['max_lat']},EPSG:4326", }, ) return parse_litho(body) def fetch_names(lon: float, lat: float) -> list[dict[str, Any]]: """점 하나의 암상 칸(한글 정상) — 그 점을 덮는 면들의 속성.""" d = 1e-4 body = _get( WMS_URL, { "SERVICE": "WMS", "VERSION": "1.1.1", "REQUEST": "GetFeatureInfo", "LAYERS": MAP_LAYER, "QUERY_LAYERS": MAP_LAYER, "SRS": "EPSG:4326", "BBOX": f"{lon - d},{lat - d},{lon + d},{lat + d}", "WIDTH": 11, "HEIGHT": 11, "X": 5, "Y": 5, "INFO_FORMAT": "application/json", "FORMAT": "image/png", "STYLES": "", }, ) got = _plain(_json(body)) return [f.get("properties") or {} for f in (got or {}).get("features", [])] def _names_for(symnum: str, symbol: str, lon: float, lat: float) -> dict[str, str]: """기호 하나의 이름 칸 — 못 받으면 기호만.""" try: rows = fetch_names(lon, lat) except RuntimeError as error: logger.warning("%s", error) rows = [] row = next((r for r in rows if str(r.get("symnum")) == symnum), rows[0] if rows else {}) return { "symbol": str(row.get("기호") or symbol), "litho": str(row.get("대표암상") or ""), "strata": str(row.get("지층") or ""), "age": str(row.get("시대") or ""), "sheet": str(row.get("mapname") or ""), } def build_layer( raw: list[dict[str, Any]], box: Box, epsg: int, names: Any = _names_for ) -> dict[str, Any]: """받은 면 → 받은 범위로 자른 겹(경위도 지물 · 기호별 요약 · 작업 좌표계 면).""" from pyproj import Transformer from shapely.geometry import box as shape_box from shapely.geometry import mapping, shape from shapely.ops import transform 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 frame = shape_box(*box) by_sym: dict[str, dict[str, Any]] = {} parts: list[tuple[str, Any]] = [] for item in raw: try: geom = transform(to_m, shape(item["geometry"])).buffer(0).intersection(frame) except (ValueError, TypeError, AttributeError): continue if geom.is_empty or geom.area <= 0: continue parts.append((item["symnum"], geom)) unit = by_sym.setdefault( item["symnum"], {"symnum": item["symnum"], "symbol": item["symbol"], "big": geom} ) if geom.area > unit["big"].area: unit["big"] = geom summary = [] for symnum, unit in by_sym.items(): lon, lat = to_ll(*unit["big"].representative_point().coords[0]) info = names(symnum, unit["symbol"], lon, lat) area = sum(g.area for s, g in parts if s == symnum) count = sum(1 for s, _ in parts if s == symnum) summary.append( { "symnum": symnum, **info, "group": rock_group(info["litho"]), "count": count, "area_m2": round(area, 1), } ) summary.sort(key=lambda s: -s["area_m2"]) index = {s["symnum"]: i for i, s in enumerate(summary)} features = [] for symnum, geom in parts: unit = summary[index[symnum]] features.append( { "type": "Feature", "properties": { "i": index[symnum], "symbol": unit["symbol"], "litho": unit["litho"], "group": unit["group"], "area_m2": round(geom.area, 1), }, "geometry": mapping(transform(to_ll, geom.simplify(1.0))), } ) return { "version": VERSION, "box_m": [float(v) for v in box], "source": SOURCE, "summary": summary, "features": features, "parts_m": [[index[s], mapping(g)] for s, g in parts], } def ensure_geology(processed: Path, box: Box, epsg: int) -> dict[str, Any]: """암상 겹 — 저장 상자가 덮으면 그대로 · 아니면 받아 지음.""" from D03_PlanAnalysis.D03_PlanAnalysis_Layers import _covers, _lonlat_bounds path = processed / FILE try: saved = json.loads(path.read_text(encoding="utf-8")) except (OSError, ValueError): saved = None if ( isinstance(saved, dict) and saved.get("version") == VERSION and _covers(saved.get("box_m"), box) ): return saved started = time.perf_counter() raw = fetch_litho(_lonlat_bounds(box, epsg)) value = build_layer(raw, box, epsg) processed.mkdir(parents=True, exist_ok=True) atomic_write_json(path, value) logger.info( "D 암상 지질도 받기: 면 %d · 기호 %d · %.2f s", len(raw), len(value["summary"]), time.perf_counter() - started, ) return value def route_segments( parts_m: list[list[Any]], points: list[list[float]], step: float = STEP_M ) -> list[dict[str, Any]]: """노선(작업 좌표계 꺾은선) 따라 암상 구간 — [{from_m, to_m, i}] · i = 요약 차례 · -1 = 자료 없음. ponytail: 1 m 걸음 표본 — 경계 오차 ≤ 1 m(1:50,000 지질 경계 어긋남 수십 m 보다 작음).""" from shapely import contains_xy from shapely.geometry import LineString, shape if len(points) < 2: return [] line = LineString([(float(p[0]), float(p[1])) for p in points]) length = line.length if length <= 0: return [] count = max(1, int(np.ceil(length / step))) edges = np.minimum(np.arange(count + 1) * step, length) mids = (edges[:-1] + edges[1:]) / 2 xy = np.array([line.interpolate(float(m)).coords[0] for m in mids]) owner = np.full(count, -1) for i, geometry in parts_m: hit = (owner < 0) & contains_xy(shape(geometry), xy[:, 0], xy[:, 1]) owner[hit] = int(i) segments: list[dict[str, Any]] = [] for k in range(count): if segments and segments[-1]["i"] == owner[k]: segments[-1]["to_m"] = round(float(edges[k + 1]), 1) else: segments.append( { "from_m": round(float(edges[k]), 1), "to_m": round(float(edges[k + 1]), 1), "i": int(owner[k]), } ) return segments def project_routes(root: Path) -> list[dict[str, Any]]: """프로젝트 노선 — 예정노선 = 세트마다 이은 선 · 그 밖 = 계획노선(수정본 → 초기 폴리라인).""" from common_util.common_util_route_geometry import ( planned_route_initial_path, planned_route_working_path, read_planned_route_csv, ) from D03_PlanAnalysis import D03_PlanAnalysis_Store as store plan = store.read_plan(root) if plan and plan.get("lines"): try: _, units = store.plan_parts(plan) except Exception as error: # 그물 오류 — 구간 표만 빠짐 logger.warning("암상 구간 — 예정노선을 못 읽음: %s", error) return [] return [ {"name": f"{u['line_name']} {u['id']}", "points": u["planned"]} for u in units if len(u["planned"]) >= 2 ] for path in (planned_route_working_path(root), planned_route_initial_path(root)): if path.is_file(): route = read_planned_route_csv(path) if route is not None and len(route.vertices) >= 2: return [{"name": "계획노선", "points": [[v.x, v.y] for v in route.vertices]}] return [] def cover_box(box: Box, routes: list[dict[str, Any]], pad: float = 100.0) -> Box: """받을 상자 = 받은 범위 ∪ 노선(+ 여유) — 노선 끝이 지형 범위 밖이어도 구간이 빈칸 안 됨.""" xs = [float(p[0]) for r in routes for p in r["points"]] ys = [float(p[1]) for r in routes for p in r["points"]] if not xs: return box return ( min(box[0], min(xs) - pad), min(box[1], min(ys) - pad), max(box[2], max(xs) + pad), max(box[3], max(ys) + pad), ) def layer_response(value: dict[str, Any], routes: list[dict[str, Any]]) -> dict[str, Any]: """화면 답 — 경위도 지물 · 요약 · 노선마다 구간(작업 좌표계 면은 안 내보냄).""" parts = value.get("parts_m", []) return { "source": value.get("source", SOURCE), "summary": value["summary"], "groups": ROCK_GROUP_NAMES, "type": "FeatureCollection", "features": value["features"], "routes": [ { "name": r["name"], "segments": route_segments(parts, r["points"]), } for r in routes ], }