ecef2geodetic: test/fix float32

Co-authored-by: unjambonakap <[email protected]>

Author: scivision

README.md

@@ -76,6 +76,9 @@ where tuple `lla` is comprised of scalar or N-D arrays `(lat,lon,alt)`.
 
 Example scripts are in the [examples](./Examples) directory.
 
+Native Python float is typically [64 bit](https://docs.python.org/3/library/stdtypes.html#typesnumeric).
+Numpy can select real precision bits: 32, 64, 128, etc.
+
 ### Functions
 
 Popular mapping toolbox functions ported to Python include the

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pymap3d
-version = 2.7.2
+version = 2.7.3
 author = Michael Hirsch, Ph.D.
 author_email = [email protected]
 description = pure Python (no prereqs) coordinate conversions, following convention of several popular Matlab routines.

src/pymap3d/ecef.py

@@ -3,7 +3,19 @@
 import typing
 
 try:
-    from numpy import radians, sin, cos, tan, arctan as atan, hypot, degrees, arctan2 as atan2, sqrt
+    from numpy import (
+        radians,
+        sin,
+        cos,
+        tan,
+        arctan as atan,
+        hypot,
+        degrees,
+        arctan2 as atan2,
+        sqrt,
+        finfo,
+        where,
+    )
     from .eci import eci2ecef, ecef2eci
 except ImportError:
     from math import radians, sin, cos, tan, atan, hypot, degrees, atan2, sqrt  # type: ignore
@@ -146,18 +158,36 @@ def ecef2geodetic(
             Beta = -pi / 2
 
     # eqn. 13
-    eps = ((ell.semiminor_axis * u - ell.semimajor_axis * huE + E ** 2) * sin(Beta)) / (
+    dBeta = ((ell.semiminor_axis * u - ell.semimajor_axis * huE + E ** 2) * sin(Beta)) / (
         ell.semimajor_axis * huE * 1 / cos(Beta) - E ** 2 * cos(Beta)
     )
 
-    Beta += eps
+    Beta += dBeta
+
+    # eqn. 4c
     # %% final output
     lat = atan(ell.semimajor_axis / ell.semiminor_axis * tan(Beta))
 
+    try:
+        # patch latitude for float32 precision loss
+        lim_pi2 = pi / 2 - finfo(dBeta.dtype).eps
+        lat = where(Beta >= lim_pi2, pi / 2, lat)
+        lat = where(Beta <= -lim_pi2, -pi / 2, lat)
+    except NameError:
+        pass
+
     lon = atan2(y, x)
 
     # eqn. 7
-    alt = hypot(z - ell.semiminor_axis * sin(Beta), Q - ell.semimajor_axis * cos(Beta))
+    cosBeta = cos(Beta)
+    try:
+        # patch altitude for float32 precision loss
+        cosBeta = where(Beta >= lim_pi2, 0, cosBeta)
+        cosBeta = where(Beta <= -lim_pi2, 0, cosBeta)
+    except NameError:
+        pass
+
+    alt = hypot(z - ell.semiminor_axis * sin(Beta), Q - ell.semimajor_axis * cosBeta)
 
     # inside ellipsoid?
     inside = (

src/pymap3d/tests/test_geodetic.py

@@ -13,6 +13,28 @@
 A = ELL.semimajor_axis
 B = ELL.semiminor_axis
 
+xyzlla = [
+    ((A - 1, 0, 0), (0, 0, -1)),
+    ((0, A - 1, 0), (0, 90, -1)),
+    ((0, 0, B - 1), (90, 0, -1)),
+    ((0, 0, B - 1), (89.999999, 0, -1)),
+    ((0, 0, B - 1), (89.99999, 0, -1)),
+    ((0, 0, -B + 1), (-90, 0, -1)),
+    ((0, 0, -B + 1), (-89.999999, 0, -1)),
+    ((0, 0, -B + 1), (-89.99999, 0, -1)),
+    ((-A + 1, 0, 0), (0, 180, -1)),
+]
+
+llaxyz = [
+    ((0, 0, -1), (A - 1, 0, 0)),
+    ((0, 90, -1), (0, A - 1, 0)),
+    ((0, -90, -1), (0, -A + 1, 0)),
+    ((90, 0, -1), (0, 0, B - 1)),
+    ((90, 15, -1), (0, 0, B - 1)),
+    ((-90, 0, -1), (0, 0, -B + 1)),
+]
+
+
 atol_dist = 1e-6  # 1 micrometer
 
 
@@ -94,7 +116,7 @@ def test_xarray():
     xyz = pm.geodetic2ecef(*xr_lla)
 
     assert xyz == approx(xyz0)
-    # %%
+
     xr_xyz = xarray.DataArray(list(xyz0))
 
     lla = pm.ecef2geodetic(*xr_xyz)
@@ -136,31 +158,12 @@ def test_ecef():
     assert pm.ecef2geodetic((A - 1) / sqrt(2), (A - 1) / sqrt(2), 0) == approx([0, 45, -1])
 
 
-@pytest.mark.parametrize(
-    "lla, xyz",
-    [
-        ((0, 0, -1), (A - 1, 0, 0)),
-        ((0, 90, -1), (0, A - 1, 0)),
-        ((0, -90, -1), (0, -A + 1, 0)),
-        ((90, 0, -1), (0, 0, B - 1)),
-        ((90, 15, -1), (0, 0, B - 1)),
-        ((-90, 0, -1), (0, 0, -B + 1)),
-    ],
-)
+@pytest.mark.parametrize("lla, xyz", llaxyz)
 def test_geodetic2ecef(lla, xyz):
     assert pm.geodetic2ecef(*lla) == approx(xyz, abs=atol_dist)
 
 
-@pytest.mark.parametrize(
-    "xyz, lla",
-    [
-        ((A - 1, 0, 0), (0, 0, -1)),
-        ((0, A - 1, 0), (0, 90, -1)),
-        ((0, 0, B - 1), (90, 0, -1)),
-        ((0, 0, -B + 1), (-90, 0, -1)),
-        ((-A + 1, 0, 0), (0, 180, -1)),
-    ],
-)
+@pytest.mark.parametrize("xyz, lla", xyzlla)
 def test_ecef2geodetic(xyz, lla):
     lat, lon, alt = pm.ecef2geodetic(*xyz)
     assert lat == approx(lla[0])
@@ -221,3 +224,19 @@ def test_somenan():
 
     lat, lon, alt = pm.ecef2geodetic(xyz[:, 0], xyz[:, 1], xyz[:, 2])
     assert (lat[0], lon[0], alt[0]) == approx(lla0)
+
+
+@pytest.mark.parametrize("xyz, lla", xyzlla)
+def test_numpy_ecef2geodetic(xyz, lla):
+    np = pytest.importorskip("numpy")
+    lat, lon, alt = pm.ecef2geodetic(
+        *np.array(
+            [
+                [xyz],
+            ],
+            dtype=np.float32,
+        ).T
+    )
+    assert lat[0] == approx(lla[0])
+    assert lon[0] == approx(lla[1])
+    assert alt[0] == approx(lla[2])