山行時のGPSログをPythonで分析する(1) gpxファイルの読込み・経緯度からの距離計算

GARMIN eTrex30xを山行毎に使うようになって半年ばかりが過ぎました。現在地を把握できる絶対的な安心感に加えて、行動の時刻歴も分かりますし、ヤマレコなどにアップすればグラフィカルに表示されるしで、GPSはぼくにとって欠かせないツールになりました。

eTrex30xからはカシミール3D経由でgpxファイルに保存できます。gpxファイルはXMLという規格のファイル形式で、例えば、2020年末に行った黄蓮谷右俣のgpxファイルの中身は以下のようになっています。

<?xml version="1.0" encoding="UTF-8"?>
<gpx
version="1.1"
creator="Kashmir3D 9.360 - http://www.kashmir3d.com"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns="http://www.topografix.com/GPX/1/1"
xmlns:kashmir3d="http://www.kashmir3d.com/namespace/kashmir3d"
xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd
	http://www.kashmir3d.com/namespace/kashmir3d http://www.kashmir3d.com/namespace/kashmir3d.xsd">
<trk>
 <name>2020-12-26~28_黄蓮谷右俣</name>
 <number>1</number>
 <extensions>
 <kashmir3d:line_color>0000ff</kashmir3d:line_color>
 <kashmir3d:line_size>2</kashmir3d:line_size>
 <kashmir3d:line_style>1</kashmir3d:line_style>
 <kashmir3d:icon>901001</kashmir3d:icon>
 </extensions>
<trkseg>
<trkpt lat="35.797310" lon="138.298132">
 <ele>749.490000</ele>
 <time>2020-12-25T21:15:43Z</time>
</trkpt>
<trkpt lat="35.797305" lon="138.298173">
 <ele>750.930000</ele>
 <time>2020-12-25T21:15:49Z</time>
</trkpt>
<trkpt lat="35.797293" lon="138.298165">
 <ele>756.210000</ele>
 <time>2020-12-25T21:15:52Z</time>
</trkpt>
<trkpt lat="35.797289" lon="138.298156">
 <ele>759.580000</ele>
 <time>2020-12-25T21:15:55Z</time>
</trkpt>
...

<trkpt>タグがトラックポイントの内容で、latが緯度、lonが経度です。入れ子になっている<ele>タグの内容が標高、<time>タグの内容が時刻です。このデータを汎用プログラミング言語であるPythonで読込めば、調理できそうです。

Pythonには、gpxpyというライブラリがあり、これを利用すれば簡単にgpxファイルを扱えるようです。僕はanacondaを使用しているので、ここの説明に従ってcondaコマンドからgpxpyをインストールしました。

conda install -c conda-forge gpxpy

先程のログ「orenmigi.gpx」を読み込んで、何も考えずに緯度・経度でプロットしてみます。カシミール3Dの出力には日本語が含まれるので、open時にencoding=’utf-8’を指定する必要があるようです。

import gpxpy
import matplotlib.pyplot as plt

filename = './orenmigi.gpx'
lat = []
lon = []

# gpxファイルを読む
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:
            lat.append(point.latitude)
            lon.append(point.longitude)

plt.plot(lat, lon)
plt.show()

実行すると、以下のグラフが表示されます。横軸が経度(東経)、縦軸が緯度(北緯)です。右上の竹宇駒ヶ岳神社から、黒戸尾根を登り、六丈沢から下降して、黄蓮谷を詰め、甲斐駒ヶ岳山頂から再び黒戸尾根を下山している様子です。きちんと軌跡が表示できて、ちょっと感動します。

今度は、標高を表示してみます。横軸は移動距離にしたいので、緯度と経度から地球表面上の距離を計算します。計算にはヒュベニの式^[1]https://www.trail-note.net/tech/calc_distance/^[2]https://dasonhouse.home.blog/2019/03/17/gpx-to-graph-ver-0-1/なるものを利用させていただきました。

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

filename = './orenmigi.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

# 度からラジアンに換算
rlat = np.radians(lat)
rlon = np.radians(lon)
    
# ヒュベニの式
def hyubeni(rlat_a, rlat_b, rlon_a, rlon_b): 
    A = 6378137.000  # WGS84測地系の楕円体の長半径a
    B = 6356752.314  # WGS84測地系の楕円体の短半径b
    E = np.sqrt((A**2 - B**2) / A**2)  # 離心率
    Dy = rlat_a - rlat_b  # 2点の緯度(latitude)の差 [rad]
    Dx = rlon_a - rlon_b  # 2点の経度(longitude)の差 [rad]
    P = (rlat_a + rlat_b) / 2  # 2点の緯度(latitude)の平均 [rad]
    W = np.sqrt(1 - E**2 * np.sin(P)**2)
    M = A * (1 - E**2) / W**3  # 子午線曲率半径
    N = A / W  # 卯酉線曲線半径
    D = np.sqrt((Dy * M)**2 + (Dx * N * np.cos(P))**2)
    return D  # 2点間の距離 [m]

# 2点間の距離を計算(ヒュベニの式)
dis = []
for i in range(len(rlat)):
    if i == 0:
        dis.append(0)
    else:
        rlat_a = rlat[i-1]
        rlat_b = rlat[i]
        rlon_a = rlon[i-1]
        rlon_b = rlon[i]
        dis.append(hyubeni(rlat_a, rlat_b, rlon_a, rlon_b))
dis = np.array(dis)
dis_cum = np.cumsum(dis)

# 標高
ele = df_gpx['elevation'].values

# 出力
print('累積距離 =', format(dis_cum[-1], '.0f'), '(m)')
plt.plot(dis_cum, ele, label='elevation', color='k')
plt.xlabel('distance (m)')
plt.ylabel('elevation (m)')
plt.legend(loc=0)
plt.show()

実行すると、以下の数字とグラフが出力されます。

累積距離 = 27025 (m)