GPSデータ解析① | MATLAB事例4

前回は複数ファイルの一括処理方法について紹介しました。

今回は緯度経度情報を活用する事例について解説します。実際に計測したデータ用いた解析方法についてご紹介いたしますので、実務で役立つ情報を加えながら説明していきます。

目標の確認

今回はアップルウォッチで計測したGPSデータを取得して、MATLABで下記のような走行時のデータを確認するアプリケーションを作りたいと思います。

具体的には下記５項目に分けて説明しようと思います。
１、アップルウォッチのデータをMATLABで読み込む
２、緯度経度情報をメートル座標系に変換する
３、速度に応じて色分けした軌跡を表示する
４、スライダーを設置して値の更新を確認する
５、矢印を用いて走行時の向きと速度を可視化する

１、アップルウォッチのデータをMATLABで読み込む

アップルウォッチのワークアウトアプリで記録したデータは、接続しているiPhoneのヘルスケアアプリで書き出すことができます。ワークアウトデータの書き出し方法は下記のとおりです。

ヘルスケアアプリを開き、右上のアカウントが表示された部分をタップ。

一番下までスクロールして、「すべてのヘルスケアデータを書き出す」をタップ

「書き出す」をタップでデータの書き出しが始まります。保存先はお好みですが、iPhoneのファイルアプリに書き出しを行い、クラウド経由でPCに送付することができます。

書き出したデータはzip形式となっており、解凍すると「apple_health_export」というフォルダが出来上がります。

「***.\apple_health_export\workout-routes」に「route_2023-04-08_9.52am.gpx」のようなgpxファイルがありますので、このデータを使用します。新規スクリプトファイルを作成し下記のmatlabコードでファイルを読み込み、データをtable形式に変換します。

xmlObj = readstruct('route_2023-04-08_9.52am.gpx', 'FileType', 'xml');
tbl = struct2table(xmlObj.trk.trkseg.trkpt);

readstruct関数でgpxデータをxml形式で読み込み、xmlObjに代入します。xmlObj.trk.trkseg.trkptの中にGPS情報が保存されていましたので、このデータをテーブルに変換して、tblに代入しています。

読み込んだデータを確認するために下記コードを実行します。

plot(tbl.lonAttribute,tbl.latAttribute)

実行すると下図のように横軸に経度、縦軸に緯度を表すグラフが作成できます。下のデータは右下の位置がスタート地点で、北に進み、大通りに出た後、北西に向かって進んでいます。折り返し地点に到着したのち同じ道を帰ってきたデータになります。（横軸、縦軸の情報は削除しました）

以上の手順でアップルウォッチで取得したヘルスケアデータからGPS情報を表示することができました。

２、緯度経度情報をメートル座標系に変換する

GPSデータは緯度と経度の情報からなりますが、走行距離などがわかりにくいためメートルの座標系に変換します。https://www.trail-note.net/tech/calc_distance/を参考にヒュベニの公式を用いて緯度経度情報から距離を計算します。下記のような関数を作成しました。ヒュベニの公式では２点の緯度経度情報をもとに距離を計算していますが、緯度方向の移動と経度方向の移動をそれぞれ別で計算するようにしました。つまり、２点間の距離を計算するとき下記２つの手順に分けます。

①緯度差を0にして経度差から横方向の移動距離を計算
②経度差を0にして緯度差から縦方向の移動距離を計算

また、移動距離がマイナスになることはありませんが、x-y座標系でplotする都合上、経度差が負、あるいは緯度差が負の時はx,y方向距離が負になるように処理を追加しています。

function d = hubeny(Pos1,Pos2)
    %Pos1 = [lat lon]
    %Pos2 = [lat lon]

    Dy = deg2rad(Pos2(1) - Pos1(1));
    Dx = deg2rad(Pos2(2) - Pos1(2));
    P = deg2rad(mean([Pos1(1) Pos2(1)]));

    Rx = 6378137.0;%[m]
    Ry = 6356752.314245;%[m]
    E2 = (Rx^2 - Ry^2)/Rx^2;
    W = sqrt(1-E2 * sin(P)^2);
    N = Rx/W;
    M = Rx*(1-E2) / W^3;
    
    if Dx==0 && Dy<0
        d = -sqrt((Dy*M)^2 + (Dx*N*cos(P))^2);
    elseif Dx<0 && Dy==0
        d = -sqrt((Dy*M)^2 + (Dx*N*cos(P))^2);
    else
        d = sqrt((Dy*M)^2 + (Dx*N*cos(P))^2);
    end
            
end

上記の関数を別ファイルで保存するか、スクリプトファイルの最後尾に記述します。関数の作成方法がわからない場合関数の作成 | MATLAB基礎18を参考にしてください。

スクリプトファイルに下記の記述を追加します。

len = length(tbl.latAttribute);
Lat = zeros(len,1);
Lon = zeros(len,1);
Dist = zeros(len,1);
for i=1:len
    Lat(i) = hubeny([tbl.latAttribute(1) tbl.lonAttribute(1)],[tbl.latAttribute(i) tbl.lonAttribute(1)]);
    Lon(i) = hubeny([tbl.latAttribute(1) tbl.lonAttribute(1)],[tbl.latAttribute(1) tbl.lonAttribute(i)]);
    if i>1
        Dist(i) = hubeny([tbl.latAttribute(i-1) tbl.lonAttribute(i-1)],[tbl.latAttribute(i) tbl.lonAttribute(i)]);
    end
end

そして、下記コードでグラフを作成し結果を確認します。

plot(Lon,Lat)

計測開始時点を原点として、北方向を縦軸の正方向、東方向を横軸の正方向に正しくplotできていることがわかります。

まとめ

今回は「MATLABでGPSデータを処理する」方法について紹介しています。最終的に作るアプリケーションの確認を行いました。最初の手順として下記２項目について実施しました。
１、アップルウォッチのデータをMATLABで読み込む
２、緯度経度情報をメートル座標系に変換する

次回は下記の項目について説明していきます。

３、速度に応じて色分けした軌跡を表示する
４、スライダーを設置して値の更新を確認する
５、矢印を用いて走行時の向きと速度を可視化する