山行時のGPSログをPythonで分析する(3) マルチピッチクライミング時のピッチ情報を推定

マルチピッチのクライミングをしたときに、どのピッチでどれだけ時間がかかったかは、後々の参考になるので記録したい情報です。しかし、現場ではクライミングにもビレイにも必死なので、ノートを取出して鉛筆でメモを取るなんてできません。

僕は、記録を残したいときにはポケットに入れておいたデジカメで写真を撮影し、撮影時刻として記録していますが、写真を撮る余裕がなかったり、撮り忘れたり、撮ってもどの場所だったか分からなくなったりします。

そこで、GPSログを分析することで、クライミング時にどのピッチでどれだけ時間を要したのか推定する方法を考えてみます。

対象は、2020年11月に登った七面山南壁のLong Hope 12p 5.11-にしてみます。当日のGPSログから、クライミング部分のみを抜き出します。こういう作業は地図上でGUI操作するのが早いので、カシミール3Dで編集して、longhope.gpxという名前のgpxファイルに保存します。

ここからはPythonで分析します。例によって、gpxpyでパースしてDataFrameに保存します。

import gpxpy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

filename = './longhope.gpx'

# gpxファイルを読む
gpx_data = []
gpx_file = open(filename, 'r', encoding='utf-8')
gpx = gpxpy.parse(gpx_file)
for track in gpx.tracks:
    for segment in track.segments:
        for point in segment.points:
            gpx_data.append([point.latitude, point.longitude, point.elevation, point.time])

# DataFrameに変換
colname = ['latitude', 'longitude', 'elevation', 'time']
df_gpx = pd.DataFrame(gpx_data, columns=colname)

# 緯度と経度をndarrayに保存
lat = df_gpx['latitude'].values  # 緯度
lon = df_gpx['longitude'].values  # 経度
ele = df_gpx['elevation'].values  # 標高

時間と標高に着目してグラフにしてみます。

t = df_gpx['time']

def time_diff(time):
    time_s = time[0]
    diff = []
    for i in range(len(time)):
        diff.append((time[i] - time_s).seconds)
    return diff

t_diff = np.array(time_diff(time))
t_diff_hr = t_diff / 3600

plt.plot(t_diff_hr, ele, label='GPS data', c='r')
plt.xlabel('time (hr)')
plt.ylabel('elevation (m)')
plt.legend(loc=0)
plt.show()

時間とともに、徐々に高度が上がっていることは分かりますが、もうちょっと情報を引き出したいところです。

時間に伴って変化する信号とみなせるので、信号処理による特徴抽出ができそうです。捉えたいのは、信号(標高)が大きく変化する部分なので、ハイパスフィルタを適用してみます。ハイパス(high pass)の名が示すように、変化が早い部分のみを通過させるフィルタです。詳しくは理解できておりませんが、参考書の例^[1]Cyrille Rossant, IPythonデータサイエンスクックブック, 2015, p. 355を真似てみました。

# バターワースフィルタを適用
N = 4  # 次数
Wn = 0.01  # カットオフ周波数
b, a = sg.butter(N, Wn, 'high') 
ele_hp = sg.filtfilt(b, a, ele)

# 描画
plt.plot(t_diff_hr, ele_hp, label='high path', c='c')
plt.xlabel('time (hr)')
plt.legend(loc=0)
plt.show()

時間の経過とともに高度変化があることが見えてきました。ほとんど変換がない部分はビレイしている時間、徐々に高度があがってる所は自分がクライミングしている時間に見えます。このような変化が起きる位置での経過時刻を読んでやれば、ピッチと時間の関係を得られそうです。手作業で読んでも良いのですが、もう少し自動処理したい感じです。

ハイパスフィルタをかけた後の信号から、変動のピークを読み取って、ピーク位置での経過時間が分かれば目的を達成できるでしょう。ピークが正負の両方にあってややこしそうなので、ひとまず2乗して、すべて正の数にします。

ele_hp2 = ele_hp ** 2
plt.plot(t_diff_hr, ele_hp2, label='high path', c='c')
plt.xlabel('time (hr)')
plt.legend(loc=0)
plt.show()

ピーク位置のインデックスを読むには、scipy.signal.argrelmaxが使えるようです^[2]https://watlab-blog.com/2020/09/26/fft-find-peaks/。引数orderの値で感度が変わります。数が少ないと細かなピークも検出し、大きくすると検出感度が鈍くなります。今回はorderを20にすると、丁度良い頃合いのピークを検出できました。

