LoginSignup
7
6

More than 3 years have passed since last update.

@toran

Python3.7によるGPXファイルの読み込みと移動距離の算出

Last updated at Posted at 2020-01-02

目的

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

参考

7
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
6