diff --git a/src/pymap3d/aer.py b/src/pymap3d/aer.py index 631345b..c8ad144 100644 --- a/src/pymap3d/aer.py +++ b/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( diff --git a/src/pymap3d/ecef.py b/src/pymap3d/ecef.py index 88423a0..b6b543d 100644 --- a/src/pymap3d/ecef.py +++ b/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 ) diff --git a/src/pymap3d/eci.py b/src/pymap3d/eci.py index 7c552c8..4687393 100644 --- a/src/pymap3d/eci.py +++ b/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]] + ) diff --git a/src/pymap3d/enu.py b/src/pymap3d/enu.py index a91764f..a33919f 100644 --- a/src/pymap3d/enu.py +++ b/src/pymap3d/enu.py @@ -1,4 +1,5 @@ """ transforms involving ENU East North Up """ + from __future__ import annotations from math import tau diff --git a/src/pymap3d/spherical.py b/src/pymap3d/spherical.py index ecfc4ed..031d0d0 100644 --- a/src/pymap3d/spherical.py +++ b/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 diff --git a/src/pymap3d/tests/test_eci.py b/src/pymap3d/tests/test_eci.py index cc8b472..4e4836d 100644 --- a/src/pymap3d/tests/test_eci.py +++ b/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]