eci2ecef: scalar output for scalar input with Numpy backend

add test for Numpy and AstroPy backends

fixes #98

Author: scivision

src/pymap3d/aer.py

@@ -199,12 +199,9 @@ def eci2aer(x, y, z, lat0, lon0, h0, t: datetime, *, deg: bool = True) -> tuple:
          slant range [meters]
     """
 
-    try:
-        xecef, yecef, zecef = eci2ecef(x, y, z, t)
-    except NameError:
-        raise ImportError("pip install numpy")
+    xe, ye, ze = eci2ecef(x, y, z, t)
 
-    return ecef2aer(xecef, yecef, zecef, lat0, lon0, h0, deg=deg)
+    return ecef2aer(xe, ye, ze, lat0, lon0, h0, deg=deg)
 
 
 def aer2eci(

src/pymap3d/ecef.py

@@ -75,9 +75,7 @@ def geodetic2ecef(
         lon = radians(lon)
 
     # radius of curvature of the prime vertical section
-    N = ell.semimajor_axis**2 / hypot(
-        ell.semimajor_axis * cos(lat), ell.semiminor_axis * sin(lat)
-    )
+    N = ell.semimajor_axis**2 / hypot(ell.semimajor_axis * cos(lat), ell.semiminor_axis * sin(lat))
     # Compute cartesian (geocentric) coordinates given (curvilinear) geodetic coordinates.
     x = (N + alt) * cos(lat) * cos(lon)
     y = (N + alt) * cos(lat) * sin(lon)
@@ -202,9 +200,7 @@ def ecef2geodetic(
 
     # inside ellipsoid?
     inside = (
-        x**2 / ell.semimajor_axis**2
-        + y**2 / ell.semimajor_axis**2
-        + z**2 / ell.semiminor_axis**2
+        x**2 / ell.semimajor_axis**2 + y**2 / ell.semimajor_axis**2 + z**2 / ell.semiminor_axis**2
         < 1
     )
 

src/pymap3d/eci.py

@@ -4,7 +4,7 @@
 
 from datetime import datetime
 
-import numpy as np
+import numpy
 
 try:
     import astropy.units as u
@@ -45,36 +45,97 @@ def eci2ecef(x, y, z, time: datetime) -> tuple:
     """
 
     try:
-        gcrs = GCRS(CartesianRepresentation(x * u.m, y * u.m, z * u.m), obstime=time)
-        itrs = gcrs.transform_to(ITRS(obstime=time))
-
-        x_ecef = itrs.x.value
-        y_ecef = itrs.y.value
-        z_ecef = itrs.z.value
+        return eci2ecef_astropy(x, y, z, time)
     except NameError:
-        x = np.atleast_1d(x)
-        y = np.atleast_1d(y)
-        z = np.atleast_1d(z)
-        gst = np.atleast_1d(greenwichsrt(juliandate(time)))
-        assert (
-            x.shape == y.shape == z.shape
-        ), f"shape mismatch: x: ${x.shape}  y: {y.shape}  z: {z.shape}"
-        if gst.size == 1 and x.size != 1:
-            gst = np.broadcast_to(gst, x.shape[0])
-        assert x.size == gst.size, f"shape mismatch: x: {x.shape}  gst: {gst.shape}"
-
-        eci = np.column_stack((x.ravel(), y.ravel(), z.ravel()))
-        ecef = np.empty((x.size, 3))
-        for i in range(eci.shape[0]):
-            ecef[i, :] = R3(gst[i]) @ eci[i, :].T
-
-        x_ecef = ecef[:, 0].reshape(x.shape)
-        y_ecef = ecef[:, 1].reshape(y.shape)
-        z_ecef = ecef[:, 2].reshape(z.shape)
+        return eci2ecef_numpy(x, y, z, time)
+
+
+def eci2ecef_astropy(x, y, z, t: datetime) -> tuple:
+    """
+    eci2ecef using Astropy
+
+    Parameters
+    ----------
+    x : float
+        ECI x-location [meters]
+    y : float
+        ECI y-location [meters]
+    z : float
+        ECI z-location [meters]
+    t : datetime.datetime
+        time of obsevation (UTC)
+
+    Results
+    -------
+    x_ecef : float
+        x ECEF coordinate
+    y_ecef : float
+        y ECEF coordinate
+    z_ecef : float
+        z ECEF coordinate
+    """
+
+    gcrs = GCRS(CartesianRepresentation(x * u.m, y * u.m, z * u.m), obstime=t)
+    itrs = gcrs.transform_to(ITRS(obstime=t))
+
+    x_ecef = itrs.x.value
+    y_ecef = itrs.y.value
+    z_ecef = itrs.z.value
 
     return x_ecef, y_ecef, z_ecef
 
 
+def eci2ecef_numpy(x, y, z, t: datetime) -> tuple:
+    """
+    eci2ecef using Numpy
+
+    less accurate than astropy, but may be good enough.
+
+    Parameters
+    ----------
+    x : float
+        ECI x-location [meters]
+    y : float
+        ECI y-location [meters]
+    z : float
+        ECI z-location [meters]
+    t : datetime.datetime
+        time of obsevation (UTC)
+
+    Results
+    -------
+    x_ecef : float
+        x ECEF coordinate
+    y_ecef : float
+        y ECEF coordinate
+    z_ecef : float
+        z ECEF coordinate
+    """
+
+    x = numpy.atleast_1d(x)
+    y = numpy.atleast_1d(y)
+    z = numpy.atleast_1d(z)
+    gst = numpy.atleast_1d(greenwichsrt(juliandate(t)))
+    assert (
+        x.shape == y.shape == z.shape
+    ), f"shape mismatch: x: ${x.shape}  y: {y.shape}  z: {z.shape}"
+
+    if gst.size == 1 and x.size != 1:
+        gst = numpy.broadcast_to(gst, x.shape[0])
+    assert x.size == gst.size, f"shape mismatch: x: {x.shape}  gst: {gst.shape}"
+
+    eci = numpy.column_stack((x.ravel(), y.ravel(), z.ravel()))
+    ecef = numpy.empty((x.size, 3))
+    for i in range(eci.shape[0]):
+        ecef[i, :] = R3(gst[i]) @ eci[i, :].T
+
+    x_ecef = ecef[:, 0].reshape(x.shape)
+    y_ecef = ecef[:, 1].reshape(y.shape)
+    z_ecef = ecef[:, 2].reshape(z.shape)
+
+    return x_ecef.squeeze()[()], y_ecef.squeeze()[()], z_ecef.squeeze()[()]
+
+
 def ecef2eci(x, y, z, time: datetime) -> tuple:
     """
     Point => Point   ECEF => ECI
@@ -112,15 +173,15 @@ def ecef2eci(x, y, z, time: datetime) -> tuple:
         y_eci = eci.y.value
         z_eci = eci.z.value
     except NameError:
-        x = np.atleast_1d(x)
-        y = np.atleast_1d(y)
-        z = np.atleast_1d(z)
-        gst = np.atleast_1d(greenwichsrt(juliandate(time)))
+        x = numpy.atleast_1d(x)
+        y = numpy.atleast_1d(y)
+        z = numpy.atleast_1d(z)
+        gst = numpy.atleast_1d(greenwichsrt(juliandate(time)))
         assert x.shape == y.shape == z.shape
         assert x.size == gst.size
 
-        ecef = np.column_stack((x.ravel(), y.ravel(), z.ravel()))
-        eci = np.empty((x.size, 3))
+        ecef = numpy.column_stack((x.ravel(), y.ravel(), z.ravel()))
+        eci = numpy.empty((x.size, 3))
         for i in range(x.size):
             eci[i, :] = R3(gst[i]).T @ ecef[i, :]
 
@@ -133,4 +194,6 @@ def ecef2eci(x, y, z, time: datetime) -> tuple:
 
 def R3(x: float):
     """Rotation matrix for ECI"""
-    return np.array([[np.cos(x), np.sin(x), 0], [-np.sin(x), np.cos(x), 0], [0, 0, 1]])
+    return numpy.array(
+        [[numpy.cos(x), numpy.sin(x), 0], [-numpy.sin(x), numpy.cos(x), 0], [0, 0, 1]]
+    )

src/pymap3d/enu.py

@@ -1,4 +1,5 @@
 """ transforms involving ENU East North Up """
+
 from __future__ import annotations
 
 from math import tau

src/pymap3d/spherical.py

@@ -3,6 +3,7 @@
 longitude, height) and geocentric spherical (spherical latitude, longitude,
 radius).
 """
