Add return_back_azimuth: bool = False to allow compatibility betwee…

…n the azimuth output of `fwd` and `fwd_intermediate`

docs/history.rst

@@ -8,6 +8,7 @@ Latest
 -----
 - BUG: Changed so that the setup.cfg depends on the version code in the __init__.py instead of the other way around (issuue #1155)
 - REF: Use upper case EPSG code when creating CRS (pull #1162)
+- ENH: Add `return_back_azimuth: bool = False` to allow compatibility between the azimuth output of `fwd` and `fwd_intermediate` (issue #1163)
 
 3.4.0
 -----

pyproj/_geod.pyi

@@ -43,6 +43,7 @@ class Geod:
         out_lons: Any,
         out_lats: Any,
         out_azis: Any,
+        return_back_azimuth: bool,
     ) -> GeodIntermediateReturn: ...
     def _line_length(self, lons: Any, lats: Any, radians: bool = False) -> float: ...
     def _polygon_area_perimeter(

pyproj/_geod.pyx

@@ -221,6 +221,7 @@ cdef class Geod:
         object out_lons,
         object out_lats,
         object out_azis,
+        bint return_back_azimuth,
     ) -> GeodIntermediateReturn:
         """
         .. versionadded:: 3.1.0
@@ -312,6 +313,11 @@ cdef class Geod:
                 lats_buff.data[iii] = plat2
                 lons_buff.data[iii] = plon2
                 if store_az:
+                    if return_back_azimuth:
+                        if pazi2 > 0:
+                            pazi2 = pazi2 - 180.
+                        elif pazi2 <= 0:
+                            pazi2 = pazi2 + 180.
                     azis_buff.data[iii] = pazi2
 
         return GeodIntermediateReturn(

pyproj/geod.py

@@ -469,6 +469,7 @@ def npts(
             out_lons=None,
             out_lats=None,
             out_azis=None,
+            return_back_azimuth=False,
         )
         return list(zip(res.lons, res.lats))
 
@@ -487,6 +488,7 @@ def inv_intermediate(
         out_lons: Any = None,
         out_lats: Any = None,
         out_azis: Any = None,
+        return_back_azimuth: bool = False,
     ) -> GeodIntermediateReturn:
         """
         .. versionadded:: 3.1.0
@@ -598,6 +600,9 @@ def inv_intermediate(
             az12(s) of the intermediate point(s)
             If None then buffers would be allocated internnaly
             unless requested otherwise by the flags
+        return_back_azimuth: bool, default=False
+            if True, out_azis will contain back azimuth,
+            Otherwise, out_azis will contain forward azimuth.
 
         Returns
         -------
@@ -618,6 +623,7 @@ def inv_intermediate(
             out_lons=out_lons,
             out_lats=out_lats,
             out_azis=out_azis,
+            return_back_azimuth=return_back_azimuth,
         )
 
     def fwd_intermediate(
@@ -634,6 +640,7 @@ def fwd_intermediate(
         out_lons: Any = None,
         out_lats: Any = None,
         out_azis: Any = None,
+        return_back_azimuth: bool = False,
     ) -> GeodIntermediateReturn:
         """
         .. versionadded:: 3.1.0
@@ -730,6 +737,9 @@ def fwd_intermediate(
             az12(s) of the intermediate point(s)
             If None then buffers would be allocated internnaly
             unless requested otherwise by the flags
+        return_back_azimuth: bool, default=False
+            if True, out_azis will contain back azimuth,
+            Otherwise, out_azis will contain forward azimuth.
 
         Returns
         -------
@@ -750,6 +760,7 @@ def fwd_intermediate(
             out_lons=out_lons,
             out_lats=out_lats,
             out_azis=out_azis,
+            return_back_azimuth=return_back_azimuth,
         )
 
     def line_length(self, lons: Any, lats: Any, radians: bool = False) -> float:

test/test_geod.py

@@ -46,6 +46,15 @@ def temporary_directory():
         shutil.rmtree(temp_dir)
 
 
+def reverse_azimuth(azi: float) -> float:
+    return azi - 180. if azi > 0 else azi + 180.
+
+
+def reverse_azimuth_np(arr: np.ndarray):
+    for i in range(len(arr)):
+        arr[i] = reverse_azimuth(arr[i])
+
+
 def test_geod_inverse_transform():
     gg = Geod(ellps="clrk66")
     lat1pt = 42.0 + (15.0 / 60.0)
@@ -230,27 +239,35 @@ def test_geod_inverse_transform():
             assert_almost_equal(res.azis, azis_a)
 
     print("test fwd_intermediate")
-    res = gg.fwd_intermediate(
-        out_lons=lons_b,
-        out_lats=lats_b,
-        out_azis=azis_b,
-        lon1=lon1pt,
-        lat1=lat1pt,
-        azi1=true_az12,
-        npts=point_count,
-        del_s=del_s,
-        initial_idx=0,
-        terminus_idx=0,
-    )
-    assert res.npts == point_count
-    assert_almost_equal(res.del_s, del_s)
-    assert_almost_equal(res.dist, dist)
-    assert_almost_equal(res.lons, lons_a)
-    assert_almost_equal(res.lats, lats_a)
-    assert_almost_equal(res.azis, azis_a)
-    assert res.lons is lons_b
-    assert res.lats is lats_b
-    assert res.azis is azis_b
+    for return_back_azimuth in [False, True]:
+        lons_b = np.empty(point_count)
+        lats_b = np.empty(point_count)
+        azis_b = np.empty(point_count)
+
+        res = gg.fwd_intermediate(
+            out_lons=lons_b,
+            out_lats=lats_b,
+            out_azis=azis_b,
+            lon1=lon1pt,
+            lat1=lat1pt,
+            azi1=true_az12,
+            npts=point_count,
+            del_s=del_s,
+            initial_idx=0,
+            terminus_idx=0,
+            return_back_azimuth=return_back_azimuth,
+        )
+        assert res.npts == point_count
+        assert_almost_equal(res.del_s, del_s)
+        assert_almost_equal(res.dist, dist)
+        assert_almost_equal(res.lons, lons_a)
+        assert_almost_equal(res.lats, lats_a)
+        if return_back_azimuth:
+            reverse_azimuth_np(azis_b)
+        assert_almost_equal(res.azis, azis_a)
+        assert res.lons is lons_b
+        assert res.lats is lats_b
+        assert res.azis is azis_b
 
     print("test inv_intermediate (by del_s)")
     for del_s_fact, flags in (