LoginSignup
0
0
@cartooh

標高とカーブ半径（推定）のグラフ化

Last updated at Posted at 2023-07-01

目的

地図上で経路を引いたときにカーブ半径と高低差を見たかった。
（そういうのが気になるお年頃ｗ）

断面図｜地理院地図の使い方を見ていて、高低差は見れるかなと思って、

を使っていたところ、グラフの保存でCSVで保存できるではないですか。

なら一緒に見れるものを作ろう。
（カーブ半径の推定は、どうしても変なところとかにフィッティングしてしまうので、経路図を見ながら見たいのもあったので）

特徴

  • ipyleaflet で地図はインタラクティブに確認
  • 気になる標高、カーブ半径（グラフ）をクリックするとそこを表示

実装

関数定義等

は説明省略

%matplotlib widget 
import ipywidgets as widgets
import matplotlib.pyplot as plt
import japanize_matplotlib
import pandas as pd
import io
import ipyleaflet
import math
import numpy as np
from matplotlib.colors import to_rgb
from IPython.display import display
from collections.abc import Iterable


def calc_xy(phi_deg, lambda_deg, phi0_deg, lambda0_deg):
    """ 緯度経度を平面直角座標に変換する
    - input:
        (phi_deg, lambda_deg): 変換したい緯度・経度[度]（分・秒でなく小数であることに注意）
        (phi0_deg, lambda0_deg): 平面直角座標系原点の緯度・経度[度]（分・秒でなく小数であることに注意）
    - output:
        x: 変換後の平面直角座標[m]
        y: 変換後の平面直角座標[m]
    """
    # 緯度経度・平面直角座標系原点をラジアンに直す
    phi_rad     = np.deg2rad(phi_deg)
    lambda_rad  = np.deg2rad(lambda_deg)
    phi0_rad    = np.deg2rad(phi0_deg)
    lambda0_rad = np.deg2rad(lambda0_deg)

    # 補助関数
    def A_array(n):
        A0 = 1 + (n**2)/4. + (n**4)/64.
        A1 = -     (3./2)*( n - (n**3)/8. - (n**5)/64. ) 
        A2 =     (15./16)*( n**2 - (n**4)/4. )
        A3 = -   (35./48)*( n**3 - (5./16)*(n**5) )
        A4 =   (315./512)*( n**4 )
        A5 = -(693./1280)*( n**5 )
        return np.array([A0, A1, A2, A3, A4, A5])

    def alpha_array(n):
        a0 = np.nan # dummy
        a1 = (1./2)*n - (2./3)*(n**2) + (5./16)*(n**3) + (41./180)*(n**4) - (127./288)*(n**5)
        a2 = (13./48)*(n**2) - (3./5)*(n**3) + (557./1440)*(n**4) + (281./630)*(n**5)
        a3 = (61./240)*(n**3) - (103./140)*(n**4) + (15061./26880)*(n**5)
        a4 = (49561./161280)*(n**4) - (179./168)*(n**5)
        a5 = (34729./80640)*(n**5)
        return np.array([a0, a1, a2, a3, a4, a5])

    # 定数 (a, F: 世界測地系-測地基準系1980（GRS80）楕円体)
    m0 = 0.9999 
    a = 6378137.
    F = 298.257222101

    # (1) n, A_i, alpha_iの計算
    n = 1. / (2*F - 1)
    A_array = A_array(n)
    alpha_array = alpha_array(n)

    # (2), S, Aの計算
    A_ = ( (m0*a)/(1.+n) )*A_array[0] # [m]
    S_ = ( (m0*a)/(1.+n) )*( A_array[0]*phi0_rad + np.dot(A_array[1:], np.sin(2*phi0_rad*np.arange(1,6))) ) # [m]

    # (3) lambda_c, lambda_sの計算
    lambda_c = np.cos(lambda_rad - lambda0_rad)
    lambda_s = np.sin(lambda_rad - lambda0_rad)

    # (4) t, t_の計算
    t = np.sinh( np.arctanh(np.sin(phi_rad)) - ((2*np.sqrt(n)) / (1+n))*np.arctanh(((2*np.sqrt(n)) / (1+n)) * np.sin(phi_rad)) )
    t_ = np.sqrt(1 + t*t)

    # (5) xi', eta'の計算
    xi2  = np.arctan(t / lambda_c) # [rad]
    eta2 = np.arctanh(lambda_s / t_)

    # (6) x, yの計算
    x = A_ * (xi2 + np.sum(alpha_array[1:].reshape(-1,1) *
                           np.sin(2*xi2*np.arange(1,6).reshape(-1,1)) *
                           np.cosh(2*eta2*np.arange(1,6).reshape(-1,1)), axis=0)) - S_ # [m]
    y = A_ * (eta2 + np.sum(alpha_array[1:].reshape(-1,1) *
                            np.cos(2*xi2*np.arange(1,6).reshape(-1,1)) *
                            np.sinh(2*eta2*np.arange(1,6).reshape(-1,1)), axis=0)) # [m]
    # return
    return x, y # [m]


