Find distance along a path (GPS coordinates)

Question

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.

Answer 1

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]]

Answer 2

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.

EDIT (to add full deployment)

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"]];

thread = Prepend[Thread[{ids, distances}], {"ID", "Distance"}];

Export["distances.csv", thread, "CSV"]