eomsangdonandClaude Opus 5.5
ebdb0610a6
- 균형 계획선 함수(서버 단독) — 종단기울기 상한 · 종단곡선(대수차) 지키며 B06 횡단 단면적으로 절 · 성토 균형 - 토공 최소 · 공사비 어림 점수를 균형 종단 토공(절 + 성 + 반입 + 사토)으로 - 추천 창에 균형 종단 한 줄 · 종단 그림 · 셈한 후보 표 · 지난 초 표시 - 균형 종단은 캐시에 안 적음 · 확정 체인 종단 규칙은 그대로 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EYJJRCfHfdCGNhmXTjQKV7
379 lines
17 KiB
Python
379 lines
17 KiB
Python
"""노선 추천 후보의 균형 종단(PLAN 63-4 · 63-2 안 A).
|
|
|
|
후보 노선 하나에 종단을 세워 절 · 성토를 **횡단 ㎥** 로 맞춘 계획선과 그 토량을 낸다.
|
|
편집 창 [노선 추천]의 점수 · 미리 보기에만 쓴다 — 파일 · 캐시 · initial_snapshot · 확정 종단을
|
|
바꾸지 않는다. 확정 체인의 종단 규칙(1차 계획고 = 지반고 · 2026-09-02 사용자 확정)은 그대로다.
|
|
|
|
계산 자리 = 서버 단독(③) — 화면은 받은 측점 · 지반고 · 계획고를 그리기만 한다.
|
|
|
|
방법:
|
|
1. 측점(20 m)마다 지반 횡단(±20 m · 0.5 m)을 DEM 에서 뜬다 — B05 측점 생성과 같은 꼴
|
|
(좌 = +offset · 상단측 = 좌우 평균 표고가 높은 쪽 · 그쪽 절토).
|
|
2. 측점마다 계획고 ±Δ 의 절 · 성 단면적 곡선 — B06 `compute_cross_design` 그대로(지반유형
|
|
리핑암 = B06 기본값) · `AREA_STEP_M` 간격으로 재고 사이는 직선 보간.
|
|
3. 동적계획(앞 두 측점 상태)으로 측점 계획고를 고른다 — 목표 = 절토(다짐 환산) + 성토 ㎥ +
|
|
λ × (성토 − 절토). λ 를 이분해 |절 − 성| / 큰 쪽이 허용(10 %) 안인 가장 작은 λ.
|
|
· 종단기울기 상한 — 넘지 못함.
|
|
· 종단곡선 — 측점마다 변화점(20 m)이라 곡선 반쪽이 직선의 0.45 를 못 넘어 L ≤ 18 m <
|
|
법정 최소 L(20 m 이상). 그래서 대수차가 곡선 생략 한도(비포장 5 %) 안이어야 함 — 넘으면
|
|
벌점(지형상 못 지킬 때만 넘음) · 넘은 변화점 수를 낸다.
|
|
4. 배관 덮개 — 제약은 끈 채(확정 체인과 같음) 관 자리 계획고 − 지반고 < 덮개 여유인 관 수만.
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import time
|
|
from pathlib import Path
|
|
from typing import TYPE_CHECKING, Any
|
|
|
|
import numpy as np
|
|
from scipy.interpolate import PchipInterpolator, RegularGridInterpolator
|
|
from shapely.geometry import LineString, Point
|
|
|
|
from B05_Profile.B05_Profile_Engine_RouteInitial_Metrics import COST_EARTH_M3
|
|
from config.config_system_design import (
|
|
EARTHWORK_CONVERSION_FACTORS,
|
|
FOREST_ROAD_PROFILE_ALIGNMENT,
|
|
FOREST_ROAD_PROFILE_CRITERIA,
|
|
)
|
|
|
|
if TYPE_CHECKING:
|
|
from B05_Profile.B05_Profile_Engine_RouteInitial_Search import Terrain
|
|
|
|
STATION_M = 20.0 # 측점 간격 — 추천 비교값(`..._Metrics.STATION_M`)과 같음
|
|
TAIL_MERGE_M = 5.0 # 종점과 이보다 가까운 정규 측점은 뺌(짧은 끝 구간의 기울기 튐 방지)
|
|
CROSS_HALF_M = 20.0 # 지반 횡단 반폭 · 간격 — B06 횡단 파일과 같음
|
|
CROSS_STEP_M = 0.5
|
|
LEVEL_RANGE_M = 6.0 # 계획고 후보 = 상한 안 기준선 ± 이 폭
|
|
LEVEL_STEP_M = 0.25
|
|
AREA_STEP_M = 2.0 # 단면적 곡선을 재는 계획고 간격(사이는 PCHIP — 1 m 간격과 3 % 안)
|
|
GROUND_TYPE = "ripping_rock" # B06 기본 지반유형(`attach_default_designs`)
|
|
CURVE_PENALTY_M3 = 500.0 # 대수차가 생략 한도를 1 %p 넘을 때 벌점(㎥)
|
|
LAMBDA_MAX = 0.95
|
|
LAMBDA_STEPS = 14
|
|
PIPE_SNAP_M = 5.0 # 관 자리가 노선에서 이 안이면 그 노선의 관으로 봄
|
|
GRADE_EPS = 1e-6
|
|
|
|
|
|
def route_sections(
|
|
vertices: list[tuple[float, float]], terrain: Terrain
|
|
) -> tuple[np.ndarray, np.ndarray, list[list[dict[str, Any]]], list[str]]:
|
|
"""측점 누가거리 · 중심 지반고 · 측점별 지반 횡단 샘플 · 단면유형(상단측 절토)."""
|
|
xy = np.asarray(vertices, dtype=np.float64)
|
|
keep = np.r_[True, np.hypot(*np.diff(xy, axis=0).T) > 1e-9]
|
|
xy = xy[keep]
|
|
seg = np.hypot(*np.diff(xy, axis=0).T)
|
|
chain = np.concatenate([[0.0], np.cumsum(seg)])
|
|
total = float(chain[-1])
|
|
at = np.arange(0.0, total, STATION_M)
|
|
if len(at) > 1 and total - at[-1] < TAIL_MERGE_M:
|
|
at = at[:-1]
|
|
at = np.append(at, total)
|
|
index = np.clip(np.searchsorted(chain, at, side="right") - 1, 0, len(seg) - 1)
|
|
tangent = np.diff(xy, axis=0)[index] / seg[index][:, None]
|
|
left = np.column_stack([-tangent[:, 1], tangent[:, 0]])
|
|
center = np.column_stack([np.interp(at, chain, xy[:, 0]), np.interp(at, chain, xy[:, 1])])
|
|
offsets = np.arange(-CROSS_HALF_M, CROSS_HALF_M + 1e-9, CROSS_STEP_M)
|
|
points = center[:, None, :] + left[:, None, :] * offsets[None, :, None]
|
|
lookup = points.reshape(-1, 2)[:, ::-1]
|
|
grid = (terrain.y, terrain.x)
|
|
z = RegularGridInterpolator(grid, terrain.z, bounds_error=False, fill_value=np.nan)(lookup)
|
|
valid = RegularGridInterpolator(
|
|
grid, terrain.valid.astype(np.float64), method="nearest", bounds_error=False, fill_value=0
|
|
)(lookup)
|
|
z = z.reshape(len(at), len(offsets))
|
|
ok = (valid.reshape(z.shape) > 0.5) & np.isfinite(z)
|
|
middle = int(np.argmin(np.abs(offsets)))
|
|
ground = np.where(ok[:, middle], z[:, middle], np.nan)
|
|
if np.isnan(ground).all():
|
|
raise ValueError("노선 위 지반고가 없어 종단을 세울 수 없습니다.")
|
|
fine = ~np.isnan(ground)
|
|
ground = np.interp(at, at[fine], ground[fine])
|
|
sections, modes = [], []
|
|
for row, flags in zip(z, ok):
|
|
sections.append(
|
|
[
|
|
{"offset_m": float(o), "elevation_m": float(e) if f else None, "valid": bool(f)}
|
|
for o, e, f in zip(offsets, row, flags)
|
|
]
|
|
)
|
|
lefts, rights = row[(offsets > 0) & flags], row[(offsets < 0) & flags]
|
|
up = "right" if lefts.size and rights.size and rights.mean() > lefts.mean() else "left"
|
|
modes.append(f"{up}_cut")
|
|
return at, ground, sections, modes
|
|
|
|
|
|
def _area_curves(
|
|
sections: list[list[dict[str, Any]]], modes: list[str], ground: np.ndarray, levels: np.ndarray
|
|
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
|
"""측점 · 계획고 후보마다 (절토 자연 ㎡, 절토 다짐 환산 ㎡, 성토 ㎡)."""
|
|
from B06_Section.B06_Section_Engine_Design import compute_cross_design
|
|
|
|
factors = EARTHWORK_CONVERSION_FACTORS
|
|
shape = levels.shape
|
|
cut, cut_c, fill = np.zeros(shape), np.zeros(shape), np.zeros(shape)
|
|
for i, (samples, mode) in enumerate(zip(sections, modes)):
|
|
delta = levels[i] - ground[i]
|
|
lo = np.floor(delta.min() / AREA_STEP_M) * AREA_STEP_M
|
|
at = np.arange(lo, delta.max() + AREA_STEP_M, AREA_STEP_M)
|
|
rows = []
|
|
for d in at:
|
|
try:
|
|
r = compute_cross_design(
|
|
samples, float(ground[i] + d), ground_type=GROUND_TYPE, section_mode=mode
|
|
)
|
|
except ValueError: # 지반 샘플이 모자란 측점 — 토량 모름(0)
|
|
rows.append((0.0, 0.0, 0.0))
|
|
continue
|
|
rock = factors.get(r.get("cut_rock_kind") or GROUND_TYPE, factors[GROUND_TYPE])
|
|
rows.append(
|
|
(
|
|
r["cut_area_m2"],
|
|
r["cut_soil_area_m2"] * factors["soil"]["compacted"]
|
|
+ r["cut_rock_area_m2"] * rock["compacted"],
|
|
r["fill_area_m2"],
|
|
)
|
|
)
|
|
table = np.asarray(rows)
|
|
# 단면적은 계획고에 대해 볼록(대략 2차)이라 직선 보간은 크게 잡힌다 — 단조 3차(PCHIP)
|
|
if len(at) > 1:
|
|
curve = PchipInterpolator(at, table, axis=0)(delta)
|
|
else:
|
|
curve = np.repeat(table, len(delta), axis=0)
|
|
cut[i], cut_c[i], fill[i] = curve[:, 0], curve[:, 1], curve[:, 2]
|
|
return cut, cut_c, fill
|
|
|
|
|
|
def _reference(ground: np.ndarray, spans: np.ndarray, max_grade: float) -> np.ndarray:
|
|
"""상한 안 기준선 — 앞 · 뒤로 기울기를 누른 두 선의 평균(`..._Metrics.grade_line` 꼴)."""
|
|
forward, backward = ground.copy(), ground.copy()
|
|
for i in range(1, len(ground)):
|
|
rise = max_grade * spans[i - 1]
|
|
forward[i] = min(max(ground[i], forward[i - 1] - rise), forward[i - 1] + rise)
|
|
for i in range(len(ground) - 2, -1, -1):
|
|
rise = max_grade * spans[i]
|
|
backward[i] = min(max(ground[i], backward[i + 1] - rise), backward[i + 1] + rise)
|
|
return 0.5 * (forward + backward)
|
|
|
|
|
|
def _dp(
|
|
levels: np.ndarray, cost: np.ndarray, spans: np.ndarray, max_grade: float, delta_max: float
|
|
) -> np.ndarray | None:
|
|
"""측점별 계획고 후보 번호 — 기울기 상한은 막고 대수차 넘침은 벌점. 길이 없으면 None."""
|
|
n = len(levels)
|
|
if n == 2:
|
|
steep = np.abs(levels[1][None, :] - levels[0][:, None]) / spans[0] > max_grade + GRADE_EPS
|
|
table = cost[0][:, None] + cost[1][None, :] + np.where(steep, np.inf, 0.0)
|
|
if not np.isfinite(table.min()):
|
|
return None
|
|
return np.array(np.unravel_index(int(table.argmin()), table.shape))
|
|
steep = np.abs(levels[1][None, :] - levels[0][:, None]) / spans[0] > max_grade + GRADE_EPS
|
|
table = cost[0][:, None] + cost[1][None, :] + np.where(steep, np.inf, 0.0) # [앞앞, 앞]
|
|
back = []
|
|
for i in range(1, n - 1):
|
|
g_in = (levels[i][None, :] - levels[i - 1][:, None]) / spans[i - 1] # [a, b]
|
|
g_out = (levels[i + 1][None, :] - levels[i][:, None]) / spans[i] # [b, c]
|
|
over = np.maximum(np.abs(g_out[None, :, :] - g_in[:, :, None]) - delta_max, 0.0)
|
|
total = table[:, :, None] + CURVE_PENALTY_M3 * 100.0 * over # [a, b, c]
|
|
arg = total.argmin(axis=0)
|
|
table = np.take_along_axis(total, arg[None], axis=0)[0]
|
|
table += np.where(np.abs(g_out) > max_grade + GRADE_EPS, np.inf, 0.0) + cost[i + 1][None, :]
|
|
back.append(arg)
|
|
if not np.isfinite(table.min()):
|
|
return None
|
|
b, c = np.unravel_index(int(table.argmin()), table.shape)
|
|
path = [int(c), int(b)]
|
|
for arg in reversed(back):
|
|
path.append(int(arg[path[-1], path[-2]]))
|
|
return np.array(path[::-1])
|
|
|
|
|
|
def balanced_profile(
|
|
chainage: np.ndarray,
|
|
ground: np.ndarray,
|
|
sections: list[list[dict[str, Any]]],
|
|
modes: list[str],
|
|
*,
|
|
max_grade: float,
|
|
design_speed: int = 20,
|
|
) -> dict[str, Any]:
|
|
"""상한(종단기울기 · 종단곡선)을 지키며 절 · 성토를 횡단 ㎥ 로 맞춘 계획선과 토량.
|
|
|
|
시점 · 종점은 지반고(확정 체인 기본 오프셋 0 과 같음) — 그 길이 없을 때만 풀어 다시 찾는다.
|
|
"""
|
|
started = time.perf_counter()
|
|
spans = np.diff(chainage)
|
|
weight = np.zeros(len(chainage))
|
|
weight[:-1] += spans / 2.0
|
|
weight[1:] += spans / 2.0
|
|
offsets = np.arange(-LEVEL_RANGE_M, LEVEL_RANGE_M + 1e-9, LEVEL_STEP_M)
|
|
levels = _reference(ground, spans, max_grade)[:, None] + offsets[None, :]
|
|
cut, cut_c, fill = (weight[:, None] * a for a in _area_curves(sections, modes, ground, levels))
|
|
criteria = (
|
|
FOREST_ROAD_PROFILE_CRITERIA["design_speed"].get(design_speed)
|
|
or (FOREST_ROAD_PROFILE_CRITERIA["design_speed"][20])
|
|
)
|
|
skip = FOREST_ROAD_PROFILE_CRITERIA["vertical_curve_skip_delta_pct"] / 100.0
|
|
room = 2.0 * FOREST_ROAD_PROFILE_ALIGNMENT["curve_tangent_max_ratio"] * float(spans.min())
|
|
fits = criteria["min_curve_length_m"] <= room
|
|
delta_max = max(skip, room / criteria["min_vertical_radius_m"]) if fits else skip
|
|
tolerance = FOREST_ROAD_PROFILE_ALIGNMENT["balance_tolerance_percent"] / 100.0
|
|
pinned = np.full(levels.shape, np.inf)
|
|
for end in (0, -1):
|
|
pinned[end, int(np.argmin(np.abs(levels[end] - ground[end])))] = 0.0
|
|
pinned[1:-1] = 0.0
|
|
|
|
def solve(lam: float) -> dict[str, Any] | None:
|
|
base = cut_c + fill + lam * (fill - cut_c)
|
|
path = _dp(levels, base + pinned, spans, max_grade, delta_max)
|
|
if path is None:
|
|
path = _dp(levels, base, spans, max_grade, delta_max)
|
|
if path is None:
|
|
return None
|
|
pick = np.arange(len(path))
|
|
c, f = float(cut_c[pick, path].sum()), float(fill[pick, path].sum())
|
|
return {
|
|
"lambda": lam,
|
|
"path": path,
|
|
"cut": float(cut[pick, path].sum()),
|
|
"cut_c": c,
|
|
"fill": f,
|
|
"imbalance": (f - c) / max(c, f, 1e-9),
|
|
}
|
|
|
|
tried = [solve(0.0)]
|
|
if tried[0] is None:
|
|
raise ValueError("종단기울기 상한 안에서 계획고를 고를 수 없습니다.")
|
|
first = tried[0]
|
|
if abs(first["imbalance"]) > tolerance:
|
|
sign = 1.0 if first["imbalance"] > 0 else -1.0
|
|
low, high = 0.0, LAMBDA_MAX
|
|
for _ in range(LAMBDA_STEPS):
|
|
mid = 0.5 * (low + high)
|
|
result = solve(sign * mid)
|
|
if result is None:
|
|
break
|
|
tried.append(result)
|
|
if abs(result["imbalance"]) <= tolerance or result["imbalance"] * sign < 0:
|
|
high = mid
|
|
else:
|
|
low = mid
|
|
within = [r for r in tried if abs(r["imbalance"]) <= tolerance]
|
|
best = (
|
|
min(within, key=lambda r: abs(r["lambda"]))
|
|
if within
|
|
else min(tried, key=lambda r: abs(r["imbalance"]))
|
|
)
|
|
plan = levels[np.arange(len(chainage)), best["path"]]
|
|
grades = np.diff(plan) / spans
|
|
changes = np.abs(np.diff(grades))
|
|
return {
|
|
"stations_m": [round(float(v), 2) for v in chainage],
|
|
"ground_m": [round(float(v), 3) for v in ground],
|
|
"plan_m": [round(float(v), 3) for v in plan],
|
|
"cut_m3": round(best["cut"], 1), # 자연 상태
|
|
"cut_compacted_m3": round(best["cut_c"], 1),
|
|
"fill_m3": round(best["fill"], 1),
|
|
"borrow_m3": round(max(best["fill"] - best["cut_c"], 0.0), 1),
|
|
"spoil_m3": round(max(best["cut_c"] - best["fill"], 0.0), 1),
|
|
"imbalance_pct": round(best["imbalance"] * 100.0, 1),
|
|
"balanced": abs(best["imbalance"]) <= tolerance,
|
|
"tolerance_pct": round(tolerance * 100.0, 1),
|
|
"max_grade_pct": round(float(np.abs(grades).max()) * 100.0 if grades.size else 0.0, 2),
|
|
"grade_limit_pct": round(max_grade * 100.0, 2),
|
|
"curve_over": int((changes > delta_max + GRADE_EPS).sum()),
|
|
"curve_delta_max_pct": round(delta_max * 100.0, 2),
|
|
"max_cut_depth_m": round(float(np.max(ground - plan)), 2),
|
|
"max_fill_height_m": round(float(np.max(plan - ground)), 2),
|
|
"seconds": round(time.perf_counter() - started, 2),
|
|
}
|
|
|
|
|
|
def route_profile(
|
|
vertices: list[tuple[float, float]], terrain: Terrain, *, max_grade: float, design_speed: int
|
|
) -> dict[str, Any]:
|
|
"""후보 노선 하나의 균형 종단 — 측점 횡단을 뜨고 `balanced_profile`."""
|
|
chainage, ground, sections, modes = route_sections(vertices, terrain)
|
|
if len(chainage) < 2:
|
|
raise ValueError("노선이 짧아 종단을 세울 수 없습니다.")
|
|
return balanced_profile(
|
|
chainage, ground, sections, modes, max_grade=max_grade, design_speed=design_speed
|
|
)
|
|
|
|
|
|
def cover_shortfall(
|
|
profile: dict[str, Any],
|
|
vertices: list[tuple[float, float]],
|
|
pipes: list[tuple[float, float, float]],
|
|
) -> dict[str, int]:
|
|
"""관(x, y, 덮개 여유 m) 중 이 노선 위(`PIPE_SNAP_M` 안)인 것 · 그중 덮개가 모자란 수.
|
|
|
|
덮개 = 관 자리 계획고 − 지반고(측점 사이 직선 보간). 제약으로 쓰지 않고 세기만 한다."""
|
|
line = LineString(vertices)
|
|
stations = np.asarray(profile["stations_m"])
|
|
rise = np.asarray(profile["plan_m"]) - np.asarray(profile["ground_m"])
|
|
on, short = 0, 0
|
|
for x, y, clearance in pipes:
|
|
point = Point(x, y)
|
|
if line.distance(point) > PIPE_SNAP_M:
|
|
continue
|
|
on += 1
|
|
if float(np.interp(line.project(point), stations, rise)) < clearance:
|
|
short += 1
|
|
return {"pipes_on_route": on, "cover_short": short}
|
|
|
|
|
|
def project_pipes(project_root: Path) -> list[tuple[float, float, float]]:
|
|
"""사업지에 놓인 관(x, y, 덮개 여유 m) — 배수유역 편집분을 읽기만 한다(좌표 없는 관은 뺌)."""
|
|
from common_util.common_util_drainage_pipes import facility_clearance_m, read_pipe_points_file
|
|
from config.config_system import (
|
|
DRAINAGE_CACHE_DIRNAME,
|
|
DRAINAGE_EDITS_DIRNAME,
|
|
DRAINAGE_PIPE_POINTS_FILENAME,
|
|
)
|
|
|
|
path = (
|
|
project_root
|
|
/ "B04_PreProcess"
|
|
/ DRAINAGE_CACHE_DIRNAME
|
|
/ DRAINAGE_EDITS_DIRNAME
|
|
/ DRAINAGE_PIPE_POINTS_FILENAME
|
|
)
|
|
return [
|
|
(float(p.x), float(p.y), facility_clearance_m(p.facility, p.options))
|
|
for p in read_pipe_points_file(path)
|
|
if p.x is not None and p.y is not None
|
|
]
|
|
|
|
|
|
def attach_profile(
|
|
metrics: dict[str, Any],
|
|
vertices: list[tuple[float, float]],
|
|
terrain: Terrain,
|
|
*,
|
|
max_grade: float,
|
|
design_speed: int,
|
|
) -> None:
|
|
"""후보 비교값에 균형 종단을 붙이고 토공 점수를 그것으로 바꾼다(한 번만 · 묶음 캐시에 남음).
|
|
|
|
점수 = 절토(자연) + 성토 + 반입 + 사토 ㎥ — 가늠값(폭 4 m 띠 · 비탈면 뺌)은
|
|
`earth_estimate_m3` 로 남기고, 공사비 어림의 토공 몫도 이 점수로 바꾼다."""
|
|
if "profile" in metrics:
|
|
return
|
|
try:
|
|
profile = route_profile(vertices, terrain, max_grade=max_grade, design_speed=design_speed)
|
|
except ValueError as error:
|
|
metrics["profile"] = {"error": str(error)}
|
|
return
|
|
estimate = metrics["cut_m3"] + metrics["fill_m3"]
|
|
metrics["profile"] = profile
|
|
metrics["earth_estimate_m3"] = round(estimate, 1)
|
|
metrics["cost_index"] = round(
|
|
metrics["cost_index"] + COST_EARTH_M3 * (profile_earth(profile) - estimate)
|
|
)
|
|
|
|
|
|
def profile_earth(profile: dict[str, Any]) -> float:
|
|
"""균형 종단 토공 점수 — 절토(자연) + 성토 + 반입 + 사토 ㎥."""
|
|
return profile["cut_m3"] + profile["fill_m3"] + profile["borrow_m3"] + profile["spoil_m3"]
|