Files
Aislo/B05_Profile/B05_Profile_Engine_RouteInitial.py
T
eomsangdonandClaude Opus 5.5 55f1eeeb7d feat(B05): 초기 계획노선 갈래 — 추종 · 토공 최소 · 기울기 순한 · 비중(POST /route/initial)
- 갈래 인자 하나로 초기노선 생성 · follow 는 지금 로직 그대로(예상노선 폴리라인화)
- 탐색 = 예상노선 둘레 100 m 복도 · 4 m 격자 · 16 방향 + 진행 방향 기억 Dijkstra(꺾을 때 R·θ 벌점)
- 비중 = 예상 이격 · 토공 · 기울기 세 목표를 복도 평균으로 정규화해 가중합(합 100 · 기본 50/25/25)
- 공통 제약(종단기울기 상한 · 최소 곡선반지름 · 계류 이격) — 탐색 벌점 · 원호 끼우기 · 초과 구간 수
- 비교값: 연장 · 최대/평균 종단기울기 · 절성토량 · 기준 초과 구간 수 · 예상노선 이격
- 파일 안 씀(캐시용) · 기준값은 지식DB 후보 임시값(사용자 확정 대기)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01TULoa94ZFL26KU6ZqVpjkF
2026-09-29 09:59:24 +09:00

465 lines
22 KiB
Python

