Files
Aislo/common_util/common_util_route_geometry.py
T
eomsangdonandClaude Opus 5 e83171b2b3 refactor(공통): 작업 좌표계 폴백을 창구 하나로 통합 + 트림 실패 진단
같은 폴백 사다리가 네 곳에 서로 다른 모양으로 흩어져 있었고, 그중 배수유역 라우터
`_route_center_lonlat` 은 노선 CSV 의 `crs_epsg` **라벨**을 그대로 변환에 썼다(라벨과
실좌표계가 다른 사례 실측 — 2026-09-01 용화 라벨 5179 / 실제 5176).

`common_util_crs.resolve_project_crs()` 신설 — 사다리는 ① 노선 crs_input ② 파일 라벨
(원본 좌표를 읽는 자리만) ③ 지형 PRJ(작업 좌표계 정본) ④ DB epsg ⑤ EPSG:5186.
작업 좌표계는 **지형 PRJ 우선으로 확정**(2026-09-03 사용자 결정 — 서피스 격자가 모델좌표의
주인이라 현행 유지). `project_epsg_from_prj()` 도 이 창구로 위임.

더해 지표면 트림이 노선을 통째로 지울 때 노선·지표면 bbox 를 함께 로그에 남긴다 —
"겹치지 않습니다" 만으로는 좌표계 문제인지 측량 범위 문제인지 갈리지 않았다.

검증 — `tmp/tests/test_project_crs_resolution.py` 5건 신설(사다리 4·진단 1),
전체 375 passed·17 skipped, ruff format 무변경.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-09-03 13:25:20 +09:00

