目的

GPSデータ形式の一つであるGPXファイルをPythonで読み込み，CSVファイルに変換する．
基本的にデータごとの緯度経度しかないため，移動距離を算出し，それを新たに列に加える．

環境

下記のAnacondaをインストールしたWindows10のPCにて，PATHを通し，Cygwin環境から実行している．

Windows10 Pro 64bit
Anaconda3-2018.12-Windows-x86_64 / Python 3.7
gpxpy 1.3.5

確認した際が以下．

toran@localhost ~
$ python -V
Python 3.7.1

toran@localhost ~
$ conda -V
conda 4.8.0

toran@localhost ~
$ conda list gpxpy
# packages in environment at C:\Users\toran\Anaconda3:
#
# Name                    Version                   Build  Channel
gpxpy                     1.3.5                      py_0    conda-forge

0. gpxpyのインストール

基本的には以下でインストール可．

toran@localhost ~
$ conda install gpxpy -c conda-forge

1. GPXファイルの読み込み

input.gpxというGPXファイルを読み込む場合，以下のtest1.pyのような書き方で可能．
ファイルを読み込んでいき，各配列にデータを追加している．

test1.py

import gpxpy
import gpxpy.gpx
from pytz import timezone

# タイムゾーン
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)

2. 移動距離の算出

計算速度の観点から，Numpyを使用する．

import numpy as np

次に，test1.pyの各配列をnumpy配列に変換する．

DateTime = np.array(DateTime, dtype=np.datetime64)
Lat      = np.array(Lat, dtype=np.float64)
Lng      = np.array(Lng, dtype=np.float64)
Alt      = np.array(Alt, dtype=np.float64)

2地点間の距離は，赤道半径$r=6378.137$kmを半径とする球体として水平距離$d_\mathrm{2D}$を求め，
標高差から斜距離を求める（簡易式であり，GPXデータの取得時間間隔が広いと誤差が大きく点に注意）．

d_\mathrm{2D} = r \cos^{-1} \left(\sin y_1 \sin y_2 + \cos y_1 \cos y_2 \cos \Delta x \right) \\
\Delta x = x_2 - x_1 \\
d_\mathrm{3D} = \sqrt{d_\mathrm{2D}^2 + (z_2 - z_1)^2}

ここから，地点1の経度および緯度を$x_1$および$y_1$，地点2の経度および緯度を$x_2$および$y_2$とし，
各々の標高を$z_1$および$z_2$とすると，二地点間の距離は以下で求められる．

# 地球の半径 [m]
r = 6378.137 * 1e3
# 二地点間の距離
dx    = x2 - x1
d_div = r * np.arccos(np.sin(y1)*np.sin(y2) + np.cos(y1)*np.cos(y2)*np.cos(dx))
d_div = np.sqrt((z2-z1)**2 + d_div**2)

GPXデータからはすべて配列であるため，$(X[i], Y[i], Z[i])$および$(X[i+1], Y[i+1], Z[i+1])$間の
地点間距離を求め，その累積和を取る形とすることで，データごとの移動距離を一度に求めることができる．
なお，GPXデータによっては地点間の緯度経度が同一である場合があるので，マスクにより除外する点に注意する．
以下が，緯度経度および標高の配列から，全体の移動距離を求める関数calc_moving_distanceである．

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

参考にしたページ

https://keisan.casio.jp/exec/system/1257670779

3. CSVへの出力

シンプルで使いやすいので，Pandasを用いる．

import pandas as pd

以下のような形でCSV出力が可能．

# 移動距離の算出
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')

ファイル全体

以下に，今までのコードを一つにしたものと，クラスにした二通りをまとめる．

その1：一つにしたもの

test2.py

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)
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')

その2：クラスにしたもの

import numpy as np
import pandas as pd
import gpxpy
import gpxpy.gpx
from pytz import timezone

class GPXModel(object):
    # タイムゾーン
    dt_tz = timezone('Asia/Tokyo')
    # 日時文字列形式
    dt_fmt = '%Y-%m-%d %H:%M:%S'

    def __init__(self, fpath):
        # 配列の初期化
        DateTime = []
        Lat      = []
        Lng      = []
        Alt      = []
        # GPXファイルの読み込み
        gpx_file_r = open(fpath, '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(__class__.dt_tz).strftime(__class__.dt_fmt)
                    lat      = point.latitude
                    lng      = point.longitude
                    alt      = point.elevation
                    # データ代入
                    DateTime.append(datetime)
                    Lat.append(lat)
                    Lng.append(lng)
                    Alt.append(alt)
        # numpy配列への変換
        self.DateTime = np.array(DateTime, dtype=np.datetime64)
        self.Lat      = np.array(Lat, dtype=np.float64)
        self.Lng      = np.array(Lng, dtype=np.float64)
        self.Alt      = np.array(Alt, dtype=np.float64)
        # 移動距離の算出
        self.calc_moving_distance()
    
    def calc_moving_distance(self):
        # 地球の半径 [m]
        r = 6378.137 * 1e3
        # データ数
        N = self.DateTime.shape[0]
        # データ抽出
        i1 = np.arange(0, N-1)
        i2 = np.arange(1, N)
        X1 = self.Lng[i1]
        Y1 = self.Lat[i1]
        Z1 = self.Alt[i1]
        X2 = self.Lng[i2]
        Y2 = self.Lat[i2]
        Z2 = self.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)
        # 移動距離
        self.D_mov = np.cumsum(D_div)

    def to_csv(self, fpath):
        # CSV出力
        df = pd.DataFrame(np.c_[self.Lat, self.Lng, self.Alt, self.D_mov],
                          columns=['Lat', 'Lng', 'Alt', 'D_mov'],
                          index=self.DateTime)
        df.to_csv(fpath)

model = GPXModel(fpath='input.gpx')
model.to_csv(fpath='output.csv')

参考