前回はアップルウォッチのデータからGPSデータを読み込みメートル座標系に変換する方法について紹介しました。今回は続きの３から説明していきます。
１、アップルウォッチのデータをMATLABで読み込む２、緯度経度情報をメートル座標系に変換する
３、速度に応じて色分けした軌跡を表示する
４、スライダーを設置して値の更新を確認する
５、矢印を用いて走行時の向きと速度を可視化する
３、速度に応じて色分けした軌跡を表示する
GPSデータの中には緯度経度以外にも速度のデータが含まれています。ここでは速度に応じて段階別に色分けしたグラフを作成します。
速度情報をグラフ化
まずは下記コードで速度情報をSpeedという変数に格納します。前回作成したtblという変数の中のextensionsにspeedという値が保存されています。この値を抽出しSpeedという変数に代入していきます。保存されている速度は[m/s]単位なので好みに応じて3.6を掛け算し[km/h]に変換することも可能です。
len = length(tbl.time);
Speed = zeros(len,1);
for i=1:len
Speed(i) = tbl.extensions(i).speed * 3.6;
end得られた速度情報の確認のため、下記コードでグラフを作成します。
plot(Speed)
横軸にサンプル数、縦軸に速度[km/h]が表示されています。1km約6分で走行するペースで走ったデータを使用しているので、10km/h程度の値が表示されておりおおよそ正しいことがわかります。また、途中信号待ちなどで停止した際の速度も正しく記録されていそうです。
タイムスタンプを抽出
後々使用しますのでタイムスタンプを抽出しておきます。tbl.timeに「2023-04-08T00:16:20Z」のようにタイムスタンプが保存されています。これらを抽出してタイムスタンプを作成しましたが、アップルウォッチでのGPSデータは1sec間隔のデータですので最初のデータを0secとして、1secずつ加算していく方法と変化ありませんでした。下記１行を追加しました。
len = length(tbl.time);
Speed = zeros(len,1);
for i=1:len
tbl.timeoffset(i) = i-1;
Speed(i) = tbl.extensions(i).speed * 3.6;
endtblにtimeoffsetを追加して、i-1を代入しています。
念のため下記コードでplotして確認しておきます。
plot(tbl.timeoffset)
速度に応じてカテゴリ分け
せっかく速度データが取得できていますので、速度に応じて色分けしたグラフを作成します。
MATLABのStatics and Machine Learning Toolboxにはgscatter関数がありますので速度域に応じたcategorical変数を作成することで簡単に色分けした散布図が作成できます。
数行のコードで同様のグラフが作成できますので、今回はgscatterを使用しない方法をご紹介します。
下記コードのように速度をSpd1~Spd4を定義して速度を５段階に分割するためのIdx1~5を作成します。Spd1~4は好みに設定して問題ありませんが、[km/h]に変換するため、3.6を掛けています。行列から条件に一致する値を抽出する方法は、配列から条件を満たす値を抽出で紹介しておりますので事前にご確認ください。
それぞれの速度領域ごとにIdxを作成し、plotしていきます。注意点としては、デフォルトのplotだと２点間を線でつないでしまうので点でplotするように「’.’」を追加します。
また、カテゴリごとに色分けするため、カラーコードで点のカラーを指定しました。
Spd1 = 1*3.6;
Spd2 = 2.8*3.6;
Spd3 = 3*3.6;
Spd4 = 3.2*3.6;
Idx1 = Speed < Spd1;
Idx2 = Speed >= Spd1 & Speed < Spd2;
Idx3 = Speed >= Spd2 & Speed < Spd3;
Idx4 = Speed >= Spd3 & Speed < Spd4;
Idx5 = Speed >= Spd4;
f2 = figure(2);
hold on;
grid on;
grid minor;
plot(Lon(Idx1),Lat(Idx1),'.','Color','#0000ff');
plot(Lon(Idx2),Lat(Idx2),'.','Color','#00ffff');
plot(Lon(Idx3),Lat(Idx3),'.','Color','#00ff00');
plot(Lon(Idx4),Lat(Idx4),'.','Color','#ff8000');
plot(Lon(Idx5),Lat(Idx5),'.','Color','#ff0000');
速度領域ごとに色分けしたグラフが作成できました。アプリ等でもよくあるような見た目に変えることができました。
グラフの見た目を整える
最後にグラフの見た目を整えます。速度域ごとに色が異なるため凡例を追加しましょう。Spd1~4は後々調整する可能性があるため、直入力するのではなくsprintfを使用します。また、axis equalで横軸と縦軸の縮尺を1:1にし、xlim,ylimで描画範囲を整えています。
legend(sprintf('< %1.1f',Spd1),...
sprintf('< %1.1f',Spd2),...
sprintf('< %1.1f',Spd3),...
sprintf('< %1.1f',Spd4),...
sprintf('> %1.1f',Spd4));
axis equal
xlim([-1700 1000])
ylim([-100 2200])
まとめ
今回はGPSデータ解析①の続きについて説明しました。GPSデータから速度データを抽出し、速度領域ごとに色分けしたグラフを作成する方法を紹介しております。Statics and Machine Learning Toolboxのgscatter関数を使用せず、自分で色分けする方法を解説しました。また、最後にこれまで作成したコードを掲載しておきますのでご活用ください。
１、アップルウォッチのデータをMATLABで読み込む２、緯度経度情報をメートル座標系に変換する３、速度に応じて色分けした軌跡を表示する
次回以降は下記について説明します。
４、スライダーを設置して値の更新を確認する
５、矢印を用いて走行時の向きと速度を可視化する
close all;clear
xmlObj = readstruct('route_2023-04-08_9.52am.gpx', 'FileType', 'xml');
tbl = struct2table(xmlObj.trk.trkseg.trkpt);
len = length(tbl.time);
Speed = zeros(len,1);
for i=1:len
tbl.timeoffset(i) = i-1;
Speed(i) = tbl.extensions(i).speed * 3.6;
end
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
Spd1 = 1*3.6;
Spd2 = 2.8*3.6;
Spd3 = 3*3.6;
Spd4 = 3.2*3.6;
Idx1 = Speed < Spd1;
Idx2 = Speed >= Spd1 & Speed < Spd2;
Idx3 = Speed >= Spd2 & Speed < Spd3;
Idx4 = Speed >= Spd3 & Speed < Spd4;
Idx5 = Speed >= Spd4;
f1 = figure(1);
plot(Lon,Lat);
axis equal
f2 = figure(2);
bx = subplot(10,1,[1 9]);
hold on;
grid on;
grid minor;
plot(Lon(Idx1),Lat(Idx1),'.','Color','#0000ff');
plot(Lon(Idx2),Lat(Idx2),'.','Color','#00ffff');
plot(Lon(Idx3),Lat(Idx3),'.','Color','#00ff00');
plot(Lon(Idx4),Lat(Idx4),'.','Color','#ff8000');
plot(Lon(Idx5),Lat(Idx5),'.','Color','#ff0000');
legend(sprintf('< %1.1f',Spd1),...
sprintf('< %1.1f',Spd2),...
sprintf('< %1.1f',Spd3),...
sprintf('< %1.1f',Spd4),...
sprintf('> %1.1f',Spd4));
axis equal
xlim([-1700 1000])
ylim([-100 2200])
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
コメント