bug: LocalCartesian.reverse with lon = 1 radian

[ThomasGreenhill]

There appears to be a bug in the LocalCartesian.reverse method when lat is pi/2 rad and lon is 1 rad.
PyGeodesy 23.7.12
Python 3.9.7

import numpy as np
from pygeodesy.ltp import LocalCartesian
from pygeodesy.ellipsoidalVincenty import LatLon
from pygeodesy.ltpTuples import Ned

origin_np = np.array([np.pi / 2, 1.0, 0])

origin = LatLon(np.degrees(origin_np[0]), np.degrees(origin_np[1]), height=origin_np[2])
lc = LocalCartesian(origin)
point_relative_ned_position = np.zeros(3)
point_lla = lc.reverse(Ned(point_relative_ned_position[0], point_relative_ned_position[1], point_relative_ned_position[2]))

print(point_lla.toLatLon())

The result is:
(90.0, 0.0, 0.0, Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84))

Whereas it should be:
(90.0, 57.3, 0.0, Datum(name='WGS84', ellipsoid=Ellipsoids.WGS84, transform=Transforms.WGS84))

Without looking into the details of the implementation, I'd guess python is interpreting 1.0 as True somewhere.

[mrJean1]

A LocalCartesian.reverse call with x and y (and z all) zero always returns lon=0 because there is no single, unique solution (especially with the LocalCartesian positioned at one of the poles).

1. See the Note as the end of the documentation of the EcefKarney.reverse method (that is the Ecef class to convert between local and geodetic/geocentric coordinates inside LocalCartesian).

2. See the private function _norm3 inside EcefKarney.reverse on line 509 and lines 517 and 518 in pygeodesy.ecef.py.

The same applies to the Ltp class usable with other Ecef classes. Since returning lon=0 is arbitrary, perhaps it should be possible to specify the lon value to be returned in this particular case.

PS) Python considers 1.0 (any non-zero value) as True, but that is not related to the cause of this issue.

[mrJean1]

Proposed enhancement: a new keyword argument named lon00 for the methods LocalCartesian.reverse and EcefKarney.reverse but with different default values, LocalCartesian.reverse(…, lon00=self.lon0) and EcefKarney.reverse(…, lon00=0). Similarly for the Ltp and other Ecef classes.

The lon00 value is returned as longitude instead of the current, arbitrary zero when both x and y of the reverse call are zero. That will break backward compatibility for the former, but not for the latter.

The next release of PyGeodesy will include this enhancement for the LocalCartesian, Ltp and all Ecef classes.

[ThomasGreenhill]

Thank you for the response and the proposed enhancement.

I'm surprised by your statement regarding non-uniqueness of the result of reverse for 0, 0, 0 local cartesian displacement from an origin LLA (which is not at one of the poles). Shouldn't the result be identical to the origin regardless of the coordinates of the origin? The library does indeed behave as I expected (and as I described here) for points outside of the poles.

I know this is a bit of contrived/niche issue, since it's obvious that any point with latitude at the poles is going to be non-unique, and I'm grateful for the detailed explanation and the improvement!

[mrJean1]

Currently, converting a geocentric (ECEF) polar location (x = y = 0) results in geodetic latitude +/-90 and longitude 0, always and in all Ecef classes. But that 0 is not the unique longitude for a polar point, any other value would be correct.

The LocalCartesian and Ltp classes depend on an Ecef class to convert between local cartesian and geodetic coordinates. The local x, y and z coordinates are first multiplied with the rotation matrix and then converted to goedetic with the Ecef class (see line 161 and line 50 where pygeodesy is transcoded from).

As a result, a LocalCartesian or Ltp at one of the poles suffers from the same arbitrary 0 geodetic longitude for local point (0, 0). But that does not occur for a LocalCartesian or Ltp located elsewhere because the rotated x and y value are never both 0.

As you pointed out, the 0 longitude does not makes sense, agreed. As you also suggested, the longitude of the LocalCartesian's or Ltp's geodetic origin is expected and that should be returned.

Breaking backward compatibility is not of an concern, since retaining the old, 0 value would be wrong.

One other point. Using the LocalCartesian or Ltp origin’s longitude for the “polar case” as default for a new keyword argument like lon00 is the straightforward way to resolve this issue.

An other option might be the following: for a “polar” LocalCartesian or Ltp, adjust the rotation matrix such that it produces a geocentric x, y and z set for reverse locations such that final geodetic (lat- and) longitude are the origin of the “polar” LocalCartesian or Ltp. Not feasible.

[mrJean1]

PyGeodesy 23.08.08 includes the new lon00 keyword argument for all Ecef, LocalCartesian and Ltp classes (and a few other things ;).

Please verify that his resolves the original issue and if so, please close this issue. Again, thank you very much for reporting the issue and for the example.