Message Boards Message Boards

GPS Mountainbike analysis

Here I will show a way to do some analysis of GPS data using new functionality present in the Wolfram Language.

I just moved to Lyon (France), where they organised a mountainbike tour last week (Lyon Free VTT). I opted for the 'expert' trail of roughly 60 km with 1000+ meters ascend. I recorded my trip using the Garmin GPSmap 62s and I'm wearing a heart rate sensor from Garmin as well.

We start by setting the current directory to the notebookdirectoy, and look for the GPX file of the day of the event:

SetDirectory[NotebookDirectory[]];
fn = FileNames["*2015-09-13*.gpx"][[1]]

"Track_2015-09-13 121636.gpx"

So we found the file. The GPX file is an XML formatted file. So we can read the data easily:

data = Import[fn, "XML"];

The data is stored in "track points" which are abbreviated as trkpt. We can find the parts of the track by using Cases.

FirstCase[data, XMLElement["trkpt", ___], Missing[], \[Infinity]]

which looks like:

XMLElement[trkpt,{lat->45.7242092583,lon->4.8268832266},
{XMLElement[ele,{},{173.17}],
 XMLElement[time,{},{2015-09-13T06:09:03Z}],
 XMLElement[extensions,{},{
  XMLElement[{http://www.garmin.com/xmlschemas/TrackPointExtension/v1,TrackPointExtension},{},{
  XMLElement[{http://www.garmin.com/xmlschemas/TrackPointExtension/v1,hr},{},{108}]}
]}
]}
]

Wow! Looks quite complicated, but the important thing is that it always starts with the Latitude and Longitude coordinates as attributes -> values, and then followed by some more data (in this case the time, elevation, and heart rate).

So lets get all the track points:

data = Cases[data, XMLElement["trkpt", {"lat" -> lat_, "lon" -> lon_}, other_] :> {ToExpression /@ {lat, lon}, other}, \[Infinity]]

Note that I immediately detect the latitude, longitude and the other data and put them as a replacement rule right in Cases.

I started recording when I left my house to go to the Start, and I also continued recording after I cycled back from the finish to my house. So we remove this data:

data = Drop[Drop[data, 246], -128];

The trail is now stored as the first element of each element of data:

trail = GeoPosition /@ data[[All, 1]];

We can now easily visualise this using the new GeoGraphics function (V10 and up):

GeoGraphics[{Thick, Line[trail]}, ImageSize -> 400]

Giving:

enter image description here

That looks pretty good! Now let's retrieve the time, elevation and heart rate from the data. We overwrite data, because we already extracted the latitude and longitude from it:

data = data[[All,2]];

First, let's try the time, we use the new FirstCase command (V10) as we only expect one time per track point entry. We make a pure function and map that over all the data:

absolutetime = time = Map[
FirstCase[#,XMLElement["time",{},{time_}]:>AbsoluteTime[time]+$TimeZone*3600,Missing[],\[Infinity]]&
,
data
] 
time -= First[time]

Note that I also add some seconds so that I have the time in my own timezone (the default is in Zulu (universal) time). In case the time is somehow not recorded it will returning Missing[]. I will also save a variable which is the time starting when I start called time (second line).

We can do something similar for the elevation:

elevation=Map[FirstCase[#,XMLElement["ele",{},{ele_}]:>ToExpression[ele],Missing[],\[Infinity]]&,data]

and for the heart rate:

hr=Map[FirstCase[#,XMLElement[{"http://www.garmin.com/xmlschemas/TrackPointExtension/v1","hr"},{},{hr_}]:>ToExpression[hr],Missing[],\[Infinity]]&,data];

Now that we have all the data, let's explore a little bit further. Let's first compute the distance travelled:

distances=BlockMap[GeoDistance@@#&,trail,2,1];
distances=Prepend[QuantityMagnitude[distances,"Meters"],0.0];
distance=Accumulate[distances];

Here I used the new BlockMap function (V10.2) to compute the distance between each pair of points. (2 to denote pairs, and 1 to have overlapping pairs). I convert the results to just number (without units), otherwise the Accumulate is very slow (why is this slow Wolfram?). Lastly I keep track of the accumulated distance using the command Accumulate.

Now we can plot this as a function of time:

DateListPlot[{absolutetime,distance}\[Transpose],FrameLabel->"Meters",FrameStyle->12,ImageSize->400]

giving:

enter image description here

The total distance travelled is:

Last[distance]
57106.

This answer is in meters because I took the magnitude in "Meters". We can convert it to miles quite easily:

UnitConvert[Quantity[%, "Meters"], "Miles"]
35.484mi

About 35.5 miles.

Let's calculate the speed:

speed=Differences[distance]/Differences[time];
speed=Join[{speed[[1]]},MovingAverage[speed,2],{speed[[-1]]}];

The first order difference will give me the speed, I'm doing a moving average so that I have the average speed at an instant, I append and prepend the first and last speed, respectively, such as to have the length of speed and the length of time to be the same.

Now we can plot the speed (which is in meter/second, so I multiply by 3.6 to get km/h):

DateListPlot[Transpose[{absolutetime,3.6speed}],FrameLabel->"Kilometer/Hour",FrameStyle->12,ImageSize->1000,AspectRatio->1/5,PlotRange->All]

enter image description here

My speed varies a lot as you can see. This is partly because this is an urban trail with climbs and descends, but also sharp corners, and up and down-going stairs. But also because of some GPS errors you have some fluctuations. This is apparent in the above plot, I'm quite fast, but not 80 km/h (50 miles/hour) fast!! So let's smooth the data over a few points:

slen=3;
smoothspeed=Join[Accumulate[speed[[;;slen]]]/Range[slen],MovingAverage[speed,2slen+1],Reverse[Accumulate[Reverse[speed[[-slen;;]]]]/Range[slen]]];

Which looks more realistic:

DateListPlot[Transpose[{absolutetime,3.6smoothspeed}],FrameLabel->"Kilometer/Hour",FrameStyle->12,ImageSize->1000,AspectRatio->1/5,PlotRange->All]

enter image description here

You can clearly see that I had a flat tire around 10:15 that took me 7-8 minutes to repair...

Let's calculate the average speed (up to some time t):

avgspeeds=Prepend[Rest[distance]/Rest[time],0];

And showing:

DateListPlot[Transpose[{absolutetime,3.6avgspeeds}],FrameLabel->"Kilometer/Hour",FrameStyle->12,ImageSize->1000,AspectRatio->1/5,PlotRange->All]

enter image description here

So I ended with a speed of 3.6 Last[avgspeeds] = 14.5 km/h which is not too bad considering it includes also a flat tire.

Plotting the elevation is also quite easy now:

DateListPlot[Transpose[{absolutetime, elevation}], FrameLabel -> "Meters", FrameStyle -> 12, ImageSize -> 1000, AspectRatio -> 1/5, PlotRange -> All]

enter image description here

There were some very steep hills in this track, and pieces where I had to change to my lowest gear to climb them.

Also plotting the heart rate is quite easy:

DateListPlot[Transpose[{absolutetime,hr}],FrameLabel->"BPM",FrameStyle->12,ImageSize->1000,AspectRatio->1/5,PlotRange->All]

enter image description here

Let's compute the total number of heart beats and the average heart rate:

hb=Interpolation[DeleteMissing[{time,hr/60.0}\[Transpose],1,2],InterpolationOrder->1];
hb=Round[Integrate[hb[t],{t,Min[time],Max[time]}]]
avghr=60hb/(Max[time]-Min[time])

38443 (* beats *)
162.8 (* average beats / minute *)

That is quite high for my age (28 years)!

Let's compute the elevation gain:

elevationgain = Total[Select[Differences[elevation], Positive]]

1193.05 (* meters *)

That explains partly the low speeds and the high heart rates! I've climbed almost 1200 meters.

QuantityMagnitude[Quantity[elevationgain,"Meters"],"Stories"]

or about 400-500 stories!! (depending on your metric of a storey)

Let's make a visual plot where the climbs are by creating a map colored by the height:

colordata = Rescale[elevation];
color = Blend[{{0, Darker@Green}, {0.5, Yellow}, {1, Red}}, #] & /@ colordata;
GeoGraphics[{AbsoluteThickness[8], JoinForm["Round"], Line[trail, VertexColors -> color]}, ImageSize -> 400]

enter image description here

Green means low, and red means high altitude. Using the same colors we can create also the accompanying plot:

DateListPlot[Transpose[{absolutetime,elevation}],FrameLabel->"Meters",FrameStyle->12,ImageSize->1000,AspectRatio->1/5,PlotRange->All,ColorFunction->(Blend[{{0,Darker@Green},{0.5,Yellow},{1,Red}},#2]&),PlotStyle->Directive[AbsoluteThickness[5],JoinForm["Round"]]]

enter image description here

Or the same but instead colored by the speed:

colordata=Rescale[smoothspeed];
cf=Blend[{{0.2,Red},{0.45,Yellow},{0.8,Darker@Green}},#]&;
color=cf/@colordata;
g1=GeoGraphics[{AbsoluteThickness[8],JoinForm["Round"],Line[trail,VertexColors->color]},ImageSize->400]
DateListPlot[Transpose[{absolutetime,3.6smoothspeed}],FrameLabel->"km/h",FrameStyle->12,ImageSize->1000,AspectRatio->1/5,PlotRange->All,ColorFunction->(cf[#2]&),PlotStyle->Directive[AbsoluteThickness[5],JoinForm["Round"]]]

enter image description here enter image description here

In the corners and in the climbs it is slow (red), while for the flat, down, and straight pieces it is fast (green). It looks like there is way more green, then red on the map, while for the bottom plot it is roughly equal, this is because for a low speed you travel only a little so just a bit of red on the map. While for high speed you travel more distance so it looks more green...

This is part I, below I will post a second part of the analysis. With some nice general functionality. PS. Also have a look at the 2009 post at the wolfram blog when GeoGraphics (and other functions) were not introduced yet. Wolfram Blog about GPS
POSTED BY: Sander Huisman
18 Replies
Sort By:

Using the code above I made an animation: https://www.youtube.com/watch?v=zf6xPomRZzM

enter image description here

enter image description here
POSTED BY: Sander Huisman

Nice ! I added some small teaser-animation.

POSTED BY: Vitaliy Kaurov

Perfect! I couldn't find an easy way to generate a gif from an MP4. I had already deleted the original frames... What is your modus operandi?
POSTED BY: Sander Huisman

When I need something quick I just capture from screen with LICEcap.
POSTED BY: Vitaliy Kaurov

During quite some years of biking I have a big collection of GPS tracks. In order to get a quick idea of any specific track I wrote some code which generates a graphic file (.png) showing the significant info. It looks like this:

enter image description here

This file is then located in the same directory and with the same filename as the GPX file (but of course with a different extension!). Well, most comments and descriptions are in German, and certainly there is much more room for improvement (e.g. when reading in more than four GPX files at a time the kernel quits - maybe somebody finds the reason). The nice feature is that the altitude is coded in color. I am using a Garmin Oregon 600t which has access to barometric pressure; the altitude data are contained in GPX file.

Regards -- Henrik

Attachments:
POSTED BY: Henrik Schachner

Thanks for sharing your code, many similarities! great job!
POSTED BY: Sander Huisman

Very well done!!! Thanks for sharing your code! :-)

It would be interesting the create a service that given the GPS data it returns the map with the colored track by height.

It reminds me a Demonstration I did from car gps data:Gaps In GPS Data
POSTED BY: Daniel Carvalho

In case your GPS device does not store the elevation or is bad at it (some devices don't have a pressure sensor, so it solely relies on GPS for the height, which is generally of bad quality) you can upload it to gpsies.com, and then download it again. It will add the terrain-elevation to it. (a trick i just found out this Sunday)

I'm not sure if the wolfram cloud allows for uploading files to the cloud. But if so, then it should not be too difficult right?
POSTED BY: Sander Huisman

In addition you could query GeoElevationData with a list of all your points. I'm not sure what the resolution is of the data, and won't work in tunnels/caves...

POSTED BY: Sander Huisman

Good idea! I gonna try and publish a new Demonstration about it! :-)
POSTED BY: Daniel Carvalho

Welcome to Lyon :-)
POSTED BY: Szabolcs Horvát

Thanks! You're also in Lyon?! We should meet up some time, I've always wondered how to pronounce your name :-P
POSTED BY: Sander Huisman

Thanks, Sander, beautifully done! So do you have a way of estimating how many calories you've burned?
POSTED BY: Vitaliy Kaurov

Computing how many calories I have burned is very tricky to be reliable. Of course I could do something like:

WolframAlpha[StringTemplate["cycling `1` minutes `2` kilometer 70 kilogram male"][240, 60], {{"MetabolicProperties", 1}, "ComputableData"}]

and get the "energy expenditure" value from the array it returns. But this does not really include the elevation climb nor does it account for surface (was it muddy, sandy, stones, or clean smooth asphalt), nor the weather conditions (temperature, humidity, as well as wind).

There are roughly 3 types of power one has to provide, 1 for the friction with the road, 1 for the friction with the air, 1 for climbing.

$P(t) = \mu(t)mv(t)+C_D v(t)^3 + mgv(t)Tan[\theta(t)]$

where $\mu$ is the friction coefficient as a function of time (depending on the type of terrain), $C_D$ the drag coefficient with the area of myself and other prefactors included already, $g$ acceleration of gravity, and $m$ my mass, and $\theta$ the slope. The total power would then be the integral of that quantity over time. The you still have to include that your muscles are not very efficient (20% or so). So you'll burn 5x the energy... For me, there is too much guessing to give a reliable answer...

Note that the last term cancels out if the track is a closed contour (arrive and depart from the same point) (there is as much ascend as there is descend).
POSTED BY: Sander Huisman

Hi Sander,

thank you very much for sharing! For quite a while I am doing analysis of my GPS tracks as well (I am using a Garmin Oregon 600t), but not on a that elaborated level as you do. So here comes just a little remark: I myself found it more useful/convenient to import the gpx-file like so:

gpx = Block[{$TimeZone = 0}, Flatten@Import[gpxFileName, "Data"]];

You are right: .gpx is a XML Format, but nevertheless one does not need to import it as such.

Regards -- Henrik
POSTED BY: Henrik Schachner

Indeed that also works quite nice! Nice trick for the TimeZone, I didn't think about it! Thanks for sharing as well!
POSTED BY: Sander Huisman

Now that I have this Import Function I can easily make more plots:

I can read my data in now like this:

SetDirectory[NotebookDirectory[]];
fns=FileNames["*2015-09-13*.gpx"];
data=ImportGarminGPX[fns[[1]],{246,128},2];

Now data will be an association with all kinds of keys. We can now make a graph to display the time for each kilometer:

plotdata=Table[data["InverseDistanceFunction"][x],{x,0,data["TotalDistance"],1000}];
AppendTo[plotdata,data["InverseDistanceFunction"][data["TotalDistance"]]];
plotdata=Differences[plotdata];
plotdata={Range[Length[plotdata]]-1,plotdata}\[Transpose];
plotdata=Append[plotdata,Last[plotdata]+{1,0}];
ListPlot[plotdata,ImageSize->500,Joined->True,Frame->True,Filling->Axis,PlotRange->{{0,All},{0,All}},PlotRangePadding->None,GridLines->{None,Automatic},ClippingStyle->Directive[Dotted,Red],FrameLabel->{"km #","Time [sec]"},InterpolationOrder->0]

Giving:

enter image description here

So my quickest kilometer (the 10th kilometer) was roughly in 100 seconds; not bad! But I don't know where that is on the map! So let's update our map with km-marks in it!

tr=data[["Trail"]];
ele=data[["Elevation"]];
colors=Blend[{{0,Darker@Green},{0.5,Yellow},{1,Red}},#]&/@Rescale[ele];
pts=Table[data["TrailFunction"][data["InverseDistanceFunction"][x]],{x,0,data["TotalDistance"],1000}];
pts=MapIndexed[Text[Style[First[#2]-1,16],#1]&,pts];
g1=GeoGraphics[{AbsoluteThickness[12],Line[tr,VertexColors->colors],Black,pts},ImageSize->600]

enter image description here

The 10th kilometer was quite easy; it was straight on smooth asphalt in a tunnel...
POSTED BY: Sander Huisman

Here is an Import function I made for my GPX files that computes various quantities:

ClearAll[ImportGarminGPX,DuplicatePositions]
DuplicatePositions[x_List]:=Module[{nums},
    nums=Select[Tally[x],Last[#]>1&][[All,1]];
    Flatten[Rest/@(Position[x,#]&/@nums)]
]
ImportGarminGPX[fn_String,ignorepts:{n_Integer,m_Integer}:{0,0},timezone_:2]:=Module[{data,hr,trail,rawtrail,elevation,elevationunits,abstime,time,date,dups,distances,distance,speed,slen,smoothspeed,avgspeeds,avgspeed,laf,lof,trailfunction,elevationgain,maxhr,hb,avghr,runningmin,runningmax,biggestclimb,biggestdrop},
data=Import[fn,"XML"];

data=Cases[data,XMLElement["trkpt",{"lat"->lat_,"lon"->lon_},other_]:>{ToExpression/@{lat,lon},other},\[Infinity]];
data=Drop[Drop[data,n],-m];
trail=rawtrail=data[[All,1]];
trail=GeoPosition/@trail;

data=data[[All,2]];

abstime=time=Map[FirstCase[#,XMLElement["time",{},{time_}]:>AbsoluteTime[time]+timezone 3600,Missing[],\[Infinity]]&,data];
date=DateObject[First[time]];
time-=First[time];

dups=List/@DuplicatePositions[time]; (* delete duplicate times *)
If[dups=!={},
    Print["duplicate time detected removing indices:",dups];
    time=Delete[time,dups];
    abstime=Delete[abstime,dups];
    trail=Delete[trail,dups];
    rawtrail=Delete[rawtrail,dups];
    data=Delete[data,dups];
];

elevationunits=elevation=Map[FirstCase[#,XMLElement["ele",{},{ele_}]:>ToExpression[ele],Missing[],\[Infinity]]&,data];
elevationunits=Quantity[elevation,"Meters"];

runningmin=FoldList[Min,First[elevation],elevation];
runningmax=Reverse[FoldList[Max,First[Reverse[elevation]],Reverse[elevation]]];
biggestclimb=Max[runningmax-runningmin];

runningmin=Reverse[FoldList[Min,First[Reverse[elevation]],Reverse[elevation]]];
runningmax=FoldList[Max,First[elevation],elevation];
biggestdrop=Max[runningmax-runningmin];

hr=Map[FirstCase[#,XMLElement[{"http://www.garmin.com/xmlschemas/TrackPointExtension/v1","hr"},{},{hr_}]:>ToExpression[hr],Missing[],\[Infinity]]&,data];

distances=BlockMap[GeoDistance@@#&,trail,2,1];
distances=Prepend[QuantityMagnitude[distances,"Meters"],0.0];
distance=Accumulate[distances];

speed=Differences[distance]/Differences[time];
speed=Join[{speed[[1]]},MovingAverage[speed,2],{speed[[-1]]}];
slen=3;
smoothspeed=Join[Accumulate[speed[[;;slen]]]/Range[slen],MovingAverage[speed,2slen+1],Reverse[Accumulate[Reverse[speed[[-slen;;]]]]/Range[slen]]];
avgspeeds=Prepend[Rest[distance]/Rest[time],0];
avgspeed=Last[avgspeeds];


laf=Interpolation[{time,rawtrail[[All,1]]}\[Transpose],InterpolationOrder->1];
lof=Interpolation[{time,rawtrail[[All,2]]}\[Transpose],InterpolationOrder->1];
trailfunction=Function[t,GeoPosition[{laf[t],lof[t]}]];
elevationgain=Total[Select[Differences[elevation],Positive]];

If[Length[DeleteMissing[hr]]>0,
    maxhr=Max[DeleteMissing[hr]];
    hb=Interpolation[DeleteMissing[{time,hr/60.0}\[Transpose],1,2],InterpolationOrder->1];
    hb=Integrate[hb[t],{t,Min[time],Max[time]}];
    avghr=60hb/(Max[time]-Min[time]);
,
    maxhr=Missing[];
    hb=Missing[];
    avghr=Missing[];
];
<|
"AbsoluteTime"->abstime,
"AverageHeartRate"->Round[avghr,0.1],
"AverageSpeed"->avgspeed,
"AverageSpeeds"->avgspeeds,
"BiggestClimb"->biggestclimb,
"BiggestDrop"->biggestdrop,
"Date"->date,
"DateString"->DateString[date,{"Day"," ","MonthName"," ","Year"}],
"Distance"->distance,
"DistanceFunction"->Interpolation[{time,distance}\[Transpose],InterpolationOrder->1],
"Distances"->distances,
"Elevation"->elevation,
"ElevationGain"->elevationgain,
"ElevationUnits"->elevationunits,
"Heartbeats"->hb,
"HeartRate"->hr,
"InverseDistanceFunction"->Interpolation[DeleteDuplicatesBy[{distance,time}\[Transpose],First],InterpolationOrder->1],
"MaximumHeartRate"->maxhr,
"MaxTime"->Max[time],
"RawTrail"->rawtrail,
"SmoothSpeed"->smoothspeed,
"Speed"->speed,
"SpeedFunction"->Interpolation[{time,speed}\[Transpose],InterpolationOrder->1],
"Speedkmh"->3.6 speed,
"Time"->time,
"TotalDistance"->Last[distance],
"Trail"->trail,
"TrailFunction"->trailfunction
|>
]

Here it also calculates various other properties:

AbsoluteTime
AverageHeartRate
AverageSpeed
AverageSpeeds
BiggestClimb
BiggestDrop
Date
DateString
Distance
DistanceFunction
Distances
Elevation
ElevationGain
ElevationUnits
Heartbeats
HeartRate
InverseDistanceFunction
MaximumHeartRate
MaxTime
RawTrail
SmoothSpeed
Speed
SpeedFunction
Speedkmh
Time
TotalDistance
Trail
TrailFunction

Note that the GPX format has a lot of variations and extensions and I can't guarantee it works on your particular GPX file. I have several GPX files from different devices or applications (from e.g. a smart phone), and they are all structurally a little bit different. I tried using FirstCase[....Missing[]..] everywhere, and I built in some safe-guards but it not completely hooligan-proof! You might have to change the code slightly to adapt it to your files.
POSTED BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract