From c99e74e4440f4d40cae9a409dc1ff8f3ea179453 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Sun, 12 May 2024 18:38:53 -0400 Subject: [PATCH] Numpy ver of length and angle in rad. Don't skip tests --- spherical_geometry/great_circle_arc.py | 33 +++++++++++--------------- spherical_geometry/polygon.py | 5 ++-- spherical_geometry/tests/test_basic.py | 16 +++++-------- spherical_geometry/tests/test_union.py | 7 ++++-- 4 files changed, 27 insertions(+), 34 deletions(-) diff --git a/spherical_geometry/great_circle_arc.py b/spherical_geometry/great_circle_arc.py index 7776e8c..469d6c1 100644 --- a/spherical_geometry/great_circle_arc.py +++ b/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: diff --git a/spherical_geometry/polygon.py b/spherical_geometry/polygon.py index cc3b452..9a070fd 100644 --- a/spherical_geometry/polygon.py +++ b/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) diff --git a/spherical_geometry/tests/test_basic.py b/spherical_geometry/tests/test_basic.py index 5349198..892580c 100644 --- a/spherical_geometry/tests/test_basic.py +++ b/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,7 +531,6 @@ 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 = [ @@ -542,7 +538,7 @@ def test_math_util_angle_nearly_coplanar_vec(): 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) diff --git a/spherical_geometry/tests/test_union.py b/spherical_geometry/tests/test_union.py index 7575daa..c66c7af 100644 --- a/spherical_geometry/tests/test_union.py +++ b/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 = []