PERF: Add scalar optimized Geod fwd/inv functions (#1206)

This adds a fast-path for scalar inputs and falls back to
the previous implementation if any inputs can't be cast to double.

pyproj/_geod.pyi

@@ -31,6 +31,15 @@ class Geod:
         radians: bool = False,
         return_back_azimuth: bool = True,
     ) -> None: ...
+    def _fwd_point(
+        self,
+        lons: float,
+        lats: float,
+        az: float,
+        dist: float,
+        radians: bool = False,
+        return_back_azimuth: bool = True,
+    ) -> Tuple[float, float, float]: ...
     def _inv(
         self,
         lons1: Any,
@@ -40,6 +49,15 @@ class Geod:
         radians: bool = False,
         return_back_azimuth: bool = False,
     ) -> None: ...
+    def _inv_point(
+        self,
+        lons1: float,
+        lats1: float,
+        lons2: float,
+        lats2: float,
+        radians: bool = False,
+        return_back_azimuth: bool = False,
+    ) -> Tuple[float, float, float]: ...
     def _inv_or_fwd_intermediate(
         self,
         lon1: float,

pyproj/_geod.pyx

@@ -166,6 +166,63 @@ cdef class Geod:
                     latbuff.data[iii] = _DG2RAD * plat2
                     azbuff.data[iii] = _DG2RAD * pazi2
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _fwd_point(
+        self,
+        object lon1in,
+        object lat1in,
+        object az1in,
+        object s12in,
+        bint radians=False,
+        bint return_back_azimuth=True,
+    ):
+        """
+        Scalar optimized function
+        forward transformation - determine longitude, latitude and back azimuth
+        of a terminus point given an initial point longitude and latitude, plus
+        forward azimuth and distance.
+        if radians=True, lons/lats are radians instead of degrees.
+        """
+        cdef double plon2, plat2, pazi2
+
+        # We do the type-checking internally here rather than declaring double as the
+        # type into the function due to automatically casting length-1 arrays to float
+        # that we don't want to return scalar for.
+        # Ex: float(np.array([0])) works and we don't want to accept numpy arrays
+        cdef double lon1 = lon1in
+        cdef double lat1 = lat1in
+        cdef double az1 = az1in
+        cdef double s12 = s12in
+        for x in (lon1in, lat1in, az1in, s12in):
+            if not isinstance(x, (float, int)):
+                raise TypeError("Scalar input is required for point based functions")
+
+        with nogil:
+            if radians:
+                lon1 = _RAD2DG * lon1
+                lat1 = _RAD2DG * lat1
+                az1 = _RAD2DG * az1
+            geod_direct(
+                &self._geod_geodesic,
+                lat1,
+                lon1,
+                az1,
+                s12,
+                &plat2,
+                &plon2,
+                &pazi2,
+            )
+            # back azimuth needs to be flipped 180 degrees
+            # to match what PROJ geod utility produces.
+            if return_back_azimuth:
+                pazi2 =_reverse_azimuth(pazi2, factor=180)
+            if radians:
+                plon2 = _DG2RAD * plon2
+                plat2 = _DG2RAD * plat2
+                pazi2 = _DG2RAD * pazi2
+        return plon2, plat2, pazi2
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _inv(
@@ -227,6 +284,56 @@ cdef class Geod:
                 # write azimuth data into lon2 buffer
                 lon2buff.data[iii] = ps12
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _inv_point(
+        self,
+        object lon1in,
+        object lat1in,
+        object lon2in,
+        object lat2in,
+        bint radians=False,
+        bint return_back_azimuth=True,
+    ):
+        """
+        Scalar optimized function
+        inverse transformation - return forward and back azimuth, plus distance
+        between an initial and terminus lat/lon pair.
+        if radians=True, lons/lats are radians instead of degree
+        """
+        cdef double pazi1, pazi2, ps12
+        # We do the type-checking internally here rather than declaring double as the
+        # type into the function due to automatically casting length-1 arrays to float
+        # that we don't want to return scalar for.
+        # Ex: float(np.array([0])) works and we don't want to accept numpy arrays
+        cdef double lon1 = lon1in
+        cdef double lat1 = lat1in
+        cdef double lon2 = lon2in
+        cdef double lat2 = lat2in
+        for x in (lon1in, lat1in, lon2in, lat2in):
+            if not isinstance(x, (float, int)):
+                raise TypeError("Scalar input is required for point based functions")
+
+        with nogil:
+            if radians:
+                lon1 = _RAD2DG * lon1
+                lat1 = _RAD2DG * lat1
+                lon2 = _RAD2DG * lon2
+                lat2 = _RAD2DG * lat2
+            geod_inverse(
+                &self._geod_geodesic,
+                lat1, lon1, lat2, lon2,
+                &ps12, &pazi1, &pazi2,
+            )
+            # back azimuth needs to be flipped 180 degrees
+            # to match what proj4 geod utility produces.
+            if return_back_azimuth:
+                pazi2 =_reverse_azimuth(pazi2, factor=180)
+            if radians:
+                pazi1 = _DG2RAD * pazi1
+                pazi2 = _DG2RAD * pazi2
+        return pazi1, pazi2, ps12
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _inv_or_fwd_intermediate(

pyproj/geod.py

@@ -294,7 +294,7 @@ def fwd(  # pylint: disable=invalid-name
             instead of returning a new array. This will fail if the input
             is not an array in C order with the double data type.
         return_back_azimuth: bool, default=True
-            if True, the third return value will be the back azimuth,
+            If True, the third return value will be the back azimuth,
             Otherwise, it will be the forward azimuth.
 
         Returns
@@ -306,6 +306,20 @@ def fwd(  # pylint: disable=invalid-name
         scalar or array:
             Back azimuth(s) or Forward azimuth(s)
         """
+        try:
+            # Fast-path for scalar input, will raise if invalid types are input
+            # and we can fallback below
+            return self._fwd_point(
+                lons,
+                lats,
+                az,
+                dist,
+                radians=radians,
+                return_back_azimuth=return_back_azimuth,
+            )
+        except TypeError:
+            pass
+
         # process inputs, making copies that support buffer API.
         inx, x_data_type = _copytobuffer(lons, inplace=inplace)
         iny, y_data_type = _copytobuffer(lats, inplace=inplace)
@@ -371,8 +385,8 @@ def inv(
             instead of returning a new array. This will fail if the input
             is not an array in C order with the double data type.
         return_back_azimuth: bool, default=True
-            if True, the second return value (azi21) will be the back azimuth
-            (flipped 180 degrees), Otherwise, it will be also a forward azimuth.
+            If True, the second return value (azi21) will be the back azimuth
+            (flipped 180 degrees), Otherwise, it will also be a forward azimuth.
 
         Returns
         -------
@@ -384,6 +398,20 @@ def inv(
             Distance(s) between initial and terminus point(s)
             in meters
         """
+        try:
+            # Fast-path for scalar input, will raise if invalid types are input
+            # and we can fallback below
+            return self._inv_point(
+                lons1,
+                lats1,
+                lons2,
+                lats2,
+                radians=radians,
+                return_back_azimuth=return_back_azimuth,
+            )
+        except TypeError:
+            pass
+
         # process inputs, making copies that support buffer API.
         inx, x_data_type = _copytobuffer(lons1, inplace=inplace)
         iny, y_data_type = _copytobuffer(lats1, inplace=inplace)

test/test_geod.py

@@ -678,6 +678,51 @@ def test_geod_fwd_honours_input_types(lon, lat, az):
     assert isinstance(outz, type(az))
 
 
+def test_geod_fwd_radians():
+    g = Geod(ellps="clrk66")
+    lon1 = 1
+    lat1 = 1
+    az1 = 1
+    dist = 1
+    assert_almost_equal(
+        np.rad2deg(g.fwd(lon1, lat1, az1, dist, radians=True)),
+        g.fwd(lon1 * 180 / np.pi, lat1 * 180 / np.pi, az1 * 180 / np.pi, dist),
+    )
+
+
+def test_geod_inv_radians():
+    g = Geod(ellps="clrk66")
+    lon1 = 0
+    lat1 = 0
+    lon2 = 1
+    lat2 = 1
+    # the third output is in distance, so we don't want to change from deg-rad there
+    out_rad = list(g.inv(lon1, lat1, lon2, lat2, radians=True))
+    out_rad[0] *= 180 / np.pi
+    out_rad[1] *= 180 / np.pi
+    assert_almost_equal(
+        out_rad,
+        g.inv(
+            lon1 * 180 / np.pi,
+            lat1 * 180 / np.pi,
+            lon2 * 180 / np.pi,
+            lat2 * 180 / np.pi,
+        ),
+    )
+
+
+@pytest.mark.parametrize("func_name", ("fwd", "inv"))
+@pytest.mark.parametrize("radians", (True, False))
+def test_geod_scalar_array(func_name, radians):
+    # verify two singlepoint calculations match an array of length two
+    g = Geod(ellps="clrk66")
+    func = getattr(g, func_name)
+    assert_almost_equal(
+        np.transpose([func(0, 0, 1, 1, radians=radians) for i in range(2)]),
+        func([0, 0], [0, 0], [1, 1], [1, 1], radians=radians),
+    )
+
+
 @pytest.mark.parametrize(
     "lons1,lats1,lons2", permutations([10.0, [10.0], (10.0,)])
 )  # 6 test cases