Numpy ver of length and angle in rad. Don't skip tests

spherical_geometry/great_circle_arc.py

@@ -19,6 +19,7 @@
 # the python versions are a fallback if the C cannot be used
 try:
     from spherical_geometry import math_util
+    math_util = None
     HAS_C_UFUNCS = False
 except ImportError:
     HAS_C_UFUNCS = False
@@ -183,7 +184,7 @@ def intersection(A, B, C, D):
     return np.where(equals, np.nan, cross)
 
 
-def length(A, B, degrees=True):
+def length(A, B):
     r"""
     Returns the angular distance between two points (in vector space)
     on the unit sphere.
@@ -193,10 +194,6 @@ def length(A, B, degrees=True):
     A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
        The endpoints of the great circle arc, in vector space.
 
-    degrees : bool, optional
-        If `True` (default) the result is returned in decimal degrees,
-        otherwise radians.
-
     Returns
     -------
     length : scalar or array of scalars
@@ -213,7 +210,6 @@ def length(A, B, degrees=True):
     if HAS_C_UFUNCS:
         result = math_util.length(A, B)
     else:
-        raise AssertionError("C version should have been used")
         approx1 = 1 + 3 * np.finfo(float).eps
         A = np.asanyarray(A)
         B = np.asanyarray(B)
@@ -240,10 +236,7 @@ def length(A, B, degrees=True):
         with np.errstate(invalid='ignore'):
             result = np.arccos(dot)
 
-    if degrees:
-        return np.rad2deg(result)
-    else:
-        return result
+    return result
 
 
 def intersects(A, B, C, D):
@@ -300,10 +293,10 @@ def intersects_point(A, B, C):
 
     length_diff = np.abs((left_length + right_length) - total_length)
 
-    return length_diff < 1e-10
+    return length_diff < 3e-11
 
 
-def angle(A, B, C, degrees=True):
+def angle(A, B, C):
     """
     Returns the angle at *B* between *AB* and *BC*.
 
@@ -312,10 +305,6 @@ def angle(A, B, C, degrees=True):
     A, B, C : (*x*, *y*, *z*) triples or Nx3 arrays of triples
         Points on sphere.
 
-    degrees : bool, optional
-        If `True` (default) the result is returned in decimal degrees,
-        otherwise radians.
-
     Returns
     -------
     angle : float or array of floats
@@ -339,14 +328,20 @@ def angle(A, B, C, degrees=True):
         ABX = _cross_and_normalize(A, B)
         BCX = _cross_and_normalize(C, B)
         X = _cross_and_normalize(ABX, BCX)
+        m = np.logical_or(
+                np.linalg.norm(ABX, axis=-1) == 0.0,
+                np.linalg.norm(BCX, axis=-1) == 0.0
+            )
+
         diff = inner1d(B, X)
         inner = inner1d(ABX, BCX)
         with np.errstate(invalid='ignore'):
+            inner = np.clip(inner, -1.0, 1.0)  # needed due to accuracy loss
             angle = np.arccos(inner)
+
         angle = np.where(diff < 0.0, (2.0 * np.pi) - angle, angle)
 
-    if degrees:
-        angle = np.rad2deg(angle)
+        angle[m] = np.nan
 
     return angle
 
@@ -405,7 +400,7 @@ def interpolate(A, B, steps=50):
     steps = int(max(steps, 2))
     t = np.linspace(0.0, 1.0, steps, endpoint=True).reshape((steps, 1))
 
-    omega = length(A, B, degrees=False)
+    omega = length(A, B)
     if omega == 0.0:
         offsets = t
     else:

spherical_geometry/polygon.py

@@ -546,8 +546,7 @@ def area(self):
             return np.array(0.0)
 
         points = np.vstack((self._points, self._points[1]))
-        angles = great_circle_arc.angle(points[:-2], points[1:-1],
-                                        points[2:], degrees=False)
+        angles = great_circle_arc.angle(points[:-2], points[1:-1], points[2:])
 
         return np.sum(angles) - (len(angles) - 2) * np.pi
 
@@ -704,7 +703,7 @@ def draw(self, m, **plot_args):
 
         alpha = 1.0
         for A, B in zip(points[0:-1], points[1:]):
