los/lox, math and numpy.vectorize

Author: scivision

pymap3d/__init__.py

@@ -37,13 +37,12 @@
 from .ned import ned2ecef, ned2geodetic, geodetic2ned, ecef2nedv, ned2aer, aer2ned, ecef2ned
 from .ecef import geodetic2ecef, ecef2geodetic, eci2geodetic, ecef2enuv, enu2ecef, ecef2enu, enu2uvw, uvw2enu
 from .sidereal import datetime2sidereal
+from .lox import isometric, meridian_dist, loxodrome_inverse
+from .los import lookAtSpheroid
+from .timeconv import str2dt
 
 try:
-    from .timeconv import str2dt
     from .azelradec import radec2azel, azel2radec
     from .eci import eci2ecef, ecef2eci
-
-    from .los import lookAtSpheroid
-    from .lox import isometric, meridian_dist, loxodrome_inverse
 except ImportError:
     from .vallado import radec2azel, azel2radec

pymap3d/los.py

@@ -1,16 +1,29 @@
 """ Line of sight intersection of space observer to ellipsoid """
 from typing import Tuple
-import numpy as np
-from math import pi
+from math import pi, nan, sqrt
+
 from .aer import aer2enu
 from .ecef import enu2uvw, geodetic2ecef, ecef2geodetic
 from .ellipsoid import Ellipsoid
+try:
+    import numpy
+except ImportError:
+    numpy = None
 
 __all__ = ['lookAtSpheroid']
 
 
 def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float,
                    ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
+    if numpy is not None:
+        fun = numpy.vectorize(lookAtSpheroid_point)
+        return fun(lat0, lon0, h0, az, tilt, ell, deg)
+    else:
+        return lookAtSpheroid_point(lat0, lon0, h0, az, tilt, ell, deg)
+
+
+def lookAtSpheroid_point(lat0: float, lon0: float, h0: float, az: float, tilt: float,
+                         ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
     Calculates line-of-sight intersection with Earth (or other ellipsoid) surface from above surface / orbit
 
@@ -47,14 +60,12 @@ def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float,
     Algorithm based on https://medium.com/@stephenhartzell/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6 Stephen Hartzell
     """
 
-    if (np.asarray(h0) < 0).any():
+    if h0 < 0:
         raise ValueError('Intersection calculation requires altitude  [0, Infinity)')
 
     if ell is None:
         ell = Ellipsoid()
 
-    tilt = np.asarray(tilt)
-
     a = ell.semimajor_axis
     b = ell.semimajor_axis
     c = ell.semiminor_axis
@@ -73,11 +84,13 @@ def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float,
     magnitude = a**2 * b**2 * w**2 + a**2 * c**2 * v**2 + b**2 * c**2 * u**2
 
 # %%   Return nan if radical < 0 or d < 0 because LOS vector does not point towards Earth
-    with np.errstate(invalid='ignore'):
-        d = np.where(radical > 0,
-                     (value - a * b * c * np.sqrt(radical)) / magnitude,
-                     np.nan)
-        d[d < 0] = np.nan
+    if radical > 0:
+        d = (value - a * b * c * sqrt(radical)) / magnitude
+    else:
+        d = nan
+
+    if d < 0:
+        d = nan
 # %% cartesian to ellipsodal
     lat, lon, _ = ecef2geodetic(x + d * u, y + d * v, z + d * w, deg=deg)
 

pymap3d/lox.py

@@ -1,5 +1,5 @@
 """ isometric latitude, meridian distance """
-import numpy as np
+from math import radians, degrees, sin, cos, atan2, tan, log, sqrt, pi
 from .ellipsoid import Ellipsoid
 
 
@@ -40,20 +40,20 @@ def isometric(lat: float, ell: Ellipsoid = None, deg: bool = True):
     f = ell.flattening  # flattening of ellipsoid
 
     if deg is True:
-        lat = np.deg2rad(lat)
+        lat = radians(lat)
 
     e2 = f * (2 - f)  # eccentricity-squared
-    e = np.sqrt(e2)  # eccentricity of ellipsoid
+    e = sqrt(e2)  # eccentricity of ellipsoid
 
-    x = e * np.sin(lat)
+    x = e * sin(lat)
     y = (1 - x) / (1 + x)
-    z = np.pi / 4 + lat / 2
+    z = pi / 4 + lat / 2
 
 #   calculate the isometric latitude
-    isolat = np.log(np.tan(z) * (y**(e / 2)))
+    isolat = log(tan(z) * (y**(e / 2)))
 
     if deg is True:
-        isolat = np.degrees(isolat)
+        isolat = degrees(isolat)
 
     return isolat
 
@@ -90,7 +90,7 @@ def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True):
     """
 
     if deg is True:
-        lat = np.radians(lat)
+        lat = radians(lat)
 
     #   set ellipsoid parameters
     if ell is None:
@@ -116,11 +116,11 @@ def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True):
     F = (693 / 131072) * e10
 
     term1 = A * lat
-    term2 = (B / 2) * np.sin(2 * lat)
-    term3 = (C / 4) * np.sin(4 * lat)
-    term4 = (D / 6) * np.sin(6 * lat)
-    term5 = (E / 8) * np.sin(8 * lat)
-    term6 = (F / 10) * np.sin(10 * lat)
+    term2 = (B / 2) * sin(2 * lat)
+    term3 = (C / 4) * sin(4 * lat)
+    term4 = (D / 6) * sin(6 * lat)
+    term5 = (E / 8) * sin(8 * lat)
+    term6 = (F / 10) * sin(10 * lat)
 
     mdist = a * (1 - e2) * (term1 - term2 + term3 - term4 + term5 - term6)
 
@@ -177,7 +177,7 @@ def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
         ell = Ellipsoid()
 
     if deg is True:
-        lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
+        lat1, lon1, lat2, lon2 = radians(lat1), radians(lon1), radians(lat2), radians(lon2)
 
     # compute isometric latitude of P1 and P2
     isolat1 = isometric(lat1, deg=False, ell=ell)
@@ -188,15 +188,15 @@ def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
     dlon = lon2 - lon1
 
     # compute azimuth
-    az12 = np.arctan2(dlon, disolat)
+    az12 = atan2(dlon, disolat)
 
     # compute distance along loxodromic curve
     m1 = meridian_dist(lat1, deg=False, ell=ell)
     m2 = meridian_dist(lat2, deg=False, ell=ell)
     dm = m2 - m1
-    lox_s = dm / np.cos(az12)
+    lox_s = dm / cos(az12)
 
     if deg is True:
-        az12 = np.degrees(az12) % 360.
+        az12 = degrees(az12) % 360.
 
     return lox_s, az12