Make great circle metric more robust (#15)

This change makes the calculation of the great circle distance more robust against rounding issues.

CITATION.cff

@@ -5,6 +5,6 @@ authors:
   given-names: "Nis"
   orcid: "https://orcid.org/0000-0002-4712-9579"
 title: "GeoRDPy"
-version: 3.0
+version: 3.1
 date-released: 2023-10-06
 url: "https://github.com/avitase/geordpy"

geordpy/great_circle.py

@@ -34,8 +34,12 @@ def cos_distance_segment(lat, lon, *, lat1, lon1, lat2, lon2, eps=1e-5):
     )
 
     n = np.cross(a, b)
-    cos_gamma = np.dot(a, b)
-    n /= max(np.sqrt((1.0 - cos_gamma) * (1.0 + cos_gamma)), np.nextafter(0.0, 1.0))
+    cos_gamma = np.clip(np.dot(a, b), a_min=-1.0, a_max=1.0)
+
+    # if |cos(gamma)| = 1, then |n| = 0. However, in case |n| is slightly larger due
+    # to numeric issues, division by 1e-100 is a decent approximation. Even smaller
+    # values for the denominator should be avoided as this can trigger overflows.
+    n /= max(np.sqrt((1.0 - cos_gamma) * (1.0 + cos_gamma)), 1e-100)
 
     s = x @ n
 

tests/rdp_test.py

@@ -66,6 +66,24 @@ def test_great_circle():
     assert all(mask == [True, False, True, True, True, True, False, True])
 
 
+def test_degenerated_great_circle():
+    points = [
+        (0.0, 10.0),  # point 0
+        (0.0, 0.0),  # point 1
+        (0.0, 0.0),  # point 2
+        (0.0, 10.0),  # point 3
+    ]
+
+    radius = 1000
+    _filter = functools.partial(rdp_filter, radius=radius)
+
+    mask = _filter(points, threshold=1)
+    assert all(mask == [True, True, False, True])
+
+    mask = _filter(points, threshold=10 * radius)
+    assert all(mask == [True, False, False, True])
+
+
 def test_rhumb_line():
     def _gd(x):
         return np.rad2deg(utils.gd(np.deg2rad(x)))