maxid = sg.argrelmax(ele_hp2, order=20) #最大値
plt.plot(t_diff_hr, ele_hp2, label='high path', c='c')
plt.plot(t_diff_hr[maxid],ele_hp2[maxid],'ro',label='peak')
plt.xlabel('time (hr)')
plt.legend(loc=0)
plt.show()

ピーク位置から取得した経過時刻を標高のグラフに重ね書きします。標高が変化している場所 = クライミング中のピッチ変化に対応している可能性が高い箇所、の自動検出に一応は成功しました。併せて、その箇所での経過時刻も出力して、csvにも保存しておきます。

plt.vlines(t_diff_hr[maxid], min(ele), max(ele), linestyle='dashed', label='changing points', linewidth=0.5, color='k')
plt.plot(t_diff_hr, ele, label='GPS data', c='r')
plt.xlabel('time (hr)')
plt.ylabel('elevation (m)')
plt.legend(loc=0)
plt.show()

df_res = pd.DataFrame(t_diff_hr[maxid])
print(df_res)
df_res.to_csv('time.csv', index=True)

0
0 0.045833
1 0.444444
2 0.751944
3 1.003611
4 1.275278
5 1.543611
6 1.978889
7 2.258056
8 2.557500
9 2.944722
10 3.159722
11 3.707500
12 4.046389
13 4.383056
14 4.524722
15 4.831667
16 5.362222
17 5.654444
18 5.940000
19 6.587222
20 7.015833
21 7.555278
22 7.881944
23 7.983333
24 8.099444

最後は、手作業で表計算ソフトを使って集計してみます。…が、やはり煩雑なので、途中断念してしまいました。高度が順調に上がっている部分が、自分が(リード、フォロー問わず)登っている部分であることは確かなので、その時間だけをカウントすると全8.1時間中の3.2時間でした。

今回のコードをまとめると以下のようになります。

import gpxpy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sg

filename = './longhope.gpx'

# gpxファイルを読む
gpx_data = []
gpx_file = open(filename, 'r', encoding='utf-8')
gpx = gpxpy.parse(gpx_file)
for track in gpx.tracks:
    for segment in track.segments:
        for point in segment.points:
            gpx_data.append([point.latitude, point.longitude, point.elevation, point.time])

# DataFrameに変換
colname = ['latitude', 'longitude', 'elevation', 'time']
df_gpx = pd.DataFrame(gpx_data, columns=colname)

# 緯度と経度をndarrayに保存
lat = df_gpx['latitude'].values  # 緯度
lon = df_gpx['longitude'].values  # 経度
ele = df_gpx['elevation'].values  # 標高

#%%
t = df_gpx['time']
def time_diff(time):
    time_s = time[0]
    diff = []
    for i in range(len(time)):
        diff.append((time[i] - time_s).seconds)
    return diff

t_diff = np.array(time_diff(t))
t_diff_hr = t_diff / 3600

#%%
plt.plot(t_diff_hr, ele, label='GPS data', c='r')
plt.xlabel('time (hr)')
plt.ylabel('elevation (m)')
plt.legend(loc=0)
plt.show()

#%%
# バターワースフィルタを適用
N = 4  # 次数
Wn = 0.01  # カットオフ周波数
b, a = sg.butter(N, Wn, 'high') 
ele_hp = sg.filtfilt(b, a, ele)

# 描画
plt.plot(t_diff_hr, ele_hp, label='high path', c='c')
plt.xlabel('time (hr)')
plt.legend(loc=0)
plt.show()

#%%
ele_hp2 = ele_hp ** 2
plt.plot(t_diff_hr, ele_hp2, label='high path', c='c')
plt.xlabel('time (hr)')
plt.legend(loc=0)
plt.show()

#%%
maxid = sg.argrelmax(ele_hp2, order=20)
plt.plot(t_diff_hr, ele_hp2, label='high path', c='c')
plt.plot(t_diff_hr[maxid],ele_hp2[maxid],'ro',label='peak')
plt.xlabel('time (hr)')
plt.legend(loc=0)
plt.show()

#%%
print(t_diff_hr[maxid])
plt.vlines(t_diff_hr[maxid], min(ele), max(ele), linestyle='dashed', label='changing points', linewidth=0.5, color='k')
plt.plot(t_diff_hr, ele, label='GPS data', c='r')
plt.xlabel('time (hr)')
plt.ylabel('elevation (m)')
plt.legend(loc=0)
plt.show()

df_res = pd.DataFrame(t_diff_hr[maxid])
print(df_res)
df_res.to_csv('time.csv', index=True)

References[+]