# Find distance along a path (GPS coordinates)

I have an imported version of a GPX file plotting the route of a trail. I also have a list of coordinates for waypoints along the route. I am trying to calculate how far along the route each waypoint is.

The cleaned up route data can be here: https://pastebin.com/8KDcvMex and can be imported directly:

ClearAll[data];
ToExpression@Import["https://pastebin.com/raw/8KDcvMex", "String"];


Using that data I then used the following code to plot the trail:

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


Then I bring in the waypoint I want to find the location of:

waypoint = GeoPosition[{34.56544401, -77.90251801}]


I can find the first point on the route with this:

start = First[trail]


And find the straight line distance with this:

GeoDistance[waypoint, start]


But how do I have Mathematica measure that distance along the path instead of straight distance?

Additionally, I have found I can find the closest route point with this:

GeoNearest[trail, waypoint]


And I can find out far the waypoint is from one of the route points with this:

Min[GeoDistance[waypoint, trail]]


In looking for possible solutions, I found this example on the Wolfram Community and he does some cool stuff using Accumulate as a way of calculating overall distance. I could use that function combined with finding the nearest route point and have the answer, but that won't always work. You see, there are gaps in some of the GPX route coordinates where at times there is up to half a mile between coordinates. I have cases where waypoints are returning the same route point as closest because they really are the closest in the dataset, even though the waypoints themselves are separated by some distance. In cases like that, my distance calculations could be off by a pretty substantial bit.

Is there another way to do what I am trying to do, or am I stuck with either the drawbacks of finding the distances to the nearest route point or manually entering the waypoint coordinates into the appropriate spot within the route points?

EDIT:

Using techniques described in the Wolfram Community post above, I was able to approximate the result, but it still doesn't account for the distance between the waypoint and the nearest route point.

Here is what I did:

First I calculated the cumulative distance of each route point from the starting point (I stripped units for simplicity).

distances =
Accumulate[Prepend[QuantityMagnitude[GeoDistanceList[trail]], 0]];


Then I combined the two datasets:

combined = Thread[{trail, distances}]


To extract the relevent data, I used:

