GPSデータ解析③ | MATLAB事例6

前回はアップルウォッチのデータからGPSデータから速度データを抽出し、速度領域ごとに色分けしたグラフを作成する方法について紹介しました。今回は「４、スライダーを設置して値の更新を確認する」と「５、矢印を用いて走行時の向きと速度を可視化する」について説明していきます。

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

４、スライダーを設置して値の更新を確認する

グラフにスライダーを設置して値を更新する方法については、スライダーを使ってfigureを更新で紹介しましたので詳細はこちらをご確認ください。

まずはじめに値を更新するためのplotを追加します。

MATLAB 
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');
h = plot(Lon(1),Lat(1),'ko');

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

% figureにテキストを設置
d = uicontrol(f2,'Style','text');
d.Position = [250 30 200 20];
d.String = 'Time 0';

% figureにスライダーを設置
c = uicontrol(f2,'Style','slider');
c.Position = [20 30 200 20];
c.Value = 1;
c.Min = 1;
c.Max = len;
c.SliderStep = [1/len 10/len];

2行目でsubplotを使いスライダーを設置するスペースを確保しています
11行目でLon,Latの１点目のみplotして「〇」で表示しています。

出来上がったグラフは下記のとおりです。

スライダーを設置したグラフ

 

これにグラフ更新用の関数と、スライダーを変更した際に呼び出される関数を定義します。

MATLAB 
c.Callback = @(src,event)updatePosition(h,c,d,Lat,Lon,tbl);

function [] = updatePosition(h,c,d,Lat,Lon,tbl)
    i = uint32(c.Value);
    d.String = sprintf('Time %4d',tbl.timeoffset(i));
    h.XData = Lon(i);
    h.YData = Lat(i);
end

スライダーが更新されたときに呼ばれる関数updatePositionでさきほどplotした「〇」の位置を更新しています。ついでに　グラフしたのTimeも更新しています。

スライダーでグラフを更新

 

現段階で上記のようなものができていると思います。現状のコードも下記に掲載します。これだけでは以前ご紹介した内容のデータが変わっただけというものになってしまいますので、さらに情報を追加したいと思います。

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');
h = plot(Lon(1),Lat(1),'ko');

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

% figureにテキストを設置
d = uicontrol(f2,'Style','text');
d.Position = [250 30 200 20];
d.String = 'Time 0';

% figureにスライダーを設置
c = uicontrol(f2,'Style','slider');
c.Position = [20 30 200 20];
c.Value = 1;
c.Min = 1;
c.Max = len;
c.SliderStep = [1/len 10/len];

% 6,figure更新
c.Callback = @(src,event)updatePosition(h,c,d,Lat,Lon,tbl);



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

function [] = updatePosition(h,c,d,Lat,Lon,tbl)
    i = uint32(c.Value);
    d.String = sprintf('Time %4d',tbl.timeoffset(i));
    h.XData = Lon(i);
    h.YData = Lat(i);
end

５、矢印を用いて走行時の向きと速度を可視化する

続いて、GPSデータに含まれていたCource情報を使用したいとおもいます。

MATLAB 
tbl = struct2table(xmlObj.trk.trkseg.trkpt);
len = length(tbl.time);
Speed = zeros(len,1);
Cource = zeros(len,1);

for i=1:len
    tbl.timeoffset(i) = i-1;
    Speed(i) = tbl.extensions(i).speed * 3.6;
    Cource(i) = tbl.extensions(i).course;
end

4, 9行目のように記述してtbl.extensions.courceに含まれていた値をCourceに代入しました。中身の確認のため、plotして確認しました。

MATLAB 
plot(Cource)
courceをplot

グラフを確認すると0~360までの値が出力されていることや走行ルートの向きから、北向きが0°で右回転を正とする方位角情報を表していることがわかります。

ここからは方位角情報を用いて、各走行時間における向きと、速度を矢印で表現します。

下記のようにArrowという配列を定義しました。矢印の大きさは速度に比例し、向きはCourceの値を使用します。

MATLAB 
tbl = struct2table(xmlObj.trk.trkseg.trkpt);
len = length(tbl.time);
Speed = zeros(len,1);
Cource = zeros(len,1);
Arrow = zeros(len,2);

for i=1:len
    tbl.timeoffset(i) = i-1;
    Speed(i) = tbl.extensions(i).speed * 3.6;
    Cource(i) = tbl.extensions(i).course;
    Arrow(i,1) = Speed(i)*sind(Cource(i))*30;
    Arrow(i,2) = Speed(i)*cosd(Cource(i))*30;
end

11,12行目は各時刻における矢印のベクトルを生成しています。矢印はSpeedに比例し、[0,1]を時計回りに回転させた向きになります。また、グラフ表示の都合上矢印が小さくなりすぎてしまうので大きさを30倍しました。

つづいて、先ほど「〇」で表現した現在位置を矢印に変更します。

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

h = quiver(Lon(1),Lat(1),Arrow(1,1),Arrow(1,2),'k','LineWidth',3);
h.MaxHeadSize = 5;

