GPS Mountainbike analysis

GROUPS:
 Sander Huisman 19 Votes 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: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: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] 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] 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] 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] 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] 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] 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"]]] 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"]]]  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 Answer 2 years ago 18 Replies  Sander Huisman 8 Votes 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. Answer 2 years ago  Sander Huisman 2 Votes 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: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] The 10th kilometer was quite easy; it was straight on smooth asphalt in a tunnel... Answer 2 years ago  Henrik Schachner 3 Votes 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
2 years ago
 Indeed that also works quite nice! Nice trick for the TimeZone, I didn't think about it! Thanks for sharing as well!
2 years ago
 Vitaliy Kaurov 1 Vote Thanks, Sander, beautifully done! So do you have a way of estimating how many calories you've burned?
2 years ago
 Sander Huisman 1 Vote 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).
2 years ago
 Szabolcs Horvát 2 Votes Welcome to Lyon :-)
2 years ago
 Thanks! You're also in Lyon?! We should meet up some time, I've always wondered how to pronounce your name :-P
2 years ago
 Daniel Carvalho 3 Votes 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
2 years ago
 Sander Huisman 2 Votes 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?
2 years ago
 Sander Huisman 1 Vote 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...
2 years ago
 Daniel Carvalho 1 Vote Good idea! I gonna try and publish a new Demonstration about it! :-)
2 years ago
 Henrik Schachner 3 Votes 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: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:
2 years ago
 Thanks for sharing your code, many similarities! great job!
2 years ago