Singularity at south pole not handled in geodetic calculation

[drsteve]

The conversion to geodetic (GDZ) coordinates blows up when evaluated at the south pole. This was reported by a SpacePy user who encountered the issue in our coordinates module, using IRBEM as the backend for the conversions.

Minimal example to reproduce issue:

To demonstrate using the IRBEM Python package (rather than SpacePy):

import IRBEM
import datetime as dt
import numpy as np

cc = IRBEM.Coords()
tt = dt.datetime.now()  # time is irrelevant for a GEO spherical to GDZ conversion
pos = np.array([[1.063871355909377, -90.0, 360.0], [1.063871355909377, -90.0, 0.0], [1.063871355909377, 90.0, 180.0]])
cc.coords_transform([tt]*3, pos, 'SPH', 'GDZ')

See spacepy/spacepy#577 for the initial report, a reproduction through both our Coords objects and our low-level wrapping, and additional discussion.

The returned geodetic (altitude, latitude, longitude), which should have units of [km, deg, deg], is:

array([[-1.31195190e+04, -9.00000000e+01, -1.40334186e-14],
       [-1.31195190e+04, -9.00000000e+01,  0.00000000e+00],
       [ 4.06014389e+02,  9.00000000e+01,  1.80000000e+02]])

The altitude returned for the first two positions (both at the southern pole) are clearly incorrect. The altitudes should all be the same, as noted in spacepy/spacepy#296 (scroll to the bottom, there's a ton of discussion), this can be handled correctly in the geodetic algorithm.

OS, IRBEM version, and compiler version information:

(One of my test setups listed here, but this is an algorithm/implementation issue and is platform-independent) Ubuntu 20.04 Latest commit: a4759c0 (as used in last release of SpacePy); 1ceaca6 (problem reproducible on latest commit) gfortran 9.3.0

Source of IRBEM library

Git repo & via SpacePy

Suggested resolution

The easiest fix is probably adding an explicit trap in the geo_gdz subroutine (starts on line 1529 of init_nouveau.f). I will note that the iterative algorithm may not be either the fastest or the most accurate these days, so it may also be worth considering a change to a more modern algorithm.

[AntoineBrunet]

Hi Steve, I'll look into this. What do you mean avout a more modern algorithm? From David Eberly's book on robust geometric computing, he seems to suggest that bissection is the more robust approach here. I guess we could implement the algorithm from Listing 2 in this document.

[drsteve]

There are exact, closed-form solutions for the conversion from cartesian to elliptical coordinates (i.e., Cartesian GEO to geodetic). See, e.g., https://doi.org/10.1109/7.303772 for a comparison of a few methods, including a number that are closed-form. Some require special handling of singularity points like poles or the equator, some don't. I used Heikkinen's exact method with one modification for numerical stability and a trap for the polar singularity, but Zhu's method is singularity-free.

The most accurate, efficient, and robust of the approximate algorithms, that I'm aware of, is Olson's 1996 method (not his older '88 method). This method should also be singularity-free.

[AntoineBrunet]

Thanks for the references. I'll look into it and submit a PR to fix this issue ASAP.