eomsangdonandClaude Opus 5.5
f6cfaf446c
- 서버: 지오빅데이터(KIGAM) 수치지질도 1:50,000 — WFS 암상 면(깨진 한글 칸은 기호별 점 조회로 채움) · 받은 범위 ∪ 노선으로 자름 · 저장본 · 노선 따라 1 m 표본 구간 · 키는 서버만(.env GEOBIGDATA_API_KEY · 쿼리 key) - 그림 조각 대리에 geology(키 쿼리 이름 칸 key_param · 옛 겹 serviceKey 그대로) - 화면(D03 · B05): 「제한·재해」에 「암상 지질도」 — 그림 조각 + 갈래 색 면 · 「기호 · 암상」 글 · 범례 · 노선 구간 표(참고 — 토사 · 암 판정 아님) - 켠 밖 겹마다 범례 맨 아래 출처 줄(브이월드 · 국토지리정보원 · 산림청 · 환경부 · 생활안전지도 공공누리 · 한국지질자원연구원 CC BY-NC) - 범례 칸 높이 한도 · 밀기(낮은 화면) - 시험: test_760_1 새로 · test_634_2 · test_553_4 겹 차례 표에 geology 한 줄 더함(새 겹 · 낡은 표) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
401 lines
15 KiB
Python
401 lines
15 KiB
Python
"""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
|
|
],
|
|
}
|