-            length = great_circle_arc.length(A, B, degrees=True)
+            length = np.rad2deg(great_circle_arc.length(A, B))
             if not np.isfinite(length):
                 length = 2
             interpolated = great_circle_arc.interpolate(A, B, length * 4)

spherical_geometry/tests/test_basic.py

@@ -161,8 +161,6 @@ def test_interpolate():
         assert abs(length - first_length) < 1.0e-10
 
 
-@pytest.mark.xfail(
-    math_util is None, reason="math_util C-ext is missing, numpy gives different results")
 def test_overlap():
     def build_polygon(offset):
         points = []
@@ -380,26 +378,26 @@ def test_great_circle_arc_length():
     B = [-90, 0]
     A = vector.lonlat_to_vector(*A)
     B = vector.lonlat_to_vector(*B)
-    assert great_circle_arc.length(A, B) == 180.0
+    assert_almost_equal(great_circle_arc.length(A, B), np.pi)
 
     A = [135, 0]
     B = [-90, 0]
     A = vector.lonlat_to_vector(*A)
     B = vector.lonlat_to_vector(*B)
-    assert_almost_equal(great_circle_arc.length(A, B), 135.0)
+    assert_almost_equal(great_circle_arc.length(A, B), (3.0 / 4.0) * np.pi)
 
     A = [0, 0]
     B = [0, 90]
     A = vector.lonlat_to_vector(*A)
     B = vector.lonlat_to_vector(*B)
-    assert_almost_equal(great_circle_arc.length(A, B), 90.0)
+    assert_almost_equal(great_circle_arc.length(A, B), np.pi / 2.0)
 
 
 def test_great_circle_arc_angle():
     A = [1, 0, 0]
     B = [0, 1, 0]
     C = [0, 0, 1]
-    assert great_circle_arc.angle(A, B, C) == 270.0
+    assert great_circle_arc.angle(A, B, C) == (3.0 / 2.0) * np.pi
 
     # TODO: More angle tests
 
@@ -522,10 +520,9 @@ def test_convex_hull(repeat_pts):
         assert b == r, "Polygon boundary has correct points"
 
 
-@pytest.mark.skipif(math_util is None, reason="math_util C-ext is missing")
 def test_math_util_angle_domain():
     assert not np.isfinite(
-        math_util.angle([[0, 0, 0]], [[0, 0, 0]], [[0, 0, 0]])[0]
+        great_circle_arc.angle([[0, 0, 0]], [[0, 0, 0]], [[0, 0, 0]])[0]
     )
 
 
@@ -534,15 +531,14 @@ def test_math_util_length_domain():
         great_circle_arc.length([[np.nan, 0, 0]], [[0, 0, np.inf]])
 
 
-@pytest.mark.skipif(math_util is None, reason="math_util C-ext is missing")
 def test_math_util_angle_nearly_coplanar_vec():
     # test from issue #222 + extra values
     vectors = [
         5 * [[1.0, 1.0, 1.0]],
         5 * [[1, 0.9999999, 1]],
         [[1, 0.5, 1], [1, 0.15, 1], [1, 0.001, 1], [1, 0.15, 1], [-1, 0.1, -1]]
     ]
-    angles = math_util.angle(*vectors)
+    angles = great_circle_arc.angle(*vectors)
 
     assert_allclose(angles[:-1], np.pi, rtol=0, atol=1e-16)
     assert_allclose(angles[-1], 0, rtol=0, atol=1e-32)

spherical_geometry/tests/test_union.py

@@ -18,6 +18,7 @@
 
 try:
     from spherical_geometry import math_util
+    math_util = None
 except ImportError:
     math_util = None
 
@@ -332,9 +333,11 @@ def test_edge_crossings():
 
 
 @pytest.mark.skipif(
-    math_util is None, reason="math_util C-ext is missing, double accuracy leads to crash")
+    math_util is None,
+    reason="math_util C-ext is missing, double accuracy leads to crash"
+)
 def test_almost_identical_polygons_multi_union():
-    filename = resolve_imagename(ROOT_DIR,'almost_same_polygons.npz')
+    filename = resolve_imagename(ROOT_DIR, 'almost_same_polygons.npz')
     polygon_data = np.load(filename)
 
     polygons = []