def CircleFitting(x, y):
    """Circle Fitting with least squared
        input: point x-y positions  
        output  cxe x center position
                cye y center position
                re  radius of circle 
    """
    x = np.array(x)
    y = np.array(y)
    sumx = x.sum()
    sumy = y.sum()
    sumx2 = (x ** 2).sum()
    sumy2 = (y ** 2).sum()
    sumxy = (x * y).sum()
    
    F = np.array([[sumx2, sumxy,   sumx],
                  [sumxy, sumy2,   sumy],
                  [sumx,  sumy,  len(x)]])
    x2y2 = x ** 2 + y ** 2;
    G = np.array([[-(x2y2 * x).sum()],
                  [-(x2y2 * y).sum()],
                  [-x2y2.sum()]])
    
    try:
        T = np.linalg.solve(F, G).squeeze()
    except np.linalg.LinAlgError:
        return 0, 0, float("inf")

    cxe = T[0] / -2
    cye = T[1] / -2

    try:
        re = np.sqrt(cxe ** 2 + cye ** 2 - T[2])
    except:
        return cxe, cye, float("inf")
    return cxe, cye, re


def calc_range(idx, n_data, npo):
    """
    npo: the number of points using
    ex) npo=1: using 3 point
        npo=2: using 5 point
        npo=3: using 7 point
    """
    lind = idx - npo
    hind = idx + npo + 1

    if lind < 0:
        lind = 0
    if hind >= n_data:
        hind = n_data

    return range(lind, hind)


def calc_sign(xs, ys):
    # sign evaluation
    c_index = int((len(xs) - 1) / 2.0)
    xc = xs[c_index]
    yc = ys[c_index]
    
    sign = (xs[0] - xc) * (ys[-1] - yc) - (
            ys[0] - yc) * (xs[-1] - xc)
    
    return sign


def check_straight_line(xs, ys):
    c_index = int((len(xs) - 1) / 2.0)
    # check straight line
    a = np.array([xs[0] - xs[c_index], ys[0] - ys[c_index]])
    b = np.array([xs[-1] - xs[c_index], ys[-1] - ys[c_index]])
    na = np.linalg.norm(a)
    nb = np.linalg.norm(b)
    
    if na == 0 or nb == 0:
        return True
    
    theta = np.degrees(np.arccos(np.dot(a, b) / (na * nb)))

    return theta == 180.0


def calc_curvature_circle_fitting(x, y, npo=1):
    """
    Calc curvature
    x,y: x-y position list
    npo: the number of points using Calculation curvature
    ex) npo=1: using 3 point
        npo=2: using 5 point
        npo=3: using 7 point
    """

    cv = np.zeros_like(x)
    n_data = len(x)

    for i in range(n_data):
        r = calc_range(i, n_data, npo)

        xs = x[r]
        ys = y[r]
        cxe, cye, re = CircleFitting(xs, ys)

        if len(xs) < 3:
            continue

        # check straight line
        is_straight = check_straight_line(xs, ys)
        if is_straight:
            # straight line
            continue
        
        # sign evaluation
        sign = calc_sign(xs, ys)

        if sign > 0:
            cv[i] = 1.0 / -re
        else:
            cv[i] = 1.0 / re
            

    return cv


def set_alpha(lines, line_alpha, marker_alpha):
    if not isinstance(lines, Iterable):
        lines = [lines]
    for line in lines:
        line.set_color((*to_rgb(line.get_color()), line_alpha))
        line.set_markerfacecolor((*to_rgb(line.get_markerfacecolor()), marker_alpha))


def get_ax_size(ax):
    fig = ax.figure
    bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    width, height = bbox.width, bbox.height
    width *= fig.dpi
    height *= fig.dpi
    return width, height


