add geodetic2eci

Author: scivision

README.md

@@ -88,7 +88,7 @@ converted to the desired coordinate system:
 
     aer2ecef  aer2enu  aer2geodetic  aer2ned
     ecef2aer  ecef2enu  ecef2enuv  ecef2geodetic  ecef2ned  ecef2nedv
-    ecef2eci  eci2ecef  eci2aer  aer2eci
+    ecef2eci  eci2ecef eci2aer aer2eci geodetic2eci eci2geodetic
     enu2aer  enu2ecef   enu2geodetic
     geodetic2aer  geodetic2ecef  geodetic2enu  geodetic2ned
     ned2aer  ned2ecef   ned2geodetic
@@ -135,7 +135,8 @@ As noted above, use list comprehension if you need vector data without Numpy.
 
 As compared to [PyProj](https://github.com/jswhit/pyproj):
 
-* PyMap3D does not require anything beyond pure Python -- not even Numpy is required except for ECI (let us know if this is an issue).
+* PyMap3D does not require anything beyond pure Python for most transforms
+* Astronomical conversions are done using (optional) AstroPy for established accuracy
 * PyMap3D API is similar to Matlab Mapping Toolbox, while PyProj's interface is quite distinct
 * PyMap3D intrinsically handles local coordinate systems such as ENU,
   while PyProj ENU requires some [additional effort](https://github.com/jswhit/pyproj/issues/105).

src/pymap3d/__init__.py

@@ -35,7 +35,7 @@
 from .ellipsoid import Ellipsoid
 from .enu import enu2geodetic, geodetic2enu, aer2enu, enu2aer
 from .ned import ned2ecef, ned2geodetic, geodetic2ned, ecef2nedv, ned2aer, aer2ned, ecef2ned
-from .ecef import geodetic2ecef, ecef2geodetic, eci2geodetic, ecef2enuv, enu2ecef, ecef2enu, enu2uvw, uvw2enu
+from .ecef import geodetic2ecef, ecef2geodetic, eci2geodetic, geodetic2eci, ecef2enuv, enu2ecef, ecef2enu, enu2uvw, uvw2enu
 from .sidereal import datetime2sidereal
 from .latitude import (
     geod2geoc,

src/pymap3d/ecef.py

@@ -1,7 +1,7 @@
 """ Transforms involving ECEF: earth-centered, earth-fixed frame """
 try:
     from numpy import radians, sin, cos, tan, arctan as atan, hypot, degrees, arctan2 as atan2, sqrt, pi, vectorize
-    from .eci import eci2ecef
+    from .eci import eci2ecef, ecef2eci
 except ImportError:
     from math import radians, sin, cos, tan, atan, hypot, degrees, atan2, sqrt, pi
 
@@ -15,7 +15,17 @@
 # py < 3.6 compatible
 tau = 2 * pi
 
-__all__ = ["geodetic2ecef", "ecef2geodetic", "ecef2enuv", "ecef2enu", "enu2uvw", "uvw2enu", "eci2geodetic", "enu2ecef"]
+__all__ = [
+    "geodetic2ecef",
+    "ecef2geodetic",
+    "ecef2enuv",
+    "ecef2enu",
+    "enu2uvw",
+    "uvw2enu",
+    "eci2geodetic",
+    "geodetic2eci",
+    "enu2ecef",
+]
 
 if typing.TYPE_CHECKING:
     from numpy import ndarray
@@ -320,10 +330,14 @@ def uvw2enu(
     return East, North, Up
 
 
-def eci2geodetic(x: "ndarray", y: "ndarray", z: "ndarray", t: datetime) -> typing.Tuple["ndarray", "ndarray", "ndarray"]:
+def eci2geodetic(
+    x: "ndarray", y: "ndarray", z: "ndarray", t: datetime, ell: Ellipsoid = None, deg: bool = True
+) -> typing.Tuple["ndarray", "ndarray", "ndarray"]:
     """
     convert Earth Centered Internal ECI to geodetic coordinates
 
+    J2000 time
+
     Parameters
     ----------
     x : "ndarray"
@@ -333,7 +347,11 @@ def eci2geodetic(x: "ndarray", y: "ndarray", z: "ndarray", t: datetime) -> typin
     z : "ndarray"
         ECI z-location [meters]
     t : datetime.datetime, "ndarray"
-          length N vector of datetime OR greenwich sidereal time angle [radians].
+        UTC time
+    ell : Ellipsoid (optional)
+        planet ellipsoid model
+    deg : bool (optional)
+        if True, degrees. if False, radians
 
     Results
     -------
@@ -344,21 +362,56 @@ def eci2geodetic(x: "ndarray", y: "ndarray", z: "ndarray", t: datetime) -> typin
     alt : "ndarray"
           altitude above ellipsoid  (meters)
 
-    Notes
-    -----
-
-    Conversion is idealized: doesn't consider nutations, perterbations,
-    etc. like the IAU-76/FK5 or IAU-2000/2006 model-based conversions
-    from ECI to ECEF
-
     eci2geodetic() a.k.a. eci2lla()
     """
     if eci2ecef is None:
         raise ImportError("pip install astropy")
 
     xecef, yecef, zecef = eci2ecef(x, y, z, t)
 
-    return ecef2geodetic(xecef, yecef, zecef)
+    return ecef2geodetic(xecef, yecef, zecef, ell, deg)
+
+
+def geodetic2eci(
+    lat: "ndarray", lon: "ndarray", alt: "ndarray", t: datetime, ell: Ellipsoid = None, deg: bool = True
+) -> typing.Tuple["ndarray", "ndarray", "ndarray"]:
+    """
+    convert geodetic coordinates to Earth Centered Internal ECI
+
+    J2000 frame
+
+    Parameters
+    ----------
+    lat : "ndarray"
+        geodetic latitude
+    lon : "ndarray"
+        geodetic longitude
+    alt : "ndarray"
+        altitude above ellipsoid  (meters)
+    t : datetime.datetime, "ndarray"
+        UTC time
+    ell : Ellipsoid (optional)
+        planet ellipsoid model
+    deg : bool (optional)
+        if True, degrees. if False, radians
+
+    Results
+    -------
+    x : "ndarray"
+        ECI x-location [meters]
+    y : "ndarray"
+        ECI y-location [meters]
+    z : "ndarray"
+        ECI z-location [meters]
+
+    geodetic2eci() a.k.a lla2eci()
+    """
+    if ecef2eci is None:
+        raise ImportError("pip install astropy")
+
+    x, y, z = geodetic2ecef(lat, lon, alt, ell, deg)
+
+    return ecef2eci(x, y, z, t)
 
 
 def enu2ecef(
@@ -384,9 +437,9 @@ def enu2ecef(
     u1 : "ndarray"
         target up ENU coordinate (meters)
     lat0 : "ndarray"
-           Observer geodetic latitude
+        Observer geodetic latitude
     lon0 : "ndarray"
-           Observer geodetic longitude
+        Observer geodetic longitude
     h0 : "ndarray"
          observer altitude above geodetic ellipsoid (meters)
     ell : Ellipsoid, optional

src/pymap3d/eci.py

@@ -35,12 +35,12 @@ def eci2ecef(x: "ndarray", y: "ndarray", z: "ndarray", time: datetime) -> typing
 
     Results
     -------
-    x : "ndarray"
-        target x ECEF coordinate
-    y : "ndarray"
-        target y ECEF coordinate
-    z : "ndarray"
-        target z ECEF coordinate
+    x_ecef : "ndarray"
+        x ECEF coordinate
+    y_ecef : "ndarray"
+        y ECEF coordinate
+    z_ecef : "ndarray"
+        z ECEF coordinate
     """
 
     if Time is None:
@@ -72,12 +72,12 @@ def ecef2eci(x: "ndarray", y: "ndarray", z: "ndarray", time: datetime) -> typing
 
     Results
     -------
-    x : "ndarray"
-        target x ECI coordinate
-    y : "ndarray"
-        target y ECI coordinate
-    z : "ndarray"
-        target z ECI coordinate
+    x_eci : "ndarray"
+        x ECI coordinate
+    y_eci : "ndarray"
+        y ECI coordinate
+    z_eci : "ndarray"
+        z ECI coordinate
     """
 
     if Time is None:

tests/test_astropy.py

@@ -3,16 +3,14 @@
 from pytest import approx
 from math import radians
 from datetime import datetime
-import pymap3d as pm
+
 import pymap3d.sidereal as pmd
 import pymap3d.haversine as pmh
 
 lon = -148
 t0 = datetime(2014, 4, 6, 8)
-lla0 = (42, -82, 200)
 sra = 2.90658
 ha = 45.482789587392013
-eci0 = (-3.977913815668146e6, -2.582332196263046e6, 4.250818828152067e6)
 
 
 @pytest.mark.parametrize("time", [t0, [t0]], ids=("scalar", "list"))
@@ -44,25 +42,5 @@ def test_anglesep_meeus():
     assert pmh.anglesep_meeus(35, 23, 84, 20) == approx(ha)
 
 
-def test_eci_geodetic():
-    pytest.importorskip("astropy")
-    t = "2013-01-15T12:00:05"
-    lla = pm.eci2geodetic(*eci0, t)
-    assert lla == approx(lla0, rel=0.001)
-
-
-def test_eci_aer():
-    pytest.importorskip("astropy")
-    t = "2013-01-15T12:00:05"
-
-    aer1 = pm.eci2aer(*eci0, 42, -100, 0, t)
-    assert aer1 == approx([83.9511, -6.66787, 1485123.89], rel=0.001)
-
-    assert pm.aer2eci(*aer1, 42, -100, 0, t) == approx(eci0, rel=0.001)
-
-    with pytest.raises(ValueError):
-        pm.aer2eci(aer1[0], aer1[1], -1, 42, -100, 0, t)
-
-
 if __name__ == "__main__":
     pytest.main([__file__])

tests/test_eci.py

@@ -25,5 +25,49 @@ def test_ecef2eci():
     assert eci == approx([-2.9818e6, 5.2070e6, 3.1616e6], rel=0.01)
 
 
+def test_eci2geodetic():
+    pytest.importorskip("astropy")
+
+    eci = [-2981784, 5207055, 3161595]
+    utc = datetime(2019, 1, 4, 12)
+    lla = pm.eci2geodetic(*eci, utc)
+    assert lla == approx([27.881, -163.722, 408850.65], rel=0.001)
+
+
+def test_geodetic2eci():
+    pytest.importorskip("astropy")
+
+    lla = [27.881, -163.722, 408850.65]
+    utc = datetime(2019, 1, 4, 12)
+    eci = pm.geodetic2eci(*lla, utc)
+    assert eci == approx([-2981784, 5207055, 3161595], rel=0.001)
+
+
+def test_eci2aer():
+    # test coords from Matlab eci2aer
+    pytest.importorskip("astropy")
+    t = datetime(1969, 7, 20, 21, 17, 40)
+
+    eci = [-3.8454e8, -0.5099e8, -0.3255e8]
+    lla = [28.4, -80.5, 2.7]
+
+    aer = pm.eci2aer(*eci, *lla, t)
+    assert aer == approx([162.55, 55.12, 384013940.9], rel=0.001)
+
+
+def test_aer2eci():
+    # test coords from Matlab aer2eci
+    pytest.importorskip("astropy")
+
+    aer = [162.55, 55.12, 384013940.9]
+    lla = [28.4, -80.5, 2.7]
+    t = datetime(1969, 7, 20, 21, 17, 40)
+
+    assert pm.aer2eci(*aer, *lla, t) == approx([-3.8454e8, -0.5099e8, -0.3255e8], rel=0.001)
+
+    with pytest.raises(ValueError):
+        pm.aer2eci(aer[0], aer[1], -1, *lla, t)
+
+
 if __name__ == "__main__":
     pytest.main([__file__])