add return_back_azimuth: bool = True to inv()

pyproj/_geod.pyi

@@ -32,7 +32,13 @@ class Geod:
         return_back_azimuth: bool = True,
     ) -> None: ...
     def _inv(
-        self, lons1: Any, lats1: Any, lons2: Any, lats2: Any, radians: bool = False
+        self,
+        lons1: Any,
+        lats1: Any,
+        lons2: Any,
+        lats2: Any,
+        radians: bool = False,
+        return_back_azimuth: bool = False,
     ) -> None: ...
     def _inv_or_fwd_intermediate(
         self,

pyproj/_geod.pyx

@@ -175,11 +175,14 @@ cdef class Geod:
         object lons2,
         object lats2,
         bint radians=False,
+        bint return_back_azimuth=True,
     ):
         """
-        inverse transformation - return forward and back azimuths, plus distance
+        inverse transformation - return forward azimuth (azi12) and back azimuths (azi21), plus distance
         between an initial and terminus lat/lon pair.
         if radians=True, lons/lats are radians instead of degree
+        if return_back_azimuth=True, azi21 is a back azimuth (180 degrees flipped),
+        otherwise azi21 is also a forward azimuth.
         """
         cdef PyBuffWriteManager lon1buff = PyBuffWriteManager(lons1)
         cdef PyBuffWriteManager lat1buff = PyBuffWriteManager(lats1)
@@ -209,12 +212,12 @@ cdef class Geod:
                     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.
+
+                # by default (return_back_azimuth=True),
+                # forward azimuth needs to be flipped 180 degrees
+                # to match the (back azimuth) output of PROJ geod utilities.
+                if return_back_azimuth:
+                    pazi2 = _reverse_azimuth(pazi2, factor=180)
                 if radians:
                     lon1buff.data[iii] = _DG2RAD * pazi1
                     lat1buff.data[iii] = _DG2RAD * pazi2

pyproj/geod.py

@@ -327,14 +327,17 @@ def inv(
         lats2: Any,
         radians: bool = False,
         inplace: bool = False,
+        return_back_azimuth: bool = True,
     ) -> Tuple[Any, Any, Any]:
         """
+
         Inverse transformation
 
         Determine forward and back azimuths, plus distances
         between initial points and terminus points.
 
         .. versionadded:: 3.5.0 inplace
+        .. versionadded:: 3.5.0 return_back_azimuth
 
         Accepted numeric scalar or array:
 
@@ -366,13 +369,16 @@ def inv(
             If True, will attempt to write the results to the input array
             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.
 
         Returns
         -------
         scalar or array:
-            Forward azimuth(s)
+            Forward azimuth(s) (azi12)
         scalar or array:
-            Back azimuth(s)
+            Back azimuth(s) or Forward azimuth(s) (azi21)
         scalar or array:
             Distance(s) between initial and terminus point(s)
             in meters
@@ -382,7 +388,9 @@ def inv(
         iny, y_data_type = _copytobuffer(lats1, inplace=inplace)
         inz, z_data_type = _copytobuffer(lons2, inplace=inplace)
         ind = _copytobuffer(lats2, inplace=inplace)[0]
-        self._inv(inx, iny, inz, ind, radians=radians)
+        self._inv(
+            inx, iny, inz, ind, radians=radians, return_back_azimuth=return_back_azimuth
+        )
         # if inputs were lists, tuples or floats, convert back.
         outx = _convertback(x_data_type, inx)
         outy = _convertback(y_data_type, iny)

test/test_geod.py

@@ -70,10 +70,15 @@ def test_geod_inverse_transform():
     true_az12 = -66.5305947876623
     true_az21 = 75.65363415556968
     print("from pyproj.Geod.inv:")
-    az12, az21, dist = gg.inv(lon1pt, lat1pt, lon2pt, lat2pt)
-    assert_almost_equal(
-        (az12, az21, dist), (true_az12, true_az21, 4164192.708), decimal=3
-    )
+    for return_back_azimuth in (True, False):
+        az12, az21, dist = gg.inv(
+            lon1pt, lat1pt, lon2pt, lat2pt, return_back_azimuth=return_back_azimuth
+        )
+        if not return_back_azimuth:
+            az21 = reverse_azimuth(az21)
+        assert_almost_equal(
+            (az12, az21, dist), (true_az12, true_az21, 4164192.708), decimal=3
+        )
 
     print("forward transform")
     print("from proj.4 geod:")