Code covered by the BSD License

# Processing gpx file

### towardthesea (view profile)

Using to read and process some information contained in a gpx file

clc;
clear all;
close all;

%Read and show track points from GPS on a figure
Lat = p.Latitude;
Lon = p.Longitude;
geoshow(Lat,Lon,'DisplayType','line','Color','red');
grid on;

%Calculate distance between two adjacent points of track line in the unit of
%degree

limLat = [min(Lat) max(Lat)];
limLon = [min(Lon) max(Lon)];
%zone = utmzone(limLat,limLon)
%ellipsoid = utmgeoid(zone); %define the local ellipsoid to calculate the line
%distance

ellipsoid = wgs84Ellipsoid;

for i=1:length(p)-1
arclen(i) = distance(Lat(i),Lon(i),Lat(i+1),Lon(i+1));
d(i) = distance(Lat(i),Lon(i),Lat(i+1),Lon(i+1),ellipsoid);
end

%cumulative_dis = cumsum(d)
total_dis = sum(d);

Time = p.Time;

timeStr = strrep(p.Time, 'T', ' ');
timeStr = strrep(timeStr, 'Z', '');
p.DateNumber = datenum(timeStr);
day = fix(p.DateNumber(1));
p.TimeOfDay = p.DateNumber - day;

t = datevec(p.DateNumber);
%t_beg = '2012-11-5 10:43:37';
%t_end = '2012-11-5 10:45:39';
t_beg = '2012-11-5 20:24:37';
t_end = '2012-11-5 20:30:51';

t1 = datevec(t_beg)
t2 = datevec(t_end)

%Find the nearest timestamp of t1
for k=1:size(t,1)-1
for j = 1:size(t,2)
if t1(j) < t(k+1,j) & t1(1:j-1) <=t (k+1,1:j-1)
t11 = t(k,:)
t12 = t(k+1,:)
k_sto = k
k = size(t,1)
break;
end
end
if k == size(t,1)
break;
end
end

%Find the nearest timestamp of t2
for l=size(t,1)-1:-1:2
for j = 1:size(t,2)
if t2(j) > t(l-1,j) & t2(1:j-1) >= t(l-1,1:j-1)
t21 = t(l-1,:)
t22 = t(l,:)
l_sto = l
l = 1
break;
end
end
if l == 1
break;
end
end

%Find the location of timestamp t1
Lat1 = (t1(6)-t11(6))/(t12(6)-t11(6)).* (Lat(k_sto+1) - Lat(k_sto)) + Lat(k_sto)
Lon1 = (t1(6)-t11(6))/(t12(6)-t11(6)).* (Lon(k_sto+1) - Lon(k_sto)) + Lon(k_sto)

%Find the location of timestamp t2
Lat2 = (t2(6)-t12(6))/(t22(6)-t12(6)).* (Lat(l_sto) - Lat(l_sto-1)) + Lat(l_sto-1)
Lon2 = (t2(6)-t12(6))/(t22(6)-t12(6)).* (Lon(l_sto) - Lon(l_sto-1)) + Lon(l_sto-1)

%Plot two above locations on map
hold on; plot(Lon1,Lat1,'b*',Lon2,Lat2,'bd')

%Calculate the distance between two locations

for i=k_sto:l_sto-1
%arclen(i) = distance(Lat(i),Lon(i),Lat(i+1),Lon(i+1));
dis_t1t2(i) = distance(Lat(i),Lon(i),Lat(i+1),Lon(i+1),ellipsoid);
end

dis_t1_t2 = sum(dis_t1t2) + distance(Lat1,Lon1,Lat(k_sto),Lon(k_sto),ellipsoid)+...
distance(Lat(l_sto),Lon(l_sto),Lat2,Lon2,ellipsoid);

sprintf('total distance is %7.2f metres',total_dis)
sprintf('distance between t1 = %s and t2 = %s is %7.2f metres',t_beg,t_end,dis_t1_t2)