def get_closed_index(ax, line, event):
    xmin, xmax = ax.get_xlim()
    ymin, ymax = ax.get_ylim()
    x, y = line.get_data()
    x = np.array(x)
    y = np.array(y)
    w, h = get_ax_size(ax)
    x = (x - event.xdata) * w / (xmax - xmin) 
    y = (y - event.ydata) * h / (ymax - ymin) 
    return np.argmin(x ** 2 + y ** 2)

ipyleafletによる地図表示

  • 移動方向がわかるようにAntPathを使用
  • Maker群を総入れ替えできるようにMarkerClusterを使用
  • Icon で1km, 2kmとか出したかったけどできず、残念
コード 
m = ipyleaflet.Map(center=(35.693406,139.749839), zoom=15, basemap=ipyleaflet.TileLayer(
    url="https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}",
    name="Google Satellite",
    attribution="Google",
    max_zoom=22,
))

ant_path = ipyleaflet.AntPath(
    locations=[],
    dash_array=[1, 10],
    delay=1000,
    color='#7590ba',
    pulse_color='#3f6fba'
)
m.add_layer(ant_path)

marker_cluster = ipyleaflet.MarkerCluster(markers=[])
m.add_layer(marker_cluster)

cur_marker = ipyleaflet.Marker(title='cur', draggable=False)
m.add_layer(cur_marker)

