geod.py - merge reverse_azimuth and reverse_azimuth_arr

geod.py, _geod.pyx - fix `inv`, `_inv` doc string for `return_back_azimuth` parameter,
test/test_geod.py - add some tests for `return_back_azimuth` and `reverse_azimuth`

pyproj/_geod.pyx

@@ -162,10 +162,11 @@ cdef class Geod:
         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, the return azimuth will be the forward azimuth instead of the forward azimuth.
+        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)

pyproj/geod.py

@@ -14,6 +14,7 @@
     "geodesic_version_str",
     "GeodIntermediateFlag",
     "GeodIntermediateReturn",
+    "reverse_azimuth",
 ]
 
 import math
@@ -369,15 +370,15 @@ 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 third return value will be the back azimuth,
-            Otherwise, it will be the forward azimuth.
+            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) or Forward azimuth(s)
+            Back azimuth(s) or Forward azimuth(s) (azi21)
         scalar or array:
             Distance(s) between initial and terminus point(s)
             in meters
@@ -1116,38 +1117,47 @@ def __eq__(self, other: Any) -> bool:
         return self.__repr__() == other.__repr__()
 
 
-def reverse_azimuth(azi: float, radians: bool = False) -> float:
-    """
-    Reverses the given azimuth (forward <-> backwards)
+try:
+    import numpy as np  # type: ignore  # pylint: disable=import-outside-toplevel
 
-    Parameters
-    ----------
-    azi: float
-        azimuth
-    radians: bool, default=False
-        If True, the input data is assumed to be in radians.
-        Otherwise, the data is assumed to be in degrees.
+    def reverse_azimuth(azi: Any, radians: bool = False) -> Any:
+        """
+        Reverses the given azimuth (forward <-> backwards)
 
-    Returns
-    -------
-    float
-        the reversed azimuth (forward <-> backwards)
-    """
-    factor = math.pi if radians else 180.0
-    return azi - factor if azi > 0 else azi + factor
+        Parameters
+        ----------
+        azi: float or array, :class:`numpy.ndarray`
+            azimuth
+        radians: bool, default=False
+            If True, the input data is assumed to be in radians.
+            Otherwise, the data is assumed to be in degrees.
+
+        Returns
+        -------
+        float or array, :class:`numpy.ndarray`
+            the reversed azimuth (forward <-> backwards)
+        """
+        factor = math.pi if radians else 180.0
+        return azi - np.sign(azi) * factor
 
+except ImportError:
 
-def reverse_azimuth_arr(arr: Any, radians: bool = False):
-    """
-    Reverses the given azimuth (forward <-> backwards) inplace for the given array
+    def reverse_azimuth(azi: Any, radians: bool = False) -> Any:
+        """
+        Reverses the given azimuth (forward <-> backwards)
 
-    Parameters
-    ----------
-    arr: array, :class:`numpy.ndarray`
-        azimuth
-    radians: bool, default=False
-        If True, the input data is assumed to be in radians.
-        Otherwise, the data is assumed to be in degrees.
-    """
-    for idx, item in enumerate(arr):
-        arr[idx] = reverse_azimuth(item, radians=radians)
+        Parameters
+        ----------
+        azi: float
+            azimuth
+        radians: bool, default=False
+            If True, the input data is assumed to be in radians.
+            Otherwise, the data is assumed to be in degrees.
+
+        Returns
+        -------
+        float
+            the reversed azimuth (forward <-> backwards)
+        """
+        factor = math.pi if radians else 180.0
+        return azi - factor if azi > 0 else azi + factor

test/test_geod.py

@@ -4,15 +4,15 @@
 import pickle
 import shutil
 import tempfile
-from contextlib import contextmanager
+from contextlib import contextmanager, nullcontext
 from itertools import permutations
 
 import numpy as np
 import pytest
 from numpy.testing import assert_almost_equal
 
 from pyproj import Geod
