ecef2geodetic 100x speedup from not using vectorize

Author: scivision

scripts/benchmark_ecef2geo.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+"""
+benchmark ecef2geodetic
+"""
+import time
+from pymap3d.ecef import ecef2geodetic
+import numpy as np
+import argparse
+
+ll0 = (42.0, 82.0)
+
+
+def bench(N: int) -> float:
+
+    x = np.random.random(N)
+    y = np.random.random(N)
+    z = np.random.random(N)
+
+    tic = time.monotonic()
+    lat, lon, alt = ecef2geodetic(x, y, z)
+
+    return time.monotonic() - tic
+
+
+if __name__ == "__main__":
+    p = argparse.ArgumentParser()
+    p.add_argument("N", type=int)
+    p = p.parse_args()
+    N = p.N
+
+    print(f"ecef2geodetic: {bench(N):.3f} seconds")

src/pymap3d/aer.py

@@ -20,15 +20,15 @@
 
 
 def ecef2aer(
-    x: float | ndarray,
-    y: float | ndarray,
-    z: float | ndarray,
+    x: ndarray,
+    y: ndarray,
+    z: ndarray,
     lat0: float,
     lon0: float,
     h0: float,
     ell: Ellipsoid = None,
     deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     compute azimuth, elevation and slant range from an Observer to a Point with ECEF coordinates.
 
@@ -69,15 +69,15 @@ def ecef2aer(
 
 
 def geodetic2aer(
-    lat: float,
-    lon: float,
-    h: float,
+    lat: ndarray,
+    lon: ndarray,
+    h: ndarray,
     lat0: float,
     lon0: float,
     h0: float,
     ell: Ellipsoid = None,
     deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     gives azimuth, elevation and slant range from an Observer to a Point with geodetic coordinates.
 
@@ -177,7 +177,7 @@ def eci2aer(
     *,
     deg: bool = True,
     use_astropy: bool = True
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     takes Earth Centered Inertial x,y,z ECI coordinates of point and gives az, el, slant range from Observer
 

src/pymap3d/ecef.py

@@ -15,12 +15,10 @@
         arctan2 as atan2,
         sqrt,
         pi,
-        vectorize,
     )
 except ImportError:
     from math import radians, sin, cos, tan, atan, hypot, degrees, atan2, sqrt, pi  # type: ignore
 
-    vectorize = None
     ndarray = typing.Any  # type: ignore
 
 from datetime import datetime
@@ -110,21 +108,6 @@ def ecef2geodetic(
     z: float | ndarray,
     ell: Ellipsoid = None,
     deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
-    if vectorize is not None:
-        fun = vectorize(ecef2geodetic_point)
-        lat, lon, alt = fun(x, y, z, ell, deg)
-        return lat[()], lon[()], alt[()]
-    else:
-        return ecef2geodetic_point(x, y, z, ell, deg)
-
-
-def ecef2geodetic_point(
-    x: float | ndarray,
-    y: float | ndarray,
-    z: float | ndarray,
-    ell: Ellipsoid = None,
-    deg: bool = True,
 ) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
     """
     convert ECEF (meters) to geodetic coordinates
@@ -155,6 +138,7 @@ def ecef2geodetic_point(
     You, Rey-Jer. (2000). Transformation of Cartesian to Geodetic Coordinates without Iterations.
     Journal of Surveying Engineering. doi: 10.1061/(ASCE)0733-9453
     """
+
     if ell is None:
         ell = Ellipsoid()
 
@@ -199,8 +183,14 @@ def ecef2geodetic_point(
         + z ** 2 / ell.semiminor_axis ** 2
         < 1
     )
-    if inside:
-        alt = -alt
+
+    try:
+        if inside.any():  # type: ignore
+            # avoid all false assignment bug
+            alt[inside] = -alt
+    except (TypeError, AttributeError):
+        if inside:
+            alt = -alt
 
     if deg:
         lat = degrees(lat)
@@ -255,15 +245,15 @@ def ecef2enuv(
 
 
 def ecef2enu(
-    x: float | ndarray,
-    y: float | ndarray,
-    z: float | ndarray,
+    x: ndarray,
+    y: ndarray,
+    z: ndarray,
     lat0: float,
     lon0: float,
     h0: float,
     ell: Ellipsoid = None,
     deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     from observer to target, ECEF => ENU
 
@@ -343,7 +333,7 @@ def enu2uvw(
 
 def uvw2enu(
     u: float, v: float, w: float, lat0: float, lon0: float, deg: bool = True
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     Parameters
     ----------

src/pymap3d/enu.py

@@ -3,11 +3,10 @@
 import typing
 
 try:
-    from numpy import radians, sin, cos, hypot, arctan2 as atan2, degrees, pi, vectorize, ndarray
+    from numpy import asarray, radians, sin, cos, hypot, arctan2 as atan2, degrees, pi, ndarray
 except ImportError:
     from math import radians, sin, cos, hypot, atan2, degrees, pi  # type: ignore
 
-    vectorize = None
     ndarray = typing.Any  # type: ignore
 
 from .ecef import geodetic2ecef, ecef2geodetic, enu2ecef, uvw2enu
@@ -20,8 +19,8 @@
 
 
 def enu2aer(
-    e: float | ndarray, n: float | ndarray, u: float | ndarray, deg: bool = True
-) -> tuple[float, float, float]:
+    e: ndarray, n: ndarray, u: ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     ENU to Azimuth, Elevation, Range
 
@@ -48,26 +47,19 @@ def enu2aer(
         slant range [meters]
     """
 
-    if vectorize is not None:
-        fun = vectorize(enu2aer_point)
-        az, el, rng = fun(e, n, u, deg)
-        return az[()], el[()], rng[()]
-    else:
-        return enu2aer_point(e, n, u, deg)
-
-
-def enu2aer_point(
-    e: float | ndarray, n: float | ndarray, u: float | ndarray, deg: bool = True
-) -> tuple[float, float, float]:
-
     # 1 millimeter precision for singularity
 
-    if abs(e) < 1e-3:
-        e = 0.0
-    if abs(n) < 1e-3:
-        n = 0.0
-    if abs(u) < 1e-3:
-        u = 0.0
+    try:
+        e[abs(e) < 1e-3] = 0.0
+        n[abs(n) < 1e-3] = 0.0
+        u[abs(u) < 1e-3] = 0.0
+    except TypeError:
+        if abs(e) < 1e-3:
+            e = 0.0  # type: ignore
+        if abs(n) < 1e-3:
+            n = 0.0  # type: ignore
+        if abs(u) < 1e-3:
+            u = 0.0  # type: ignore
 
     r = hypot(e, n)
     slantRange = hypot(r, u)
@@ -82,17 +74,6 @@ def enu2aer_point(
 
 
 def aer2enu(az: float, el: float, srange: float, deg: bool = True) -> tuple[float, float, float]:
-    if vectorize is not None:
-        fun = vectorize(aer2enu_point)
-        e, n, u = fun(az, el, srange, deg)
-        return e[()], n[()], u[()]
-    else:
-        return aer2enu_point(az, el, srange, deg)
-
-
-def aer2enu_point(
-    az: float, el: float, srange: float, deg: bool = True
-) -> tuple[float, float, float]:
     """
     Azimuth, Elevation, Slant range to target to East, North, Up
 
@@ -120,8 +101,12 @@ def aer2enu_point(
         el = radians(el)
         az = radians(az)
 
-    if srange < 0:
-        raise ValueError("Slant range  [0, Infinity)")
+    try:
+        if (asarray(srange) < 0).any():
+            raise ValueError("Slant range  [0, Infinity)")
+    except UnboundLocalError:
+        if srange < 0:
+            raise ValueError("Slant range  [0, Infinity)")
 
     r = srange * cos(el)
 
@@ -177,15 +162,15 @@ def enu2geodetic(
 
 
 def geodetic2enu(
-    lat: float,
-    lon: float,
-    h: float,
+    lat: ndarray,
+    lon: ndarray,
+    h: ndarray,
     lat0: float,
     lon0: float,
     h0: float,
     ell: Ellipsoid = None,
     deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     Parameters
     ----------

src/pymap3d/ned.py

@@ -45,7 +45,9 @@ def aer2ned(
     return n, e, -u
 
 
-def ned2aer(n: float, e: float, d: float, deg: bool = True) -> tuple[float, float, float]:
+def ned2aer(
+    n: ndarray, e: ndarray, d: ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
     """
     converts North, East, Down to azimuth, elevation, range
 
@@ -170,9 +172,9 @@ def ned2ecef(
 
 
 def ecef2ned(
-    x: float,
-    y: float,
-    z: float,
+    x: ndarray,
+    y: ndarray,
+    z: ndarray,
     lat0: float,
     lon0: float,
     h0: float,
@@ -219,9 +221,9 @@ def ecef2ned(
 
 
 def geodetic2ned(
-    lat: float,
-    lon: float,
-    h: float,
+    lat: ndarray,
+    lon: ndarray,
+    h: ndarray,
     lat0: float,
     lon0: float,
     h0: float,

src/pymap3d/tests/test_geodetic.py

@@ -37,17 +37,11 @@ def test_3d_geodetic2ecef():
     assert (x0, y0, z0) == approx(xyz0)
 
 
-@pytest.mark.parametrize(
-    "xyz", [(xyz0[0], xyz0[1], xyz0[2]), ([xyz0[0]], [xyz0[1]], [xyz0[2]])], ids=("scalar", "list")
-)
-def test_scalar_ecef2geodetic(xyz):
+def test_scalar_ecef2geodetic():
     """
     verify we can handle the wide variety of input data type users might use
     """
-    if isinstance(xyz[0], list):
-        pytest.importorskip("numpy")
-
-    lat, lon, alt = pm.ecef2geodetic(*xyz)
+    lat, lon, alt = pm.ecef2geodetic(xyz0[0], xyz0[1], xyz0[2])
 
     assert [lat, lon, alt] == approx(lla0, rel=1e-4)