5
$\begingroup$

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

enter image description here

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

enter image description here

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.

$\endgroup$

2 Answers 2

3
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ Excellent - thanks much for explaining the GeoPath/Line discrepancy. $\endgroup$ Jan 24 at 19:56
0
$\begingroup$

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]

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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