Building/calc/fig-laba2-plan.py

# -*- coding: utf-8 -*-
"""
Лаба 2, рисунок «План трассы». Горизонтали векторизуются из скана варианта 6
(.work/laba2-contours-paths.pkl, координаты в пикселях кадра карты),
затем на карту наносится трасса А–Б: тангенсы с круговыми кривыми,
пикетаж, тоннель у водораздела.

Схема трассы задаётся списком вершин углов поворота ВУ (пиксели кадра) и
радиусами кривых; всё остальное (пикетаж, пересечения горизонталей, элементы
кривых) вычисляется. Данные профиля сохраняются в .work/laba2-profile-data.pkl.
"""
import numpy as np, pickle, json, math

with open('.work/laba2-contours-paths.pkl', 'rb') as f:
    raw = pickle.load(f)

# --- 1. Склейка разорванных горизонталей и упорядочивание -------------------
# (пути 11+10 и 12+13 разорваны маркером пункта А)
contours = [
    np.vstack([raw[11], raw[10]]),   # 250 (через пункт А)
    np.vstack([raw[12], raw[13]]),   # 240 (нижняя)
] + [raw[i] for i in range(10)]      # остальные, с юга на север: p9..p0

# юг -> север по положению на западном крае (x < 170)
def west_y(p):
    m = (p[:, 0] >= 20) & (p[:, 0] <= 170)
    return np.mean(p[m, 1]) if m.any() else p[:, 1].mean()

contours.sort(key=west_y, reverse=True)          # южные первыми (y вниз)
# южный склон 240->290, водораздел 300 (седьмая линия), северный склон 290->250
ELEV = [240, 250, 260, 270, 280, 290, 300, 290, 280, 270, 260, 250]

# проверка унимодальности отметок на вертикальных разрезах
for xq in (60, 200, 400, 600, 760):
    seq = []
    for p, z in zip(contours, ELEV):
        m = np.abs(p[:, 0] - xq) < 30
        if m.any():
            seq.append((np.mean(p[m, 1]), z))
    seq.sort()
    zs = [z for _, z in seq]
    peak = zs.index(max(zs))
    ok = zs[:peak+1] == sorted(zs[:peak+1]) and zs[peak:] == sorted(zs[peak:], reverse=True)
    print(f'x={xq}: {zs} {"OK" if ok else "!! НЕ унимодально"}')

# --- 2. Масштаб --------------------------------------------------------------
PX_PER_KM = 262.7           # кадр карты 788 px = 3.0 км (масштаб 1:10 000)
M_PER_PX = 1000.0 / PX_PER_KM

