Aislo/B05_Profile/B05_Profile_Engine_RouteInitial_Metrics.py

"""초기 계획노선 갈래의 비교값 · 제약 못 지킨 구간 목록(PLAN 20장).

`B05_Profile_Engine_RouteInitial` 이 만든 폴리라인 하나를 받아 연장 · 종단기울기 · 절성토량 ·
예상노선 이격을 재고, 공통 제약(종단기울기 상한 · 최소 곡선반지름 · 계류 이격)을 못 지킨
구간을 측점 · 값 · 사유로 모은다. `over_limit` 수는 이 목록 수와 같다.
"""

from __future__ import annotations

import math
from typing import TYPE_CHECKING, Any

import numpy as np
import shapely
from scipy.interpolate import RegularGridInterpolator
from shapely.geometry import LineString, MultiLineString, Point

if TYPE_CHECKING:
    from B05_Profile.B05_Profile_Engine_RouteInitial import Terrain
    from common_util.common_util_route_polyline import RouteCurve

SAMPLE_M = 2.0  # 측정 간격(m)
# 종단기울기를 잴 지반 = DEM 을 이 폭(가우스 σ, m)으로 고른 것 — 도상(1/5,000) 측정 가늠. 지반의
# 잔굴곡(수 m 턱 · 절개면)은 노면 계획고가 흡수하므로 노선 기울기에서 뺀다(사용자 확정 대기).
GRADE_SMOOTH_M = 5.0
STATION_M = 20.0  # 종단기울기 측정 간격(m) — 측점 체계 · 타당성평가 20 m 간격
GRADE_EPS = 1e-6  # 상한과 같은 기울기는 지킨 것으로 봄

REASON_TERRAIN = "지형상 불가능 — 복도 안 어떤 길도 {bound:.1f}% 이상 기울기를 지나야 이어짐"
REASON_GRADE_FIT = "곡선 맞춤 뒤 지반 차 — 측점 사이 지반이 상한을 넘음"
REASON_RADIUS = "곡선 둘 자리 모자람 — 앞뒤 직선이 짧아 최소 반지름 원호가 안 들어감"
REASON_STREAM = "계류 버퍼 안 — 건너기 · 기점 · 종점 봐줌을 넘음"
REASON_STREAM_TERRAIN = "지형상 불가능 — 기울기 · 반지름을 지키며 계류 버퍼를 비켜 잇는 길 없음"


def smoothed_ground(terrain: Terrain, grade_smooth_m: float = GRADE_SMOOTH_M) -> np.ndarray:
    """종단기울기를 잴 지반(`grade_smooth_m` 로 고른 DEM) — 탐색과 비교값이 같은 것을 본다.

    폭은 화면에서 바꿀 수 있다(PLAN 20장 확정 · 기본값 `GRADE_SMOOTH_M`)."""
    from scipy.ndimage import gaussian_filter

    spacing = 0.5 * (_spacing(terrain.x) + _spacing(terrain.y))
    return gaussian_filter(np.asarray(terrain.z, dtype=np.float64), grade_smooth_m / spacing)