"""초기 계획노선 갈래 — 예상노선 추종 · 토공 최소 · 기울기 순한 노선 · 비중(PLAN 20장).
계산 자리 = 서버 단독(③). 결과는 화면 캐시용 — 이 모듈은 파일을 쓰지 않는다. 적용은 화면이
받은 노드로 기존 `POST /route/replan`([확인])을 부른다.
갈래(`mode`):
· follow — 지금 로직 그대로: 예상노선 점 묶음을 폴리라인화(`build_planned_polyline`).
`B05_Profile_Router_Replan._ensure_planned_initial` 과 같은 입력 · 같은 결과.
· earthwork — 절 · 성토 합 최소(DEM 옆경사로 가늠한 단면 토량 + 상한 넘는 지반 기울기).
· gentle — 종단기울기 최소(지반 기울기 제곱 합).
· weighted — 예상노선 이격 · 토공 · 기울기 세 목표를 **정규화해 가중합**(합 100).
정규화 = 목표마다 복도 안 링크 평균값으로 나눔 — 단위가 다른 세 값을 1 m 당
「보통 크기 1」로 맞춘 뒤 비중을 곱한다.
탐색 = 예상노선 둘레 복도(`CORRIDOR_M`) 안 `SEARCH_CELL_M` 격자의 8-연결 Dijkstra(scipy csgraph).
탐색 격자 경로는 follow 와 같은 폴리라인화(꺾임점 뽑기 + 원호 끼우기)를 지난다.
공통 제약(모든 갈래) — 종단기울기 상한 · 최소 곡선반지름 · 계류 이격:
· 탐색 비용 — 상한 넘는 지반 기울기 링크 · 계류 버퍼 안 링크에 벌점(막지 않음 — 길이 없으면
노선이 끊기므로).
· 폴리라인화 — 최소 곡선반지름으로 원호를 끼우고 못 끼우면 위반 표시.
· 결과 — `metrics.over_limit` 에 제약별 초과 구간 수.
⚠ 사용자 확정 대기 — 아래 임시값과 호출부가 넘기는 종단기울기 상한 · 최소 곡선반지름은 지식DB
(`resources/knowledge/technical_info/01_임도`) 후보일 뿐이다(PLAN 20장 기준값 후보 표).
확정되면 이 머리의 상수만 고친다.
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from pathlib import Path
from typing import Any
import numpy as np
import shapely
from scipy import ndimage
from scipy.interpolate import RegularGridInterpolator
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import dijkstra
from shapely.geometry import LineString, MultiLineString
from common_util.common_util_route_polyline import build_planned_polyline
# ── 사용자 확정 대기(지식DB 후보 · PLAN 20장) ─────────────────────────────────────
# 계류 이격(m) — 법정 이격 수치 없음. 타당성평가 「산지계류 100 m 이내 구간 비율」(별표1 비고,
# 1/25,000 주요 수계)을 후보로 씀.
STREAM_OFFSET_M = 100.0
# 이격을 따질 계류 — 도엽 하천중심선 `구분`. 「주요 수계」 가늠으로 세류 · 연결선은 뺌.
STREAM_CLASSES = ("소하천", "지방하천", "국가하천")
# 토공 가늠 폭(m) = 유효너비 3 + 길어깨 0.5 × 2(별표2 간선).
ROAD_WIDTH_M = 4.0
# 탐색 복도 반폭(m) — 예상노선에서 이보다 멀리는 가지 않음(법령 값 아님 · 알고리즘 값).
CORRIDOR_M = 100.0
PENDING = ("max_grade_pct", "min_radius_m", "stream_offset_m", "road_width_m", "corridor_m")
# ── 탐색 · 측정 값 ───────────────────────────────────────────────────────────────
MODES = ("follow", "earthwork", "gentle", "weighted")
DEFAULT_WEIGHTS = {"expected": 50.0, "earthwork": 25.0, "grade": 25.0}
MODE_WEIGHTS = {
"follow": {"expected": 100.0, "earthwork": 0.0, "grade": 0.0},
"earthwork": {"expected": 0.0, "earthwork": 100.0, "grade": 0.0},
"gentle": {"expected": 0.0, "earthwork": 0.0, "grade": 100.0},
}
SEARCH_CELL_M = 4.0 # 탐색 격자 칸(m) — 비용면(2 m)을 솎아 씀
SAMPLE_M = 2.0 # 결과 측정 간격(m)
STATION_M = 20.0 # 종단기울기 측정 간격(m) — 타당성평가 20 m 간격 · 측점 체계
LENGTH_COST = 0.05 # 1 m 당 거리 비용 — 목표가 평평한 곳에서 헤매지 않게
# 지반 기울기가 상한을 넘는 링크 1 m 당 벌점 = 이 값 × (넘은 기울기 / 상한). 세게 두면(20)
# 제약만 남고 비중이 안 먹음(랩탑_보조 실측: 비중 100/0/0 · 80/10/10 · 50/25/25 결과 같음).
# 2 에서 비중을 올릴수록 이격 · 토공이 한 방향으로 바뀜(이격 4 → 16 m · 토공 4.9 → 4.0천 m³).
GRADE_OVER_COST = 2.0
STREAM_COST = 3.0 # 계류 버퍼 안 링크 1 m 당 벌점 — 건너기는 되고 나란히 달리기는 피함
TURN_COST = 1.0 # 꺾을 때 벌점 = 최소 곡선반지름 원호 길이(R·θ) × 이 값
MAX_TURN_STEPS = 2 # 한 걸음에 꺾는 방향 칸 수 한도(16 방향 중 ±2 = 약 45°)
# 16 방향(행, 열) — +x 에서 반시계. 나이트 걸음(1, 2) 등을 넣어 22.5° 안팎 방향을 곧게 달림.
DIRECTIONS = (
(0, 1),
(1, 2),
(1, 1),
(2, 1),
(1, 0),
(2, -1),
(1, -1),
(1, -2),
(0, -1),
(-1, -2),
(-1, -1),
(-2, -1),
(-1, 0),
(-2, 1),
(-1, 1),
(-1, 2),
)
# 계류 버퍼 안 한 구간이 건너기로 인정되는 길이 — 45° 이상으로 가로지를 때의 버퍼 안 길이.
STREAM_CROSS_MAX_M = 2.0 * STREAM_OFFSET_M * math.sqrt(2.0)
@dataclass
class Terrain:
"""지반 격자 — 탐색 비용면(`B05_Profile_Engine_Solver._load_or_build_cost_surface`) 그대로."""
x: np.ndarray # 열 좌표(오름차순)
y: np.ndarray # 행 좌표(오름차순)
z: np.ndarray # [행, 열] 표고
valid: np.ndarray # [행, 열] 지형 있음
def load_terrain_and_streams(
project_root: Path, selection: dict[str, Any]
) -> tuple[Terrain, list[list[tuple[float, float]]]]:
"""지반 격자(탐색 비용면 캐시 그대로)와 이격을 따질 계류(사업지 좌표계)를 읽는다."""
from pyproj import Transformer
from B04_PreProcess.B04_PreProcess_Router_Watershed import (
STREAM_FILE,
_load_features,
_reproject_features,
)
from B05_Profile.B05_Profile_Engine_Solver import _MODELS_SUBDIR, _load_or_build_cost_surface
from common_util.common_util_crs import resolve_project_crs
x, y, z, valid, *_ = _load_or_build_cost_surface(
project_root,
project_root / _MODELS_SUBDIR,
str(selection.get("source_filter")),
str(selection.get("method") or "dtm"),
bool(selection.get("smooth")),
)
to_metric = Transformer.from_crs("EPSG:4326", resolve_project_crs(project_root), always_xy=True)
features = _load_features(project_root / "B04_PreProcess" / "processed", STREAM_FILE)
streams = stream_lines(_reproject_features(features, to_metric))
return Terrain(np.asarray(x), np.asarray(y), np.asarray(z), np.asarray(valid)), streams
def resolve_weights(mode: str, weights: dict[str, float] | None) -> dict[str, float]:
"""갈래의 비중(합 100). weighted 가 아니면 갈래가 정한 값."""
if mode != "weighted":
return dict(MODE_WEIGHTS[mode])
merged = {**DEFAULT_WEIGHTS, **(weights or {})}
return {key: float(merged[key]) for key in DEFAULT_WEIGHTS}
def stream_lines(features: list[dict[str, Any]]) -> list[list[tuple[float, float]]]:
"""이격을 따질 계류만 좌표 목록으로(사업지 좌표계로 이미 옮긴 피처)."""
lines: list[list[tuple[float, float]]] = []
for feature in features:
if (feature.get("properties") or {}).get("구분") not in STREAM_CLASSES:
continue
geometry = feature.get("geometry") or {}
parts = geometry.get("coordinates") or []
if geometry.get("type") == "LineString":
parts = [parts]
elif geometry.get("type") != "MultiLineString":
continue
for part in parts:
points = [(float(p[0]), float(p[1])) for p in part if len(p) >= 2]
if len(points) >= 2:
lines.append(points)
return lines
def generate_initial_route(
expected: list[tuple[float, float]],
terrain: Terrain,
streams: list[list[tuple[float, float]]],
*,
mode: str,
weights: dict[str, float] | None,
max_grade: float,
min_radius_m: float,
) -> dict[str, Any]:
"""갈래 하나로 초기 계획노선을 만들고 비교값을 붙여 돌려준다. `max_grade` 는 비율(0.14)."""
if mode not in MODES:
raise ValueError(f"모르는 갈래: {mode}")
if len(expected) < 2:
raise ValueError("예상노선 정점이 2개 미만입니다.")
used = resolve_weights(mode, weights)
points = (
list(expected)
if mode == "follow"
else _search_route(expected, terrain, streams, used, max_grade, min_radius_m)
)
outline = build_planned_polyline(points, min_radius_m=min_radius_m, simplify=True)
metrics = route_metrics(
outline.vertices,
expected,
terrain,
streams,
max_grade=max_grade,
radius_runs=sum(1 for curve in outline.curves if curve.violations),
)
return {
"mode": mode,
"weights": used,
"planned": [[round(x, 4), round(y, 4)] for x, y in outline.vertices],
"nodes": [node.as_dict() for node in outline.nodes],
"curves": [curve.as_dict() for curve in outline.curves],
"metrics": metrics,
"criteria": {
"max_grade_pct": round(max_grade * 100.0, 2),
"min_radius_m": round(min_radius_m, 2),
"stream_offset_m": STREAM_OFFSET_M,
"road_width_m": ROAD_WIDTH_M,
"corridor_m": CORRIDOR_M,
"pending": list(PENDING),
},
}
# ── 탐색 ─────────────────────────────────────────────────────────────────────────
def _cell_size(coords: np.ndarray) -> float:
return float((coords[-1] - coords[0]) / (len(coords) - 1)) if len(coords) > 1 else 1.0
def _distance_grid(
lines: list[list[tuple[float, float]]], xs: np.ndarray, ys: np.ndarray, cell: float
) -> np.ndarray:
"""격자 칸마다 선까지 거리(m). 선이 격자에 안 닿으면 모두 무한."""
marked = np.zeros((len(ys), len(xs)), dtype=bool)
for line in lines:
for (x0, y0), (x1, y1) in zip(line, line[1:]):
count = max(2, int(math.hypot(x1 - x0, y1 - y0) / (cell * 0.5)) + 1)
px = np.linspace(x0, x1, count)
py = np.linspace(y0, y1, count)
cols = np.rint((px - xs[0]) / cell).astype(int)
rows = np.rint((py - ys[0]) / cell).astype(int)
inside = (cols >= 0) & (cols < len(xs)) & (rows >= 0) & (rows < len(ys))
marked[rows[inside], cols[inside]] = True
if not marked.any():
return np.full(marked.shape, np.inf)
return ndimage.distance_transform_edt(~marked) * cell
def _search_route(
expected: list[tuple[float, float]],
terrain: Terrain,
streams: list[list[tuple[float, float]]],
weights: dict[str, float],
max_grade: float,
min_radius_m: float,
) -> list[tuple[float, float]]:
"""복도 안 격자에서 진행 방향을 기억하는 Dijkstra — 시점 · 종점은 예상노선 양 끝 그대로.
상태 = (칸, 들어온 방향 16 가지). 한 걸음에 22.5° · 45° 까지만 꺾고, 꺾을 때마다 최소
곡선반지름 원호 길이만큼 벌점을 매겨 곧게 달리게 한다. 16 방향(체스 나이트 걸음 포함)이라
8 방향 격자가 급한 옆경사에서 등고선을 못 따라 계단처럼 꺾이던 것을 피한다.
"""
step = max(1, int(round(SEARCH_CELL_M / _cell_size(terrain.x))))
xs, ys = terrain.x[::step], terrain.y[::step]
z = np.asarray(terrain.z[::step, ::step], dtype=np.float64)
cell = _cell_size(xs)
d_expected = _distance_grid([expected], xs, ys, cell)
d_stream = _distance_grid(streams, xs, ys, cell)
corridor = (d_expected <= CORRIDOR_M) & np.asarray(terrain.valid[::step, ::step], dtype=bool)
grad_y, grad_x = np.gradient(z, cell)
rows, cols = z.shape
index = np.full((rows, cols), -1, dtype=np.int64)
cells = int(corridor.sum())
index[corridor] = np.arange(cells)
# 방향마다 링크(칸 → 칸)와 세 목표 · 벌점(1 m 당)
links: list[tuple[np.ndarray, np.ndarray, float, np.ndarray, np.ndarray]] = []
for dr, dc in DIRECTIONS:
first = (slice(max(0, -dr), rows - max(0, dr)), slice(max(0, -dc), cols - max(0, dc)))
second = (slice(max(0, dr), rows - max(0, -dr)), slice(max(0, dc), cols - max(0, -dc)))
a, b = index[first], index[second]
ok = (a >= 0) & (b >= 0)
norm = math.hypot(dr, dc)
grade = np.abs(z[second][ok] - z[first][ok]) / (cell * norm)
over = np.maximum(0.0, grade - max_grade)
# 옆경사 = 진행 방향에 수직인 지반 기울기(두 칸 평균)
cross = np.abs(
(grad_x[first][ok] + grad_x[second][ok]) * -dr
+ (grad_y[first][ok] + grad_y[second][ok]) * dc
) / (2.0 * norm)
# 토공 가늠(1 m 당 m³) — 옆경사 위 반절 · 반성 단면 W²·s/4 + 상한 넘는 기울기가 측점
# 반 간격 동안 쌓는 높이 차 단면
earth = ROAD_WIDTH_M**2 * cross / 4.0 + ROAD_WIDTH_M * over * STATION_M / 2.0
offset = 0.5 * (d_expected[first][ok] + d_expected[second][ok])
in_buffer = 0.5 * (
(d_stream[first][ok] < STREAM_OFFSET_M).astype(float)
+ (d_stream[second][ok] < STREAM_OFFSET_M)
)
penalty = GRADE_OVER_COST * over / max_grade + STREAM_COST * in_buffer
links.append(
(
a[ok],
b[ok],
cell * norm,
np.stack([offset, earth, (grade / max_grade) ** 2]),
penalty,
)
)
if not any(link[0].size for link in links):
raise ValueError("예상노선 둘레에 탐색할 지형이 없습니다.")
# 세 목표 정규화 — 복도 안 링크 평균이 1 이 되게(길이 가중). 평균이 0 이면 그 목표는 0.
lengths = np.concatenate([np.full(link[0].size, link[2]) for link in links])
stacked = np.hstack([link[3] for link in links])
scale = (stacked * lengths).sum(axis=1) / max(lengths.sum(), 1e-9)
order = ("expected", "earthwork", "grade")
share = np.array([weights[key] for key in order]) / max(sum(weights.values()), 1e-9)
factor = np.divide(share, scale, out=np.zeros_like(share), where=scale > 1e-12)
count = len(DIRECTIONS)
angles = [math.atan2(dr, dc) for dr, dc in DIRECTIONS]
heads, tails, costs = [], [], []
for out, (a, b, length, term, penalty) in enumerate(links):
move = length * (LENGTH_COST + factor @ term + penalty)
for turn in range(-MAX_TURN_STEPS, MAX_TURN_STEPS + 1):
came = (out + turn) % count
bend = abs((angles[out] - angles[came] + math.pi) % (2.0 * math.pi) - math.pi)
heads.append(a * count + came)
tails.append(b * count + out)
costs.append(move + TURN_COST * min_radius_m * bend)
graph = csr_matrix(
(np.concatenate(costs), (np.concatenate(heads), np.concatenate(tails))),
shape=(cells * count, cells * count),
)
cell_rows, cell_cols = np.nonzero(corridor)
centres = np.column_stack([xs[cell_cols], ys[cell_rows]])
def nearest(point: tuple[float, float]) -> int:
return int(np.argmin(np.hypot(centres[:, 0] - point[0], centres[:, 1] - point[1])))
start, end = nearest(expected[0]), nearest(expected[-1])
distance, predecessors, _ = dijkstra(
graph,
indices=[start * count + d for d in range(count)],
return_predecessors=True,
min_only=True,
)
state = end * count + int(np.argmin(distance[end * count : (end + 1) * count]))
if not np.isfinite(distance[state]):
raise ValueError("복도 안에서 시점과 종점을 잇는 길을 찾지 못했습니다.")
path = [state]
while predecessors[path[-1]] >= 0:
path.append(int(predecessors[path[-1]]))
points = [
(float(centres[s // count, 0]), float(centres[s // count, 1])) for s in reversed(path)
]
# 양 끝은 격자 칸 가운데가 아니라 예상노선 끝점 그대로 — 기점 · 종점은 옮기지 않는다.
points[0], points[-1] = tuple(expected[0]), tuple(expected[-1])
return points
# ── 비교값 ───────────────────────────────────────────────────────────────────────
def _resample(vertices: list[tuple[float, float]], step: float) -> tuple[np.ndarray, np.ndarray]:
"""폴리라인을 `step` 간격 점으로 — (점 [n,2], 누가거리 [n])."""
xy = np.asarray(vertices, dtype=np.float64)
seg = np.hypot(*np.diff(xy, axis=0).T)
chain = np.concatenate([[0.0], np.cumsum(seg)])
total = float(chain[-1])
at = np.append(np.arange(0.0, total, step), total) if total > 0 else np.array([0.0])
return np.column_stack([np.interp(at, chain, xy[:, 0]), np.interp(at, chain, xy[:, 1])]), at
def _runs(flags: np.ndarray) -> list[tuple[int, int]]:
"""참이 이어진 구간들의 (첫, 끝) 번호."""
runs: list[tuple[int, int]] = []
start = None
for i, flag in enumerate(flags):
if flag and start is None:
start = i
elif not flag and start is not None:
runs.append((start, i - 1))
start = None
if start is not None:
runs.append((start, len(flags) - 1))
return runs
def grade_line(ground: np.ndarray, step: float, max_grade: float) -> np.ndarray:
"""상한 안의 계획고 가늠 — 앞으로 · 뒤로 기울기를 누른 두 선의 평균(둘 다 상한 안 → 평균도)."""
rise = max_grade * step
forward, backward = ground.copy(), ground.copy()
for i in range(1, len(ground)):
forward[i] = min(max(ground[i], forward[i - 1] - rise), forward[i - 1] + rise)
for i in range(len(ground) - 2, -1, -1):
backward[i] = min(max(ground[i], backward[i + 1] - rise), backward[i + 1] + rise)
return 0.5 * (forward + backward)
def route_metrics(
vertices: list[tuple[float, float]],
expected: list[tuple[float, float]],
terrain: Terrain,
streams: list[list[tuple[float, float]]],
*,
max_grade: float,
radius_runs: int,
) -> dict[str, Any]:
"""연장 · 최대/평균 종단기울기(지반) · 절성토량(가늠) · 기준 초과 구간 수 · 예상노선 이격."""
sampler = RegularGridInterpolator(
(terrain.y, terrain.x), terrain.z, bounds_error=False, fill_value=None
)
points, chain = _resample(vertices, SAMPLE_M)
length = float(chain[-1])
ground = sampler(points[:, ::-1])
# 종단기울기 — 측점(20 m) 사이 지반 기울기
stations = np.append(np.arange(0.0, length, STATION_M), length)
station_z = np.interp(stations, chain, ground)
spans = np.diff(stations)
keep = spans > 0.5
grades = np.abs(np.diff(station_z))[keep] / spans[keep]
max_grade_run = float(grades.max()) if grades.size else 0.0
avg_grade = float((grades * spans[keep]).sum() / spans[keep].sum()) if grades.size else 0.0
# 절 · 성토(가늠) — 상한 안 계획고와 옆경사 단면을 폭 W 로 적분(비탈면 제외)
grad_y, grad_x = np.gradient(terrain.z, _cell_size(terrain.y), _cell_size(terrain.x))
gx = RegularGridInterpolator(
(terrain.y, terrain.x), grad_x, bounds_error=False, fill_value=None
)
gy = RegularGridInterpolator(
(terrain.y, terrain.x), grad_y, bounds_error=False, fill_value=None
)
tangent = np.gradient(points, axis=0)
tangent /= np.maximum(np.hypot(tangent[:, 0], tangent[:, 1]), 1e-9)[:, None]
lookup = points[:, ::-1]
cross = np.abs(gx(lookup) * -tangent[:, 1] + gy(lookup) * tangent[:, 0])
height = grade_line(ground, SAMPLE_M, max_grade) - ground
across = np.linspace(-ROAD_WIDTH_M / 2.0, ROAD_WIDTH_M / 2.0, 9)
depth = cross[:, None] * across[None, :] - height[:, None] # +: 지반이 노면 위 = 절토
weight = np.gradient(chain) if len(chain) > 1 else np.zeros_like(chain)
cut = float((np.maximum(depth, 0.0).mean(axis=1) * ROAD_WIDTH_M * weight).sum())
fill = float((np.maximum(-depth, 0.0).mean(axis=1) * ROAD_WIDTH_M * weight).sum())
# 예상노선 이격
route_points = shapely.points(points)
offsets = shapely.distance(route_points, LineString(expected))
# 계류 이격 — 버퍼 안 구간 중 봐줄 길이를 넘는 것. 봐줌 = 건너기(계류에 닿음) 한 번의
# 버퍼 안 길이 + 기점 · 종점이 버퍼 안이면 45° 로 빠져나가는 길이(끝점은 옮길 수 없음)
stream_runs = 0
if streams:
near = shapely.distance(route_points, MultiLineString(streams))
for first, last in _runs(near < STREAM_OFFSET_M):
ends = int(first == 0) + int(last == len(near) - 1)
crosses = float(near[first : last + 1].min()) <= SAMPLE_M
allowance = ends * STREAM_OFFSET_M * math.sqrt(2.0) + crosses * STREAM_CROSS_MAX_M
if chain[last] - chain[first] > allowance:
stream_runs += 1
grade_runs = len(_runs(grades > max_grade))
return {
"length_m": round(length, 2),
"max_grade_pct": round(max_grade_run * 100.0, 2),
"avg_grade_pct": round(avg_grade * 100.0, 2),
"cut_m3": round(cut, 1),
"fill_m3": round(fill, 1),
"over_limit": {
"grade": grade_runs,
"radius": int(radius_runs),
"stream": stream_runs,
"total": grade_runs + int(radius_runs) + stream_runs,
},
"offset_avg_m": round(float(offsets.mean()), 2),
"offset_max_m": round(float(offsets.max()), 2),
}