9
$\begingroup$

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?

$\endgroup$
7
$\begingroup$

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 = GeoGraphics`GeoEvaluate[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]]}
]

enter image description here

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

enter image description here

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

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.