gpxファイルをMATLABで読み込む

gpxファイルはGPSロガーが出力するファイル形式です。MATLABで読めるのでメモをしておきます。 gpxreadやgeoshowはMapping Toolboxが必要です。

ここでは各時刻での移動距離(m)、速度(km/h)、標高(m)、勾配(degree)を計算しています。

勾配も計算していますが、誤差が大きいようです(そもそもGPSの標高データ(Elevation)の誤差が大きい)。 asinで計算していますので、引数が1を超えると複素数を返してきておかしくなります(期待しているものと違う値が返ってくる)。 そこで苦肉の策ですがスムージングをしています。 スムージングのパラメタはあくまでご参考程度です。

% gpxreadにはMapping Toolboxが必要
% Calculate disntace, velocity, and slope from gps data.
%% Read gpx file
%gps = gpxread('2020-05-08-18-24-51.gpx');
gpxFile = '2020-05-08-19-01-02.gpx';
g = gpxread(gpxFile);
lat = g.Latitude;
lon = g.Longitude;
alt = g.Elevation;

gps = struct;

%% Store altitude
altS = smooth(alt, 10, 'sgolay', 5);
gps.altitude = altS;
% 表示
geoshow(lat, lon, 'DisplayType','line', 'Color', 'red');
grid on;

%% Format time data
% GPSデータの時間列は次のような形式(ISO 8601): '2020-05-08T09:24:57.000Z'
% Tは日時と時刻を分ける
% ZはUTCを意味する(TZが違えば例えば +09:00 などになる。
%
time = g.Time;
timefmt = 'yyyy-MM-dd''T''HH:mm:ss.SSS''Z''';
tt = datetime(time , 'InputFormat', timefmt);

gps.time = tt - tt(1);

%% Caclulate distance and velocity
% 楕円体はWGS84を選択する
ellipsoid = wgs84Ellipsoid;

d = double.empty(length(g) , 0); % distance
t = double.empty(length(g) , 0); % time interval
v = double.empty(length(g) , 0); % velocity
a = double.empty(length(g) , 0); % altitude difference
s = double.empty(length(g) , 0); % slope in degree

% 原点はゼロにしておく
d(1, 1) = 0.0; t(1, 1) = 1.0; v(1, 1) = 0.0;
a(1, 1) = 0.0; s(1, 1) = 0.0;

for i = 1:(length(g)-1)
    d(i+1, 1)   = distance(lat(i), lon(i), lat(i+1), lon(i+1), ellipsoid);
    t(i+1, 1)   = seconds(tt(i+1) - tt(i));
    a(i+1, 1)   = altS(i+1) - altS(i);
end
gps.distance = d; % 

v = d ./ t * 3600 / 1e3;  % m/s to km/h.
gps.velocity = v;

%% Calculation of slope
dS = smooth(d(:,1), 10,  'sgolay', 5); % smooth distnace
id= find(dS > 1); % index of non-zero d.
for i = id
%    s(i, 1) = rad2deg(asin( a(i, 1) ./ d(i, 1)  ));
    s(i, 1) = ( a(i, 1) ./ dS(i, 1)  );
end
%%% Smooth the slope
ss = smooth(s(:,1), 20, 'sgolay', 5 );
ss = asin(ss); % in radian.
ss = rad2deg(ss); % convert to degree.

gps.slope = ss;

%% Save gps data.
save("gpsData.mat", 'gps');

%% Show data in Figure
figure
ax1 = subplot(4, 1, 1);
plot(d)
ylabel("distance [m]")
title("distance at every 1 second.")
ax2 = subplot(4, 1, 2);
plot(v)
ylabel("speed [km/h]")
title("speed (km/h) at every 1 second.")
ax3 = subplot(4, 1, 3);
plot(altS)
ylabel("altitude [m]")
title("altitude at every 1 second.")
ax4 = subplot(4, 1, 4);
plot(ss)
ylabel("slope [degree]")
title("slope in degree at every 1 second.")

linkaxes([ ax1, ax2, ax3, ax4], 'x')
savefig("gpsData.fig")