-from pyproj.geod import GeodIntermediateFlag, reverse_azimuth_arr
+from pyproj.geod import GeodIntermediateFlag, reverse_azimuth
 
 try:
     from shapely.geometry import (
@@ -70,17 +70,27 @@ 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:")
-    endlon, endlat, backaz = gg.fwd(lon1pt, lat1pt, az12, dist)
-    assert_almost_equal(
-        (endlon, endlat, backaz), (lon2pt, lat2pt, true_az21), decimal=3
-    )
+    for return_back_azimuth in (True, False):
+        endlon, endlat, backaz = gg.fwd(
+            lon1pt, lat1pt, az12, dist, return_back_azimuth=return_back_azimuth
+        )
+        if not return_back_azimuth:
+            backaz = reverse_azimuth(backaz)
+        assert_almost_equal(
+            (endlon, endlat, backaz), (lon2pt, lat2pt, true_az21), decimal=3
+        )
 
     inc_exc = ["excluding", "including"]
     res_az12_az21_dists_all = [
@@ -150,33 +160,40 @@ def test_geod_inverse_transform():
             ).transpose()
 
     del_s = dist / (point_count - 1)
-    lons_a = np.empty(point_count)
-    lats_a = np.empty(point_count)
-    azis_a = np.empty(point_count)
-
-    print("test inv_intermediate (by npts) with azi output")
-    res = gg.inv_intermediate(
-        out_lons=lons_a,
-        out_lats=lats_a,
-        out_azis=azis_a,
-        lon1=lon1pt,
-        lat1=lat1pt,
-        lon2=lon2pt,
-        lat2=lat2pt,
-        npts=point_count,
-        initial_idx=0,
-        terminus_idx=0,
-        return_back_azimuth=False,
-    )
-    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)
-    assert_almost_equal(res.lats, lats)
-    assert_almost_equal(res.azis[:-1], azis12[1:])
-    assert res.lons is lons_a
-    assert res.lats is lats_a
-    assert res.azis is azis_a
+    for return_back_azimuth in [None, True, False]:
+        print(f"test inv_intermediate (by npts) with azi output {return_back_azimuth=}")
+        lons_a = np.empty(point_count)
+        lats_a = np.empty(point_count)
+        azis_a = np.empty(point_count)
+
+        with pytest.warns(
+            UserWarning
+        ) if return_back_azimuth is None else nullcontext():
+            res = gg.inv_intermediate(
+                out_lons=lons_a,
+                out_lats=lats_a,
+                out_azis=azis_a,
+                lon1=lon1pt,
+                lat1=lat1pt,
+                lon2=lon2pt,
+                lat2=lat2pt,
+                npts=point_count,
+                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)
+        assert_almost_equal(res.lats, lats)
+        out_azis = res.azis[:-1]
+        if return_back_azimuth in [True, None]:
+            out_azis = reverse_azimuth(out_azis)
+        assert_almost_equal(out_azis, azis12[1:])
+        assert res.lons is lons_a
+        assert res.lats is lats_a
+        assert res.azis is azis_a
 
     for flags in (GeodIntermediateFlag.AZIS_DISCARD, GeodIntermediateFlag.AZIS_KEEP):
         print("test inv_intermediate (by npts) without azi output, no buffers")
@@ -221,47 +238,50 @@ def test_geod_inverse_transform():
             return_back_azimuth=False,
         )
         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 res.lons is lons_b
         assert res.lats is lats_b
         if flags == GeodIntermediateFlag.AZIS_DISCARD:
             assert res.azis is None
         else:
             assert_almost_equal(res.azis, azis_a)
+        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)
 
-    print("test fwd_intermediate")
-    for return_back_azimuth in [False, True]:
+    for return_back_azimuth in [False, True, None]:
+        print(f"test fwd_intermediate with {return_back_azimuth=}")
         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,
-        )
+        with pytest.warns(
+            UserWarning
+        ) if return_back_azimuth is None else nullcontext():
+            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 res.lons is lons_b
+        assert res.lats is lats_b
+        assert res.azis is azis_b
         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_arr(azis_b)
+        if return_back_azimuth in [True, None]:
+            res.azis = reverse_azimuth(res.azis)
         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 (
@@ -702,3 +722,14 @@ def test_geod__build_kwargs(kwarg):
     assert_almost_equal(gg.b, gg2.b)
     assert_almost_equal(gg.f, gg2.f)
     assert_almost_equal(gg.es, gg2.es)
+
+
+@pytest.mark.parametrize("radians", [False, True])
+def test_geod__reverse_azimuth(radians):
+    f = math.pi / 180 if radians else 1
+    xx = np.array([10, 20, -10])
+    yy = np.array([10 - 180, 20 - 180, -10 + 180])
+    for x, y in zip(xx, yy):
+        assert_almost_equal(reverse_azimuth(x * f, radians=radians), y * f)
+
+    assert_almost_equal(reverse_azimuth(xx * f, radians=radians), yy * f)