# Finding the closest interpreted point in GeoDistance

I'd like to know how far away I am standing from the path of totality for the 2024 solar eclipse.

loc = GeoPosition[{43.215580, -77.937510}];
eclipse =
SolarEclipse[{DateObject[{2024, 4, 8, 0, 0}]},
"TotalPhaseCenterLine", EclipseType -> "Total"];
bport = Entity["City", {"Brockport", "NewYork", "UnitedStates"}];
GeoGraphics[{
GeoMarker[loc],
eclipse},
GeoCenter -> bport,
GeoRange -> Quantity[2, "Miles"]]


GeoDistance[loc, First@eclipse["Values"]]
(* 0.735... miles *)


That seems about right. Now, I would like to know the GeoPosition of that point. The problem is, the solar eclipse data is (in this case) ~ 195 points, and the closest data point is actually 19 miles away.

minpt = GeoPosition@
With[{pts = GeoDistance[loc, #] & /@ eclipse["Values"][[1, 1, 1]]},
(*Position[pts,Min@pts]*)
eclipse["Values"][[1, 1, 1, Position[pts, Min@pts][[1, 1]]]]]
GeoGraphics[GeoMarker[{loc, minpt}], PlotLabel -> GeoDistance[loc, minpt]]


It is nice that GeoDistance is interpolating the eclipse line. Do I have to do the same in order to figure out the point on the solar eclipse path that is closest to my current location? I thought perhaps one could use GeoDirection along with GeoDestination however the former does not accept the eclipse path as an argument.

We can binary search for the nearest point.

First we find our initial search bracket:

dists0 = Normal[GeoDistance[loc, eclipse["Values"][[1, 1]]]];

i1 = Ordering[QuantityMagnitude[dists0, "Miles"], {1}][[1]];
i2 = Which[
i1 === 1, 2,
i1 === Length[dists0], i1-1,
dists0[[i1 + 1]] < dists0[[i1 - 1]], i1 + 1,
True, i1 - 1
];

{p1, p2} = Thread[eclipse["Values"][[1, 1]]][[{i1, i2}]];


Next define the binary search function:

BinarySearch[f_, l_, h_, distf_, meanf_, tol_] :=
Block[{lo = l, hi = h, mid},
While[distf[lo, hi] > tol,
mid = meanf[lo, hi];
If[TrueQ[f[lo, hi, mid]], hi = mid, lo = mid]
];
lo
]


Then our custom functions:

bfunc = GeoDistance[loc, #1] < GeoDistance[loc, #2]&;

mfunc = GeoPosition[Mean[{#1[[1]], #2[[1]]}]]&;

dfunc = GeoDistance;


Then run:

pmin = BinarySearch[bfunc, p1, p2, dfunc, mfunc, Quantity[10^-9.0, "Meters"]]

GeoPosition[{43.205, -77.930}]


Visualize:

GeoGraphics[
{GeoMarker[loc], eclipse, GeoMarker[pmin, "Color" -> Blue]},
GeoCenter -> bport,
GeoRange -> Quantity[2, "Miles"]
]


If the path of totality was returned as a GeoPath instead of Line we would need to replace mfunc to ensure the midpoint remains on the path. Something like this:

mfunc2 = GeoDestination[#1, GeoDisplacement[{0.5GeoDistance[##], GeoDirection[##]}]]&;


Now, note that pmin is farther away than what GeoDistance gives us for the path:

GeoDistance[loc, eclipse["Values"][[1]]]

0.735... miles

GeoDistance[loc, pmin]

0.785... miles


However the distance reported by pmin looks correct:

GeoGraphics[{GeoMarker[loc], eclipse,
{Red, GeoCircle[loc, GeoDistance[loc, eclipse["Values"][[1]]]]},
{Blue, GeoCircle[loc, GeoDistance[loc, pmin]]}},
GeoCenter -> bport,
GeoRange -> Quantity[2, "Miles"]
]


What's happening is GeoDistance is treating the path as a GeoPath, whereas GeoGraphics is staying true to Line. If we work with geodesics across the board we can fix the discrepancy:

pmin2 = BinarySearch[bfunc, p1, p2, dfunc, mfunc2, Quantity[10^-9.0, "Meters"]]

GeoPosition[{43.2061, -77.9308}]

GeoDistance[loc, pmin2]

0.735... miles


Replacing Line with GeoPath in the graphic:

GeoGraphics[{GeoMarker[loc], GeoPath @@ eclipse["Values"][[1]],
{Blue, GeoCircle[loc, GeoDistance[loc, pmin2]]}},
GeoCenter -> bport,
GeoRange -> Quantity[2, "Miles"]
]


Lastly, a completely different approach is to numerically solve for where the relevant GeoCircle and the path of totality intersect. The advantage here is we don't need to scrutinize the 'path type' of our curve.

path = eclipse["Values"][[1]];

d = GeoDistance[loc, path];

distα[α_?NumericQ] :=
QuantityMagnitude[
GeoDistance[GeoDestination[loc, GeoDisplacement[{d, α}]], path],
"Miles"
]

αmin = α /. FindRoot[distα[α], {α, 150, 0, 360}];
pmin3 = GeoDestination[loc, GeoDisplacement[{d, αmin}]]

GeoPosition[{43.2061, -77.9307}]

GeoDistance[loc, pmin3]

0.735... miles

• Excellent - thanks much for explaining the GeoPath/Line discrepancy. Commented Jan 24 at 19:56

A possible solution is to use Projection.

near = Nearest[GeoPosition /@ eclipse["Values"][[1, 1, 1]], loc, 2] /.
GeoPosition[pt_] :> pt
ptonline =
Projection[loc["LatitudeLongitude"] - First@near,
Normalize@Flatten@Differences@near] + First@near // GeoPosition
GeoGraphics[{GeoMarker[loc], eclipse, GeoMarker[ptonline]},
GeoCenter -> bport, GeoRange -> Quantity[1, "Miles"]]
GeoDistance[loc, ptonline]


Although note that GeoDistance[loc,ptonline] gives a slightly longer distance (0.80 mi vs 0.735) which might be due to working with vectors on a sphere. That math is well beyond me.