# Finding position at given distance in a GeoPath

Is there a way to find the geoposition of a given distance from start in a GeoPath? I want to mark equidistant positions along a track, for example, a mark every 500 km along the path given by

path=GeoGraphics[
GeoPath[{
Entity["City", {"Boston", "Massachusetts", "UnitedStates"}],
Entity["City", {"Rochester", "NewYork", "UnitedStates"}],
Entity["City", {"Chicago", "Illinois", "UnitedStates"}]
}, "Rhumb"]
]


Is there a way to find the pos that gives GeoDistance[Entity["City", {"Boston", "Massachusetts", "UnitedStates"}],pos]==500 km along the path?

## 2 Answers

Here's my attempt to parametrize the path by arclength, where here arclength is GeoLength.

First I build up a function that can be used on many values:

ParametrizeGeoPath[g_GeoPath, t_] := ParametrizeGeoPath[g][t]

ParametrizeGeoPath[GeoPath[locs_, args___]] :=
Block[{line, nodes, lens, acc, nf, n1, n2, solver},
line = GeoGraphicsGeoEvaluate[GeoPath[locs, args]];
nodes = GeoPosition /@ Reverse[line[[1]], {2}];
lens = QuantityMagnitude[UnitConvert[GeoLength[GeoPath[#, args]]& /@ Partition[nodes, 2, 1], "Kilometers"]];
acc = Accumulate[lens];
nf = Nearest[acc -> {"Index", "Element"}];

GeoPathParametricFunction[acc, nodes, nf, args]
]


Given a target distance, we can invert GeoLength with FindRoot:

GeoPathParametricFunction[_, nodes_, __][d_] /; d == 0 := First[nodes]

GeoPathParametricFunction[acc_, nodes_, __][d_] /; d == Last[acc] := Last[nodes]

GeoPathParametricFunction[acc_, nodes_, nf_, args___][d_] /; 0 < d < Last[acc] :=
Block[{i, v, s, n1, n2, dist, root, t = Unique["t"]},
{i, v} = First[nf[d]];
If[v > d, i--];
s = If[i == 0, 0, acc[[i]]];
n1 = nodes[[i+1, 1]];
n2 = nodes[[i+2, 1]];

dist[t_?NumericQ] := s + QuantityMagnitude[UnitConvert[GeoLength[GeoPath[{GeoPosition[n1], GeoPosition[(1-t)n1 + t n2]}, args]], "Kilometers"]];

root = Quiet@FindRoot[dist[t] == d, {t, .5, 0, 1}];

(
GeoPosition[(1-t)n1 + t n2] /. root

) /; MatchQ[root, {t -> _Real}]
]

GeoPathParametricFunction[___][d_?NumericQ] = Indeterminate;

GeoPathParametricFunction /: MakeBoxes[expr:GeoPathParametricFunction[__], _] := InterpretationBox[RowBox[{"GeoPathParametricFunction", "[", "\"\[Ellipsis]\"", "]"}], expr]


Your example:

path = GeoPath[{
Entity["City", {"Boston", "Massachusetts", "UnitedStates"}],
Entity["City", {"Rochester", "NewYork", "UnitedStates"}],
Entity["City", {"Chicago", "Illinois", "UnitedStates"}]
}, "Rhumb"];

gpf = ParametrizeGeoPath[path];

gpf[500]

 GeoPosition[{43.0932, -77.0359}]

Manipulate[GeoGraphics[{
path,
GeoMarker[gpf[d]]
},
PlotLabel -> Row[{"Distance: ", Quantity[d, "Kilometers"]}]],
{d, 0, gpf[[1, -1]]}
]


The points returned are very close to the initial path:

ListLinePlot[
GeoDistance[path, g /@ Range[0, 1300, 100]],
TargetUnits -> "Meters",
AxesLabel -> Automatic,
DataRange -> {0, 1300}
]


• Excellent work. I will certainly be using some aspects of your answer in other geo data related work I have been doing. – kickert Nov 28 '18 at 16:01

Here is a function for finding GeoPositions between 2 cities with certain step

city1 = Entity["City", {"Boston", "Massachusetts", "UnitedStates"}];
city2 = Entity["City", {"Rochester", "NewYork", "UnitedStates"}];
city3 = Entity["City", {"Chicago", "Illinois", "UnitedStates"}];

geopath[c1_, c2_, step_] := Module[{s = First@GeoDistance[c1, c2],
a = GeoDirection[c1, c2]},
GeoPath[NestList[GeoDestination[#,{1000 step,a}]&,c1,Floor[s/step]],"Rhumb"]]


Now if you want to find positions from city1 to city2 every 100 km, type

geopath[city1, city2, 100]
`

and you will get the positions

GeoPath[{Entity["City", {"Boston", "Massachusetts", "UnitedStates"}], GeoPosition[{42.5124, -72.2106}], GeoPosition[{42.6928, -73.4044}], GeoPosition[{42.8732, -74.6017}], GeoPosition[{43.0535, -75.8026}], GeoPosition[{43.2338, -77.0069}]}, "Rhumb"]