やりたいこと
Ride with GPS等で作ったルートを自分の走力でどのぐらいの時間で走れるのか知りたい。
やったこと
pythonでGPXファイルを処理して、Excelで勾配と速度を紐づけて所要時間を算出しました。
GPXファイルを処理する
こちらを参考にしました。
Python3.7によるGPXファイルの読み込みと移動距離の算出
これで累積距離（D_mov）とその地点の高度（Alt）を含むデータをcsvに吐き出せます。
Google ColaboratoryやJupyter Notebookで実行するとよいと思います。
import numpy as np
import pandas as pd
import gpxpy
import gpxpy.gpx
from pytz import timezone
def calc_moving_distance(Lat, Lng, Alt):
# 地球の半径 [m]
r = 6378.137 * 1e3
# データ数
N = Lat.shape[0]
# データ抽出
i1 = np.arange(0, N-1)
i2 = np.arange(1, N)
X1 = Lng[i1]
Y1 = Lat[i1]
Z1 = Alt[i1]
X2 = Lng[i2]
Y2 = Lat[i2]
Z2 = Alt[i2]
# deg -> rad
X1 = np.deg2rad(X1)
Y1 = np.deg2rad(Y1)
Z1 = np.deg2rad(Z1)
X2 = np.deg2rad(X2)
Y2 = np.deg2rad(Y2)
Z2 = np.deg2rad(Z2)
# データマスク
mask = ((X1 != X2) & (Y1 != Y2))
# データ間距離
dx = X2 - X1
D_div = np.zeros(N-1)
D_div[mask] = r * np.arccos(np.sin(Y1[mask])*np.sin(Y2[mask])
+ np.cos(Y1[mask])*np.cos(Y2[mask])*np.cos(dx[mask]))
D_div = np.sqrt((Z2-Z1)**2 + D_div**2)
D_div = np.insert(D_div, 0, 0.0)
# 移動距離
D_mov = np.cumsum(D_div)
return D_mov
# タイムゾーン
dt_tz = timezone('Asia/Tokyo')
# 日時文字列形式
dt_fmt = '%Y-%m-%d %H:%M:%S'
# 配列の初期化
DateTime = []
Lat = []
Lng = []
Alt = []
# GPXファイルの読み込み
gpx_file_r = open('input.gpx', 'r')
# GPXファイルのパース
gpx = gpxpy.parse(gpx_file_r)
# GPXデータの読み込み
for track in gpx.tracks:
for segment in track.segments:
# ポイントデータリストの読み込み
points = segment.points
# ポイントデータの長さ
N = len(points)
# ポイントデータの読み込み
for i in range(N):
# ポイントデータ
point = points[i]
# データ抽出
#datetime = point.time.astimezone(dt_tz).strftime(dt_fmt)
lat = point.latitude
lng = point.longitude
alt = point.elevation
# データ代入
#DateTime.append(datetime)
Lat.append(lat)
Lng.append(lng)
Alt.append(alt)
# numpy配列への変換
#DateTime = np.array(DateTime, dtype=np.datetime64) #RidewithGPSのgpxには日付がないのでコメントアウト
Lat = np.array(Lat, dtype=np.float64)
Lng = np.array(Lng, dtype=np.float64)
Alt = np.array(Alt, dtype=np.float64)
# 移動距離の算出
D_mov = calc_moving_distance(Lat, Lng, Alt)
# CSV出力
df = pd.DataFrame(np.c_[Lat, Lng, Alt, D_mov],
columns=['Lat', 'Lng', 'Alt', 'D_mov'])
#index=DateTime)
df.to_csv('output.csv')
gpx_file_r.close()
勾配と速度の計算
参考文献：ヒルクライム計算、タイムからWを計算
ある速度v[km/h]で勾配x[%]の坂を登るために必要な出力P[W]は、登坂抵抗U[W]（位置エネルギーを稼ぐのに必要な出力）、空気抵抗$\gamma$[W]、機械損失$\eta$[W]とすると次式になります。
P=U+\gamma+\eta
各項について確認していきましょう。
登坂抵抗
単位時間あたりに得る上昇量から、単位時間あたりに得る位置エネルギーを計算しましょう。
速度をv[km/h]、勾配をx[%]、総重量（人間+機材）をm[kg]、重力加速度をg[m/$s^2$]として、次式のとおりです。
U=mg * (1000/60^2)v * (1/100)x
位置エネルギーの定義通りですね。
空気抵抗
参考文献からもってきました。あんまり神経質にならなくていいと思います。
$\gamma$=1/2×空気密度$\rho$×面積A×係数$C_d$×$(速度v)^3$とします。
空気密度$\rho$は1.205[kg/m^3]
面積Aは0.4046[m^2]
空気抵抗係数$C_d$は0.468
を採用しました。
機械損失
参考文献からざっくり13.72Wとしました。
求めたいのは速度なんじゃが？
今回の目標はgpxデータと自分の体についてのデータ（体重やPWR、FTP）から所要時間を求めることでした。
Q.求めたいのは速度なんじゃが？
A.三次方程式なんか解いていられるか！
→excelのソルバーアドインで解を探してしまいましょう。
自分の体はそんな急激には変わりませんから、再計算の手間には目をつぶります。
勾配と速度の早見表
参考までに体重65kg、機材重量10kg、PWR2.5倍くらいで計算した結果の表を載せておきます。（貧脚ですまんな）
下り勾配ではブレーキをかけて40km/hを超えない、平坦～緩勾配は力を抜く…といったことを仮定しています。
ここらへんの匙加減が一番のキモでは？と思います。
また、計算しなくとも、実績と比べてPDCA回していい塩梅になるように調整していけばいいとも思います。
| 勾配 | 出力 | 速度 |
| -10.5 | -687.284716 | 40 |
| -10 | -646.4513827 | 40 |
| -9.5 | -605.6180494 | 40 |
| -9 | -564.784716 | 40 |
| -8.5 | -523.9513827 | 40 |
| -8 | -483.1180494 | 40 |
| -7.5 | -442.284716 | 40 |
| -7 | -401.4513827 | 40 |
| -6.5 | -360.6180494 | 40 |
| -6 | -319.784716 | 40 |
| -5.5 | -278.9513827 | 40 |
| -5 | -238.1180494 | 40 |
| -4.5 | -197.284716 | 40 |
| -4 | -156.4513827 | 40 |
| -3.5 | -115.6180494 | 40 |
| -3 | -74.78471605 | 40 |
| -2.5 | -33.95138272 | 40 |
| -2 | 2.2599E-05 | 39.06779279 |
| -1.5 | -3.92232E-06 | 32.89141455 |
| -1 | 0.885189767 | 25 |
| -0.5 | 26.4060231 | 25 |
| 0 | 51.92685643 | 25 |
| 0.5 | 72.02298133 | 24 |
| 1 | 90.42955396 | 23 |
| 1.5 | 107.1319029 | 22 |
| 2 | 122.1153566 | 21 |
| 2.5 | 135.3652438 | 20 |
| 3 | 139.999968 | 18.20744278 |
| 3.5 | 144.9999495 | 16.76043236 |
| 4 | 149.9999769 | 15.55946952 |
| 4.5 | 154.99999 | 14.55650449 |
| 5 | 160.0001388 | 13.71194506 |
| 5.5 | 159.9999281 | 12.59201416 |
| 6 | 159.9999617 | 11.62743307 |
| 6.5 | 159.9999797 | 10.79112849 |
| 7 | 159.9999893 | 10.06108453 |
| 7.5 | 159.9999956 | 9.419516273 |
| 8 | 159.9999908 | 8.852073774 |
| 8.5 | 160.0000022 | 8.347152852 |
| 9 | 159.9999987 | 7.895321964 |
| 9.5 | 160.0000039 | 7.488876571 |
| 10 | 160.0000102 | 7.121479125 |
| 10.5 | 160.0000128 | 6.787883572 |
| 11 | 160.0000017 | 6.48371851 |
| 11.5 | 160.0000006 | 6.205319484 |
| 12 | 159.9999999 | 5.949593002 |
| 12.5 | 159.9999997 | 5.713913509 |
| 13 | 159.9999985 | 5.496039563 |
| 13.5 | 159.9999955 | 5.294047378 |
| 14 | 159.999991 | 5.10627738 |
| 14.5 | 159.9999929 | 4.931291195 |
| 15 | 159.9999956 | 4.767835807 |
| 15.5 | 159.9999791 | 4.614814871 |
| 16 | 159.9999671 | 4.471266911 |
| 16.5 | 159.9999597 | 4.336343592 |
| 17 | 159.9999566 | 4.209294263 |
| 17.5 | 159.9999574 | 4.089452401 |
| 18 | 159.9999615 | 3.976224226 |
| 18.5 | 159.9999686 | 3.869079109 |
| 19 | 159.9999779 | 3.767541442 |
| 19.5 | 159.9999892 | 3.671183722 |
| 20 | 159.9999999 | 3.579620606 |
GPXデータ（勾配）と速度を紐づけ
GPXデータから勾配を求め（Alt差/D_mov差）、勾配に対応する速度は先ほどの表からあてはめていきます。
前項でexcelのソルバーアドイン使った流れでExcelでやりました。
ExcelでやったことをQiitaに貼るの難しいな…
こんな感じの表を作っていきました。
| A | B | C | D | E | F | G | H | I | J |
| # | Lat | Lng | Alt | D_mov | ポイント間距離[m] | 勾配[%] | 勾配0.5丸め[%] | 勾配時速[km/h] | 所要時間[min] |
| 0 | 36.74469 | 138.36522 | 369.1 | 0 | 0 | 0 | 0 | =IF(H2<勾配と速度!$A$13,勾配と速度!$F$13,VLOOKUP(H2,勾配と速度!$A$13:$G$74,6)) | =F2/(1000*I2/60) |
| 1 | 36.74491 | 138.36501 | 368.1 | 30.8329275216206 | =E3-E2 | =IF(F3<1,0,100*(D3-D2)/(F3)) | =MAX(MIN(IF(G3=0,0,G3/ABS(G3)*MROUND(ABS(G3),0.5)),20),-20) | =IF(H2<勾配と速度!$A$13,勾配と速度!$F$13,VLOOKUP(H2,勾配と速度!$A$13:$G$74,6)) | =F2/(1000*I2/60) |
平均して30m間隔くらいでGPXデータポイントがあるので、各ポイント間の所要時間を合計していけば、目的地までのおおよその所要時間がわかるという寸法です。
チェックポイント間の所要時間を合計して、予定表に平均速度として組み込むのも良いでしょう。
予定表の例：ロングライド予定表の作成・運用法
まとめと感想
Ride with GPS等でルートを作り、自分の走力でどのぐらいの時間で走れるのかの目安を算出してみました。
Ride with GPSに課金すればいいじゃん。
途中「GPXデータから勾配計算するのに0距離ポイントとか邪魔やなあ」とか「あ～これ数値積分やなあ」とか「GPXデータのポイントが不等間隔だけど同じ重みでいいのか？」とか「移動平均とか取ってスムースにしないとほんとはダメだろうな」とか思いましたが、そもそも高精度で計算できたとして計算通りの出力を淡々と維持できるわけでもありませんから、精度を求めてもあまり意味はない気がします。
結局は脚
このくらいざっくりでも計算すれば、ゆるポタ詐欺は防げるんじゃないかなと思います（笑）