+
 from __future__ import annotations
 
 from .ellipsoid import Ellipsoid

src/pymap3d/tests/test_eci.py

@@ -1,4 +1,4 @@
-from datetime import datetime
+import datetime
 
 import pymap3d as pm
 import pytest
@@ -9,25 +9,51 @@
 except ImportError:
     astropy = None
 
+ECI = (-2981784, 5207055, 3161595)
+ECEF = [-5762640, -1682738, 3156028]
+UTC = datetime.datetime(2019, 1, 4, 12, tzinfo=datetime.timezone.utc)
+
 
 def test_eci2ecef():
     pytest.importorskip("numpy")
     # this example from Matlab eci2ecef docs
-    eci = [-2981784, 5207055, 3161595]
-    utc = datetime(2019, 1, 4, 12)
-    ecef = pm.eci2ecef(*eci, utc)
+    ecef = pm.eci2ecef(*ECI, UTC)
+
+    assert isinstance(ecef[0], float)
+    assert isinstance(ecef[1], float)
+    assert isinstance(ecef[2], float)
+
+
+def test_eci2ecef_numpy():
+    pytest.importorskip("numpy")
+
+    ecef = pm.eci2ecef(*ECI, UTC)
 
-    rel = 0.025 if astropy is None else 0.0001
+    rel = 0.025
 
-    assert ecef == approx([-5.7627e6, -1.6827e6, 3.1560e6], rel=rel)
+    assert ecef == approx(ECEF, rel=rel)
+    assert isinstance(ecef[0], float)
+    assert isinstance(ecef[1], float)
+    assert isinstance(ecef[2], float)
+
+
+def test_eci2ecef_astropy():
+    pytest.importorskip("astropy")
+
+    ecef = pm.eci2ecef(*ECI, UTC)
+
+    rel = 0.0001
+
+    assert ecef == approx(ECEF, rel=rel)
+    assert isinstance(ecef[0], float)
+    assert isinstance(ecef[1], float)
+    assert isinstance(ecef[2], float)
 
 
 def test_ecef2eci():
     pytest.importorskip("numpy")
     # this example from Matlab ecef2eci docs
-    ecef = [-5762640, -1682738, 3156028]
-    utc = datetime(2019, 1, 4, 12)
-    eci = pm.ecef2eci(*ecef, utc)
+    eci = pm.ecef2eci(*ECEF, UTC)
 
     rel = 0.01 if astropy is None else 0.0001
 
@@ -37,9 +63,7 @@ def test_ecef2eci():
 def test_eci2geodetic():
     pytest.importorskip("numpy")
 
-    eci = [-2981784, 5207055, 3161595]
-    utc = datetime(2019, 1, 4, 12)
-    lla = pm.eci2geodetic(*eci, utc)
+    lla = pm.eci2geodetic(*ECI, UTC)
 
     rel = 0.01 if astropy is None else 0.0001
 
@@ -50,8 +74,8 @@ def test_geodetic2eci():
     pytest.importorskip("numpy")
 
     lla = [27.880801, -163.722058, 408850.646]
-    utc = datetime(2019, 1, 4, 12)
-    eci = pm.geodetic2eci(*lla, utc)
+
+    eci = pm.geodetic2eci(*lla, UTC)
 
     rel = 0.01 if astropy is None else 0.0001
 
@@ -61,7 +85,7 @@ def test_geodetic2eci():
 def test_eci_aer():
     # test coords from Matlab eci2aer
     pytest.importorskip("numpy")
-    t = datetime(2022, 1, 2, 3, 4, 5)
+    t = datetime.datetime(2022, 1, 2, 3, 4, 5, tzinfo=datetime.timezone.utc)
 
     eci = [4500000, -45000000, 3000000]
     lla = [28, -80, 100]