グラフエリアの作成

  • インタラクティブにグラフを書き換えるため、枠だけ作成
  • 各グラフの左上にクリックで選んだ場所の情報が出るようにした
  • カーブ半径が飛び飛びの値で、補完されているので、
    どこにデータがあるかわかるようにマーカーとラインを描画
    (はじめクリックして変なところ選ばれると思って、だいぶはまった・・・）
  • 経路は軸の目盛を消してみた
コード 
fig, (ax_cv, ax_el, ax_xy) = plt.subplots(3,1, figsize=(10,10))

line_cv,  = ax_cv.plot([], [], ".-",            c="C0", label="カーブ半径")
point_cv, = ax_cv.plot([], [],  "x", alpha=0.5, c="C1", label="ポイント")
cur_cv,   = ax_cv.plot([], [],  "*",            c="C2", label="現在地")
label_cv  = ax_cv.text(0.01, 0.99, "", verticalalignment='top', transform=ax_cv.transAxes)
set_alpha(line_cv, 0.1, 1.0)
ax_cv.set_title("カーブ半径")
ax_cv.set_xlabel("移動距離[km]")
ax_cv.set_ylabel("カーブ半径[m]")
ax_cv.grid()
ax_cv.legend()

line_el,  = ax_el.plot([], [],       alpha=0.5, c="C0", label="標高")
point_el, = ax_el.plot([], [],  "x", alpha=0.5, c="C1", label="ポイント")
cur_el,   = ax_el.plot([], [],  "*",            c="C2", label="現在地")
label_el  = ax_el.text(0.01, 0.99, "", verticalalignment='top', transform=ax_el.transAxes)
ax_el.set_title("標高")
ax_el.set_xlabel("移動距離[km]")
ax_el.set_ylabel("標高[m]")
ax_el.grid()
ax_el.legend()

line_xy,  = ax_xy.plot([], []     , c="C0", label="経路")
point_xy, = ax_xy.plot([], [], "x", c="C1", label="ポイント")
cur_xy,   = ax_xy.plot([], [], "*", c="C2", label="現在地")
ax_xy.set_aspect("equal")
ax_xy.set_title("経路(場所選択用)")
ax_xy.grid()
ax_xy.legend()
ax_xy.tick_params(labelbottom=False, labelleft=False, labelright=False, labeltop=False)
ax_xy.tick_params(bottom=False, left=False, right=False, top=False)


fig.tight_layout()

graph = widgets.Output()
with graph:
    fig.canvas.show()

グラフをクリックしたときの処理

  • クリックしたときグラフに一番近い点を選択(get_closed_index)
    Pixelサイズを考慮して、見た目で近いところを選べてるはず。
    あたりまえだけど、x,yでスケールが全然違うと変なところが選ばれる。
  • ipyleaflet の方も更新
コード 
def onclick(event):
    line = None

    if event.inaxes == ax_xy:
        line = line_xy
    elif event.inaxes == ax_cv:
        line = line_cv
    elif event.inaxes == ax_el:
        line = line_el
    else:
        pass

    if line is None:
        return
    
    i = get_closed_index(event.inaxes, line, event)

    for cur, line in zip([cur_xy,  cur_cv,  cur_el],
                         [line_xy, line_cv, line_el]):
        x, y = line.get_data()
        cur.set_data([x[i]], [y[i]])
    
    (x_cv,), (y_cv,) = cur_cv.get_data()
    (x_el,), (y_el,) = cur_el.get_data()
    label_cv.set_text(f"{x_cv:.2f}km, {y_cv:.0f}m")
    label_el.set_text(f"{x_el:.2f}km, {y_el:.0f}m")
    
    location = ant_path.locations[i]
    m.center = location
    cur_marker.location = location
    
        
fig.canvas.mpl_connect('button_press_event', onclick)

グラフの描画処理と実際の表示

  • 国土地理院の断面のCSVをアップロードすると描画
    他のファイルは未対応
  • ipyleaflet で数字を表示するマーカがなかたので、DivIconで代用
  • 1キロ(以下)毎にマーカを追加
    最後の行を選択しているけど、確かに可読性が悪いね・・・
  • グラフエリアに用意したlineやmarkerに値をセット。これももう少し手軽に書けるといいけど
コード 
npo = widgets.IntText(
    description="カーブ半径推定に使用するサンプル数: ",
    value=4, min=2, max=50,
    style={'description_width': 'initial'},
    layout = widgets.Layout(width='350px')
)
upload_label = widgets.Label(value="CSV file: ")
upload = widgets.FileUpload(
    accept='.csv',
    multiple=False
)

def draw_graph(e):
    if len(upload.value) < 1:
        return
    
    content = list(upload.value.values())[0]["content"]
    df = pd.read_csv(io.BytesIO(content))
    lat = df['lat']
    lng = df['lng']
    dist = df['distance']
    elev = df['elevation']
    
    m.center = (lat.mean(), lng.mean())
    ant_path.locations = df[['lat', 'lng']].values.tolist()
    
    def gen_marker(idx):
        html=f'<span style=" display: block; text-align: center; border-radius: 30px; background-color: green;">{dist[idx]}</span>'
        icon = ipyleaflet.DivIcon(html=html, bg_pos=(0, 0), icon_size=[32, 32])
        
        return ipyleaflet.Marker(
            location=(lat[idx], lng[idx]), 
            title=f'{dist[idx]}',
            draggable=False,
            icon=icon
        )
    
    lst = [0]
    for i in range(math.ceil(dist.max() / 1000)):
        idx = len(dist) - (dist <= (i + 1) * 1000)[::-1].argmax() - 1
        lst.append(idx)
    
    marker_cluster.markers = [gen_marker(idx) for idx in lst]
    
    x, y = calc_xy(lat.values, lng.values, 36., 134+20./60)
    # curvature_circle_fitting
    cv = calc_curvature_circle_fitting(x, y, npo=npo.value)

    # radius_circle_fitting
    rc = np.divide(1, np.abs(cv), where=cv!=0, out=np.full_like(cv, np.nan, dtype=np.float64))
    #rc[rc > 1000] = 1000
    
    line_xy.set_data(y - y[0], x - x[0])
    point_xy.set_data(y[lst] - y[0], x[lst] - x[0])
    
    line_cv.set_data(dist/1000, rc)
    point_cv.set_data(dist[lst]/1000, rc[lst])
    
    line_el.set_data(dist/1000, elev)
    point_el.set_data(dist[lst]/1000, elev[lst])
    
    ax_xy.relim()
    ax_cv.relim()
    ax_cv.set_ylim((0, 1000))
    ax_el.relim()
    ax_xy.autoscale_view()
    ax_cv.autoscale_view()
    ax_el.autoscale_view()
    ax_xy.legend()
    ax_cv.legend()
    ax_el.legend()
    fig.tight_layout()
    fig.canvas.draw()



upload.observe(draw_graph, names='value')
npo.observe(draw_graph, names='value')
display(widgets.VBox([
    widgets.HBox([upload_label, upload, npo]),
    m,
    graph
]))

表示イメージ

image.png

おまけ

Details - 折りたたみを使ってみたけど、使いにくいね。
うまく使えば可読性あがるんだろうけど。。。

どうでもいいけど、ネストもいける。さすが

ネスト 
## xxx
- xxx
- xxx
- xxx

<details><summary>コード</summary>


```python

```
</details>
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0