ENH: Add single point optimized Geod fwd/inv functions

pyproj/_geod.pyi

@@ -25,9 +25,20 @@ class Geod:
     def _fwd(
         self, lons: Any, lats: Any, az: Any, dist: Any, radians: bool = False
     ) -> None: ...
+    def _fwd_point(
+        self, lons: float, lats: float, az: float, dist: float, radians: bool = False
+    ) -> Tuple[float, float, float]: ...
     def _inv(
         self, lons1: Any, lats1: Any, lons2: Any, lats2: Any, radians: bool = False
     ) -> None: ...
+    def _inv_point(
+        self,
+        lons1: float,
+        lats1: float,
+        lons2: float,
+        lats2: float,
+        radians: bool = False,
+    ) -> Tuple[float, float, float]: ...
     def _inv_or_fwd_intermediate(
         self,
         lon1: float,

pyproj/_geod.pyx

@@ -146,6 +146,52 @@ cdef class Geod:
                     latbuff.data[iii] = _DG2RAD * plat2
                     azbuff.data[iii] = _DG2RAD * pazi2
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _fwd_point(
+        self,
+        double lon1,
+        double lat1,
+        double az1,
+        double s12,
+        bint radians=False,
+    ):
+        """
+        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
+        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 pazi2 > 0:
+                pazi2 = pazi2 - 180.
+            elif pazi2 <= 0:
+                pazi2 = pazi2 + 180.
+            if radians:
+                plon2 = _DG2RAD * plon2
+                plat2 = _DG2RAD * plat2
+                pazi2 = _DG2RAD * pazi2
+        return plon2, plat2, pazi2
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _inv(
@@ -204,6 +250,45 @@ cdef class Geod:
                 # write azimuth data into lon2 buffer
                 lon2buff.data[iii] = ps12
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _inv_point(
+        self,
+        double lon1,
+        double lat1,
+        double lon2,
+        double lat2,
+        bint radians=False,
+    ):
+        """
+        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
+        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 pazi2 > 0:
+                pazi2 = pazi2 - 180.
+            elif pazi2 <= 0:
+                pazi2 = pazi2 + 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

@@ -298,6 +298,15 @@ def fwd(  # pylint: disable=invalid-name
         scalar or array:
             Back azimuth(s)
         """
+        if not inplace:
+            # The scalar implementation is not inplace
+            try:
+                # Fast-path for scalar input that we can try immediately and fallback
+                # to the buffer implementation below if we get any wrong input
+                return self._fwd_point(lons, lats, az, dist, radians=radians)
+            except TypeError:
+                # wrong types passed in, so fallback to creating the buffers below
+                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)
@@ -368,6 +377,15 @@ def inv(
             Distance(s) between initial and terminus point(s)
             in meters
         """
+        if not inplace:
+            # The scalar implementation is not inplace
+            try:
+                # Fast-path for scalar input that we can try immediately and fallback
+                # to the buffer implementation below if we get any wrong input
+                return self._inv_point(lons1, lats1, lons2, lats2, radians=radians)
+            except TypeError:
+                # wrong types passed in, so fallback to creating the buffers below
+                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

@@ -646,6 +646,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