Antipodes don't match when calculated numerically

[boroli]

import scipy.optimize as opt
import pyproj
g = pyproj.Geod(ellps='WGS84')

latlona = 52.5252, 13.37173, 43.85

def antipods(dst):
    # positive branch
    antipods.a1 = g.fwd(latlona[1], latlona[0], latlona[2], dst*1e3)
    # negative branch
    antipods.a2 = g.fwd(latlona[1], latlona[0], latlona[2], -dst*1e3)
    (az12, az21, dist) = g.inv(antipods.a1[0], antipods.a1[1], antipods.a2[0], antipods.a2[1])
    antipods.dist = dist*0.001
    return dist*0.001

r = opt.minimize_scalar(antipods,bounds=[15000,25000],method="bounded",tol=1e-6)

print("Great Circle Distance to Antipodes:",r['x'])
print("Distance between calculated Antipodes",antipods.dist)

print("Antipode Positive",antipods.a1[1], antipods.a1[0])
print("Antipode Negative",antipods.a2[1], antipods.a2[0])
print("Antipode Analytic",-latlona[0],",",-180+latlona[1])

Problem description

When calculating the antipodes numerically, there's quite a huge gap between the two numerically calculated Antipodes. In the specific example above it's about 25 km. That's only about 0.12% of the great circle distance, but in absolute numbers I still found it quite a big distance.

Is this "just" a numerical accuracy lack? Or is it a systematic modelling error?

Expected Output

I'm not aware of the exact modelling, but I expected that the two antipodes match within reasonable accuracy, and here I mean something like much less than 1 km to the analytical one. The error varies slightly with the azimuth, but in all cases it is still visible.

My original goal was to plot a full 360° great circle, and there I figured out that after a full circle the great circle doesn't match the start point anymore, but has a decent offset.

Environment Information

pyproj info: pyproj: 3.3.0 PROJ: 8.2.1 data dir: /usr/share/proj user_data_dir: ~/.local/share/proj PROJ DATA (recommended version): 1.8 PROJ Database: 1.2 EPSG Database: v10.041 [2021-12-03] ESRI Database: ArcMap 12.8 [2021-05-06] IGNF Database: 3.1.0 [2019-05-24]

System: python: 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0] executable: /usr/bin/python3 machine: Linux-5.15.0-107-generic-x86_64-with-glibc2.35

Python deps: certifi: 2020.06.20 pip: 22.0.2 setuptools: 59.6.0 Cython: None

Installation method

Ubuntu 22.04 python3-package

[snowman2]

@cffk, apologies for bothering you. However, I think you are best suited to answer this question. Thank you!

[cffk]

With a few exceptions, geodesics are not closed curves on an ellipsoid. If you follow a geodesic a full circuit around the world, you arrive at a point that falls short of the starting point by a distance on the order of f a where f is the flattening of the ellipsoid and a is its equatorial radius.

Naturally, there's a similar offset following a geodesic halfway around the world in opposite directions.

The simplest way to compute the antipodal point is to add ±180° to the longitude and to change the sign of the latitude. But I understand that this probably isn't your goal.

[boroli]

Ok, got it. Thanks for the quick explanation. Then I can close that issue.

[jjimenezshaw]

@cffk does it mean that the shortest distance between a point and its antipode in an ellipsoid is not a geodetic?

[cffk]

No, a shortest path is always a geodesic. However the shortest (geodesic) path between a point and its antipode is via a pole (in the case of an oblate ellipsoid). (And not in some arbitrary direction as is the case for a sphere.)

[jjimenezshaw]

Then, if I understand correctly, the problem is the azimuth of both geodesics, one going through the north pole and another going through the south pole, are not the same. In the example it is using the same azimuth, 43.85, in both directions. Am I right? (Because both geodesics must have the same length)

[cffk]

"Up to a point, Lord Copper!" For an oblate ellipsoid, the two shortest geodesics connecting antipodal points have azimuths 0deg and 180deg. These geodesics over of the poles (obviously). There are two shortest geodesics connecting some points which are nearly antipodal (the points have to be on opposite circle of latitude) and, in this case, the two azimuths are supplementary angles. (In general, these geodesics do not go over the pole.)