目的
地図上で経路を引いたときにカーブ半径と高低差を見たかった。
（そういうのが気になるお年頃ｗ）
断面図｜地理院地図の使い方を見ていて、高低差は見れるかなと思って、
を使っていたところ、グラフの保存で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
]))
表示イメージ
おまけ
Details - 折りたたみを使ってみたけど、使いにくいね。
うまく使えば可読性あがるんだろうけど。。。
どうでもいいけど、ネストもいける。さすが
ネスト
## xxx
- xxx
- xxx
- xxx
<details><summary>コード</summary>
```python
```
</details>