h = plot***をquiverを使って矢印を描画するように変更しています。最初の２つの引数が矢印の始点の位置で、3,4番目の引数が矢印のベクトルが入力となります。また適宜、見た目の調整を行っています。

この状態でスクリプトファイルを実行すると、原点付近の「〇」が小さい矢印に変更されていることがわかります。

このままだと、updatePositionでh.XData, h.YDataのみ更新しているため、矢印の位置が変わるだけで、大きさや向きが変わることがありません。ですのでupdatePositionを下記のように変更します。

MATLAB 
d.String = 'Speed 0 Time 0';

function [] = updatePosition(h,c,d,Lat,Lon,Speed,tbl,Arrow)
    i = uint32(c.Value);
    d.String = sprintf('Speed %2.1f Time %4d',round(Speed(i),1),tbl.timeoffset(i));
    h.XData = Lon(i);
    h.YData = Lat(i);
    h.UData = Arrow(i,1);
    h.VData = Arrow(i,2);

end

8,9行目のh.UData, h.VDataを使って矢印を表すベクトルの値を更新することが可能です。

1,5行目でテキストの表示にSpeedも追加しました。

updatePositionを呼び出す関数部分については入力引数にSpeedとArrowを追加したため、下記に変更しました。

MATLAB 
c.Callback = @(src,event)updatePosition(h,c,d,Lat,Lon,Speed,tbl,Arrow);

以上の変更を行い、スクリプトを実行すると下記のグラフが表示され、スライダーを使った値の更新ができます。速度と方向に応じて矢印の大きさと向きが変わることが確認できます。

GPS位置、速度、向きを表すグラフ

まとめ

これまで3記事にわたって、アップルウォッチで取得したGPSデータを用いて簡単なアプリケーションをMATLABで作成する方法をご紹介しました。業務や研究でも活用できるMATLABの使い方についても解説しました。

GPSの緯度経度情報はメートル座標系に変換することで直観的に理解できる有用なデータとなります。また、速度と向きを可視化することで、スポーツなどのデータ活用にも利用の範囲を広げることが可能です。

MATLAB/Simulinkに関する困りごとを解決するサービスをココナラで出品中です。自分一人では解決できない/誰かに相談したい等MATLABに関する困り事がありましたら、お見積もりは無料ですので、是非ご相談ください。

<!-- wp:code --><pre class="wp-block-code"><code><a class="coconala-widget" href="https://coconala.com/services/2403887" data-service_id="2403887" data-width="468" data-comment="1" data-invite="1" data-user_id="1994574">MATLAB初心者向け！環境構築からサポートします 自動運転エンジニアがMATLAB初心者のお悩みを解決！！</a><script>!function(d,s,id){var js,fjs=d.getElementsByTagName(s)[0],p=/^http:/.test(d.location)?'http':'https';if(!d.getElementById(id)){js=d.createElement(s);js.id=id;js.src=p+'://coconala.com/js/coconala_widget.js';fjs.parentNode.insertBefore(js,fjs);}}(document,'script','coconala-wjs');</script></code></pre><!-- /wp:code -->

今回作成したコードを下記に掲載しておきますので是非ご活用ください。

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;
    Cource(i) = tbl.extensions(i).course;
    Arrow(i,1) = Speed(i)*sind(Cource(i))*30;
    Arrow(i,2) = Speed(i)*cosd(Cource(i))*30;
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');%速度区間５

%矢印を描画（原点：現在位置、向き：方位角、大きさ：速度に比例）
h = quiver(Lon(1),Lat(1),Arrow(1,1),Arrow(1,2),'k','LineWidth',3);
h.MaxHeadSize = 5;

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

% figureにテキストを設置
d = uicontrol(f2,'Style','text');
d.Position = [250 30 200 20];
d.String = 'Speed 0 Time 0';

% figureにスライダーを設置
c = uicontrol(f2,'Style','slider');
c.Position = [20 30 200 20];
c.Value = 1;
c.Min = 1;
c.Max = len;
c.SliderStep = [1/len 10/len];

% figure更新
c.Callback = @(src,event)updatePosition(h,c,d,Lat,Lon,Speed,tbl,Arrow);

%２点の緯度経度から、距離を計算する関数
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

%sliderのボタンが押されたときに呼び出される関数
function [] = updatePosition(h,c,d,Lat,Lon,Speed,tbl,Arrow)
    i = uint32(c.Value);
    d.String = sprintf('Speed %2.1f Time %4d',round(Speed(i),1),tbl.timeoffset(i));
    h.XData = Lon(i);
    h.YData = Lat(i);
    h.UData = Arrow(i,1);
    h.VData = Arrow(i,2);

end
GPSデータ解析② | MATLAB事例5 前の記事

関連記事

  • MATLAB事例

複数ファイルの一括処理 | MATLAB事例紹介3

前回はスライダーを使ってfigureを更新する方法について２記事にわたって紹介しました。今回は研究や業務の現場でも非常によく使う複数のファイルを読み込んで同じ処理を実施する方法についてご紹介します。１ファイルについての処理を行うスクリプトなどを一生懸命作成したのち、同じ処理を複数のファイル

コメント

  1. Nice post! 1754776668