585 lines
22 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
"""계획 노선 기하 공용 유틸 — 정점·누가거리·세류 교차점.
배수유역 분석(B04)과 관 편집·세부유역(B05)이 같은 노선 표현을 써야 하므로 여기 한 곳에만
정의한다. 어느 한쪽 페이지 폴더에 두면 반대 방향 import가 생긴다.
노선 원천은 두 가지다.
· B03에 업로드된 **계획 노선 파일**(CSV·shapefile) — 배수유역 분석의 입력
· DB `route_points` — B05에서 탐색·확정한 노선
둘 다 같은 `RouteVertex` 목록으로 바꿔 아래 함수들이 그대로 받는다.
"""
from __future__ import annotations
import csv
import logging
import math
import struct
from dataclasses import dataclass
from pathlib import Path
from typing import Any
from shapely.geometry import LineString, Point, shape
from config.config_system import SURFACE_ROUTE_EDGE_TRIM_M
logger = logging.getLogger(__name__)
# 계획 노선 CSV 열 이름 후보. B03이 여러 형식을 받게 되므로 흔한 표기를 모두 받아 준다.
_X_KEYS = ("x", "X", "동", "easting", "EASTING")
_Y_KEYS = ("y", "Y", "북", "northing", "NORTHING")
_Z_KEYS = ("z", "Z", "표고", "elevation", "ELEV")
_ORDER_KEYS = ("sequence", "order", "seq", "no", "번호")
_EPSG_KEYS = ("crs_epsg", "epsg", "EPSG")
@dataclass
class RouteVertex:
"""노선 폴리라인의 한 점. chainage는 시점 기준 누가거리(m)."""
x: float
y: float
z: float
chainage_m: float
@dataclass
class StructureCandidate:
"""관 매설 구조물 측점 후보."""
chainage_m: float
x: float
y: float
# "stream"=세류 교차, "spacing"=최대 간격 규칙 보충, "confirmed"=사용자 확정
reason: str
stream_name: str | None = None
@dataclass
class PlannedRoute:
"""계획 노선 파일에서 읽은 노선."""
vertices: list[RouteVertex]
epsg: int | None
name: str | None
source: Path
# `Transformer.from_crs` 입력 문자열("EPSG:n" 또는 **원문 WKT**). 좌표계를 EPSG 코드로
# 가리면 라벨이 안 붙는 PRJ가 통째로 막히므로, 변환은 언제나 이 값으로 한다.
# `epsg`는 표시·로그용 라벨일 뿐이다 (2026-08-31 사용자 확정).
crs_input: str | None = None
@property
def line(self) -> LineString:
return LineString([(vertex.x, vertex.y) for vertex in self.vertices])
def read_planned_route_csv(path: Path) -> PlannedRoute | None:
"""계획 노선 CSV를 읽어 정점 목록으로 바꾼다.
열 이름은 대소문자·한글 표기를 함께 받아 준다(B03이 여러 형식을 수용할 예정).
`sequence`가 있으면 그 순서로 정렬하고, 없으면 파일에 적힌 순서를 그대로 쓴다.
"""
try:
with path.open("r", encoding="utf-8-sig", newline="") as file:
rows = list(csv.DictReader(file))
except (OSError, csv.Error, UnicodeDecodeError):
logger.warning("계획 노선 CSV를 읽지 못했습니다: %s", path)
return None
if not rows:
return None
epsg = _first_int(rows[0], _EPSG_KEYS)
name = _first_text(rows[0], ("route_name", "name", "노선명"))
parsed: list[tuple[float, float, float, float]] = [] # (정렬키, x, y, z)
for index, row in enumerate(rows):
x = _first_float(row, _X_KEYS)
y = _first_float(row, _Y_KEYS)
if x is None or y is None:
continue
order = _first_float(row, _ORDER_KEYS)
parsed.append(
(float(index) if order is None else order, x, y, _first_float(row, _Z_KEYS) or 0.0)
)
if len(parsed) < 2:
logger.warning("계획 노선 CSV에 좌표가 2점 미만입니다: %s", path)
return None
parsed.sort(key=lambda item: item[0])
vertices: list[RouteVertex] = []
cumulative = 0.0
previous: tuple[float, float] | None = None
for _, x, y, z in parsed:
if previous is not None:
cumulative += math.dist(previous, (x, y))
vertices.append(RouteVertex(x=x, y=y, z=z, chainage_m=cumulative))
previous = (x, y)
logger.info(
"계획 노선 %s: 정점 %d개, 연장 %.0fm, EPSG %s", path.name, len(vertices), cumulative, epsg
)
return PlannedRoute(
vertices=vertices,
epsg=epsg,
name=name,
source=path,
crs_input=f"EPSG:{epsg}" if epsg else None,
)
def read_planned_route_shapefile(path: Path) -> PlannedRoute | None:
"""계획 노선 shapefile을 읽어 정점 목록으로 바꾼다.
원청 자료는 **2차원 폴리라인**이라 z가 없다. 지반고는 확정 서피스에서 뽑으므로
(`build_surface_sampler`) 여기서는 0으로 채운다 — CSV의 z와 같은 취급이다.
파트가 여럿이면 가장 긴 것을 노선으로 본다.
"""
from B03_FileInput.B03_FileInput_Engine_Shapefile import (
read_shapefile_attributes,
read_shapefile_parts,
shapefile_crs_input,
shapefile_epsg_label,
)
try:
parts = read_shapefile_parts(path)
except (OSError, ValueError, struct.error) as error:
logger.warning("계획 노선 shapefile을 읽지 못했습니다: %s (%s)", path, error)
return None
if not parts:
logger.warning("계획 노선 shapefile에 폴리라인이 없습니다: %s", path)
return None
points = max(parts, key=len)
if len(points) < 2:
logger.warning("계획 노선 shapefile에 좌표가 2점 미만입니다: %s", path)
return None
vertices: list[RouteVertex] = []
cumulative = 0.0
previous: tuple[float, float] | None = None
for x, y, z in points:
if previous is not None:
cumulative += math.dist(previous, (x, y))
vertices.append(RouteVertex(x=x, y=y, z=z, chainage_m=cumulative))
previous = (x, y)
attributes = read_shapefile_attributes(path)
name = next(
(
attributes[key]
for key in ("대상지", "노선명", "route_name", "name")
if attributes.get(key)
),
path.stem,
)
epsg = shapefile_epsg_label(path)
logger.info(
"계획 노선 %s: 정점 %d개, 연장 %.0fm, EPSG %s(라벨)",
path.name,
len(vertices),
cumulative,
epsg,
)
return PlannedRoute(
vertices=vertices,
epsg=epsg,
name=name,
source=path,
crs_input=shapefile_crs_input(path),
)
def read_planned_route(path: Path) -> PlannedRoute | None:
"""계획 노선 파일을 확장자에 맞는 판독기로 읽는다(CSV·shapefile)."""
if path.suffix.lower() == ".shp":
return read_planned_route_shapefile(path)
return read_planned_route_csv(path)
def find_planned_route_file(input_dir: Path) -> Path | None:
"""B03 입력 폴더에서 계획 노선 파일을 찾는다.
shapefile을 CSV보다 우선한다 — 원청 정식 노선이 shapefile로 오고, CSV는 예전
자료거나 좌표만 뽑아 둔 보조본인 경우가 많다. 같은 종류가 여럿이면 가장 최근 것.
"""
if not input_dir.exists():
return None
for pattern in ("*.shp", "*.csv"):
candidates = sorted(
input_dir.rglob(pattern), key=lambda item: item.stat().st_mtime, reverse=True
)
if candidates:
return candidates[0]
return None
def load_design_route(
project_root: Path, surface_params: dict[str, Any] | None = None
) -> PlannedRoute | None:
"""설계가 쓸 계획노선 한 벌을 만든다 — 읽기·좌표계 변환·트림·조밀화를 여기서 끝낸다.
노선을 읽는 곳이 여럿이라(체인·배수유역·유입·도엽) 각자 읽으면 트림이 적용된 곳과
안 된 곳이 갈린다. 실제로 그렇게 갈려 관 측점이 트림 전(2,136m) 기준으로 찍히고
확정 노선(1,070m)과 어긋나 종단 계획선이 직선으로 나왔다(2026-09-01).
**설계 계통은 전부 이 함수를 지난다.**
`surface_params`(확정 필터·방식·스무딩)를 주면 지표면이 덮지 못하는 구간을 잘라 내고,
B05 격자 탐색이 계획노선을 바꾸지 않도록 정점 간격을 직결 문턱 아래로 좁힌다.
주지 않으면 읽어서 좌표계만 맞춘 원본을 돌려준다.
"""
from B04_PreProcess.B04_PreProcess_Engine_Extent import project_epsg_from_prj
from common_util.common_util_initial_snapshot import design_route_csv_path
from config.config_system import (
ROUTE_DIRECT_LINK_CELL_FACTOR,
ROUTE_GRID_RES_M,
ROUTE_PLANNED_DENSIFY_SAFETY,
)
# 체인이 남긴 CSV 정본이 있으면 그것이 설계 노선이다 — 좌표계 변환·트림·조밀화가 이미
# 끝난 값이라 다시 하지 않는다(노선이 shapefile로 왔어도 여기서는 CSV 한 벌이다).
# 지표면·노선이 바뀌면 `discard_initial_snapshot()`이 폴더째 지우므로 이 경로가 저절로
# 닫히고 원본 재판독으로 되돌아간다. 트림 **전** 원본이 필요한 호출(도엽 범위 —
# surface_params 없음)은 여기를 타지 않는다.
if surface_params:
master = design_route_csv_path(project_root)
if master.is_file():
stored = read_planned_route_csv(master)
if stored is not None and len(stored.vertices) >= 2:
return replace_vertices(
stored,
[(v.x, v.y) for v in stored.vertices],
crs_input=project_epsg_from_prj(project_root),
)
route_file = find_planned_route_file(project_root / "B03_FileInput" / "input")
planned = read_planned_route(route_file) if route_file else None
if planned is None or len(planned.vertices) < 2:
return None
target_crs = project_epsg_from_prj(project_root)
points = [(float(v.x), float(v.y)) for v in planned.vertices]
source_crs = planned.crs_input or target_crs
if source_crs.upper() != target_crs.upper():
from pyproj import Transformer
transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)
points = [transformer.transform(x, y) for x, y in points]
if surface_params:
from common_util.common_util_surface_sampler import build_surface_sampler
try:
sampler = build_surface_sampler(
project_root / "B04_PreProcess" / "models",
str(surface_params["source_filter"]),
str(surface_params["method"]),
bool(surface_params["smooth"]),
)
except (FileNotFoundError, KeyError, OSError) as exc:
logger.warning("설계 노선: 지표면을 열지 못해 트림을 건너뜁니다 — %s", exc)
else:
before = points
points = trim_route_to_surface(points, sampler)
if len(points) < 2:
_log_trim_wipeout(before, sampler, target_crs)
points = densify_route(
points,
ROUTE_DIRECT_LINK_CELL_FACTOR * ROUTE_GRID_RES_M * ROUTE_PLANNED_DENSIFY_SAFETY,
)
if len(points) < 2:
return None
return replace_vertices(planned, points, crs_input=target_crs)
def replace_vertices(
planned: PlannedRoute, points: list[tuple[float, float]], *, crs_input: str | None = None
) -> PlannedRoute:
"""XY 목록으로 노선 정점을 갈아 끼우고 누가거리를 다시 센다."""
vertices: list[RouteVertex] = []
cumulative = 0.0
previous: tuple[float, float] | None = None
for x, y in points:
if previous is not None:
cumulative += math.dist(previous, (x, y))
vertices.append(RouteVertex(x=float(x), y=float(y), z=0.0, chainage_m=cumulative))
previous = (x, y)
return PlannedRoute(
vertices=vertices,
epsg=planned.epsg,
name=planned.name,
source=planned.source,
crs_input=crs_input or planned.crs_input,
)
def _log_trim_wipeout(points: list[tuple[float, float]], sampler: Any, target_crs: str) -> None:
"""트림이 노선을 통째로 지운 이유를 **수치로** 남긴다.
노선과 지표면의 좌표계가 어긋나면 트림 결과가 빈 목록이 되는데, 그때 로그가
"겹치지 않습니다" 뿐이라 좌표계 문제인지 범위 문제인지 갈리지 않았다(2026-09-03 실사고:
도엽 서피스가 노선 원본 좌표계로 만들어져 노선이 통째로 지워짐). 두 bbox 를 함께 찍어
한 줄로 판별되게 한다 — 겹침이 0이면 좌표계, 일부 겹치면 측량 범위 문제다.
"""
xs = [x for x, _ in points]
ys = [y for _, y in points]
grid_x = getattr(sampler, "x", None)
grid_y = getattr(sampler, "y", None)
surface = (
f"[{float(min(grid_x)):.0f}~{float(max(grid_x)):.0f}, "
f"{float(min(grid_y)):.0f}~{float(max(grid_y)):.0f}]"
if grid_x is not None and grid_y is not None and len(grid_x) and len(grid_y)
else "미상"
)
logger.warning(
"설계 노선: 지표면 트림 결과가 비었습니다 — 노선 bbox [%.0f~%.0f, %.0f~%.0f] / "
"지표면 bbox %s (작업 좌표계 %s). 두 범위가 전혀 안 겹치면 좌표계 불일치입니다.",
min(xs),
max(xs),
min(ys),
max(ys),
surface,
target_crs,
)
def trim_route_to_surface(
points: list[tuple[float, float]],
sampler: Any,
edge_trim_m: float = SURFACE_ROUTE_EDGE_TRIM_M,
) -> list[tuple[float, float]]:
"""지표면이 덮는 구간만 남기고 계획노선을 자른다.
판정은 sampler가 돌려주는 `valid` — 확정 DTM의 valid_mask가 곧 불규칙한 실제
외곽이다(bounds 사각형이 아니다). 가장 긴 연속 유효 구간을 남긴다.
`edge_trim_m`은 **잘라 낸 쪽 끝에만** 적용한다. 서피스 가장자리는 점 밀도가 떨어져
지반고가 못 미덥기 때문이다. 노선 본래 끝점이 서피스 안이면 깎지 않는다 — 멀쩡한
구간을 짧게 만들 이유가 없다 (2026-09-01 사용자 확정).
전부 유효하면 입력을 그대로 돌려준다. 남는 구간이 2점 미만이면 빈 목록.
"""
import numpy as np
if len(points) < 2:
return list(points)
xy = np.asarray(points, dtype=np.float64)
try:
_, valid = sampler.sample_xy(xy)
except (ValueError, OSError) as exc:
logger.warning("노선 트림: 지표면 샘플링 실패 — %s", exc)
return list(points)
valid = np.asarray(valid, dtype=bool)
if valid.all():
return list(points)
if not valid.any():
logger.warning("노선 트림: 노선 전체가 지표면 밖입니다.")
return []
# 가장 긴 연속 유효 구간 — 가장자리에서 한두 점이 튀어도 본 구간을 잃지 않는다.
best_start = best_end = start = -1
for index, ok in enumerate([*valid.tolist(), False]):
if ok and start < 0:
start = index
elif not ok and start >= 0:
if index - start > best_end - best_start:
best_start, best_end = start, index
start = -1
kept = [(float(x), float(y)) for x, y in xy[best_start:best_end]]
trim_head = best_start > 0
trim_tail = best_end < len(xy)
kept = _trim_ends(kept, edge_trim_m if trim_head else 0.0, edge_trim_m if trim_tail else 0.0)
logger.info(
"노선 트림: 정점 %d개 → %d개 (지표면 밖 %d개, 가장자리 여유 %.0fm %s)",
len(points),
len(kept),
int((~valid).sum()),
edge_trim_m,
"앞뒤" if trim_head and trim_tail else ("앞" if trim_head else "뒤"),
)
return kept
def densify_route(
points: list[tuple[float, float]], max_spacing_m: float
) -> list[tuple[float, float]]:
"""원래 정점은 모두 남기고, 간격이 `max_spacing_m`을 넘는 구간에만 점을 끼워 넣는다.
**평면 형상은 바뀌지 않는다** — 같은 직선 위에 점을 더 찍을 뿐이다.
B05 경로 탐색은 제어점 간격이 `ROUTE_DIRECT_LINK_CELL_FACTOR × ROUTE_GRID_RES_M`
이하일 때만 격자 탐색을 건너뛰고 원좌표를 그대로 잇는다
(`B05_Profile_Engine_Solver.py:394`). 예정노선처럼 정점이 성긴 선(용화 평균 17.6m)은
그 문턱을 넘어 탐색을 타고, 급경사 현에서 "통과 경로 없음"으로 끊긴다.
간격만 좁혀 주면 계획노선이 손대지 않은 채 그대로 채택된다 (2026-09-01 사용자 확정).
"""
if len(points) < 2 or max_spacing_m <= 0:
return list(points)
dense: list[tuple[float, float]] = [points[0]]
for start, end in zip(points, points[1:]):
distance = math.dist(start, end)
steps = int(math.ceil(distance / max_spacing_m))
for step in range(1, steps):
ratio = step / steps
dense.append(
(
start[0] + (end[0] - start[0]) * ratio,
start[1] + (end[1] - start[1]) * ratio,
)
)
dense.append(end)
return dense
def _trim_ends(
points: list[tuple[float, float]], head_m: float, tail_m: float
) -> list[tuple[float, float]]:
"""폴리라인 앞뒤에서 지정 길이만큼 잘라 낸다. 남는 게 2점 미만이면 빈 목록."""
if len(points) < 2 or (head_m <= 0 and tail_m <= 0):
return points
line = LineString(points)
start = min(head_m, line.length)
end = max(start, line.length - tail_m)
if end - start <= 0:
logger.warning("노선 트림: 여유를 깎고 나니 남는 구간이 없습니다.")
return []
cumulative = 0.0
kept: list[tuple[float, float]] = [line.interpolate(start).coords[0]]
for index in range(1, len(points)):
cumulative += math.dist(points[index - 1], points[index])
if start < cumulative < end:
kept.append(points[index])
kept.append(line.interpolate(end).coords[0])
return kept if len(kept) >= 2 else []
def build_route_vertices(points: list[dict[str, Any]]) -> list[RouteVertex]:
"""DB route_points 행을 누가거리가 채워진 정점 목록으로 바꾼다."""
vertices: list[RouteVertex] = []
cumulative = 0.0
previous: tuple[float, float] | None = None
for row in points:
x = float(row["x"])
y = float(row["y"])
z = float(row.get("z") or 0.0)
if previous is not None:
cumulative += math.dist(previous, (x, y))
chainage = row.get("chainage_m")
vertices.append(
RouteVertex(
x=x, y=y, z=z, chainage_m=float(chainage) if chainage is not None else cumulative
)
)
previous = (x, y)
return vertices
def interpolate_vertex(
vertices: list[RouteVertex], chainage_m: float
) -> tuple[float, float, float]:
"""누가거리 위치의 (x, y, z)를 선형 보간한다. 범위 밖은 끝점으로 당긴다."""
if not vertices:
return (0.0, 0.0, 0.0)
if chainage_m <= vertices[0].chainage_m:
return (vertices[0].x, vertices[0].y, vertices[0].z)
for previous, current in zip(vertices, vertices[1:]):
if chainage_m <= current.chainage_m:
span = current.chainage_m - previous.chainage_m
ratio = 0.0 if span <= 0 else (chainage_m - previous.chainage_m) / span
return (
previous.x + (current.x - previous.x) * ratio,
previous.y + (current.y - previous.y) * ratio,
previous.z + (current.z - previous.z) * ratio,
)
last = vertices[-1]
return (last.x, last.y, last.z)
def is_uphill_at(vertices: list[RouteVertex], chainage_m: float, window_m: float = 20.0) -> bool:
"""해당 위치가 오르막(절토부)인지 종단 계획선의 국소 기울기 부호로 판정한다."""
_, _, back_z = interpolate_vertex(vertices, max(0.0, chainage_m - window_m))
_, _, forward_z = interpolate_vertex(vertices, chainage_m + window_m)
return forward_z >= back_z
def find_stream_crossings(
vertices: list[RouteVertex],
stream_features: list[dict[str, Any]],
) -> list[StructureCandidate]:
"""노선 평면 선형과 세류선의 교차 지점을 누가거리 순으로 찾는다."""
if len(vertices) < 2:
return []
route_line = LineString([(vertex.x, vertex.y) for vertex in vertices])
candidates: list[StructureCandidate] = []
for feature in stream_features:
geometry = feature.get("geometry")
if not geometry:
continue
try:
stream = shape(geometry)
except Exception: # noqa: BLE001 - 손상된 피처는 건너뛴다
continue
if stream.is_empty:
continue
intersection = route_line.intersection(stream)
if intersection.is_empty:
continue
name = _stream_name(feature)
for point in _collect_points(intersection):
candidates.append(
StructureCandidate(
chainage_m=route_line.project(point),
x=point.x,
y=point.y,
reason="stream",
stream_name=name,
)
)
candidates.sort(key=lambda item: item.chainage_m)
return candidates
def _stream_name(feature: dict[str, Any]) -> str | None:
properties = feature.get("properties") or {}
for key in ("명칭", "하천명", "NAME", "name"):
value = properties.get(key)
if value:
return str(value)
return None
def _collect_points(geometry: Any) -> list[Point]:
"""교차 결과(Point/MultiPoint/LineString 등)에서 대표 점들을 뽑는다."""
if geometry.geom_type == "Point":
return [geometry]
if geometry.geom_type in {"MultiPoint", "GeometryCollection"}:
points: list[Point] = []
for part in geometry.geoms:
points.extend(_collect_points(part))
return points
# 선분끼리 겹쳐 선으로 나온 경우는 중점을 대표로 쓴다.
if geometry.geom_type in {"LineString", "MultiLineString"}:
return [geometry.interpolate(0.5, normalized=True)]
return []
def _first_text(row: dict[str, Any], keys: tuple[str, ...]) -> str | None:
for key in keys:
value = row.get(key)
if value not in (None, ""):
return str(value).strip()
return None
def _first_float(row: dict[str, Any], keys: tuple[str, ...]) -> float | None:
text = _first_text(row, keys)
if text is None:
return None
try:
return float(text)
except ValueError:
return None
def _first_int(row: dict[str, Any], keys: tuple[str, ...]) -> int | None:
value = _first_float(row, keys)
return None if value is None else int(value)