nearest =
GeoNearest[combined[[All, 1]],
waypoint][[1]]; Select[combined, #[[1]] == nearest &]


Again, this gets close since the output is as expected with the geoposition of the nearest route point and the distance to that route point, but isn't accurate enough.

A couple notes on the code above. Remember that Mathematica returns GeoNearest[combined[[All, 1]], waypoint] as a single item list so you will either need to call the first item with [[1]] as I did below or use First@GeoNearest[combined[[All, 1]], waypoint]. Also, if I don't precalculate nearest and use the following code, Mathematica gets the right answer, but takes several minutes compared to near instant results:

Select[combined, #[[1]] ==
GeoNearest[combined[[All, 1]], randomblaze][[1]] &]


Definitely something in code optimization that I am missing.

• I can include the actual route data if that is helpful, but because of it's nature, it is over 2000 points. – kickert Jun 27 '18 at 18:00
• yes, please include the route and waypoints. If they are too big to fit in the question, you can share the file, or a text only Mathematica expression on pastebin or some similar service. – MarcoB Jun 27 '18 at 18:09
• @MarcoB, Updated for you. – kickert Jun 27 '18 at 18:24
• I was talking with a friend who works in GIS and he was able to take a GPX route, convert it to a line, and then export a new version that had route points spaced every 10 meters. That data source would effectively solve the problem since the margin of error would be +/- 5m instead of +/- 500m. Is there a way to redefine the datapoints in a similar manner in Mathematica? – kickert Jun 28 '18 at 2:12
• That’s an interesting question. It seems essentially equivalent to your original problem though, so if we can solve one, we can probably solve the other as well. Have you seen GeoPath? I’ve been playing with it a bit, seems promising, but not really there yet. Your edited approach is interesting too. – MarcoB Jun 28 '18 at 2:17

Given these results from the question, we want the distance following the path to the waypoint. Edit: There are 97 duplicate points: Length@trail - Length@DeleteDuplicates[trail]. The commented line removes them.

data=ToExpression@
Import["https://pastebin.com/raw/8KDcvMex", "String"];
trail = GeoPosition /@ data[[All, 1]];
(*trail = DeleteDuplicates[GeoPosition/@data[[All,1]]];*)
waypoint = GeoPosition[{34.56544401, -77.90251801}];
trailNearest = First@GeoNearest[trail, waypoint];


First, find some distances to help validate the final answer. The straight-line distance from the start of trail to trailNearest is:

GeoDistance[First@trail, trailNearest]
(* 31.1795 mi *)


The total distance along the trail is:

Total[GeoDistanceList[trail]]
(* 66.2333 mi *)


Check this result against the geodesic distance using GeoPath and GeoLength:

Total[GeoDistanceList[trail]] - GeoLength@GeoPath[trail, "Geodesic"]
(* -2.45927*10^-9 mi *)


Close enough.

We'd like know to the index of trailNearest along the points on the trail. PositionIndex[trail] assigns an index to each of the 2,128 coordinates on the trail. Then PositionIndex[trail][trailNearest] gives the index of trailNearest on the trail.

indexNearest = First@PositionIndex[trail][trailNearest]
(* 1714 *)


Now, finding the distance following the trail to the point nearest the waypoint is easy!

Total@GeoDistanceList[trail[[1 ;; indexNearest]]]
(* 54.1951 mi *)


Adding the distance from trailNearest to the waypoint gives the required distance along the trail to the waypoint.

Total@GeoDistanceList[trail[[1 ;; indexNearest]]] + Min[GeoDistance[waypoint, trail]]
(* 54.2081 mi *)


This is consistent with the straight-line distance, 31.1795 mi, and the total length of the trail, 66.2333 mi. Here's the path along the trail to the point nearest the waypoint (colored red), which looks credible.

Show[GeoListPlot[trail, PlotStyle -> Green],
GeoListPlot[trail[[1 ;; indexNearest]], GeoBackground -> None]]


Summarizing, to find the distance along the path to a waypoint:

trailNearest = First@GeoNearest[trail, waypoint];
indexNearest = First@PositionIndex[trail][trailNearest];
Total@GeoDistanceList[trail[[1 ;; indexNearest]]] + Min[GeoDistance[waypoint, trail]]

• I like what you did with PositionIndex and GeoDistanceList. It is cleaner than my use of Accumulate above. However it still doesn't address the issue of solving for gaps in the route points. If the waypoint was 400m before the nearest route point, then the actual distance would be off by +800m. In this dataset Max[GeoDistanceList[trail]] shows more than half a mile. – kickert Jun 28 '18 at 11:24
• I see. On the other hand, very few distances exceed 0.1 miles: Histogram[Normal@GeoDistanceList[trail], {.1}]. Could interpolate between points that exceed a threshold, but checking for problem waypoints might be simple considering how few there are. – creidhne Jun 28 '18 at 11:45
• In this real world dataset there were 14 times the separation exceeds 0.25 miles with the highest being 0.54. – kickert Jun 28 '18 at 11:52
• Yes, 14 out of 2,128. But considering that priorpoint is further away from the waypoint than trailNearest, does that really improve the error? I also found feet, but I conspired to avoid the problem. – creidhne Jun 28 '18 at 11:59
• Those 14 route points cover over 6% of the trail length, so it gets a bit more significant. 10% of the trail is on routepoints over 0.2 miles apart. In terms of accuracy, selecting a route point further up the trail is always going to introduce error. If I am measuring total distance to a point 200 meters in front of me and then adding my distance to that point to a distance that is already too large, then the calculation would be off by 400 meters. With this dataset, worst case scenario would give a distance over half a mile too long. When you are hiking that is significant. – kickert Jun 28 '18 at 12:24

The answer from @creidhne using IndexPosition and GeoDistanceList is elegant and works for a vast majority of cases. However, in instances where the spacing between route points is significant (up to half a mile in real world worst case scenarios), and the nearest route point is further up the trail, the reported distance will be over by twice the distance from the waypoint to the route point since you are essentially calculating the distance to double back.

Using the code from @creidhne as a base, I was able to modify it so that it finds the nearest route point prior to the waypoint and then adds the remaining distance. This should give the right answer 100% of the time since it removes any chance of reporting distance from doubling back from a far off route point up the trail.

Here is the complete code:

Bring in the data:

ClearAll[data];
ToExpression@Import["https://pastebin.com/raw/8KDcvMex", "String"];
trail = GeoPosition /@ data[[All, 1]];
waypoint = GeoPosition[{34.56544401, -77.90251801}];


Calculate the Distance

priorpoint =
Min[PositionIndex[trail] /@ GeoNearest[trail, waypoint, 2]];
Total@GeoDistanceList[trail[[1 ;; priorpoint]]] +
UnitConvert[GeoDistance[waypoint, trail[[priorpoint]]], "Miles"]


Note: I had to include UnitConvert since my answer was coming out in feet instead of miles.

As I mentioned in my initial question, I was actually looking at performing this calculation for a list of waypoints. Below is the full code deployment using an imported GPX file and an XLSX file where Waypoint IDs were in column 6, Lats were in column 12 and Longs were in Column 13. It can then be downloaded as a CSV file.

gpx = Import[
"route.gpx", "XML"];
data =
Cases[gpx,
XMLElement["trkpt", {"lat" -> lat_, "lon" -> lon_},
other_] :> {ToExpression /@ {lat, lon}}, \[Infinity]];
trail =
GeoPosition /@ data[[All, 1]];

waypointdata =
Import["waypoints.xlsx", {"Data"}][[1, All, {6, 12, 13}]] //
Rest; ids = waypointdata[[All, 1]];
lat = waypointdata[[All, 2]];
long = waypointdata[[All, 3]];
locations =
GeoPosition /@ waypointdata[[All, {2, 3}]];

distances =
QuantityMagnitude[
UnitConvert[
Total[GeoDistanceList[
trail[[1 ;;
Min[PositionIndex[trail] /@ GeoNearest[trail, #, 2]]]]]] +
GeoDistance[#,
trail[[ Min[
PositionIndex[trail] /@ GeoNearest[trail, #, 2]]]]] & /@
locations, "Miles"]];