Featured Map for 4 August 2023:
Near Antipodes

Most Featured Maps are about airlines, aircraft, or places, often trying to show off features of the site in the process. Today's Featured Map is different: it's about the "sausage making." Fasten your seatbelts because trigonometry and calculus are both on the itinerary. I won't be offended if you also strap into an ejection seat and punch out before the end.

Assuming Earth is a sphere, calculating the great circle path and distance between two points is a fairly simple calculation in spherical trigonometry. Unfortunately, a sphere is a crude approximation of Earth, so navigation usually uses a better model: an oblate spheroid—an ellipsoid which is like sphere that has been squashed at the axis of rotation, or that has a beer belly if you prefer. The most commonly used reference ellipsoid is WGS 84, the model used by the Global Positioning System (GPS). The WGS 84 ellipsoid has a polar radius of 6356.752 km and an equatorial radius of 6378.137 km.

The shortest path between two points on an ellipsoid is not really a great circle and is properly called a geodesic, though great circle is commonly used as an approximation. (Indeed, the very name of this site is guilty of that misnomer!)

Calculations of geodesics are rather more complicated than great circles. There are two different problems to solve: the direct problem is to compute the location of a point as a function of the distance and azimuth (direction) from a starting point while the inverse problem is to calculate the (geodesic) distance and azimuth between two points.

(When plotting a path, The Great Circle Mapper starts with the inverse method, then repeatedly uses the direct method with n steps from the origin in the calculated direction to generate the points along the path which are to be plotted. Well, unless you specify some other type of path such as a Bézier path.)

In 1975, Thaddeus Vincenty, a Polish-American geodesist, published a pair of iterative methods to solve both direct and inverse problems with high accuracy. Charles Karney subsequently published multiple refinements to Vincenty's work. (See FAQ: Great Circe Calculations for more detail.)

For the inverse problem, Vincenty's algorithm iterates until it gets close enough—or until it exceeds 20 iterations and gives up, which happens as the two points approach being antipodes. Antipodes are points which are directly opposite each other on the globe, such as BINP (North Pole) and NZSP (South Pole). (Er, not quite, because NZSP is defined as being at 89.999996ºS rather than 90ºS in hopes of avoiding calculation bugs. BINP is whimsy and depends on Rudolph's nose to avoid numerical navigational problems.)

For the cases where Vincenty's inverse algorithm fails to converge within a reasonable number of iterations, the Great Circle Mapper has long used a hack that offered a worst-case distance and path. This resulted in inflated distances and poor paths. Occasionally, it led to the nonsensical result that A-B was longer than A-C-B. Responses to a Featured Map earlier this year, Via Taipei: Shortest One-Stop Auckland to Vienna, finally triggered a quest to fix this problem.

Charles Karney published a way to handle these cases. Alas, he is a rocket scientist or rock star or otherwise rocks, whereas your favourite map-software developer is closer to a rock. Rocks can be stubborn in their own way, however, and 26.5 years after the challenge of "develop a web site to draw great circle maps" this became a new challenge to be stubbornly vanquished. Ideas of farming out the problem were considered but stubbornly rejected. (Stubborn is a recurring theme, yes?!)

If the inverse case fails to converge, a reasonable approach is to pick an approximate midpoint, apply Vincenty's algorithm, and then pick better "midpoints" until the result is close enough. An iterative algorithm wrapped around an iterative algorithm. As a computer scientist, that's bloody scary, especially when contemplating the "farthest" information included in all of the Airport Information pages which requires calculating the distance to many thousands of airports.

More specifically, choose a longitude midway between the endpoints and pick whichever pole is closest to the average latitude of the endpoints. The shortest path between the two points will cross the midway longitude at a latitude somewhere between the equator and that pole, so the problem is to find the latitude that results in the shortest distance from origin to midway point to destination.

One of the test cases ended up being Shanghai Pudong International Airport. Amongst the ten farthest airports from PVG, just two failed to converge: Concordia, Argentina and Salto, Uruguay. Here's a graph of distance for these two paths as a function of latitude, with the minima marked:

To find the minimum of one of these curves, start with the distance via the pole and via the equator. Pick a latitude midway between them and calculate that distance, then calculate (well, approximate) the derivative (calculus, oh no!) to decide whether a shorter distance is towards the pole or towards the equator. Now there's a 45-degree range of latitudes to consider. Repeat until the distance via the two bounds falls within a desired value: the Great Circle Mapper currrently uses 0.1 meters as the value of "good enough."

Testing revealed that the cost is minimal, partly because cases where the basic algorithm fails to converge are rare, and for those cases the two-step approach via a midpoint converges fairly quickly. Karney's approach is undoubtedly more efficient, but with the power of modern computers it's often the case that an inefficient approach is harmless. Careful benchmarking focuses optimization where it matters.

There can still be cases where the "shortest" path is longer than going via an intermediate point, but the difference will be small. For example, LMC-PGK is shorter via 83ºN 16º10'35"E but only by 4 millimeters. Yes, GCmap now offers distances in millimeters!

Today's Featured Map shows the three farthest airports from Shanghai, Shanghai's antipode, and an approximate range that bounds where Vincenty's algorithm fails to converge.

References and additional information:

• FAQ: Great Circe Calculations
• Vincenty's formulae – Wikipedia
  
  
• Featured Map: Via Taipei: Shortest One-Stop Auckland to Vienna (6 March 2023)