def station_label(chainage_m: float) -> str:
    """측점 표기 — No.(20 m 번호)+나머지(m)."""
    number = int(chainage_m // STATION_M)
    return f"No.{number}+{chainage_m - number * STATION_M:.1f}"


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 _item(kind: str, start: float, end: float, value: float, limit: float, reason: str) -> dict:
    return {
        "kind": kind,
        "from_m": round(start, 1),
        "to_m": round(end, 1),
        "station": station_label(start),
        "value": round(value, 2),
        "limit": round(limit, 2),
        "reason": reason,
    }


def route_metrics(
    vertices: list[tuple[float, float]],
    expected: list[tuple[float, float]],
    terrain: Terrain,
    streams: list[list[tuple[float, float]]],
    curves: list[RouteCurve],
    *,
    max_grade: float,
    min_radius_m: float,
    stream_offset_m: float,
    road_width_m: float,
    bound_pct: float | None,
    grade_smooth_m: float = GRADE_SMOOTH_M,
) -> dict[str, Any]:
    """비교값과 제약 못 지킨 구간 목록. 종단기울기 = 측점(20 m) 사이 고른 지반 기울기.

    `bound_pct` 는 제약을 다 지키는 길이 없을 때만 — 복도 안 어떤 길도 넘어야 하는 기울기(%).
    """
    grid = (terrain.y, terrain.x)
    sampler = RegularGridInterpolator(grid, terrain.z, bounds_error=False, fill_value=None)
    smooth = RegularGridInterpolator(
        grid, smoothed_ground(terrain, grade_smooth_m), bounds_error=False, fill_value=None
    )
    points, chain = _resample(vertices, SAMPLE_M)
    length = float(chain[-1])
    lookup = points[:, ::-1]
    ground = sampler(lookup)
    relaxed = bound_pct is not None
    terrain_reason = REASON_TERRAIN.format(bound=bound_pct or 0.0)
    violations: list[dict[str, Any]] = []

    # 종단기울기 — 측점 사이 고른 지반 기울기 · 넘은 측점 구간을 이어 한 항목
    stations = np.append(np.arange(0.0, length, STATION_M), length)
    station_z = np.interp(stations, chain, smooth(lookup))
    raw_z = np.interp(stations, chain, ground)
    spans = np.diff(stations)
    keep = spans > 0.5
    grades = np.abs(np.diff(station_z))[keep] / spans[keep]
    raw_grades = np.abs(np.diff(raw_z))[keep] / spans[keep]
    starts, ends = stations[:-1][keep], stations[1:][keep]
    for first, last in _runs(grades > max_grade + GRADE_EPS):
        violations.append(
            _item(
                "grade",
                float(starts[first]),
                float(ends[last]),
                float(grades[first : last + 1].max()) * 100.0,
                max_grade * 100.0,
                terrain_reason if relaxed else REASON_GRADE_FIT,
            )
        )

    # 최소 곡선반지름 — 원호가 하한 아래로 줄어든 곡선
    line = LineString(vertices)
    for curve in curves:
        if curve.violations:
            violations.append(
                _item(
                    "radius",
                    line.project(Point(curve.start)),
                    line.project(Point(curve.end)),
                    curve.radius_m,
                    min_radius_m,
                    REASON_RADIUS,
                )
            )

    # 계류 이격 — 버퍼 안 구간 중 봐줄 길이를 넘는 것. 봐줌 = 건너기(계류에 닿음) 한 번의
    # 버퍼 안 길이(45° 로 가로지름) + 기점 · 종점이 버퍼 안이면 45° 로 빠져나가는 길이
    route_points = shapely.points(points)
    if streams:
        near = shapely.distance(route_points, MultiLineString(streams))
        for first, last in _runs(near < stream_offset_m):
            ends_in = int(first == 0) + int(last == len(near) - 1)
            crosses = float(near[first : last + 1].min()) <= SAMPLE_M
            allowance = (ends_in + 2 * crosses) * stream_offset_m * math.sqrt(2.0)
            if chain[last] - chain[first] > allowance:
                violations.append(
                    _item(
                        "stream",
                        float(chain[first]),
                        float(chain[last]),
                        float(near[first : last + 1].min()),
                        stream_offset_m,
                        REASON_STREAM_TERRAIN if relaxed else REASON_STREAM,
                    )
                )

    # 절 · 성토(가늠) — 상한 안 계획고와 옆경사 단면을 폭 W 로 적분(비탈면 제외)
    grad_y, grad_x = np.gradient(terrain.z, _spacing(terrain.y), _spacing(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]
    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())

    offsets = shapely.distance(route_points, LineString(expected))
    counts = {
        kind: sum(1 for v in violations if v["kind"] == kind)
        for kind in ("grade", "radius", "stream")
    }
    violations.sort(key=lambda item: item["from_m"])
    return {
        "length_m": round(length, 2),
        "max_grade_pct": round(float(grades.max()) * 100.0 if grades.size else 0.0, 2),
        "avg_grade_pct": round(
            float((grades * spans[keep]).sum() / spans[keep].sum()) * 100.0 if grades.size else 0.0,
            2,
        ),
        # 참고 — 고르지 않은 DEM 측점 사이 최대 기울기(잔굴곡 포함 · 판정에 안 씀)
        "raw_max_grade_pct": round(float(raw_grades.max()) * 100.0 if raw_grades.size else 0.0, 2),
        "cut_m3": round(cut, 1),
        "fill_m3": round(fill, 1),
        "over_limit": {**counts, "total": len(violations)},
        "offset_avg_m": round(float(offsets.mean()), 2),
        "offset_max_m": round(float(offsets.max()), 2),
        "violations": violations,
    }


def _spacing(coords: np.ndarray) -> float:
    return float((coords[-1] - coords[0]) / (len(coords) - 1)) if len(coords) > 1 else 1.0