# --- 3. Сглаживание горизонталей --------------------------------------------
def smooth(p, k=7):
    if len(p) < 2 * k + 1:
        return p
    ker = np.hanning(k)
    ker /= ker.sum()
    sm = np.empty_like(p)
    for c in (0, 1):
        sm[:, c] = np.convolve(p[:, c], ker, mode='same')
    sm[:k//2] = p[:k//2]; sm[-(k//2):] = p[-(k//2):]
    return sm

contours = [smooth(p) for p in contours]

# --- 4. Трасса: тангенсы и круговые кривые ----------------------------------
# Вершины углов поворота (пиксели кадра, y вниз); НТ = начало трассы (А),
# КТ = конец трассы (Б). (A, R): None для первого/последнего участка
VUS = [
    ((260.0, 402.0), None),    # НТ = А
    ((150.0, 372.0), 300.0),   # ВУ1
    ((105.0, 255.0), 600.0),   # ВУ2
    ((55.0, 165.0), 250.0),    # ВУ3
    ((18.0, 142.0), None),     # КТ = Б
]

def build_alignment(vus):
    pts = []          # плотная полилиния оси (px)
    pieces = []       # (тип, геометрия, длина_м)
    P = [np.array(v, float) for v, _ in vus]
    R = [r for _, r in vus]
    elem = []         # элементы кривых для надписей
    # 1) тангенсы и дуги
    segs = []         # последовательность ('l', A, B) | ('a', C, A, B, R, ccw)
    cur = P[0]
    for i in range(1, len(P) - 1):
        V, Pn = P[i], P[i + 1]
        u1 = (V - cur) / np.linalg.norm(V - cur)
        u2 = (Pn - V) / np.linalg.norm(Pn - V)
        cosA = np.clip(np.dot(u1, u2), -1, 1)
        d = math.acos(cosA)                       # угол поворота
        Rp = R[i] / M_PER_PX                      # радиус в px
        T = Rp * math.tan(d / 2)
        A = V - u1 * T
        B = V + u2 * T
        cross = u1[0] * u2[1] - u1[1] * u2[0]     # >0: поворот по часовой (y вниз)
        n = np.array([-u1[1], u1[0]]) * np.sign(cross)
        C = A + n * Rp
        segs.append(('l', cur, A))
        segs.append(('a', C, A, B, Rp, cross < 0))
        elem.append({'i': i, 'V': V, 'A': A, 'B': B, 'C': C, 'R_px': Rp,
                     'd': math.degrees(d), 'T_px': T,
                     'K_m': Rp * d * M_PER_PX, 'T_m': T * M_PER_PX,
                     'R_m': R[i], 'right': cross < 0})
        cur = B
    segs.append(('l', cur, P[-1]))
    # 2) плотная выборка и длины
    total_px = 0.0
    dense = []
    for s in segs:
        if s[0] == 'l':
            _, a, b = s
            L = np.linalg.norm(b - a)
            n = max(2, int(L))
            t = np.linspace(0, 1, n)
            pts.append(a[None, :] + t[:, None] * (b - a)[None, :])
            total_px += L
        else:
            _, C, a, b, Rp, ccw = s
            a0 = math.atan2(a[1] - C[1], a[0] - C[0])
            a1 = math.atan2(b[1] - C[1], b[0] - C[0])
            da = a1 - a0
            if ccw:
                while da > 0: da -= 2 * math.pi
            else:
                while da < 0: da += 2 * math.pi
            n = max(4, int(abs(da) * Rp))
            t = np.linspace(a0, a0 + da, n)
            pts.append(C[None, :] + Rp * np.column_stack([np.cos(t), np.sin(t)]))
            total_px += abs(da) * Rp
    align = np.vstack(pts)
    s_m = np.r_[0, np.cumsum(np.hypot(*np.diff(align, axis=0).T))] * M_PER_PX
    return align, s_m, segs, elem

align, s_m, segs, curves = build_alignment(VUS)
L_total = s_m[-1]
print(f'длина трассы: {L_total:.0f} м; кривых: {len(curves)}')
for c in curves:
    print(f"  ВУ{c['i']}: У={c['d']:.1f}° R={c['R_m']} T={c['T_m']:.1f} K={c['K_m']:.1f} "
          f"{'вправо' if c['right'] else 'влево'}")

# --- 5. Пересечения трассы с горизонталями ----------------------------------
from scipy.spatial import cKDTree
trees = [cKDTree(p) for p in contours]

def elev_along(points):
    """ближайшая горизонталь -> отметка (ступенчатая), для контроля"""
    z = np.empty(len(points))
    for j, tree in enumerate(trees):
        pass
    d = np.stack([tree.query(points)[0] for tree in trees])
    return np.array(ELEV)[np.argmin(d, axis=0)]

crossings = []   # (s, z)
zone_prev = False
idx = trees_ind = np.zeros(len(align), int)
D = np.stack([tree.query(align)[0] for tree in trees])
nearest = np.argmin(D, axis=0)
near = D.min(axis=0)
inzone = near < 1.1
runs = []
i = 0
while i < len(align):
    if inzone[i]:
        j = i
        while j < len(align) and inzone[j]:
            j += 1
        runs.append((i, j))
        i = j
    else:
        i += 1
for a, b in runs:
    mid = (a + b) // 2
    crossings.append((s_m[mid], ELEV[nearest[mid]], a, b))
print('пересечения с горизонталями:')
for s, z, a, b in crossings:
    print(f'  s={s:7.1f} м  z={z} м  (px {a}..{b})')

# --- 6. Поверхность земли по пересечениям -----------------------------------
# отметки на концах трассы: А стоит на горизонтали 250, Б — между 280 и 270
S_A_Z, S_B_Z = 250.0, 277.5

def ground(s):
    cs = [0.0] + [c[0] for c in crossings] + [L_total]
    zs = [S_A_Z] + [c[1] for c in crossings] + [S_B_Z]
    return float(np.interp(s, cs, zs))

# --- 7. Проектная (красная) линия: участки уклонов ---------------------------
# (s_начало, уклон ‰): 2 выемки у порталов тоннеля 1060->1265
GRADES = [(0.0, 27.0), (400.0, 45.0), (700.0, 10.0), (1060.0, 10.0), (1265.0, -10.0)]
S_P1, S_P2 = 1060.0, 1265.0     # порталы тоннеля

def red_h(s):
    h = 250.0
    for (s1, i1), (s2, _i2) in zip(GRADES, GRADES[1:]):
        if s1 <= s < s2:
            return h + (s - s1) * i1 / 1000.0
        h += (s2 - s1) * i1 / 1000.0
    return h + (s - GRADES[-1][0]) * GRADES[-1][1] / 1000.0

# --- 8. Сохранение данных ---------------------------------------------------
with open('.work/laba2-plan-data.pkl', 'wb') as f:
    pickle.dump({
        'contours': contours, 'ELEV': ELEV, 'align': align, 's_m': s_m,
        'curves': curves, 'crossings': crossings, 'M_PER_PX': M_PER_PX,
        'A': np.array(VUS[0][0]), 'B': np.array(VUS[-1][0]),
        'GRADES': GRADES, 'S_P1': S_P1, 'S_P2': S_P2,
        'S_A_Z': S_A_Z, 'S_B_Z': S_B_Z, 'L_total': L_total,
    }, f)
print('данные сохранены: .work/laba2-plan-data.pkl')

# --- 9. Чертёж плана --------------------------------------------------------
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe
from matplotlib.patches import Circle

# Чертёжный шрифт ГОСТ 2.304-81 тип А (прямой, из системных шрифтов Windows)
import matplotlib.font_manager as fm
fm.fontManager.addfont('C:/Windows/Fonts/GOST2304A.ttf')
matplotlib.rcParams['font.family'] = 'GOST 2.304 type A'
matplotlib.rcParams['mathtext.default'] = 'regular'

RED = '#cc1f2d'
W, H_PX = 788.0, 468.0          # кадр карты, px
ML, MR, MT, MB = 12.0, 12.0, 62.0, 96.0   # поля листа вокруг рамки карты

fig = plt.figure(figsize=((W + ML + MR) / 50, (H_PX + MT + MB) / 50), dpi=200)
ax = fig.add_axes([0, 0, 1, 1])
ax.set_xlim(-ML, W + MR)
ax.set_ylim(H_PX + MB, -MT)     # y вниз, как на плане
ax.axis('off')
ax.set_facecolor('white')
fig.patches.append(plt.Rectangle((0, 0), 1, 1, transform=fig.transFigure,
                                 fill=False, ec='black', lw=1.4))

# рамка карты
ax.add_patch(plt.Rectangle((0, 0), W, H_PX, fill=False, ec='black', lw=1.1))

def m2px(pt):
    return np.asarray(pt, float)

# километровая сетка через 0,5 км
step = 500.0 / M_PER_PX
xg = np.arange(step, W, step)
yg = np.arange(step, H_PX, step)
for x in xg:
    ax.plot([x, x], [0, H_PX], color='0.55', lw=0.4, alpha=0.55, zorder=0.5)
for y in yg:
    ax.plot([0, W], [y, y], color='0.55', lw=0.4, alpha=0.55, zorder=0.5)

# горизонтали с подписями отметок
for p, z in zip(contours, ELEV):
    ax.plot(p[:, 0], p[:, 1], color='black', lw=0.75, zorder=1)
    for xq in (120, 360, 600, 740):
        i = np.argmin(np.abs(p[:, 0] - xq))
        if abs(p[i, 0] - xq) > 90 or 4 < p[i, 1] < H_PX - 4:
            pass
        ax.text(p[i, 0] + 3, p[i, 1] - 4, str(z), fontsize=9.5, color='black',
                zorder=3, path_effects=[pe.withStroke(linewidth=2.4, foreground='white')])

# трасса: участки вне тоннеля
def s_to_pts(s1, s2):
    i1 = np.searchsorted(s_m, s1)
    i2 = min(len(s_m) - 1, np.searchsorted(s_m, s2))
    return align[i1:i2 + 1]

tun = s_to_pts(S_P1, S_P2)
s_tun = s_m[np.searchsorted(s_m, S_P1):np.searchsorted(s_m, S_P2) + 1]
open1 = s_to_pts(0, S_P1)
open2 = s_to_pts(S_P2, s_m[-1])
for seg in (open1, open2):
    ax.plot(seg[:, 0], seg[:, 1], color=RED, lw=2.0, zorder=4, solid_capstyle='round')
# тоннель: две параллельные линии со штрихами
d = np.gradient(tun, axis=0)
nrm = np.column_stack([d[:, 1], -d[:, 0]])
nrm /= np.hypot(nrm[:, 0], nrm[:, 1])[:, None] + 1e-9
for sgn in (+1, -1):
    off = tun + sgn * 3.2 * nrm
    ax.plot(off[:, 0], off[:, 1], color=RED, lw=1.3, zorder=4)
for s in (S_P1, S_P2):
    k = np.searchsorted(s_tun, s)
    base = tun[min(k, len(tun) - 1)]
    ax.plot([base[0] - 6.5 * nrm[k, 0], base[0] + 6.5 * nrm[k, 0]],
            [base[1] - 6.5 * nrm[k, 1], base[1] + 6.5 * nrm[k, 1]],
            color=RED, lw=3.2, zorder=4)
ax.text(96, 180, 'Тоннель, 205 м',
        fontsize=11, color=RED, zorder=6, path_effects=[pe.withStroke(linewidth=2.6, foreground='white')])

# пикеты: each 100 м, подпись each 2
pk = 0
while pk * 100.0 <= s_m[-1]:
    s = pk * 100.0
    k = np.searchsorted(s_m, s)
    if k >= len(s_m) - 1:
        break
    base = align[k]
    tang = align[min(k + 2, len(s_m) - 1)] - align[max(k - 2, 0)]
    tang = tang / (np.hypot(*tang) + 1e-9)
    nr = np.array([-tang[1], tang[0]])   # вправо по ходу (y вниз)
    L = 7.0 if pk % 2 == 0 else 4.5
    ax.plot([base[0], base[0] + L * nr[0]], [base[1], base[1] + L * nr[1]],
            color=RED, lw=1.4, zorder=4)
    if pk % 2 == 0:
        ax.text(base[0] + (L + 2.5) * nr[0], base[1] + (L + 2.5) * nr[1], str(pk),
                fontsize=9.5, color=RED, ha='center', va='center', zorder=5,
                path_effects=[pe.withStroke(linewidth=2.2, foreground='white')])
    pk += 1

# километровые знаки
for km, s in ((1, 1000.0),):
    k = np.searchsorted(s_m, s)
    base = align[k]
    tang = align[min(k + 2, len(s_m) - 1)] - align[max(k - 2, 0)]
    tang /= np.hypot(*tang)
    nr = np.array([tang[1], -tang[0]])   # влево по ходу
    c = base + 20 * nr
    ax.add_patch(Circle(c, 8.5, fill=True, fc='white', ec=RED, lw=1.3, zorder=5))
    ax.text(c[0], c[1], str(km), fontsize=11, color=RED, ha='center', va='center', zorder=6)

# пункты А и Б
for pt, name, dx, dy in ((np.array(VUS[0][0]), 'А', 10, 14), (np.array(VUS[-1][0]), 'Б', 9, -2)):
    ax.plot(pt[0], pt[1], marker='o', ms=5, color='black', zorder=6)
    ax.text(pt[0] + dx, pt[1] + dy, name, fontsize=15, weight='normal', color='black', zorder=6,
            path_effects=[pe.withStroke(linewidth=2.4, foreground='white')])

# вершины углов поворота
for c in curves:
    ax.plot(c['V'][0], c['V'][1], marker='+', ms=7, color=RED, mew=1.1, zorder=5)
    ax.text(c['V'][0] + 5, c['V'][1] - 5, f"ВУ{c['i']}", fontsize=8.5, color=RED, zorder=6,
            path_effects=[pe.withStroke(linewidth=2.2, foreground='white')])

# выноски элементов кривых
def fmt_deg(x):
    dd = int(x); mm = int(round((x - dd) * 60))
    return f"{dd}°{mm:02d}'" if mm else f'{dd}°'

slots = {1: (240, 470), 2: (300, 290), 3: (140, 210)}
for c in curves:
    tx, ty = slots[c['i']]
    mid = c['A'] * 0.35 + c['B'] * 0.65
    # начало линии под рамкой выноски: конец скрывается белой подложкой бокса
    ax.annotate('', xy=(mid[0], mid[1]), xytext=(tx + 60, ty + 6),
                arrowprops=dict(arrowstyle='-', color=RED, lw=0.7, shrinkA=0, shrinkB=0))
    ax.text(tx, ty, f"Кривая {c['i']}: У = {fmt_deg(c['d'])}; R = {c['R_m']:.0f} м\n"
                    f"Т = {c['T_m']:.1f} м; К = {c['K_m']:.1f} м".replace('.', ','),
            fontsize=9.5, color=RED, va='top', zorder=6,
            bbox=dict(boxstyle='square,pad=0.25', fc='white', ec='0.6', lw=0.5, alpha=0.92))

# стрелка севера
nx, ny = W - 34, 34
ax.annotate('', xy=(nx, ny - 26), xytext=(nx, ny + 8),
            arrowprops=dict(arrowstyle='-|>', color='black', lw=1.4, mutation_scale=16))
ax.text(nx, ny + 20, 'С', fontsize=12, ha='center', weight='normal')

# заголовок и подписи под рамкой
ax.text((W + ML + MR) / 2 - ML, -34, 'ПЛАН ТРАССЫ АВТОМОБИЛЬНОЙ ДОРОГИ А–Б',
        fontsize=15.5, weight='normal', ha='center', va='center')
ax.text((W + ML + MR) / 2 - ML, -13, '(трассирование по плану в горизонталях, вариант 6)',
        fontsize=11, ha='center', va='center', color='0.25')

# условные обозначения (двумя рядами) и масштаб (справа)
ly = H_PX + 26
ax.text(6, ly - 16, 'Условные обозначения:', fontsize=10.5, weight='normal')
ax.text(W - 4, ly - 16, 'Масштаб 1:10 000', fontsize=11.5, ha='right', weight='normal')
# ряд 1
ax.plot([10, 52], [ly + 2, ly + 2], color=RED, lw=2.0)
ax.text(58, ly + 2, 'трасса (принятый вариант)', fontsize=10, va='center')
ax.plot([268, 298], [ly - 2, ly - 2], color=RED, lw=1.1)
ax.plot([268, 298], [ly + 4, ly + 4], color=RED, lw=1.1)
ax.text(304, ly + 2, 'тоннель', fontsize=10, va='center')
ax.plot([392, 392], [ly - 4, ly + 6], color=RED, lw=1.4)
ax.text(399, ly + 2, 'пикет', fontsize=10, va='center')
ax.add_patch(Circle((475, ly + 2), 7, fc='white', ec=RED, lw=1.2))
ax.text(475, ly + 2, '1', fontsize=9.5, color=RED, ha='center', va='center')
ax.text(487, ly + 2, 'километровый знак', fontsize=10, va='center')
# ряд 2
ax.plot([10, 60], [ly + 26, ly + 26], color='black', lw=0.75)
ax.text(66, ly + 26, 'горизонталь и её отметка, м', fontsize=10, va='center')
ax.text(W - 4, ly + 26, 'горизонтали проведены через 10 м; координатная сетка через 500 м; пикеты через 100 м',
        fontsize=10, ha='right')

fig.savefig('assets/images/laba2-plan-trassy.png', dpi=200, facecolor='white')
print('сохранено: assets/images/laba2-plan-trassy.png')