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

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

[断面図｜地理院地図の使い方](https://maps.gsi.go.jp/help/intro/kinolist/4-danmen.html)を見ていて、高低差は見れるかなと思って、

https://maps.gsi.go.jp/

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

なら一緒に見れるものを作ろう。
（[カーブ半径の推定](https://github.com/AtsushiSakai/CurvatureNumeric/tree/master)は、どうしても変なところとかにフィッティングしてしまうので、経路図を見ながら見たいのもあったので）

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

# 実装

<details><summary>関数定義等</summary>

は説明省略

```python
%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)

```
</details>

## ipyleafletによる地図表示
- 移動方向がわかるようにAntPathを使用
- Maker群を総入れ替えできるようにMarkerClusterを使用
- Icon で1km, 2kmとか出したかったけどできず、残念

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


```python
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)

```
</details>

## グラフエリアの作成
- インタラクティブにグラフを書き換えるため、枠だけ作成
- 各グラフの左上にクリックで選んだ場所の情報が出るようにした
- カーブ半径が飛び飛びの値で、補完されているので、<br>どこにデータがあるかわかるようにマーカーとラインを描画<br>(はじめクリックして変なところ選ばれると思って、だいぶはまった・・・）
- 経路は軸の目盛を消してみた

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


```python
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()

```
</details>



## グラフをクリックしたときの処理
- クリックしたときグラフに一番近い点を選択(```get_closed_index```)<br>Pixelサイズを考慮して、見た目で近いところを選べてるはず。<br>あたりまえだけど、x,yでスケールが全然違うと変なところが選ばれる。
- ipyleaflet の方も更新

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


```python
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)
```
</details>



## グラフの描画処理と実際の表示
- 国土地理院の断面のCSVをアップロードすると描画<br>他のファイルは未対応
- ipyleaflet で数字を表示するマーカがなかたので、DivIconで代用
- 1キロ(以下)毎にマーカを追加<br>[最後の行](https://ikatakos.com/pot/programming/python/packages/pandas/first_row_where_true)を選択しているけど、確かに可読性が悪いね・・・
- グラフエリアに用意したlineやmarkerに値をセット。これももう少し手軽に書けるといいけど

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


```python
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
]))
```
</details>

# 表示イメージ
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/455423/bcd6dd46-25fa-7f50-92fb-f610e3f4643f.png)


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

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

<details><summary>ネスト</summary>

``````markdown:
## xxx
- xxx
- xxx
- xxx

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


```python

```
</details>
``````

</details>