reveal aer2eci, cleanup typing, add lat/lon sanity checks

add add'l test coverage

include tests

Author: scivision

MANIFEST.in

@@ -1 +1,2 @@
 include LICENSE
+recursive-include tests *.py

README.md

@@ -20,21 +20,10 @@
 
 PyMap3D is intended for non-interactive use on massively parallel (HPC) and embedded systems.
 Includes some relevant
-[Vallado's algorithms](http://www.smad.com/vallado/fortran/fortran.html).
+[Vallado algorithms](http://www.smad.com/vallado/fortran/fortran.html).
 
 [API docs](https://www.scivision.co/pymap3d)
 
-Why not [PyProj](https://github.com/jswhit/pyproj)?
-
--   PyMap3D does not require anything beyond pure Python + Numpy.
--   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 for PyProj ENU requires some [additional
-    effort](https://github.com/jswhit/pyproj/issues/105).
--   PyProj is oriented towards points on the planet surface, while
-    PyMap3D handles points on or above the planet surface equally well,
-    particularly important for airborne vehicles and remote sensing.
-
 ## Prerequisites
 
 * Python ≥ 3.5 or PyPy3
@@ -96,7 +85,7 @@ converted to the desired coordinate system:
 
     aer2ecef  aer2enu  aer2geodetic  aer2ned
     ecef2aer  ecef2enu  ecef2enuv  ecef2geodetic  ecef2ned  ecef2nedv
-    ecef2eci  eci2ecef
+    ecef2eci  eci2ecef  eci2aer  aer2eci
     enu2aer  enu2ecef   enu2geodetic
     geodetic2aer  geodetic2ecef  geodetic2enu  geodetic2ned
     ned2aer  ned2ecef   ned2geodetic
@@ -108,7 +97,7 @@ converted to the desired coordinate system:
 
 Additional functions:
 
-* `loxodrome_inverse`: rhumb line distance and azimuth between ellipsoid points (lat,lon)  akin to Matlab `distance('rh', ...)` and `azimuth('rh', ...)`
+`loxodrome_inverse`: rhumb line distance and azimuth between ellipsoid points (lat,lon)  akin to Matlab `distance('rh', ...)` and `azimuth('rh', ...)`
 
 
 Abbreviations:
@@ -126,3 +115,16 @@ Abbreviations:
   Would need to update code to add these input parameters (just start a GitHub Issue to request).
 * Planetary perturbations and nutation etc. not fully considered.
 
+## Notes
+
+As compared to [PyProj](https://github.com/jswhit/pyproj):
+
+-   PyMap3D does not require anything beyond pure Python + Numpy.
+-   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 for PyProj ENU requires some [additional
+    effort](https://github.com/jswhit/pyproj/issues/105).
+-   PyProj is oriented towards points on the planet surface, while
+    PyMap3D handles points on or above the planet surface equally well,
+    particularly important for airborne vehicles and remote sensing.
+

angle_distance.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
-from numpy.testing import assert_allclose
 from pymap3d.haversine import anglesep, anglesep_meeus
 from argparse import ArgumentParser
+from pytest import approx
 
 
 def main():
@@ -17,7 +17,7 @@ def main():
 
     print('{:.6f} deg sep'.format(dist_deg))
 
-    assert_allclose(dist_deg, dist_deg_astropy)
+    assert dist_deg == approx(dist_deg_astropy)
 
 
 if __name__ == '__main__':

pymap3d/__init__.py

@@ -13,6 +13,6 @@
 from .ecef import geodetic2ecef, ecef2geodetic, eci2geodetic, Ellipsoid, ecef2enuv, enu2ecef, ecef2enu, enu2uvw, uvw2enu  # noqa: F401
 from .ned import ned2ecef, ned2geodetic, geodetic2ned, ecef2nedv, ned2aer, aer2ned, ecef2ned  # noqa: F401
 from .enu import enu2geodetic, geodetic2enu, aer2enu, enu2aer  # noqa: F401
-from .aer import ecef2aer, aer2ecef, geodetic2aer, aer2geodetic, eci2aer  # noqa: F401
+from .aer import ecef2aer, aer2ecef, geodetic2aer, aer2geodetic, eci2aer, aer2eci  # noqa: F401
 from .los import lookAtSpheroid  # noqa: F401
 from .lox import isometric, meridian_dist, loxodrome_inverse  # noqa: F401

pymap3d/aer.py

@@ -1,4 +1,4 @@
-from typing import Tuple, Sequence
+from typing import Tuple
 from datetime import datetime
 import numpy as np
 from .ecef import ecef2enu, geodetic2ecef, ecef2geodetic, enu2uvw
@@ -54,7 +54,8 @@ def geodetic2aer(lat: float, lon: float, h: float,
 
 
 def aer2geodetic(az: float, el: float, srange: float,
-                 lat0: float, lon0: float, h0: float, deg: bool=True) -> Tuple[float, float, float]:
+                 lat0: float, lon0: float, h0: float,
+                 deg: bool=True) -> Tuple[float, float, float]:
     """
     Input:
     -----
@@ -75,8 +76,9 @@ def aer2geodetic(az: float, el: float, srange: float,
     return ecef2geodetic(x, y, z, deg=deg)
 
 
-def eci2aer(eci: Sequence[float], lat0: float, lon0: float, h0: float,
-            t: Sequence[datetime]) -> Tuple[float, float, float]:
+def eci2aer(eci: Tuple[float, float, float],
+            lat0: float, lon0: float, h0: float,
+            t: datetime) -> Tuple[float, float, float]:
     """
     Observer => Point
 
@@ -93,7 +95,7 @@ def eci2aer(eci: Sequence[float], lat0: float, lon0: float, h0: float,
     azimuth, elevation (degrees/radians)                             [0,360),[0,90]
     slant range [meters]                                             [0,Infinity)
     """
-    ecef = eci2ecef(eci, t)
+    ecef = np.atleast_2d(eci2ecef(eci, t))
 
     return ecef2aer(ecef[:, 0], ecef[:, 1], ecef[:, 2], lat0, lon0, h0)
 

pymap3d/azelradec.py

@@ -1,34 +1,42 @@
 """
 Azimuth / elevation <==> Right ascension, declination
 """
-from typing import Tuple, Union
+from typing import Tuple
 from datetime import datetime
 import numpy as np
-from .timeconv import str2dt
+from .vallado import azel2radec as vazel2radec, radec2azel as vradec2azel
 try:
     from astropy.time import Time
     from astropy import units as u
     from astropy.coordinates import Angle, SkyCoord, EarthLocation, AltAz, ICRS
 except ImportError:
-    from .vallado import vazel2radec, vradec2azel
     Time = None
 
 
 def azel2radec(az_deg: float, el_deg: float,
                lat_deg: float, lon_deg: float,
-               time: Union[str, datetime]) -> Tuple[float, float]:
-    """convert astronomical target horizontal azimuth, elevation to
-       ecliptic right ascension, declination (degrees)
+               time: datetime, usevallado: bool=False) -> Tuple[float, float]:
+    """
+    viewing angle (az, el) to sky coordinates (ra, dec)
+
+    inputs
+    ------
+    azimuth: degrees clockwize from North
+    elevation: degrees above horizon (neglecting aberration)
+    observer latitude [-90, 90], longitude [-180, 180] (degrees)
+    time: datetime of observation
+
+    Outputs
+    -------
+    ecliptic right ascension, declination (degrees)
     """
 
-    if Time is None:  # non-AstroPy method, less accurate
+    if usevallado or Time is None:  # non-AstroPy method, less accurate
         return vazel2radec(az_deg, el_deg, lat_deg, lon_deg, time)
 
-    t = str2dt(time)
-
     obs = EarthLocation(lat=lat_deg * u.deg, lon=lon_deg * u.deg)
 
-    direc = AltAz(location=obs, obstime=Time(t),
+    direc = AltAz(location=obs, obstime=Time(time),
                   az=az_deg * u.deg, alt=el_deg * u.deg)
 
     sky = SkyCoord(direc.transform_to(ICRS()))
@@ -38,32 +46,37 @@ def azel2radec(az_deg: float, el_deg: float,
 
 def radec2azel(ra_deg: float, dec_deg: float,
                lat_deg: float, lon_deg: float,
-               time: Union[str, datetime]) -> Tuple[float, float]:
-    """convert astronomical target ecliptic right ascension, declination to
-       horizontal azimuth, eelvation (degrees)
+               time: datetime, usevallado: bool=False) -> Tuple[float, float]:
+    """
+    sky coordinates (ra, dec) to viewing angle (az, el)
+
+    inputs
+    ------
+    ecliptic right ascension, declination (degrees)
+    observer latitude [-90, 90], longitude [-180, 180] (degrees)
+    time: datetime of observation
+
+    Outputs
+    -------
+    azimuth: degrees clockwize from North
+    elevation: degrees above horizon (neglecting aberration)
     """
-    if Time is None:
+
+    if usevallado or Time is None:
         return vradec2azel(ra_deg, dec_deg, lat_deg, lon_deg, time)
 # %% input trapping
-    t = str2dt(time)
     lat = np.atleast_1d(lat_deg)
     lon = np.atleast_1d(lon_deg)
     ra = np.atleast_1d(ra_deg)
     dec = np.atleast_1d(dec_deg)
 
-    if not(lat.size == 1 & lon.size == 1):
-        raise ValueError('radec2azel is designed for one observer and one or more points (ra,dec).')
-
-    if ra.shape != dec.shape:
-        raise ValueError('ra and dec must be the same shape ndarray')
-
     obs = EarthLocation(lat=lat * u.deg,
                         lon=lon * u.deg)
 
     points = SkyCoord(Angle(ra, unit=u.deg),
                       Angle(dec, unit=u.deg),
                       equinox='J2000.0')
 
-    altaz = points.transform_to(AltAz(location=obs, obstime=Time(t)))
+    altaz = points.transform_to(AltAz(location=obs, obstime=Time(time)))
 
     return altaz.az.degree, altaz.alt.degree

pymap3d/datetime2hourangle.py

@@ -1,10 +1,8 @@
 # Copyright (c) 2014-2018 Michael Hirsch, Ph.D.
 from math import pi
-from typing import Union, List
 from datetime import datetime
+import numpy as np
 from .timeconv import str2dt
-#
-nan = float('nan')
 try:
     from astropy.time import Time
     import astropy.units as u
@@ -18,9 +16,9 @@
 """
 
 
-def datetime2sidereal(time: Union[str, datetime],
+def datetime2sidereal(time: datetime,
                       lon_radians: float,
-                      usevallado: bool=True) -> Union[float, List[float]]:
+                      usevallado: bool=True) -> float:
     """
     Convert ``datetime`` to sidereal time
 
@@ -33,34 +31,33 @@ def datetime2sidereal(time: Union[str, datetime],
     lon
         longitude in RADIANS
     """
+    usevallado = usevallado or Time is None
     if usevallado:
-        jd = datetime2julian(time)
+        jd = datetime2julian(str2dt(time))
 # %% Greenwich Sidereal time RADIANS
-        gst = julian2sidereal(jd)
+        gst = np.atleast_1d(julian2sidereal(jd))
 # %% Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME
-        if isinstance(gst, float):
-            tsr = gst + lon_radians  # type:  Union[float, List[float]]
-        else:
-            tsr = [g + lon_radians for g in gst]
-    else:  # astropy
-        if Time is not None:
-            tsr = Time(time).sidereal_time(kind='apparent',
-                                           longitude=Longitude(lon_radians, unit=u.radian)).radian
-        else:
-            return datetime2sidereal(time, lon_radians, True)
+        tsr = gst + lon_radians
+    else:
+        tsr = Time(time).sidereal_time(kind='apparent',
+                                       longitude=Longitude(lon_radians, unit=u.radian)).radian
 
     return tsr
 
 
-def datetime2julian(time: Union[str, datetime, List[datetime]]) -> Union[float, List[float]]:
+def datetime2julian(time: datetime) -> float:
     """
     Python datetime to Julian time
 
     from D.Vallado Fundamentals of Astrodynamics and Applications p.187
      and J. Meeus Astronomical Algorithms 1991 Eqn. 7.1 pg. 61
     """
 
-    def _dt2julian(t: datetime) -> float:
+    times = np.atleast_1d(time)
+    assert times.ndim == 1
+
+    jd = np.empty(times.size)
+    for i, t in enumerate(times):
         if t.month < 3:
             year = t.year - 1
             month = t.month + 12
@@ -72,22 +69,13 @@ def _dt2julian(t: datetime) -> float:
         B = 2 - A + int(A / 4.)
         C = ((t.second / 60. + t.minute) / 60. + t.hour) / 24.
 
-        jd = (int(365.25 * (year + 4716)) +
-              int(30.6001 * (month + 1)) + t.day + B - 1524.5 + C)
-
-        return jd
+        jd[i] = (int(365.25 * (year + 4716)) +
+                 int(30.6001 * (month + 1)) + t.day + B - 1524.5 + C)
 
-    time = str2dt(time)
+    return jd.squeeze()
 
-    assert isinstance(time, datetime) or isinstance(time[0], datetime)
 
-    if isinstance(time, datetime):
-        return _dt2julian(time)
-    else:
-        return [_dt2julian(t) for t in time]
-
-
-def julian2sidereal(juliandate: Union[float, List[float]]) -> Union[float, List[float]]:
+def julian2sidereal(juliandate: float) -> float:
     """
     Convert Julian time to sidereal time
 
@@ -100,7 +88,11 @@ def julian2sidereal(juliandate: Union[float, List[float]]) -> Union[float, List[
 
     """
 
-    def _jd2sr(jd: float):
+    jdate = np.atleast_1d(juliandate)
+    assert jdate.ndim == 1
+
+    tsr = np.empty(jdate.size)
+    for i, jd in enumerate(jdate):
         # %% Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1
         tUT1 = (jd - 2451545.0) / 36525.
 
@@ -109,11 +101,6 @@ def _jd2sr(jd: float):
                     tUT1 + 0.093104 * tUT1**2 - 6.2e-6 * tUT1**3)
 
         # 1/86400 and %(2*pi) implied by units of radians
-        tsr = gmst_sec * (2 * pi) / 86400. % (2 * pi)
+        tsr[i] = gmst_sec * (2 * pi) / 86400. % (2 * pi)
 
-        return tsr
-
-    if isinstance(juliandate, float):
-        return _jd2sr(juliandate)
-    else:
-        return [_jd2sr(jd) for jd in juliandate]
+    return tsr.squeeze()