From 2bb52b9088c19783e48b66a12858fb107ee55433 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sat, 14 Sep 2024 16:25:05 -0700 Subject: [PATCH 01/24] Initial PR Open --- test/test_arcs.py | 35 ++++++--------- uxarray/constants.py | 8 +++- uxarray/grid/arcs.py | 100 ++++++++++++++++++++++++++++--------------- 3 files changed, 85 insertions(+), 58 deletions(-) diff --git a/test/test_arcs.py b/test/test_arcs.py index 17d04fe57..057864c00 100644 --- a/test/test_arcs.py +++ b/test/test_arcs.py @@ -30,13 +30,6 @@ class TestIntersectionPoint(TestCase): def test_pt_within_gcr(self): # The GCR that's eexactly 180 degrees will have Value Error raised - gcr_180degree_cart = [ - _lonlat_rad_to_xyz(0.0, 0.0), - _lonlat_rad_to_xyz(np.pi, 0.0) - ] - pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) - with self.assertRaises(ValueError): - point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart)) gcr_180degree_cart = [ _lonlat_rad_to_xyz(0.0, np.pi / 2.0), @@ -45,7 +38,7 @@ def test_pt_within_gcr(self): pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) with self.assertRaises(ValueError): - point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart)) + point_within_gca(pt_same_lon_in, gcr_180degree_cart) # Test when the point and the GCA all have the same longitude gcr_same_lon_cart = [ @@ -53,19 +46,19 @@ def test_pt_within_gcr(self): _lonlat_rad_to_xyz(0.0, -1.5) ] pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0) - self.assertTrue(point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_same_lon_cart))) + self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart)) pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001) - res = point_within_gca(np.asarray(pt_same_lon_out), np.asarray(gcr_same_lon_cart)) - self.assertFalse(res) + res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart) + self.assertTrue(res) pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0) - res = point_within_gca(np.asarray(pt_same_lon_out_2), np.asarray(gcr_same_lon_cart)) + res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart) self.assertFalse(res) # And if we increase the digital place by one, it should be true again pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001) - res = point_within_gca(np.asarray(pt_same_lon_out_add_one_place), np.asarray(gcr_same_lon_cart)) + res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart) self.assertTrue(res) # Normal case @@ -76,7 +69,7 @@ def test_pt_within_gcr(self): -0.997]]) pt_cart_within = np.array( [0.25616109352676675, 0.9246590335292105, -0.010021496695000144]) - self.assertTrue(point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_2), True)) + self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True)) # Test other more complicate cases : The anti-meridian case @@ -87,16 +80,16 @@ def test_pt_within_gcr_antimeridian(self): gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]]) pt_cart = np.array( [0.9438777657502077, 0.1193199333436068, 0.922714737029319]) - self.assertTrue(point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True)) + self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True)) # If we swap the gcr, it should throw a value error since it's larger than 180 degree gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724, 0.593]]) with self.assertRaises(ValueError): - point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=True) + point_within_gca(pt_cart, gcr_cart_flip, is_directed=True) # If we flip the gcr in the undirected mode, it should still work self.assertTrue( - point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=False)) + point_within_gca(pt_cart, gcr_cart_flip, is_directed=False)) # 2nd anti-meridian case # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828], @@ -107,9 +100,9 @@ def test_pt_within_gcr_antimeridian(self): pt_cart_within = np.array( [0.6136726305712109, 0.28442243941920053, -0.365605190899831]) self.assertFalse( - point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=True)) + point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True)) self.assertFalse( - point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=False)) + point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False)) # The first case should not work and the second should work v1_rad = [0.1, 0.0] @@ -119,10 +112,10 @@ def test_pt_within_gcr_antimeridian(self): gcr_cart = np.array([v1_cart, v2_cart]) pt_cart = _lonlat_rad_to_xyz(0.01, 0.0) with self.assertRaises(ValueError): - point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True) + point_within_gca(pt_cart, gcr_cart, is_directed=True) gcr_car_flipped = np.array([v2_cart, v1_cart]) self.assertTrue( - point_within_gca(np.asarray(pt_cart), np.asarray(gcr_car_flipped), is_directed=True)) + point_within_gca(pt_cart, gcr_car_flipped, is_directed=True)) def test_pt_within_gcr_cross_pole(self): gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]]) diff --git a/uxarray/constants.py b/uxarray/constants.py index a8baa9c1c..2acc1bf93 100644 --- a/uxarray/constants.py +++ b/uxarray/constants.py @@ -10,9 +10,13 @@ # half of the working precision is already the most optimal value for the error tolerance, # more tailored values will be used in the future. -ERROR_TOLERANCE = 1.0e-8 +ERROR_TOLERANCE = np.float64(1.0e-8) -ENABLE_JIT_CACHE = True +# The below value is the machine epsilon for the float64 data type, it will be used in the most basic operations as a +# error tolerance, mainly in the intersection calculations. +MACHINE_EPSILON = np.finfo(float).eps + +ENABLE_JIT_CACHE = False ENABLE_JIT = True ENABLE_FMA = False diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 63344473e..83e2d94e7 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -3,13 +3,14 @@ # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place from uxarray.grid.coordinates import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar -from uxarray.constants import ERROR_TOLERANCE +from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON -from uxarray.utils.computing import isclose, cross, dot +from uxarray.utils.computing import isclose, cross, dot, allclose from numba import njit + def _to_list(obj): if not isinstance(obj, list): if isinstance(obj, np.ndarray): @@ -23,63 +24,92 @@ def _to_list(obj): @njit def _point_within_gca_body( - angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed + angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed ): angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) - if isclose(angle, np.pi, rtol=0.0, atol=ERROR_TOLERANCE): - raise ValueError( - "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " - "Consider breaking the Great Circle Arc" - "into two Great Circle Arcs" - ) + # See if the point is on the plane of the GCA, because we are dealing with floating point numbers with np.dot now + # just using the rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON, but consider using the more proper error tolerance + # in the future + if not allclose( + dot(np.cross(gca_cart[0], gca_cart[1]), pt), + 0, + rtol=MACHINE_EPSILON, + atol=MACHINE_EPSILON, + ): + return False if not isclose( - dot(cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1])), pt), - 0, - rtol=0.0, - atol=ERROR_TOLERANCE, + dot(cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1])), pt), + 0, + rtol=MACHINE_EPSILON, + atol=MACHINE_EPSILON, ): return False - if isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE): - # If the pt and the GCA are on the same longitude (the y coordinates are the same) - if isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE): - # Now use the latitude to determine if the pt falls between the interval - return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) - else: - # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA - return False - # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point - if isclose( - abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0.0, atol=ERROR_TOLERANCE + # Or if one of the endpoints is on the pole point, then the GCA goes through the pole point + if ( + isclose( + abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0.0, atol=MACHINE_EPSILON + ) + or isclose( + abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + ) + or isclose( + abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + ) ): # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] - if isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0.0, atol=ERROR_TOLERANCE): + # Since the point is our calculated properly, we use the atol=ERROR_TOLERANCE and rtol=ERROR_TOLERANCE + if isclose( + abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + ): pt_lonlat[0] = GCRv0_lonlat[0] + + # Special case, if one of the GCA endpoints is on the pole point, and another endpoint is not + # then we need to check if the pt is on the GCA + if isclose( + abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0 + ) or isclose(abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0): + # Identify the non-pole endpoint + non_pole_endpoint = None + if not isclose( + abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0 + ): + non_pole_endpoint = GCRv0_lonlat + elif not isclose( + abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0 + ): + non_pole_endpoint = GCRv1_lonlat + + if non_pole_endpoint is not None and not np.isclose( + non_pole_endpoint[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0 + ): + return False + if not isclose( - GCRv0_lonlat[0], pt_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE + GCRv0_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0 ) and not isclose( - GCRv1_lonlat[0], pt_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE + GCRv1_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0 ): return False else: # Determine the pole latitude and latitude extension if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or ( - GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 + GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 ): pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2 lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) - + np.pi / 2 - + abs(GCRv1_lonlat[1]) + abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + + np.pi / 2 + + abs(GCRv1_lonlat[1]) ) else: pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) - + np.pi / 2 - + abs(GCRv1_lonlat[1]) + abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + + np.pi / 2 + + abs(GCRv1_lonlat[1]) ) if is_directed and lat_extend >= np.pi: @@ -309,7 +339,7 @@ def extreme_gca_latitude(gca_cart, extreme_type): d_a_max = ( np.clip(d_a_max, 0, 1) if isclose(d_a_max, 0, atol=ERROR_TOLERANCE) - or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) + or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) else d_a_max ) From d706ce5319d17afc46cea6f85abdfa3f324a3acd Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sat, 14 Sep 2024 16:44:11 -0700 Subject: [PATCH 02/24] Finish the change of `pt_within_gcr` --- uxarray/constants.py | 2 +- uxarray/grid/arcs.py | 87 +++++++++++++++++++++++--------------------- 2 files changed, 46 insertions(+), 43 deletions(-) diff --git a/uxarray/constants.py b/uxarray/constants.py index 2acc1bf93..606b6ab2c 100644 --- a/uxarray/constants.py +++ b/uxarray/constants.py @@ -14,7 +14,7 @@ # The below value is the machine epsilon for the float64 data type, it will be used in the most basic operations as a # error tolerance, mainly in the intersection calculations. -MACHINE_EPSILON = np.finfo(float).eps +MACHINE_EPSILON = np.float64(np.finfo(float).eps) ENABLE_JIT_CACHE = False ENABLE_JIT = True diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 83e2d94e7..15640d249 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -10,7 +10,6 @@ from numba import njit - def _to_list(obj): if not isinstance(obj, list): if isinstance(obj, np.ndarray): @@ -24,92 +23,96 @@ def _to_list(obj): @njit def _point_within_gca_body( - angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed + angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed ): angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1]) + if allclose(angle, np.pi, rtol=0.0, atol=MACHINE_EPSILON): + raise ValueError( + "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. " + "Consider breaking the Great Circle Arc" + "into two Great Circle Arcs" + ) + # See if the point is on the plane of the GCA, because we are dealing with floating point numbers with np.dot now # just using the rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON, but consider using the more proper error tolerance # in the future - if not allclose( - dot(np.cross(gca_cart[0], gca_cart[1]), pt), - 0, - rtol=MACHINE_EPSILON, - atol=MACHINE_EPSILON, - ): - return False + cross_product = cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1])) - if not isclose( - dot(cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1])), pt), - 0, - rtol=MACHINE_EPSILON, - atol=MACHINE_EPSILON, + if not allclose( + dot(np.asarray(cross_product), np.asarray(pt)), # Custom dot function + 0, + rtol=MACHINE_EPSILON, + atol=MACHINE_EPSILON, ): return False # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point # Or if one of the endpoints is on the pole point, then the GCA goes through the pole point if ( - isclose( - abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0.0, atol=MACHINE_EPSILON - ) - or isclose( - abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE - ) - or isclose( - abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE - ) + isclose( + abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), + np.pi, + rtol=0.0, + atol=MACHINE_EPSILON, + ) + or isclose( + abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + ) + or isclose( + abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + ) ): # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0] # Since the point is our calculated properly, we use the atol=ERROR_TOLERANCE and rtol=ERROR_TOLERANCE if isclose( - abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE + abs(pt_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE ): pt_lonlat[0] = GCRv0_lonlat[0] # Special case, if one of the GCA endpoints is on the pole point, and another endpoint is not # then we need to check if the pt is on the GCA if isclose( - abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0 - ) or isclose(abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0): + abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 + ) or isclose(abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0): # Identify the non-pole endpoint non_pole_endpoint = None if not isclose( - abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0 + abs(GCRv0_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 ): non_pole_endpoint = GCRv0_lonlat elif not isclose( - abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0 + abs(GCRv1_lonlat[1]), np.pi / 2, rtol=ERROR_TOLERANCE, atol=0.0 ): non_pole_endpoint = GCRv1_lonlat if non_pole_endpoint is not None and not np.isclose( - non_pole_endpoint[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0 + non_pole_endpoint[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 ): return False if not isclose( - GCRv0_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0 + GCRv0_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 ) and not isclose( - GCRv1_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0 + GCRv1_lonlat[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 ): return False else: # Determine the pole latitude and latitude extension - if (GCRv0_lonlat[1] > 0 and GCRv1_lonlat[1] > 0) or ( - GCRv0_lonlat[1] < 0 and GCRv1_lonlat[1] < 0 + if (GCRv0_lonlat[1] > 0.0 and GCRv1_lonlat[1] > 0.0) or ( + GCRv0_lonlat[1] < 0.0 and GCRv1_lonlat[1] < 0.0 ): - pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0 else -np.pi / 2 + pole_lat = np.pi / 2 if GCRv0_lonlat[1] > 0.0 else -np.pi / 2 lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) - + np.pi / 2 - + abs(GCRv1_lonlat[1]) + abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + + np.pi / 2 + + abs(GCRv1_lonlat[1]) ) else: pole_lat = _decide_pole_latitude(GCRv0_lonlat[1], GCRv1_lonlat[1]) lat_extend = ( - abs(np.pi / 2 - abs(GCRv0_lonlat[1])) - + np.pi / 2 - + abs(GCRv1_lonlat[1]) + abs(np.pi / 2 - abs(GCRv0_lonlat[1])) + + np.pi / 2 + + abs(GCRv1_lonlat[1]) ) if is_directed and lat_extend >= np.pi: @@ -134,7 +137,7 @@ def _point_within_gca_body( # The necessary condition: the pt longitude is on the opposite side of the anti-meridian # Case 2: The anti-meridian case where 180 -->x0 --> 0 lon --> x1 --> 180 lon - elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0: + elif 2 * np.pi > GCRv0_lonlat[0] > np.pi > GCRv1_lonlat[0] > 0.0: return in_between( GCRv0_lonlat[0], pt_lonlat[0], 2 * np.pi ) or in_between(0, pt_lonlat[0], GCRv1_lonlat[0]) @@ -339,7 +342,7 @@ def extreme_gca_latitude(gca_cart, extreme_type): d_a_max = ( np.clip(d_a_max, 0, 1) if isclose(d_a_max, 0, atol=ERROR_TOLERANCE) - or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) + or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) else d_a_max ) From 2ce675755c6251d8a7fc25bdba9056c6870c083e Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sat, 14 Sep 2024 17:13:13 -0700 Subject: [PATCH 03/24] Finish the change of `_pole_point_inside_polygon` --- test/test_geometry.py | 3 ++ uxarray/grid/geometry.py | 103 ++++++++++++++++++++++++++++++++++----- 2 files changed, 94 insertions(+), 12 deletions(-) diff --git a/test/test_geometry.py b/test/test_geometry.py index 953942293..6041109c8 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -56,6 +56,8 @@ def test_linecollection_execution(self): lines = uxgrid.to_linecollection() + + class TestPredicate(TestCase): def test_pole_point_inside_polygon_from_vertice_north(self): @@ -161,6 +163,7 @@ def test_pole_point_inside_polygon_from_vertice_cross(self): self.assertTrue(result, "North pole should be inside the polygon") + class TestLatlonBoundUtils(TestCase): def _max_latitude_rad_iterative(self, gca_cart): diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 08229d1a0..74a6eaa3e 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -24,6 +24,57 @@ REFERENCE_POINT_EQUATOR = np.array([1.0, 0.0, 0.0]) +def _unique_points(points, tolerance=ERROR_TOLERANCE): + """Identify unique intersection points from a list of points, considering + floating point precision errors. + + Parameters + ---------- + points : list of array-like + A list containing the intersection points, where each point is an array-like structure (e.g., list, tuple, or numpy array) containing x, y, and z coordinates. + tolerance : float, optional + The distance threshold within which two points are considered identical. Default is 1e-6. + + Returns + ------- + list of np.ndarray + A list of unique snapped points in Cartesian coordinates. + + Notes + ----- + Given the nature of the mathematical equations and the spherical surface, it is more reasonable to calculate the "error radius" of the results using the following formula. + In the equation below, \(\tilde{P}_x\) and \(\tilde{P}_y\) are the calculated results, and \(P_x\) and \(P_y\) are the actual intersection points for the \(x\) and \(y\) coordinates, respectively. + The \(z\) coordinate is always the input \(z_0\), the constant latitude, so it is not included in this error calculation. + + .. math:: + \begin{aligned} + &\frac{\sqrt{(\tilde{v}_x - v_x)^2 + (\tilde{v}_y - v_y)^2 + (\tilde{v}_z - v_z)^2}}{\sqrt{v_x^2 + v_y^2 + v_z^2}}\\ + &= \sqrt{(\tilde{v}_x - v_x)^2 + (\tilde{v}_y - v_y)^2 + (\tilde{v}_z - v_z)^2} + \end{aligned} + + This method ensures that small numerical inaccuracies do not lead to multiple close points being considered different. + """ + unique_points = [] + points = [np.array(point) for point in points] # Ensure all points are numpy arrays + + def error_radius(p1, p2): + """Calculate the error radius between two points in 3D space.""" + numerator = np.sqrt( + (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2 + (p1[2] - p2[2]) ** 2 + ) + denominator = np.sqrt(p2[0] ** 2 + p2[1] ** 2 + p2[2] ** 2) + return numerator / denominator + + for point in points: + if not any( + error_radius(point, unique_point) < tolerance + for unique_point in unique_points + ): + unique_points.append(point) + + return unique_points + + # General Helpers for Polygon Viz # ---------------------------------------------------------------------------------------------------------------------- @@ -433,17 +484,27 @@ def _pole_point_inside_polygon(pole, face_edge_cart): ref_edge_south = np.array([-pole_point, REFERENCE_POINT_EQUATOR]) north_edges = face_edge_cart[np.any(face_edge_cart[:, :, 2] > 0, axis=1)] - south_edges = face_edge_cart[np.any(face_edge_cart[:, :, 2] < 0, axis=1)] - + south_edges = face_edge_cart[ + ~np.isin(face_edge_cart, north_edges).all(axis=(1, 2)) + ] + # The equator one is assigned to the south edges return ( _check_intersection(ref_edge_north, north_edges) + _check_intersection(ref_edge_south, south_edges) ) % 2 != 0 + elif ( + location == "North" + and pole == "South" + or location == "South" + and pole == "North" + ): + return False else: - warnings.warn( - "The given face should not contain both pole points.", UserWarning + raise ValueError( + "Invalid pole point query. Current location: {}, query pole point: {}".format( + location, pole + ) ) - return False def _check_intersection(ref_edge, edges): @@ -463,17 +524,34 @@ def _check_intersection(ref_edge, edges): Count of intersections. """ pole_point, ref_point = ref_edge - intersection_count = 0 + intersection_points = [] for edge in edges: intersection_point = gca_gca_intersection(ref_edge, edge) if intersection_point.size != 0: - if allclose(intersection_point, pole_point, atol=ERROR_TOLERANCE): - return True - intersection_count += 1 + # Handle both single point and multiple points case + if intersection_point.ndim == 1: # If there's only one point, make it a 2D array + intersection_point = [intersection_point] # Convert to list of points + + # for each intersection point, check if it is a pole point + for point in intersection_point: + if np.allclose(point, pole_point, atol=ERROR_TOLERANCE): + return True + intersection_points.append(point) - return intersection_count + # Only return the unique intersection points, the unique tolerance is set to ERROR_TOLERANCE + intersection_points = _unique_points(intersection_points, tolerance=ERROR_TOLERANCE) + + # If the unique intersection point is one and it is exactly one of the nodes of the face, return 0 + if len(intersection_points) == 1: + for edge in edges: + if np.allclose( + intersection_points[0], edge[0], atol=ERROR_TOLERANCE + ) or np.allclose(intersection_points[0], edge[1], atol=ERROR_TOLERANCE): + return 0 + + return len(intersection_points) def _classify_polygon_location(face_edge_cart): @@ -490,14 +568,15 @@ def _classify_polygon_location(face_edge_cart): Returns either 'North', 'South' or 'Equator' based on the polygon's location. """ z_coords = face_edge_cart[:, :, 2] - if all(z_coords > 0): + if np.all(z_coords > 0): return "North" - elif all(z_coords < 0): + elif np.all(z_coords < 0): return "South" else: return "Equator" + def _get_latlonbox_width(latlonbox_rad): """Calculate the width of a latitude-longitude box in radians. The box should be represented by a 2x2 array in radians and lon0 represent the From 6ffb63273cd275b6fe79d8887d6f4f24208125d8 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sat, 14 Sep 2024 17:52:35 -0700 Subject: [PATCH 04/24] Finish the change of `gca_constlat_intersection` --- test/test_arcs.py | 2 +- test/test_integrate.py | 418 ++++++++++++++++++++++++++++++++-- test/test_intersections.py | 24 +- uxarray/grid/arcs.py | 13 ++ uxarray/grid/geometry.py | 7 +- uxarray/grid/integrate.py | 270 +++++++++++++++++----- uxarray/grid/intersections.py | 160 ++++++++----- uxarray/grid/utils.py | 49 ++-- 8 files changed, 781 insertions(+), 162 deletions(-) diff --git a/test/test_arcs.py b/test/test_arcs.py index 057864c00..fa6bc8499 100644 --- a/test/test_arcs.py +++ b/test/test_arcs.py @@ -50,7 +50,7 @@ def test_pt_within_gcr(self): pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001) res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart) - self.assertTrue(res) + self.assertFalse(res) pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0) res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart) diff --git a/test/test_integrate.py b/test/test_integrate.py index 078f355cf..1b001818d 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -60,9 +60,169 @@ def test_multi_dim(self): class TestFaceWeights(TestCase): + gridfile_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" + dsfile_var2_ne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_var2.nc" + + def test_get_faces_constLat_intersection_info_one_intersection(self): + face_edges_cart = np.array([ + [[-5.4411371445381629e-01, -4.3910468172333759e-02, -8.3786164521844386e-01], + [-5.4463903501502697e-01, -6.6699045092185599e-17, -8.3867056794542405e-01]], + + [[-5.4463903501502697e-01, -6.6699045092185599e-17, -8.3867056794542405e-01], + [-4.9999999999999994e-01, -6.1232339957367648e-17, -8.6602540378443871e-01]], + + [[-4.9999999999999994e-01, -6.1232339957367648e-17, -8.6602540378443871e-01], + [-4.9948581138450826e-01, -4.5339793804534498e-02, -8.6513480297773349e-01]], + + [[-4.9948581138450826e-01, -4.5339793804534498e-02, -8.6513480297773349e-01], + [-5.4411371445381629e-01, -4.3910468172333759e-02, -8.3786164521844386e-01]] + ]) + + latitude_cart = -0.8660254037844386 + is_directed=False + is_latlonface=False + is_GCA_list=None + unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed) + # The expected unique_intersections length is 1 + self.assertEqual(len(unique_intersections), 1) + + def test_get_faces_constLat_intersection_info_encompass_pole(self): + face_edges_cart = np.array([ + [[0.03982285692494229, 0.00351700770436231, 0.9992005658140627], + [0.00896106681877875, 0.03896060263227105, 0.9992005658144913]], + + [[0.00896106681877875, 0.03896060263227105, 0.9992005658144913], + [-0.03428461218295055, 0.02056197086916728, 0.9992005658132106]], + + [[-0.03428461218295055, 0.02056197086916728, 0.9992005658132106], + [-0.03015012448894485, -0.02625260499902213, 0.9992005658145248]], + + [[-0.03015012448894485, -0.02625260499902213, 0.9992005658145248], + [0.01565081128889155, -0.03678697293262131, 0.9992005658167203]], + + [[0.01565081128889155, -0.03678697293262131, 0.9992005658167203], + [0.03982285692494229, 0.00351700770436231, 0.9992005658140627]] + ]) + + latitude_cart = 0.9993908270190958 + # Convert the latitude to degrees + latitude_rad = np.arcsin(latitude_cart) + latitude_deg = np.rad2deg(latitude_rad) + print(latitude_deg) + is_directed=False + is_latlonface=False + is_GCA_list=None + unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed) + # The expected unique_intersections length should be no greater than 2* n_edges + self.assertLessEqual(len(unique_intersections), 2*len(face_edges_cart)) + + def test_get_faces_constLat_intersection_info_on_pole(self): + face_edges_cart = np.array([ + [[-5.2264427688714095e-02, -5.2264427688714102e-02, -9.9726468863423734e-01], + [-5.2335956242942412e-02, -6.4093061293235361e-18, -9.9862953475457394e-01]], + + [[-5.2335956242942412e-02, -6.4093061293235361e-18, -9.9862953475457394e-01], + [6.1232339957367660e-17, 0.0000000000000000e+00, -1.0000000000000000e+00]], + + [[6.1232339957367660e-17, 0.0000000000000000e+00, -1.0000000000000000e+00], + [3.2046530646617680e-18, -5.2335956242942412e-02, -9.9862953475457394e-01]], + + [[3.2046530646617680e-18, -5.2335956242942412e-02, -9.9862953475457394e-01], + [-5.2264427688714095e-02, -5.2264427688714102e-02, -9.9726468863423734e-01]] + ]) + latitude_cart = -0.9998476951563913 + is_directed=False + is_latlonface=False + is_GCA_list=None + unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed) + # The expected unique_intersections length is 2 + self.assertEqual(len(unique_intersections), 2) + + + def test_get_faces_constLat_intersection_info_near_pole(self): + face_edges_cart = np.array([ + [[-5.1693346290592648e-02, 1.5622531297347531e-01, -9.8636780641686628e-01], + [-5.1195320928843470e-02, 2.0763904784932552e-01, -9.7686491641537532e-01]], + [[-5.1195320928843470e-02, 2.0763904784932552e-01, -9.7686491641537532e-01], + [1.2730919333264125e-17, 2.0791169081775882e-01, -9.7814760073380580e-01]], + [[1.2730919333264125e-17, 2.0791169081775882e-01, -9.7814760073380580e-01], + [9.5788483443923397e-18, 1.5643446504023048e-01, -9.8768834059513777e-01]], + [[9.5788483443923397e-18, 1.5643446504023048e-01, -9.8768834059513777e-01], + [-5.1693346290592648e-02, 1.5622531297347531e-01, -9.8636780641686628e-01]] + ]) + + latitude_cart = -0.9876883405951378 + latitude_rad = np.arcsin(latitude_cart) + latitude_deg = np.rad2deg(latitude_rad) + is_directed=False + is_latlonface=False + is_GCA_list=None + unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed) + # The expected unique_intersections length is 2 + self.assertEqual(len(unique_intersections), 1) + + + def test_get_faces_constLat_intersection_info_2(self): + """This might test the case where the calculated intersection points + might suffer from floating point errors If not handled properly, the + function might return more than 2 unique intersections.""" + + face_edges_cart = np.array([[[0.6546536707079771, -0.37796447300922714, -0.6546536707079772], + [0.6652465971273088, -0.33896007142593115, -0.6652465971273087]], + + [[0.6652465971273088, -0.33896007142593115, -0.6652465971273087], + [0.6949903639307233, -0.3541152775760984, -0.6257721344312508]], + + [[0.6949903639307233, -0.3541152775760984, -0.6257721344312508], + [0.6829382762718700, -0.39429459764546304, -0.6149203859609873]], + + [[0.6829382762718700, -0.39429459764546304, -0.6149203859609873], + [0.6546536707079771, -0.37796447300922714, -0.6546536707079772]]]) + + latitude_cart = -0.6560590289905073 + is_directed=False + is_latlonface=False + is_GCA_list=None + unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed) + # The expected unique_intersections length is 2 + self.assertEqual(len(unique_intersections), 2) + + + def test_get_faces_constLat_intersection_info_2(self): + """This might test the case where the calculated intersection points + might suffer from floating point errors If not handled properly, the + function might return more than 2 unique intersections.""" + + face_edges_cart = np.array([[[0.6546536707079771, -0.37796447300922714, -0.6546536707079772], + [0.6652465971273088, -0.33896007142593115, -0.6652465971273087]], + + [[0.6652465971273088, -0.33896007142593115, -0.6652465971273087], + [0.6949903639307233, -0.3541152775760984, -0.6257721344312508]], + + [[0.6949903639307233, -0.3541152775760984, -0.6257721344312508], + [0.6829382762718700, -0.39429459764546304, -0.6149203859609873]], + + [[0.6829382762718700, -0.39429459764546304, -0.6149203859609873], + [0.6546536707079771, -0.37796447300922714, -0.6546536707079772]]]) + + latitude_cart = -0.6560590289905073 + is_directed=False + is_latlonface=False + is_GCA_list=None + unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed) + # The expected unique_intersections length is 2 + self.assertEqual(len(unique_intersections), 2) def test_get_zonal_face_interval(self): - """Test that the zonal face weights are correct.""" + """Test the _get_zonal_face_interval function for correct interval + computation. + + This test verifies that the _get_zonal_face_interval function + accurately computes the zonal face intervals given a set of face + edge nodes and a constant latitude value (constZ). The expected + intervals are compared against the calculated intervals to + ensure correctness. + """ vertices_lonlat = [[1.6 * np.pi, 0.25 * np.pi], [1.6 * np.pi, -0.25 * np.pi], [0.4 * np.pi, -0.25 * np.pi], @@ -94,8 +254,117 @@ def test_get_zonal_face_interval(self): # Asserting almost equal arrays nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + + def test_get_zonal_face_interval_empty_interval(self): + """Test the _get_zonal_face_interval function for cases where the + interval is empty. + + This test verifies that the _get_zonal_face_interval function + correctly returns an empty interval when the latitude only + touches the face but does not intersect it. + """ + face_edges_cart = np.array([ + [[-5.4411371445381629e-01, -4.3910468172333759e-02, -8.3786164521844386e-01], + [-5.4463903501502697e-01, -6.6699045092185599e-17, -8.3867056794542405e-01]], + + [[-5.4463903501502697e-01, -6.6699045092185599e-17, -8.3867056794542405e-01], + [-4.9999999999999994e-01, -6.1232339957367648e-17, -8.6602540378443871e-01]], + + [[-4.9999999999999994e-01, -6.1232339957367648e-17, -8.6602540378443871e-01], + [-4.9948581138450826e-01, -4.5339793804534498e-02, -8.6513480297773349e-01]], + + [[-4.9948581138450826e-01, -4.5339793804534498e-02, -8.6513480297773349e-01], + [-5.4411371445381629e-01, -4.3910468172333759e-02, -8.3786164521844386e-01]] + ]) + + latitude_cart = -0.8660254037844386 + face_latlon_bounds = np.array([ + [-1.04719755, -0.99335412], + [3.14159265, 3.2321175] + ]) + + res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds, is_directed=False) + expected_res = pd.DataFrame({"start": [0.0], "end": [0.0]}) + pd.testing.assert_frame_equal(res, expected_res) + + def test_get_zonal_face_interval_encompass_pole(self): + """Test the _get_zonal_face_interval function for cases where the face + encompasses the pole inside.""" + face_edges_cart = np.array([ + [[0.03982285692494229, 0.00351700770436231, 0.9992005658140627], + [0.00896106681877875, 0.03896060263227105, 0.9992005658144913]], + + [[0.00896106681877875, 0.03896060263227105, 0.9992005658144913], + [-0.03428461218295055, 0.02056197086916728, 0.9992005658132106]], + + [[-0.03428461218295055, 0.02056197086916728, 0.9992005658132106], + [-0.03015012448894485, -0.02625260499902213, 0.9992005658145248]], + + [[-0.03015012448894485, -0.02625260499902213, 0.9992005658145248], + [0.01565081128889155, -0.03678697293262131, 0.9992005658167203]], + + [[0.01565081128889155, -0.03678697293262131, 0.9992005658167203], + [0.03982285692494229, 0.00351700770436231, 0.9992005658140627]] + ]) + + latitude_cart = 0.9993908270190958 + + face_latlon_bounds = np.array([ + [np.arcsin(0.9992005658145248), 0.5*np.pi], + [0, 2*np.pi] + ]) + # Expected result DataFrame + expected_df = pd.DataFrame({ + 'start': [0.000000, 1.101091, 2.357728, 3.614365, 4.871002, 6.127640], + 'end': [0.331721, 1.588358, 2.844995, 4.101632, 5.358270, 6.283185] + }) + + # Call the function to get the result + res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds, is_directed=False) + + # Assert the result matches the expected DataFrame + pd.testing.assert_frame_equal(res, expected_df) + + + + def test_get_zonal_face_interval_FILL_VALUE(self): + """Test the _get_zonal_face_interval function for cases where there are + dummy nodes.""" + dummy_node = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] + vertices_lonlat = [[1.6 * np.pi, 0.25 * np.pi], + [1.6 * np.pi, -0.25 * np.pi], + [0.4 * np.pi, -0.25 * np.pi], + [0.4 * np.pi, 0.25 * np.pi]] + vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat] + + face_edge_nodes = np.array([[vertices[0], vertices[1]], + [vertices[1], vertices[2]], + [vertices[2], vertices[3]], + [vertices[3], vertices[0]], + [dummy_node,dummy_node]]) + + constZ = np.sin(0.20) + # The latlon bounds for the latitude is not necessarily correct below since we don't use the latitudes bound anyway + interval_df = _get_zonal_face_interval(face_edge_nodes, constZ, + np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi, + 0.4 * np.pi]]), + is_directed=False) + expected_interval_df = pd.DataFrame({ + 'start': [1.6 * np.pi, 0.0], + 'end': [2.0 * np.pi, 00.4 * np.pi] + }) + # Sort both DataFrames by 'start' column before comparison + expected_interval_df_sorted = expected_interval_df.sort_values(by='start').reset_index(drop=True) + + # Converting the sorted DataFrames to NumPy arrays + actual_values_sorted = interval_df[['start', 'end']].to_numpy() + expected_values_sorted = expected_interval_df_sorted[['start', 'end']].to_numpy() + + # Asserting almost equal arrays + nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + + def test_get_zonal_face_interval_GCA_constLat(self): - """Test that the zonal face weights are correct.""" vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, -0.25 * np.pi], [0.4 * np.pi, -0.25 * np.pi], @@ -127,8 +396,10 @@ def test_get_zonal_face_interval_GCA_constLat(self): # Asserting almost equal arrays nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + def test_get_zonal_face_interval_equator(self): - """Test that the zonal face weights are correct.""" + """Test that the face interval is correctly computed when the latitude + is at the equator.""" vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, 0.0], [0.4 * np.pi, 0.0], [0.4 * np.pi, 0.25 * np.pi]] @@ -177,8 +448,9 @@ def test_get_zonal_face_interval_equator(self): # Asserting almost equal arrays nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) - def test_process_overlapped_intervals(self): - # Example data that has overlapping intervals and gap + + def test_process_overlapped_intervals_overlap_and_gap(self): + # Test intervals data that has overlapping intervals and gap intervals_data = [ { 'start': 0.0, @@ -226,6 +498,7 @@ def test_process_overlapped_intervals(self): nt.assert_array_equal(overlap_contributions, expected_overlap_contributions) + def test_process_overlapped_intervals_antimerdian(self): intervals_data = [ { @@ -273,6 +546,7 @@ def test_process_overlapped_intervals_antimerdian(self): nt.assert_array_equal(overlap_contributions, expected_overlap_contributions) + def test_get_zonal_faces_weight_at_constLat_equator(self): face_0 = [[1.7 * np.pi, 0.25 * np.pi], [1.7 * np.pi, 0.0], [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]] @@ -329,13 +603,11 @@ def test_get_zonal_faces_weight_at_constLat_equator(self): weight_df = _get_zonal_faces_weight_at_constLat(np.array([ face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes, face_3_edge_nodes - ]), - 0.0, - latlon_bounds, - is_directed=False) + ]), 0.0, latlon_bounds, is_directed=False) nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3) + def test_get_zonal_faces_weight_at_constLat_regular(self): face_0 = [[1.7 * np.pi, 0.25 * np.pi], [1.7 * np.pi, 0.0], [0.3 * np.pi, 0.0], [0.3 * np.pi, 0.25 * np.pi]] @@ -388,19 +660,123 @@ def test_get_zonal_faces_weight_at_constLat_regular(self): 'weight': [0.375, 0.0625, 0.3125, 0.25] }) - - # Assert the results is the same to the 3 decimal places weight_df = _get_zonal_faces_weight_at_constLat(np.array([ face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes, face_3_edge_nodes - ]), - np.sin(0.1 * np.pi), - latlon_bounds, - is_directed=False) + ]), np.sin(0.1 * np.pi), latlon_bounds, is_directed=False) nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3) + def test_get_zonal_faces_weight_at_constLat_on_pole_one_face(self): + #The face is touching the poleļ¼Œ so the weight should be 1.0 since there's only 1 face + face_edges_cart = np.array([[ + [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01], + [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]], + + [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01], + [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]], + + [[6.12323400e-17, 0.00000000e+00, -1.00000000e+00], + [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]], + + [[3.20465306e-18, -5.23359562e-02, -9.98629535e-01], + [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]] + ]]) + + # Corrected face_bounds + face_bounds = np.array([ + [-1.57079633, -1.4968158], + [3.14159265, 0.] + ]) + constLat_cart = -1 + + weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds, is_directed=False) + # Define the expected DataFrame + expected_weight_df = pd.DataFrame({"face_index": [0], "weight": [1.0]}) + + # Assert that the resulting should have weight is 1.0 + pd.testing.assert_frame_equal(weight_df, expected_weight_df) + + + def test_get_zonal_faces_weight_at_constLat_on_pole_faces(self): + #there will be 4 faces touching the pole, so the weight should be 0.25 for each face + face_edges_cart = np.array([ + [ + [[5.22644277e-02, -5.22644277e-02, 9.97264689e-01], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]], + [[5.23359562e-02, 0.00000000e+00, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]], + [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]], + [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [5.22644277e-02, -5.22644277e-02, 9.97264689e-01]] + ], + [ + [[5.23359562e-02, 0.00000000e+00, 9.98629535e-01], [5.22644277e-02, 5.22644277e-02, 9.97264689e-01]], + [[5.22644277e-02, 5.22644277e-02, 9.97264689e-01], [3.20465306e-18, 5.23359562e-02, 9.98629535e-01]], + [[3.20465306e-18, 5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]], + [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [5.23359562e-02, 0.00000000e+00, 9.98629535e-01]] + ], + [ + [[3.20465306e-18, -5.23359562e-02, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]], + [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [-5.23359562e-02, -6.40930613e-18, 9.98629535e-01]], + [[-5.23359562e-02, -6.40930613e-18, 9.98629535e-01], + [-5.22644277e-02, -5.22644277e-02, 9.97264689e-01]], + [[-5.22644277e-02, -5.22644277e-02, 9.97264689e-01], [3.20465306e-18, -5.23359562e-02, 9.98629535e-01]] + ], + [ + [[6.12323400e-17, 0.00000000e+00, 1.00000000e+00], [3.20465306e-18, 5.23359562e-02, 9.98629535e-01]], + [[3.20465306e-18, 5.23359562e-02, 9.98629535e-01], [-5.22644277e-02, 5.22644277e-02, 9.97264689e-01]], + [[-5.22644277e-02, 5.22644277e-02, 9.97264689e-01], [-5.23359562e-02, -6.40930613e-18, 9.98629535e-01]], + [[-5.23359562e-02, -6.40930613e-18, 9.98629535e-01], [6.12323400e-17, 0.00000000e+00, 1.00000000e+00]] + ] + ]) + + face_bounds = np.array([ + [[1.4968158, 1.57079633], [4.71238898, 0.0]], + [[1.4968158, 1.57079633], [0.0, 1.57079633]], + [[1.4968158, 1.57079633], [3.14159265, 0.0]], + [[1.4968158, 1.57079633], [0.0, 3.14159265]] + ]) + + constLat_cart = 1.0 + + weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds, is_directed=False) + # Define the expected DataFrame + expected_weight_df = pd.DataFrame({ + 'face_index': [0, 1, 2, 3], + 'weight': [0.25, 0.25, 0.25, 0.25] + }) + + # Assert that the DataFrame matches the expected DataFrame + pd.testing.assert_frame_equal(weight_df, expected_weight_df) + + + def test_get_zonal_face_interval_pole(self): + #The face is touching the pole + face_edges_cart = np.array([ + [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01], + [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]], + + [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01], + [6.12323400e-17, 0.00000000e+00, -1.00000000e+00]], + + [[6.12323400e-17, 0.00000000e+00, -1.00000000e+00], + [3.20465306e-18, -5.23359562e-02, -9.98629535e-01]], + + [[3.20465306e-18, -5.23359562e-02, -9.98629535e-01], + [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]] + ]) + + # Corrected face_bounds + face_bounds = np.array([ + [-1.57079633, -1.4968158], + [3.14159265, 0.] + ]) + constLat_cart = -0.9986295347545738 + + weight_df = _get_zonal_face_interval(face_edges_cart, constLat_cart, face_bounds, is_directed=False) + # No Nan values should be present in the weight_df + self.assertFalse(weight_df.isnull().values.any()) + + def test_get_zonal_faces_weight_at_constLat_latlonface(self): face_0 = [[np.deg2rad(350), np.deg2rad(40)], [np.deg2rad(350), np.deg2rad(20)], [np.deg2rad(10), np.deg2rad(20)], [np.deg2rad(10), np.deg2rad(40)]] @@ -449,10 +825,7 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self): # Assert the results is the same to the 3 decimal places weight_df = _get_zonal_faces_weight_at_constLat(np.array([ face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes - ]), - np.sin(np.deg2rad(20)), - latlon_bounds, - is_directed=False, is_latlonface=True) + ]), np.sin(np.deg2rad(20)), latlon_bounds, is_directed=False, is_latlonface=True) nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3) @@ -463,8 +836,5 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self): # It's edges are all GCA with self.assertRaises(ValueError): _get_zonal_faces_weight_at_constLat(np.array([ - face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes - ]), - np.deg2rad(20), - latlon_bounds, - is_directed=False) + face_0_edge_nodes, face_1_edge_nodes, face_2_edge_nodes + ]), np.deg2rad(20), latlon_bounds, is_directed=False) diff --git a/test/test_intersections.py b/test/test_intersections.py index b134d9029..769859b51 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -13,10 +13,6 @@ class TestGCAGCAIntersection(TestCase): def test_get_GCA_GCA_intersections_antimeridian(self): # Test the case where the two GCAs are on the antimeridian - - - - GCA1 = _lonlat_rad_to_xyz(np.deg2rad(170.0), np.deg2rad(89.99)) GCR1_cart = np.array([ _lonlat_rad_to_xyz(np.deg2rad(170.0), @@ -45,6 +41,7 @@ def test_get_GCA_GCA_intersections_antimeridian(self): ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = res_cart[0] # Test if the result is normalized self.assertTrue( @@ -70,7 +67,10 @@ def test_get_GCA_GCA_intersections_parallel(self): _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0) ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) - self.assertTrue(np.array_equal(res_cart, np.array([]))) + res_cart = res_cart[0] + expected_res = np.array(_lonlat_rad_to_xyz(0.5 * np.pi, 0.0)) + # Test if two results are equal within the error tolerance + self.assertAlmostEqual(np.linalg.norm(res_cart - expected_res), 0.0, delta=ERROR_TOLERANCE) def test_get_GCA_GCA_intersections_perpendicular(self): # Test the case where the two GCAs are perpendicular to each other @@ -85,7 +85,7 @@ def test_get_GCA_GCA_intersections_perpendicular(self): _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0]) ]) res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) - + res_cart = res_cart[0] # Test if the result is normalized self.assertTrue( np.allclose(np.linalg.norm(res_cart, axis=0), @@ -139,3 +139,15 @@ def test_GCA_constLat_intersections_two_pts(self): res = gca_constLat_intersection(GCR1_cart, np.sin(query_lat), verbose=False) self.assertTrue(res.shape[0] == 2) + + + def test_GCA_constLat_intersections_no_convege(self): + # It should return an one single point and a warning about unable to be converged should be raised + GCR1_cart = np.array([[-0.59647278, 0.59647278, -0.53706651], + [-0.61362973, 0.61362973, -0.49690755]]) + + constZ = -0.5150380749100542 + + with self.assertWarns(UserWarning): + res = gca_constLat_intersection(GCR1_cart, constZ, verbose=False) + self.assertTrue(res.shape[0] == 1) diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 15640d249..c06a938bd 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -46,6 +46,19 @@ def _point_within_gca_body( ): return False + if isclose( + GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON + ): + # If the pt and the GCA are on the same longitude (the y coordinates are the same) + if np.isclose( + GCRv0_lonlat[0], pt_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON + ): + # Now use the latitude to determine if the pt falls between the interval + return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1]) + else: + # If the pt and the GCA are not on the same longitude when the GCA is a longnitude arc, then the pt is not on the GCA + return False + # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point # Or if one of the endpoints is on the pole point, then the GCA goes through the pole point if ( diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 74a6eaa3e..1281c554d 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -24,6 +24,7 @@ REFERENCE_POINT_EQUATOR = np.array([1.0, 0.0, 0.0]) + def _unique_points(points, tolerance=ERROR_TOLERANCE): """Identify unique intersection points from a list of points, considering floating point precision errors. @@ -75,7 +76,6 @@ def error_radius(p1, p2): return unique_points - # General Helpers for Polygon Viz # ---------------------------------------------------------------------------------------------------------------------- @njit @@ -531,7 +531,9 @@ def _check_intersection(ref_edge, edges): if intersection_point.size != 0: # Handle both single point and multiple points case - if intersection_point.ndim == 1: # If there's only one point, make it a 2D array + if ( + intersection_point.ndim == 1 + ): # If there's only one point, make it a 2D array intersection_point = [intersection_point] # Convert to list of points # for each intersection point, check if it is a pole point @@ -576,7 +578,6 @@ def _classify_polygon_location(face_edge_cart): return "Equator" - def _get_latlonbox_width(latlonbox_rad): """Calculate the width of a latitude-longitude box in radians. The box should be represented by a 2x2 array in radians and lon0 represent the diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 16cc13194..74629b5d9 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -8,9 +8,9 @@ def _get_zonal_faces_weight_at_constLat( - faces_edges_cart, + faces_edges_cart_candidate, latitude_cart, - face_latlon_bound, + face_latlon_bound_candidate, is_directed=False, is_latlonface=False, is_face_GCA_list=None, @@ -20,14 +20,19 @@ def _get_zonal_faces_weight_at_constLat( Parameters ---------- - face_edges_cart : np.ndarray - A list of face polygon represented by edges in Cartesian coordinates. Shape: (n_faces, n_edges, 2, 3) + faces_edges_cart_candidate : np.ndarray + A list of the candidate face polygon represented by edges in Cartesian coordinates. + The faces must be selected in the previous step such that they will be intersected by the constant latitude. + It should have the same shape as the face_latlon_bound_candidate. + The input should not contain any 'INT_FILL_VALUE'. Shape: (n_faces(candidate), n_edges, 2, 3) latitude_cart : float The latitude in Cartesian coordinates (The normalized z coordinate) - face_latlon_bound : np.ndarray - The list of latitude and longitude bounds of faces. Shape: (n_faces,2, 2), + face_latlon_bound_candidate : np.ndarray + The list of latitude and longitude bounds of candidate faces. + It should have the same shape as the face_edges_cart_candidate. + Shape: (n_faces(candidate),,2, 2), [...,[lat_min, lat_max], [lon_min, lon_max],...] is_directed : bool, optional (default=False) @@ -52,11 +57,43 @@ def _get_zonal_faces_weight_at_constLat( - 'face_index': The index of the face (integer). - 'weight': The calculated weight of the face in radian (float). The DataFrame is indexed by the face indices, providing a mapping from each face to its corresponding weight. + + Notes + ----- + Special handling is implemented for the cases when the latitude_cart is close to 1 or -1, + which corresponds to the poles (90 and -90 degrees). In these cases, if a pole point is + inside a face, that face's value is the only value that should be considered. If the pole + point is not inside any face, it lies on the boundary of surrounding faces, and their weights + are considered evenly since they only contain points rather than intervals. + This treatment is hard-coded in the function and should be tested with appropriate test cases. """ + + # Special case if the latitude_cart is 1 or -1, meaning right at the pole + # If the latitude_cart is close to 1 or -1 (indicating the pole), handle it separately. + # The -90 and 90 treatment is hard-coded in the function, based on: + # If a pole point is inside a face, then this face's value is the only value that should be considered. + # If the pole point is not inside any face, then it's on the boundary of faces around it, so their weights are even. + if np.isclose(latitude_cart, 1, atol=ERROR_TOLERANCE) or np.isclose( + latitude_cart, -1, atol=ERROR_TOLERANCE + ): + # Now all candidate faces( the faces around the pole) are considered as the same weight + # If the face encompases the pole, then the weight is 1 + weights = { + face_index: 1 / len(faces_edges_cart_candidate) + for face_index in range(len(faces_edges_cart_candidate)) + } + # Convert weights to DataFrame + weights_df = pd.DataFrame( + list(weights.items()), columns=["face_index", "weight"] + ) + return weights_df + intervals_list = [] # Iterate through all faces and their edges - for face_index, face_edges in enumerate(faces_edges_cart): + for face_index, face_edges in enumerate(faces_edges_cart_candidate): + # Remove the Int_fill_value from the face_edges + face_edges = face_edges[np.all(face_edges != INT_FILL_VALUE, axis=(1, 2))] if is_face_GCA_list is not None: is_GCA_list = is_face_GCA_list[face_index] else: @@ -64,28 +101,55 @@ def _get_zonal_faces_weight_at_constLat( face_interval_df = _get_zonal_face_interval( face_edges, latitude_cart, - face_latlon_bound[face_index], + face_latlon_bound_candidate[face_index], is_directed=is_directed, is_latlonface=is_latlonface, is_GCA_list=is_GCA_list, ) - for _, row in face_interval_df.iterrows(): - intervals_list.append( - {"start": row["start"], "end": row["end"], "face_index": face_index} - ) + # If any end of the interval is NaN + if face_interval_df.isnull().values.any(): + # Skip this face as it is just being touched by the constant latitude + continue + # Check if the DataFrame is empty (start and end are both 0) + if (face_interval_df["start"] == 0).all() and ( + face_interval_df["end"] == 0 + ).all(): + # Skip this face as it is just being touched by the constant latitude + continue + else: + for _, row in face_interval_df.iterrows(): + intervals_list.append( + {"start": row["start"], "end": row["end"], "face_index": face_index} + ) intervals_df = pd.DataFrame(intervals_list) - overlap_contributions, total_length = _process_overlapped_intervals(intervals_df) + try: + overlap_contributions, total_length = _process_overlapped_intervals( + intervals_df + ) - # Calculate weights for each face - weights = { - face_index: overlap_contributions.get(face_index, 0.0) / total_length - for face_index in range(len(faces_edges_cart)) - } + # Calculate weights for each face + weights = { + face_index: overlap_contributions.get(face_index, 0.0) / total_length + for face_index in range(len(faces_edges_cart_candidate)) + } - # Convert weights to DataFrame - weights_df = pd.DataFrame(list(weights.items()), columns=["face_index", "weight"]) - return weights_df + # Convert weights to DataFrame + weights_df = pd.DataFrame( + list(weights.items()), columns=["face_index", "weight"] + ) + return weights_df + + except ValueError: + print(f"Face index: {face_index}") + print(f"Face edges information: {face_edges}") + print(f"Constant z0: {latitude_cart}") + print( + f"Face latlon bound information: {face_latlon_bound_candidate[face_index]}" + ) + print(f"Face interval information: {face_interval_df}") + # Handle the exception or propagate it further if necessary + raise def _is_edge_gca(is_GCA_list, is_latlonface, edges_z): @@ -149,8 +213,8 @@ def _get_faces_constLat_intersection_info( tuple A tuple containing: - intersections_pts_list_cart (list): A list of intersection points where each point is where an edge intersects with the latitude. - - overlap_flag (bool): A boolean indicating if any overlap with the latitude was detected. - - overlap_edge (np.ndarray or None): The edge that overlaps with the latitude, if any; otherwise, None. + - pt_lon_min (float): The min longnitude of the interseted intercal in radian if any; otherwise, None.. + - pt_lon_max (float): The max longnitude of the interseted intercal in radian, if any; otherwise, None. """ valid_edges_mask = ~(np.any(face_edges_cart == DUMMY_EDGE_VALUE, axis=(1, 2))) @@ -181,6 +245,7 @@ def _get_faces_constLat_intersection_info( intersections = gca_constLat_intersection( edge, latitude_cart, is_directed=is_directed ) + if intersections.size == 0: continue elif intersections.shape[0] == 2: @@ -188,8 +253,9 @@ def _get_faces_constLat_intersection_info( else: intersections_pts_list_cart.append(intersections[0]) - # Handle unique intersections and check for convex or concave cases + # Find the unique intersection points unique_intersections = np.unique(intersections_pts_list_cart, axis=0) + if len(unique_intersections) == 2: # TODO: vectorize? unique_intersection_lonlat = np.array( @@ -199,10 +265,28 @@ def _get_faces_constLat_intersection_info( sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) pt_lon_min, pt_lon_max = sorted_lonlat[:, 0] return unique_intersections, pt_lon_min, pt_lon_max - elif len(unique_intersections) != 0: - raise ValueError( - "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" - ) + elif len(unique_intersections) == 1: + return unique_intersections, None, None + elif len(unique_intersections) != 0 and len(unique_intersections) != 1: + # If the unique intersections numbers is larger than n_edges * 2, then it means the face is concave + if len(unique_intersections) > len(valid_edges) * 2: + raise ValueError( + "UXarray doesn't support concave face with intersections points as currently, please modify your grids accordingly" + ) + else: + # Now return all the intersections points and the pt_lon_min, pt_lon_max + unique_intersection_lonlat = np.array( + [_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]) for pt in unique_intersections] + ) + + sorted_lonlat = np.sort(unique_intersection_lonlat, axis=0) + # Extract the minimum and maximum longitudes + pt_lon_min, pt_lon_max = ( + np.min(sorted_lonlat[:, 0]), + np.max(sorted_lonlat[:, 0]), + ) + + return unique_intersections, pt_lon_min, pt_lon_max elif len(unique_intersections) == 0: raise ValueError( "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" @@ -257,40 +341,87 @@ def _get_zonal_face_interval( face_lon_bound_left, face_lon_bound_right = face_latlon_bound[1] # Call the vectorized function to process all edges - ( - unique_intersections, - pt_lon_min, - pt_lon_max, - ) = _get_faces_constLat_intersection_info( - face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed - ) + try: + # Call the vectorized function to process all edges + unique_intersections, pt_lon_min, pt_lon_max = ( + _get_faces_constLat_intersection_info( + face_edges_cart, latitude_cart, is_GCA_list, is_latlonface, is_directed + ) + ) - # Convert intersection points to longitude-latitude - longitudes = np.array( - [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections] - ) + # Handle the special case where the unique_intersections is 1, which means the face is just being touched + if len(unique_intersections) == 1: + # If the face is just being touched, then just return the empty DataFrame + return pd.DataFrame({"start": [0.0], "end": [0.0]}, index=[0]) + + # Convert intersection points to longitude-latitude + longitudes = np.array( + [_xyz_to_lonlat_rad(*pt.tolist())[0] for pt in unique_intersections] + ) - # Handle special wrap-around cases by checking the face bounds - if face_lon_bound_left >= face_lon_bound_right: - if not ( - (pt_lon_max >= np.pi and pt_lon_min >= np.pi) - or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) + # Handle special wrap-around cases by checking the face bounds + if face_lon_bound_left >= face_lon_bound_right or ( + face_lon_bound_left == 0 and face_lon_bound_right == 2 * np.pi + ): + if not ( + (pt_lon_max >= np.pi and pt_lon_min >= np.pi) + or (0 <= pt_lon_max <= np.pi and 0 <= pt_lon_min <= np.pi) + ): + # If the anti-meridian is crossed, instead of just being touched,add the wrap-around points + if pt_lon_max != 2 * np.pi and pt_lon_min != 0: + # They're at different sides of the 0-lon, adding wrap-around points + longitudes = np.append(longitudes, [0.0, 2 * np.pi]) + elif pt_lon_max >= np.pi and pt_lon_min == 0: + # That means the face is actually from pt_lon_max to 2*pi. + # Replace the 0 in longnitude with 2*pi + longitudes[longitudes == 0] = 2 * np.pi + + # Ensure longitudes are sorted + longitudes = np.unique(longitudes) + longitudes.sort() + + # Split the sorted longitudes into start and end points of intervals + starts = longitudes[::2] # Start points + ends = longitudes[1::2] # End points + + # Create the intervals DataFrame + Intervals_df = pd.DataFrame({"start": starts, "end": ends}) + # For consistency, sort the intervals by start + interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) + return interval_df_sorted + + except ValueError as e: + default_print_options = np.get_printoptions() + if ( + str(e) + == "No intersections are found for the face, please make sure the build_latlon_box generates the correct results" ): - # They're at different sides of the 0-lon, adding wrap-around points - longitudes = np.append(longitudes, [0.0, 2 * np.pi]) + # Set print options for full precision + np.set_printoptions(precision=16, suppress=False) - # Ensure longitudes are sorted - longitudes.sort() + print( + "ValueError: No intersections are found for the face, please make sure the build_latlon_box generates the correct results" + ) + print(f"Face edges information: {face_edges_cart}") + print(f"Constant z_0: {latitude_cart}") + print(f"Face latlon bound information: {face_latlon_bound}") - # Split the sorted longitudes into start and end points of intervals - starts = longitudes[::2] # Start points - ends = longitudes[1::2] # End points + # Reset print options to default + np.set_printoptions(**default_print_options) + + raise + else: + # Set print options for full precision + np.set_printoptions(precision=17, suppress=False) - # Create the intervals DataFrame - Intervals_df = pd.DataFrame({"start": starts, "end": ends}) - # For consistency, sort the intervals by start - interval_df_sorted = Intervals_df.sort_values(by="start").reset_index(drop=True) - return interval_df_sorted + print(f"Face edges information: {face_edges_cart}") + print(f"Constant z_0: {latitude_cart}") + print(f"Face latlon bound information: {face_latlon_bound}") + + # Reset print options to default + np.set_printoptions(**default_print_options) + + raise # Re-raise the exception if it's not the expected ValueError def _process_overlapped_intervals(intervals_df): @@ -327,7 +458,7 @@ def _process_overlapped_intervals(intervals_df): events.append((row["start"], "start", row["face_index"])) events.append((row["end"], "end", row["face_index"])) - events.sort() # Sort by position and then by start/end + events.sort(key=lambda x: (x[0], x[1])) active_faces = set() last_position = None @@ -335,6 +466,8 @@ def _process_overlapped_intervals(intervals_df): overlap_contributions = {} for position, event_type, face_idx in events: + if face_idx == 51: + pass if last_position is not None and active_faces: segment_length = position - last_position segment_weight = segment_length / len(active_faces) if active_faces else 0 @@ -345,9 +478,28 @@ def _process_overlapped_intervals(intervals_df): total_length += segment_length if event_type == "start": - active_faces.add(face_idx) + # use try catch to handle the case where the face_idx is not be able to be added + try: + active_faces.add(face_idx) + except Exception as e: + print(f"An error occurred: {e}") + print(f"Face index: {face_idx}") + print(f"Position: {position}") + print(f"Event type: {event_type}") + print(f"Active faces: {active_faces}") + print(f"Last position: {last_position}") + print(f"Total length: {total_length}") + print(f"Overlap contributions: {overlap_contributions}") + print(f"Intervals data: {intervals_df}") + raise + elif event_type == "end": - active_faces.remove(face_idx) + if face_idx in active_faces: + active_faces.remove(face_idx) + else: + raise ValueError( + f"Error: Trying to remove face_idx {face_idx} which is not in active_faces" + ) last_position = position diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index d57969db5..7324338c3 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -1,15 +1,13 @@ import numpy as np -from uxarray.constants import ERROR_TOLERANCE +from uxarray.constants import MACHINE_EPSILON, ERROR_TOLERANCE from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat -from uxarray.grid.arcs import point_within_gca +from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform import warnings -from uxarray.utils.computing import cross_fma, allclose, cross, dot +from uxarray.utils.computing import cross_fma -import uxarray.constants - -def gca_gca_intersection(gca1_cart, gca2_cart): +def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False): """Calculate the intersection point(s) of two Great Circle Arcs (GCAs) in a Cartesian coordinate system. @@ -29,21 +27,13 @@ def gca_gca_intersection(gca1_cart, gca2_cart): Returns ------- np.ndarray - Cartesian coordinates of the intersection point(s). + Cartesian coordinates of the intersection point(s). Returns an empty array if no intersections, + a 2D array with one row if one intersection, and a 2D array with two rows if two intersections. Raises ------ ValueError If the input GCAs are not in the cartesian [x, y, z] format. - - - - Warning - ------- - If the current input data cannot be computed accurately using floating-point arithmetic. Use with care - - If running on the Windows system with fma_disabled=False since the C/C++ implementation of FMA in MS Windows - is fundamentally broken. (bug report: https://bugs.python.org/msg312480) """ # Support lists as an input @@ -56,13 +46,12 @@ def gca_gca_intersection(gca1_cart, gca2_cart): w0, w1 = gca1_cart v0, v1 = gca2_cart - if not uxarray.constants.ENABLE_FMA: - # Compute normals and orthogonal bases without FMA - w0w1_norm = cross(w0, w1) - v0v1_norm = cross(v0, v1) - cross_norms = cross(w0w1_norm, v0v1_norm) + # Compute normals and orthogonal bases using FMA + if fma_disabled: + w0w1_norm = np.cross(w0, w1) + v0v1_norm = np.cross(v0, v1) + cross_norms = np.cross(w0w1_norm, v0v1_norm) else: - # Compute normals and orthogonal bases using FMA w0w1_norm = cross_fma(w0, w1) v0v1_norm = cross_fma(v0, v1) cross_norms = cross_fma(w0w1_norm, v0v1_norm) @@ -76,46 +65,54 @@ def gca_gca_intersection(gca1_cart, gca2_cart): ) # Check perpendicularity conditions and floating-point arithmetic limitations - if not allclose(dot(w0w1_norm, w0), 0, atol=ERROR_TOLERANCE) or not allclose( - dot(w0w1_norm, w1), 0, atol=ERROR_TOLERANCE - ): + if not np.allclose( + np.dot(w0w1_norm, w0), 0, atol=MACHINE_EPSILON + ) or not np.allclose(np.dot(w0w1_norm, w1), 0, atol=MACHINE_EPSILON): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution." ) - if not allclose(dot(v0v1_norm, v0), 0, atol=ERROR_TOLERANCE) or not allclose( - dot(v0v1_norm, v1), 0, atol=ERROR_TOLERANCE - ): + if not np.allclose( + np.dot(v0v1_norm, v0), 0, atol=MACHINE_EPSILON + ) or not np.allclose(np.dot(v0v1_norm, v1), 0, atol=MACHINE_EPSILON): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " ) - if not allclose( - dot(cross_norms, v0v1_norm), 0, atol=ERROR_TOLERANCE - ) or not allclose(dot(cross_norms, w0w1_norm), 0, atol=ERROR_TOLERANCE): + if not np.allclose( + np.dot(cross_norms, v0v1_norm), 0, atol=MACHINE_EPSILON + ) or not np.allclose(np.dot(cross_norms, w0w1_norm), 0, atol=MACHINE_EPSILON): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " ) # If the cross_norms is zero, the two GCAs are parallel - if allclose(cross_norms, 0, atol=ERROR_TOLERANCE): - return np.array([]) + if np.allclose(cross_norms, 0, atol=MACHINE_EPSILON): + res = [] + # Check if the two GCAs are overlapping + if point_within_gca(v0, [w0, w1]): + # The two GCAs are overlapping, return both end points + res.append(v0) + + if point_within_gca(v1, [w0, w1]): + res.append(v1) + return np.array(res) # Normalize the cross_norms cross_norms = cross_norms / np.linalg.norm(cross_norms) x1 = cross_norms x2 = -x1 - res = np.array([]) + res = [] # Determine which intersection point is within the GCAs range if point_within_gca(x1, [w0, w1]) and point_within_gca(x1, [v0, v1]): - res = np.append(res, x1) + res.append(x1) - elif point_within_gca(x2, [w0, w1]) and point_within_gca(x2, [v0, v1]): - res = np.append(res, x2) + if point_within_gca(x2, [w0, w1]) and point_within_gca(x2, [v0, v1]): + res.append(x2) - return res + return np.array(res) def gca_constLat_intersection( @@ -155,11 +152,39 @@ def gca_constLat_intersection( ------- If running on the Windows system with fma_disabled=False since the C/C++ implementation of FMA in MS Windows is fundamentally broken. (bug report: https://bugs.python.org/msg312480) + + If the intersection point cannot be converged using the Newton-Raphson method, the initial guess intersection + point is used instead, proceed with caution. """ x1, x2 = gca_cart + # Check if the constant latitude has the same latitude as the GCA endpoints + # We are using the relative tolerance and ERROR_TOLERANCE since the constZ is calculated from np.sin, which + # may have some floating-point error. + res = None + if np.isclose(x1[2], constZ, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE): + res = np.array([x1]) if res is None else np.vstack((res, x1)) + + if np.isclose(x2[2], constZ, rtol=ERROR_TOLERANCE, atol=ERROR_TOLERANCE): + res = np.array([x2]) if res is None else np.vstack((res, x2)) + + if res is not None: + return res + + # If the constant latitude is not the same as the GCA endpoints, calculate the intersection point + lat_min = extreme_gca_latitude(gca_cart, extreme_type="min") + lat_max = extreme_gca_latitude(gca_cart, extreme_type="max") + + constLat_rad = np.arcsin(constZ) + + # Check if the constant latitude is within the GCA range + # Because the constant latitude is calculated from np.sin, which may have some floating-point error, + if not in_between(lat_min, constLat_rad, lat_max): + pass + return np.array([]) + if fma_disabled: - n = cross(x1, x2) + n = np.cross(x1, x2) else: # Raise a warning for Windows users @@ -173,7 +198,7 @@ def gca_constLat_intersection( nx, ny, nz = n - s_tilde = np.sqrt(nx**2 + ny**2 - np.linalg.norm(n) ** 2 * constZ**2) + s_tilde = np.sqrt(nx**2 + ny**2 - (nx**2 + ny**2 + nz**2) * constZ**2) p1_x = -(1.0 / (nx**2 + ny**2)) * (constZ * nx * nz + s_tilde * ny) p2_x = -(1.0 / (nx**2 + ny**2)) * (constZ * nx * nz - s_tilde * ny) p1_y = -(1.0 / (nx**2 + ny**2)) * (constZ * ny * nz - s_tilde * nx) @@ -186,19 +211,50 @@ def gca_constLat_intersection( # Now test which intersection point is within the GCA range if point_within_gca(p1, gca_cart, is_directed=is_directed): - converged_pt = _newton_raphson_solver_for_gca_constLat( - p1, gca_cart, verbose=verbose - ) - res = ( - np.array([converged_pt]) if res is None else np.vstack((res, converged_pt)) - ) + try: + converged_pt = _newton_raphson_solver_for_gca_constLat( + p1, gca_cart, verbose=verbose + ) + + if converged_pt is None: + # The point is not be able to be converged using the jacobi method, raise a warning and continue with p2 + warnings.warn( + "The intersection point cannot be converged using the Newton-Raphson method. " + "The initial guess intersection point is used instead, procced with caution." + ) + res = np.array([p1]) if res is None else np.vstack((res, p1)) + else: + res = ( + np.array([converged_pt]) + if res is None + else np.vstack((res, converged_pt)) + ) + except RuntimeError: + print(f"Error encountered with initial guess: {p1}") + print(f"gca_cart: {gca_cart}") + raise if point_within_gca(p2, gca_cart, is_directed=is_directed): - converged_pt = _newton_raphson_solver_for_gca_constLat( - p2, gca_cart, verbose=verbose - ) - res = ( - np.array([converged_pt]) if res is None else np.vstack((res, converged_pt)) - ) + try: + converged_pt = _newton_raphson_solver_for_gca_constLat( + p2, gca_cart, verbose=verbose + ) + if converged_pt is None: + # The point is not be able to be converged using the jacobi method, raise a warning and continue with p2 + warnings.warn( + "The intersection point cannot be converged using the Newton-Raphson method. " + "The initial guess intersection point is used instead, procced with caution." + ) + res = np.array([p2]) if res is None else np.vstack((res, p2)) + else: + res = ( + np.array([converged_pt]) + if res is None + else np.vstack((res, converged_pt)) + ) + except RuntimeError: + print(f"Error encountered with initial guess: {p2}") + print(f"gca_cart: {gca_cart}") + raise return res if res is not None else np.array([]) diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py index 020aae1d9..33f9f75da 100644 --- a/uxarray/grid/utils.py +++ b/uxarray/grid/utils.py @@ -1,5 +1,5 @@ import numpy as np -from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE +from uxarray.constants import INT_FILL_VALUE, MACHINE_EPSILON import warnings import uxarray.utils.computing as ac_utils @@ -116,11 +116,22 @@ def _inv_jacobian(x0, x1, y0, y1, z0, z1, x_i_old, y_i_old): [ac_utils._fmms(y0, z1, z0, y1), ac_utils._fmms(x0, z1, z0, x1)], [2 * x_i_old, 2 * y_i_old], ] + + # First check if the Jacobian matrix is singular + if np.linalg.matrix_rank(jacobian) < 2: + warnings.warn("The Jacobian matrix is singular.") + return None + try: inverse_jacobian = np.linalg.inv(jacobian) - except np.linalg.LinAlgError: - raise warnings("Warning: Singular Jacobian matrix encountered.") - return None + except np.linalg.LinAlgError as e: + # Print out the error message + + cond_number = np.linalg.cond(jacobian) + print(f"Condition number: {cond_number}") + print(f"Jacobian matrix:\n{jacobian}") + print(f"An error occurred: {e}") + raise return inverse_jacobian @@ -141,7 +152,7 @@ def _newton_raphson_solver_for_gca_constLat( Returns: np.ndarray or None: The intersection point or None if the solver fails to converge. """ - tolerance = ERROR_TOLERANCE + tolerance = MACHINE_EPSILON * 100 w0_cart, w1_cart = gca_cart error = float("inf") constZ = init_cart[2] @@ -164,19 +175,23 @@ def _newton_raphson_solver_for_gca_constLat( ] ) - j_inv = _inv_jacobian( - w0_cart[0], - w1_cart[0], - w0_cart[1], - w1_cart[1], - w0_cart[2], - w1_cart[2], - y_guess[0], - y_guess[1], - ) + try: + j_inv = _inv_jacobian( + w0_cart[0], + w1_cart[0], + w0_cart[1], + w1_cart[1], + w0_cart[2], + w1_cart[2], + y_guess[0], + y_guess[1], + ) - if j_inv is None: - return None + if j_inv is None: + return None + except RuntimeError as e: + print(f"Encountered an error: {e}") + raise y_new = y_guess - np.matmul(j_inv, f_vector) error = np.max(np.abs(y_guess - y_new)) From 94952cade4b94f92b612d8a65ccfa7c605aa79b9 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sat, 14 Sep 2024 18:13:31 -0700 Subject: [PATCH 05/24] fix CI fail --- test/test_integrate.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_integrate.py b/test/test_integrate.py index 1b001818d..8fa5b2eb8 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -8,8 +8,10 @@ import numpy.testing as nt import uxarray as ux +from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE from uxarray.grid.coordinates import _lonlat_rad_to_xyz -from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, _get_zonal_faces_weight_at_constLat +from uxarray.grid.integrate import _get_zonal_face_interval, _process_overlapped_intervals, _get_zonal_faces_weight_at_constLat,_get_faces_constLat_intersection_info + current_path = Path(os.path.dirname(os.path.realpath(__file__))) From ffad2907e7f51820fcecf469ce42c937d1e72167 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 11:47:27 -0700 Subject: [PATCH 06/24] fix facebounds fail --- test/test_geometry.py | 115 ++++++++++++++++++++++++++++++++++++++- uxarray/grid/geometry.py | 71 ++++++++++++------------ 2 files changed, 148 insertions(+), 38 deletions(-) diff --git a/test/test_geometry.py b/test/test_geometry.py index 6041109c8..f9aba5db8 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -163,7 +163,6 @@ def test_pole_point_inside_polygon_from_vertice_cross(self): self.assertTrue(result, "North pole should be inside the polygon") - class TestLatlonBoundUtils(TestCase): def _max_latitude_rad_iterative(self, gca_cart): @@ -369,6 +368,19 @@ def test_extreme_gca_latitude_max(self): expected_max_latitude, delta=ERROR_TOLERANCE) + def test_extreme_gca_latitude_max_short(self): + # Define a great circle arc in 3D space that has a small span + gca_cart = np.array( [[ 0.65465367, -0.37796447, -0.65465367], [ 0.6652466, -0.33896007, -0.6652466 ]]) + + # Calculate the maximum latitude + max_latitude = ux.grid.arcs.extreme_gca_latitude(gca_cart, 'max') + + # Check if the maximum latitude is correct + expected_max_latitude = self._max_latitude_rad_iterative(gca_cart) + self.assertAlmostEqual(max_latitude, + expected_max_latitude, + delta=ERROR_TOLERANCE) + def test_extreme_gca_latitude_min(self): # Define a great circle arc that is symmetrical around 0 degrees longitude gca_cart = np.array([ @@ -624,6 +636,103 @@ def test_populate_bounds_antimeridian(self): bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat) nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_equator(self): + # the face is touching the equator + face_edges_cart = np.array([ + [[0.99726469, -0.05226443, -0.05226443], [0.99862953, 0.0, -0.05233596]], + [[0.99862953, 0.0, -0.05233596], [1.0, 0.0, 0.0]], + [[1.0, 0.0, 0.0], [0.99862953, -0.05233596, 0.0]], + [[0.99862953, -0.05233596, 0.0], [0.99726469, -0.05226443, -0.05226443]] + ] + ) + # Apply the inverse transformation to get the lat lon coordinates + face_edges_lonlat = np.array( + [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart]) + + bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat) + expected_bounds = np.array([[-0.05235988, 0], [6.23082543, 0]]) + nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + + def test_populate_bounds_southSphere(self): + # The face is near the south pole but doesn't contains the pole + face_edges_cart = np.array([ + [[-1.04386773e-01, -5.20500333e-02, -9.93173799e-01], [-1.04528463e-01, -1.28010448e-17, -9.94521895e-01]], + [[-1.04528463e-01, -1.28010448e-17, -9.94521895e-01], [-5.23359562e-02, -6.40930613e-18, -9.98629535e-01]], + [[-5.23359562e-02, -6.40930613e-18, -9.98629535e-01], [-5.22644277e-02, -5.22644277e-02, -9.97264689e-01]], + [[-5.22644277e-02, -5.22644277e-02, -9.97264689e-01], [-1.04386773e-01, -5.20500333e-02, -9.93173799e-01]] + ]) + + # Apply the inverse transformation to get the lat lon coordinates + face_edges_lonlat = np.array( + [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart]) + + bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat) + expected_bounds = np.array([[-1.51843645, -1.45388627], [3.14159265, 3.92699082]]) + nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + + def test_populate_bounds_near_pole(self): + # The face is near the south pole but doesn't contains the pole + face_edges_cart = np.array([ + [[3.58367950e-01, 0.00000000e+00, -9.33580426e-01], [3.57939780e-01, 4.88684203e-02, -9.32465008e-01]], + [[3.57939780e-01, 4.88684203e-02, -9.32465008e-01], [4.06271283e-01, 4.78221112e-02, -9.12500241e-01]], + [[4.06271283e-01, 4.78221112e-02, -9.12500241e-01], [4.06736643e-01, 2.01762691e-16, -9.13545458e-01]], + [[4.06736643e-01, 2.01762691e-16, -9.13545458e-01], [3.58367950e-01, 0.00000000e+00, -9.33580426e-01]] + ]) + + # Apply the inverse transformation to get the lat lon coordinates + face_edges_lonlat = np.array( + [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart]) + + bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat) + expected_bounds = np.array([[-1.20427718, -1.14935491], [0,0.13568803]]) + nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + + def test_populate_bounds_near_pole2(self): + # The face is near the south pole but doesn't contains the pole + face_edges_cart = np.array([ + [[3.57939780e-01, -4.88684203e-02, -9.32465008e-01], [3.58367950e-01, 0.00000000e+00, -9.33580426e-01]], + [[3.58367950e-01, 0.00000000e+00, -9.33580426e-01], [4.06736643e-01, 2.01762691e-16, -9.13545458e-01]], + [[4.06736643e-01, 2.01762691e-16, -9.13545458e-01], [4.06271283e-01, -4.78221112e-02, -9.12500241e-01]], + [[4.06271283e-01, -4.78221112e-02, -9.12500241e-01], [3.57939780e-01, -4.88684203e-02, -9.32465008e-01]] + ]) + + + # Apply the inverse transformation to get the lat lon coordinates + face_edges_lonlat = np.array( + [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart]) + + bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat) + expected_bounds = np.array([[-1.20427718, -1.14935491], [6.147497,4.960524e-16]]) + nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + + def test_populate_bounds_long_face(self): + """Test case where one of the face edges is a longitude GCA.""" + face_edges_cart = np.array([ + [[9.9999946355819702e-01, -6.7040475551038980e-04, 8.0396590055897832e-04], + [9.9999439716339111e-01, -3.2541253603994846e-03, -8.0110825365409255e-04]], + [[9.9999439716339111e-01, -3.2541253603994846e-03, -8.0110825365409255e-04], + [9.9998968839645386e-01, -3.1763643492013216e-03, -3.2474612817168236e-03]], + [[9.9998968839645386e-01, -3.1763643492013216e-03, -3.2474612817168236e-03], + [9.9998861551284790e-01, -8.2993711112067103e-04, -4.7004125081002712e-03]], + [[9.9998861551284790e-01, -8.2993711112067103e-04, -4.7004125081002712e-03], + [9.9999368190765381e-01, 1.7522916896268725e-03, -3.0944822356104851e-03]], + [[9.9999368190765381e-01, 1.7522916896268725e-03, -3.0944822356104851e-03], + [9.9999833106994629e-01, 1.6786820488050580e-03, -6.4892979571595788e-04]], + [[9.9999833106994629e-01, 1.6786820488050580e-03, -6.4892979571595788e-04], + [9.9999946355819702e-01, -6.7040475551038980e-04, 8.0396590055897832e-04]] + ]) + + face_edges_lonlat = np.array( + [[_xyz_to_lonlat_rad(*edge[0]), _xyz_to_lonlat_rad(*edge[1])] for edge in face_edges_cart]) + + bounds = _populate_face_latlon_bound(face_edges_cart, face_edges_lonlat) + + # The expected bounds should not contains the south pole [0,-0.5*np.pi] + self.assertTrue(bounds[1][0] != 0.0) + + + + def test_populate_bounds_node_on_pole(self): # Generate a normal face that is crossing the antimeridian vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] @@ -1119,6 +1228,8 @@ def test_populate_bounds_normal(self): is_GCA_list=[True, False, True, False]) nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + + def test_populate_bounds_antimeridian(self): # Generate a normal face that is crossing the antimeridian vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] @@ -1176,7 +1287,7 @@ def test_populate_bounds_node_on_pole(self): nt.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) def test_populate_bounds_edge_over_pole(self): - # Generate a normal face that is crossing the antimeridian + # Generate a normal face that is around the north pole vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 36831940b..1af7ed73f 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -666,7 +666,7 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): >>> _insert_pt_in_latlonbox(np.array([[1.0, 2.0], [3.0, 4.0]]),np.array([1.5, 3.5])) array([[1.0, 2.0], [3.0, 4.0]]) """ - if all(new_pt == INT_FILL_VALUE): + if np.all(new_pt == INT_FILL_VALUE): return old_box latlon_box = np.copy(old_box) # Create a copy of the old box @@ -694,17 +694,18 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): # Check for pole points and update latitudes is_pole_point = ( lon_pt == INT_FILL_VALUE - and isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE) - or isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE) + and np.isclose( + new_pt[0], [0.5 * np.pi, -0.5 * np.pi], atol=ERROR_TOLERANCE + ).any() ) if is_pole_point: # Check if the new point is close to the North Pole - if isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE): + if np.isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE): latlon_box[0][1] = 0.5 * np.pi # Check if the new point is close to the South Pole - elif isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE): + elif np.isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE): latlon_box[0][0] = -0.5 * np.pi return latlon_box @@ -850,7 +851,9 @@ def _populate_face_latlon_bound( ) # Check if the node matches the pole point or if the pole point is within the edge_cart - if allclose(n1_cart, pole_point, atol=ERROR_TOLERANCE) or point_within_gca( + if np.allclose( + n1_cart, pole_point, atol=ERROR_TOLERANCE + ) or point_within_gca( pole_point, np.array([n1_cart, n2_cart]), is_directed=False ): is_center_pole = False @@ -929,16 +932,16 @@ def _populate_face_latlon_bound( ) # Insert extreme latitude points into the latlonbox if they differ from the node latitudes - if not isclose( + if not np.isclose( node1_lat_rad, lat_max, atol=ERROR_TOLERANCE - ) and not isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE): + ) and not np.isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE): # Insert the maximum latitude face_latlon_array = _insert_pt_in_latlonbox( face_latlon_array, np.array([lat_max, node1_lon_rad]) ) - elif not isclose( + elif not np.isclose( node1_lat_rad, lat_min, atol=ERROR_TOLERANCE - ) and not isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE): + ) and not np.isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE): # Insert the minimum latitude face_latlon_array = _insert_pt_in_latlonbox( face_latlon_array, np.array([lat_min, node1_lon_rad]) @@ -1018,7 +1021,7 @@ def _populate_bounds( intervals_tuple_list = [] intervals_name_list = [] - face_edges_cartesian = _get_cartesian_face_edge_nodes( + faces_edges_cartesian = _get_cartesian_face_edge_nodes( grid.face_node_connectivity.values, grid.n_face, grid.n_max_face_edges, @@ -1027,7 +1030,7 @@ def _populate_bounds( grid.node_z.values, ) - face_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( + faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes( grid.face_node_connectivity.values, grid.n_face, grid.n_max_face_edges, @@ -1035,42 +1038,37 @@ def _populate_bounds( grid.node_lat.values, ) - face_node_connectivity = grid.face_node_connectivity.values - - # # TODO: vectorize dummy face check - s = face_edges_cartesian.shape - dummy_face_face_edges_cart = np.any( - face_edges_cartesian.reshape((s[0], s[1] * s[2] * s[3])) == INT_FILL_VALUE, - axis=1, - ) + for face_idx, face_nodes in enumerate(grid.face_node_connectivity): + face_edges_cartesian = faces_edges_cartesian[face_idx] - s = face_edges_lonlat_rad.shape - dummy_face_face_edges_latlon = np.any( - face_edges_lonlat_rad.reshape((s[0], s[1] * s[2] * s[3])) == INT_FILL_VALUE, - axis=1, - ) + # Remove the edge in the face that contains the fill value + face_edges_cartesian = face_edges_cartesian[ + np.all(face_edges_cartesian != INT_FILL_VALUE, axis=(1, 2)) + ] - dummy_faces = dummy_face_face_edges_cart | dummy_face_face_edges_latlon + face_edges_lonlat_rad = faces_edges_lonlat_rad[face_idx] - for face_idx, face_nodes in enumerate(face_node_connectivity): - if dummy_faces[face_idx]: - # skip dummy faces - continue + # Remove the edge in the face that contains the fill value + face_edges_lonlat_rad = face_edges_lonlat_rad[ + np.all(face_edges_lonlat_rad != INT_FILL_VALUE, axis=(1, 2)) + ] is_GCA_list = ( is_face_GCA_list[face_idx] if is_face_GCA_list is not None else None ) temp_latlon_array[face_idx] = _populate_face_latlon_bound( - face_edges_cartesian[face_idx], - face_edges_lonlat_rad[face_idx], + face_edges_cartesian, + face_edges_lonlat_rad, is_latlonface=is_latlonface, is_GCA_list=is_GCA_list, ) - # # do we need this ? - # assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1] - # assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1] + if temp_latlon_array[face_idx][0][0] == temp_latlon_array[face_idx][0][1]: + pass + + assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1] + assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1] lat_array = temp_latlon_array[face_idx][0] # Now store the latitude intervals in the tuples @@ -1090,7 +1088,7 @@ def _populate_bounds( attrs={ "cf_role": "face_latlon_bounds", "_FillValue": INT_FILL_VALUE, - "long_name": "Latitude and Longitude bounds for each face in radians.", + "long_name": "Provides the latitude and longitude bounds for each face in radians.", "start_index": INT_DTYPE(0), "latitude_intervalsIndex": intervalsIndex, "latitude_intervals_name_map": df_intervals_map, @@ -1103,6 +1101,7 @@ def _populate_bounds( grid._ds["bounds"] = bounds + def _construct_hole_edge_indices(edge_face_connectivity): """Index the missing edges on a partial grid with holes, that is a region of the grid that is not covered by any geometry.""" From 8e53754ae0d1af611760f5f8649cc5000fe2db87 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 12:55:42 -0700 Subject: [PATCH 07/24] fix `extreme_gca` --- test/test_geometry.py | 31 +++++++++++++++++++++++++++++++ uxarray/grid/arcs.py | 17 +++++++++-------- uxarray/grid/geometry.py | 2 -- 3 files changed, 40 insertions(+), 10 deletions(-) diff --git a/test/test_geometry.py b/test/test_geometry.py index f9aba5db8..21b85cbe5 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -27,6 +27,13 @@ grid_files = [gridfile_CSne8, gridfile_geoflow] data_files = [datafile_CSne30, datafile_geoflow] +grid_quad_hex = current_path/ "meshfiles" / "ugrid" / "quad-hexagon" / "grid.nc" +grid_geoflow = current_path/ "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc" +grid_mpas = current_path/ "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc" + + +# List of grid files to test +grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas] class TestAntimeridian(TestCase): @@ -1402,6 +1409,30 @@ def test_populate_bounds_GCAList_mix(self): nt.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) +# Test class +class TestLatlonBoundsFiles: + + def test_face_bounds(self): + """Test to ensure ``Grid.face_bounds`` works correctly for all grid files.""" + for grid_path in grid_files_latlonBound: + try: + # Open the grid file + self.uxgrid = ux.open_grid(grid_path) + + # Test: Ensure the bounds are obtained + bounds = self.uxgrid.bounds + assert bounds is not None, f"Grid.face_bounds should not be None for {grid_path}" + + except Exception as e: + # Print the failing grid file and re-raise the exception + print(f"Test failed for grid file: {grid_path}") + raise e + + finally: + # Clean up the grid object + del self.uxgrid + + class TestGeoDataFrame(TestCase): def test_to_gdf(self): diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index c06a938bd..41f3dc15e 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -2,7 +2,7 @@ # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place -from uxarray.grid.coordinates import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar +from uxarray.grid.coordinates import _xyz_to_lonlat_rad, _normalize_xyz, _xyz_to_lonlat_rad_no_norm,_normalize_xyz_scalar from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON from uxarray.utils.computing import isclose, cross, dot, allclose @@ -50,7 +50,7 @@ def _point_within_gca_body( GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON ): # If the pt and the GCA are on the same longitude (the y coordinates are the same) - if np.isclose( + if isclose( GCRv0_lonlat[0], pt_lonlat[0], rtol=MACHINE_EPSILON, atol=MACHINE_EPSILON ): # Now use the latitude to determine if the pt falls between the interval @@ -347,20 +347,20 @@ def extreme_gca_latitude(gca_cart, extreme_type): raise ValueError("extreme_type must be either 'max' or 'min'") n1, n2 = gca_cart - - dot_n1_n2 = dot(np.asarray(n1), np.asarray(n2)) + dot_n1_n2 = np.dot(n1, n2) denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom d_a_max = ( np.clip(d_a_max, 0, 1) - if isclose(d_a_max, 0, atol=ERROR_TOLERANCE) - or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) + if np.isclose(d_a_max, [0, 1], atol=MACHINE_EPSILON).any() else d_a_max ) - _, lat_n1 = _xyz_to_lonlat_rad_no_norm(n1[0], n1[1], n1[2]) - _, lat_n2 = _xyz_to_lonlat_rad_no_norm(n2[0], n2[1], n2[2]) + # Before we make sure the grid coordinates are normalized, do not try to skip the normalization steps! + _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) + _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) + if 0 < d_a_max < 1: node3 = (1 - d_a_max) * n1 + d_a_max * n2 @@ -374,3 +374,4 @@ def extreme_gca_latitude(gca_cart, extreme_type): ) else: return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) + diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 1af7ed73f..5bb808097 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -1064,8 +1064,6 @@ def _populate_bounds( is_GCA_list=is_GCA_list, ) - if temp_latlon_array[face_idx][0][0] == temp_latlon_array[face_idx][0][1]: - pass assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1] assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1] From 0111b0e101c55559a7fab87ff9c47d797f97a3c9 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 12:56:50 -0700 Subject: [PATCH 08/24] pre-commit --- test/test_geometry.py | 3 ++- uxarray/grid/arcs.py | 8 +++++--- uxarray/grid/geometry.py | 3 --- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/test/test_geometry.py b/test/test_geometry.py index 21b85cbe5..0f25f6f07 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -1413,7 +1413,8 @@ def test_populate_bounds_GCAList_mix(self): class TestLatlonBoundsFiles: def test_face_bounds(self): - """Test to ensure ``Grid.face_bounds`` works correctly for all grid files.""" + """Test to ensure ``Grid.face_bounds`` works correctly for all grid + files.""" for grid_path in grid_files_latlonBound: try: # Open the grid file diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 41f3dc15e..94f9b0fd2 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -2,7 +2,11 @@ # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place -from uxarray.grid.coordinates import _xyz_to_lonlat_rad, _normalize_xyz, _xyz_to_lonlat_rad_no_norm,_normalize_xyz_scalar +from uxarray.grid.coordinates import ( + _xyz_to_lonlat_rad, + _xyz_to_lonlat_rad_no_norm, + _normalize_xyz_scalar, +) from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON from uxarray.utils.computing import isclose, cross, dot, allclose @@ -361,7 +365,6 @@ def extreme_gca_latitude(gca_cart, extreme_type): _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) - if 0 < d_a_max < 1: node3 = (1 - d_a_max) * n1 + d_a_max * n2 node3 = np.array(_normalize_xyz_scalar(node3[0], node3[1], node3[2])) @@ -374,4 +377,3 @@ def extreme_gca_latitude(gca_cart, extreme_type): ) else: return max(lat_n1, lat_n2) if extreme_type == "max" else min(lat_n1, lat_n2) - diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 5bb808097..00757d68e 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -7,7 +7,6 @@ _get_lonlat_rad_face_edge_nodes, ) -from uxarray.utils.computing import isclose, allclose, all import warnings import pandas as pd import xarray as xr @@ -1064,7 +1063,6 @@ def _populate_bounds( is_GCA_list=is_GCA_list, ) - assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1] assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1] lat_array = temp_latlon_array[face_idx][0] @@ -1099,7 +1097,6 @@ def _populate_bounds( grid._ds["bounds"] = bounds - def _construct_hole_edge_indices(edge_face_connectivity): """Index the missing edges on a partial grid with holes, that is a region of the grid that is not covered by any geometry.""" From fd26c4298079e92e3078c88371c41b64bc1b4d40 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Sun, 15 Sep 2024 16:24:28 -0500 Subject: [PATCH 09/24] Update asv.conf.json --- benchmarks/asv.conf.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/asv.conf.json b/benchmarks/asv.conf.json index 31a43921d..a049ddf43 100644 --- a/benchmarks/asv.conf.json +++ b/benchmarks/asv.conf.json @@ -94,7 +94,7 @@ "setuptools_scm": [""], "xarray": [""], "netcdf4": [""], - "pip+pyfma": [""] + "pip+git+https://github.com/philipc2/pyfma.git: [""] }, From d6e9416c60cf5dd7f2f93f3d13eff998ae434ed7 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Sun, 15 Sep 2024 16:25:39 -0500 Subject: [PATCH 10/24] Update asv.conf.json --- benchmarks/asv.conf.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/asv.conf.json b/benchmarks/asv.conf.json index a049ddf43..b108ec396 100644 --- a/benchmarks/asv.conf.json +++ b/benchmarks/asv.conf.json @@ -94,7 +94,7 @@ "setuptools_scm": [""], "xarray": [""], "netcdf4": [""], - "pip+git+https://github.com/philipc2/pyfma.git: [""] + "pip+git+https://github.com/philipc2/pyfma.git": [""] }, From bde08545e159299b056e3887396932cb7af0fc6e Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 14:28:01 -0700 Subject: [PATCH 11/24] Update uxarray/grid/intersections.py Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> --- uxarray/grid/intersections.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 7324338c3..afcca1fcd 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -7,7 +7,7 @@ from uxarray.utils.computing import cross_fma -def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False): +def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): """Calculate the intersection point(s) of two Great Circle Arcs (GCAs) in a Cartesian coordinate system. From e2f60f75ba54e1584e0b6db466c39e418b2449bf Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 20:40:35 -0700 Subject: [PATCH 12/24] enable `JIT_CACHE` --- uxarray/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/constants.py b/uxarray/constants.py index 606b6ab2c..d55e25b49 100644 --- a/uxarray/constants.py +++ b/uxarray/constants.py @@ -16,7 +16,7 @@ # error tolerance, mainly in the intersection calculations. MACHINE_EPSILON = np.float64(np.finfo(float).eps) -ENABLE_JIT_CACHE = False +ENABLE_JIT_CACHE = True ENABLE_JIT = True ENABLE_FMA = False From 3808824eecccc87d16d8ffef47a335c45eb9b745 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 22:11:52 -0700 Subject: [PATCH 13/24] replace using the numba decorate function --- uxarray/grid/arcs.py | 9 +++++---- uxarray/grid/geometry.py | 26 +++++++++++++------------- uxarray/grid/intersections.py | 22 +++++++++++----------- 3 files changed, 29 insertions(+), 28 deletions(-) diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index 94f9b0fd2..bb161bf7b 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -102,7 +102,7 @@ def _point_within_gca_body( ): non_pole_endpoint = GCRv1_lonlat - if non_pole_endpoint is not None and not np.isclose( + if non_pole_endpoint is not None and not isclose( non_pole_endpoint[0], pt_lonlat[0], rtol=ERROR_TOLERANCE, atol=0.0 ): return False @@ -351,16 +351,17 @@ def extreme_gca_latitude(gca_cart, extreme_type): raise ValueError("extreme_type must be either 'max' or 'min'") n1, n2 = gca_cart - dot_n1_n2 = np.dot(n1, n2) + + dot_n1_n2 = dot(np.asarray(n1), np.asarray(n2)) denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0) d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom d_a_max = ( np.clip(d_a_max, 0, 1) - if np.isclose(d_a_max, [0, 1], atol=MACHINE_EPSILON).any() + if isclose(d_a_max, 0, atol=ERROR_TOLERANCE) + or isclose(d_a_max, 1, atol=ERROR_TOLERANCE) else d_a_max ) - # Before we make sure the grid coordinates are normalized, do not try to skip the normalization steps! _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index 00757d68e..dc6b5665e 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -6,7 +6,7 @@ _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes, ) - +from uxarray.utils.computing import allclose, isclose import warnings import pandas as pd import xarray as xr @@ -537,7 +537,7 @@ def _check_intersection(ref_edge, edges): # for each intersection point, check if it is a pole point for point in intersection_point: - if np.allclose(point, pole_point, atol=ERROR_TOLERANCE): + if allclose(point, pole_point, atol=ERROR_TOLERANCE): return True intersection_points.append(point) @@ -547,9 +547,9 @@ def _check_intersection(ref_edge, edges): # If the unique intersection point is one and it is exactly one of the nodes of the face, return 0 if len(intersection_points) == 1: for edge in edges: - if np.allclose( + if allclose( intersection_points[0], edge[0], atol=ERROR_TOLERANCE - ) or np.allclose(intersection_points[0], edge[1], atol=ERROR_TOLERANCE): + ) or allclose(intersection_points[0], edge[1], atol=ERROR_TOLERANCE): return 0 return len(intersection_points) @@ -693,18 +693,18 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True): # Check for pole points and update latitudes is_pole_point = ( lon_pt == INT_FILL_VALUE - and np.isclose( - new_pt[0], [0.5 * np.pi, -0.5 * np.pi], atol=ERROR_TOLERANCE + and isclose( + new_pt[0], np.asarray([0.5 * np.pi, -0.5 * np.pi]), atol=ERROR_TOLERANCE ).any() ) if is_pole_point: # Check if the new point is close to the North Pole - if np.isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE): + if isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE): latlon_box[0][1] = 0.5 * np.pi # Check if the new point is close to the South Pole - elif np.isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE): + elif isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE): latlon_box[0][0] = -0.5 * np.pi return latlon_box @@ -850,7 +850,7 @@ def _populate_face_latlon_bound( ) # Check if the node matches the pole point or if the pole point is within the edge_cart - if np.allclose( + if allclose( n1_cart, pole_point, atol=ERROR_TOLERANCE ) or point_within_gca( pole_point, np.array([n1_cart, n2_cart]), is_directed=False @@ -931,16 +931,16 @@ def _populate_face_latlon_bound( ) # Insert extreme latitude points into the latlonbox if they differ from the node latitudes - if not np.isclose( + if not isclose( node1_lat_rad, lat_max, atol=ERROR_TOLERANCE - ) and not np.isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE): + ) and not isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE): # Insert the maximum latitude face_latlon_array = _insert_pt_in_latlonbox( face_latlon_array, np.array([lat_max, node1_lon_rad]) ) - elif not np.isclose( + elif not isclose( node1_lat_rad, lat_min, atol=ERROR_TOLERANCE - ) and not np.isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE): + ) and not isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE): # Insert the minimum latitude face_latlon_array = _insert_pt_in_latlonbox( face_latlon_array, np.array([lat_min, node1_lon_rad]) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index afcca1fcd..7a5625829 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -4,7 +4,7 @@ from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform import warnings -from uxarray.utils.computing import cross_fma +from uxarray.utils.computing import cross_fma, allclose, cross, dot def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): @@ -65,29 +65,29 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): ) # Check perpendicularity conditions and floating-point arithmetic limitations - if not np.allclose( - np.dot(w0w1_norm, w0), 0, atol=MACHINE_EPSILON - ) or not np.allclose(np.dot(w0w1_norm, w1), 0, atol=MACHINE_EPSILON): + if not allclose( + dot(w0w1_norm, w0), 0.0, atol=MACHINE_EPSILON + ) or not allclose(dot(w0w1_norm, w1), 0.0, atol=MACHINE_EPSILON): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution." ) - if not np.allclose( - np.dot(v0v1_norm, v0), 0, atol=MACHINE_EPSILON - ) or not np.allclose(np.dot(v0v1_norm, v1), 0, atol=MACHINE_EPSILON): + if not allclose( + dot(v0v1_norm, v0), 0.0,0.0, atol=MACHINE_EPSILON + ) or not allclose(dot(v0v1_norm, v1), 0.0, atol=MACHINE_EPSILON): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " ) - if not np.allclose( - np.dot(cross_norms, v0v1_norm), 0, atol=MACHINE_EPSILON - ) or not np.allclose(np.dot(cross_norms, w0w1_norm), 0, atol=MACHINE_EPSILON): + if not allclose( + dot(cross_norms, v0v1_norm), 0.0, atol=MACHINE_EPSILON + ) or not allclose(dot(cross_norms, w0w1_norm), 0.0, atol=MACHINE_EPSILON): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " ) # If the cross_norms is zero, the two GCAs are parallel - if np.allclose(cross_norms, 0, atol=MACHINE_EPSILON): + if allclose(cross_norms, 0.0, atol=MACHINE_EPSILON): res = [] # Check if the two GCAs are overlapping if point_within_gca(v0, [w0, w1]): From 5531e3d00df21e0d58277ccae6af0713906789b0 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Sun, 15 Sep 2024 22:12:34 -0700 Subject: [PATCH 14/24] pre-commit --- uxarray/grid/geometry.py | 4 +--- uxarray/grid/intersections.py | 14 +++++++------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index dc6b5665e..ac3624e38 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -850,9 +850,7 @@ def _populate_face_latlon_bound( ) # Check if the node matches the pole point or if the pole point is within the edge_cart - if allclose( - n1_cart, pole_point, atol=ERROR_TOLERANCE - ) or point_within_gca( + if allclose(n1_cart, pole_point, atol=ERROR_TOLERANCE) or point_within_gca( pole_point, np.array([n1_cart, n2_cart]), is_directed=False ): is_center_pole = False diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 7a5625829..b1ff096eb 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -4,7 +4,7 @@ from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform import warnings -from uxarray.utils.computing import cross_fma, allclose, cross, dot +from uxarray.utils.computing import cross_fma, allclose, dot def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): @@ -65,16 +65,16 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): ) # Check perpendicularity conditions and floating-point arithmetic limitations - if not allclose( - dot(w0w1_norm, w0), 0.0, atol=MACHINE_EPSILON - ) or not allclose(dot(w0w1_norm, w1), 0.0, atol=MACHINE_EPSILON): + if not allclose(dot(w0w1_norm, w0), 0.0, atol=MACHINE_EPSILON) or not allclose( + dot(w0w1_norm, w1), 0.0, atol=MACHINE_EPSILON + ): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution." ) - if not allclose( - dot(v0v1_norm, v0), 0.0,0.0, atol=MACHINE_EPSILON - ) or not allclose(dot(v0v1_norm, v1), 0.0, atol=MACHINE_EPSILON): + if not allclose(dot(v0v1_norm, v0), 0.0, 0.0, atol=MACHINE_EPSILON) or not allclose( + dot(v0v1_norm, v1), 0.0, atol=MACHINE_EPSILON + ): warnings.warn( "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. " ) From 79d035c589f78e488df593a6659a46139c9b56dd Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 12:34:06 -0500 Subject: [PATCH 15/24] Update uxarray/grid/integrate.py --- uxarray/grid/integrate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 74629b5d9..f329a95be 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -93,7 +93,7 @@ def _get_zonal_faces_weight_at_constLat( # Iterate through all faces and their edges for face_index, face_edges in enumerate(faces_edges_cart_candidate): # Remove the Int_fill_value from the face_edges - face_edges = face_edges[np.all(face_edges != INT_FILL_VALUE, axis=(1, 2))] + face_edges = face_edges[all(face_edges != INT_FILL_VALUE, axis=(1, 2))] if is_face_GCA_list is not None: is_GCA_list = is_face_GCA_list[face_index] else: From e01453b94ee09e5c545dfe36e4d2fb31eab5b611 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 12:34:14 -0500 Subject: [PATCH 16/24] Update uxarray/grid/intersections.py --- uxarray/grid/intersections.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index b1ff096eb..31309572a 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -48,9 +48,9 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): # Compute normals and orthogonal bases using FMA if fma_disabled: - w0w1_norm = np.cross(w0, w1) - v0v1_norm = np.cross(v0, v1) - cross_norms = np.cross(w0w1_norm, v0v1_norm) + w0w1_norm = cross(w0, w1) + v0v1_norm = cross(v0, v1) + cross_norms = cross(w0w1_norm, v0v1_norm) else: w0w1_norm = cross_fma(w0, w1) v0v1_norm = cross_fma(v0, v1) From b203dbe8739d1e70c062480ef063467547a65a02 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 12:34:19 -0500 Subject: [PATCH 17/24] Update uxarray/grid/intersections.py --- uxarray/grid/intersections.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 31309572a..9370fc40e 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -184,7 +184,7 @@ def gca_constLat_intersection( return np.array([]) if fma_disabled: - n = np.cross(x1, x2) + n = cross(x1, x2) else: # Raise a warning for Windows users From a9e823abf7041d5ff588becba7101359eb7504eb Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 12:34:28 -0500 Subject: [PATCH 18/24] Update uxarray/grid/integrate.py --- uxarray/grid/integrate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index f329a95be..0748f2fa9 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -73,7 +73,7 @@ def _get_zonal_faces_weight_at_constLat( # The -90 and 90 treatment is hard-coded in the function, based on: # If a pole point is inside a face, then this face's value is the only value that should be considered. # If the pole point is not inside any face, then it's on the boundary of faces around it, so their weights are even. - if np.isclose(latitude_cart, 1, atol=ERROR_TOLERANCE) or np.isclose( + if isclose(latitude_cart, 1, atol=ERROR_TOLERANCE) or isclose( latitude_cart, -1, atol=ERROR_TOLERANCE ): # Now all candidate faces( the faces around the pole) are considered as the same weight From b53a6b3f90471e37a0361e712fb19d0ec3486f67 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 12:55:15 -0500 Subject: [PATCH 19/24] add cross to intersections.py --- benchmarks/face_bounds.py | 13 +++++++++++++ uxarray/grid/intersections.py | 2 +- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/benchmarks/face_bounds.py b/benchmarks/face_bounds.py index b249e7b99..122ab5901 100644 --- a/benchmarks/face_bounds.py +++ b/benchmarks/face_bounds.py @@ -28,3 +28,16 @@ def time_face_bounds(self, grid_path): def peakmem_face_bounds(self, grid_path): """Peak memory usage obtain ``Grid.face_bounds.""" face_bounds = self.uxgrid.bounds + +from asv_runner.benchmarks.mark import skip_benchmark_if, timeout_class_at +@timeout_class_at(1000) +class Bounds: + def setup(self): + self.uxgrid = ux.open_grid(r"C:\Users\chmie\PycharmProjects\ncar-uxarray\uxarray-hongyu\benchmarks\oQU120.grid.nc") + + def teardown(self): + del self.uxgrid + + def time_face_bounds(self): + """Time to obtain ``Grid.face_bounds``""" + self.uxgrid.bounds diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 9370fc40e..dadf8caf3 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -4,7 +4,7 @@ from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform import warnings -from uxarray.utils.computing import cross_fma, allclose, dot +from uxarray.utils.computing import cross_fma, allclose, dot, cross def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): From bcb11e7b1ca8035878a84204450254cc3ca562c6 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 13:00:18 -0500 Subject: [PATCH 20/24] add isclose to integrate.py --- uxarray/grid/integrate.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index 0748f2fa9..d138f5325 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -4,6 +4,8 @@ from uxarray.grid.coordinates import _xyz_to_lonlat_rad import pandas as pd +from uxarray.utils.computing import isclose + DUMMY_EDGE_VALUE = [INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE] From 28f1acb8b5c072e4b0a0633eebcf14becab31afd Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 13:10:14 -0500 Subject: [PATCH 21/24] fix all() failing with params --- uxarray/grid/integrate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uxarray/grid/integrate.py b/uxarray/grid/integrate.py index d138f5325..a562fc2ca 100644 --- a/uxarray/grid/integrate.py +++ b/uxarray/grid/integrate.py @@ -95,7 +95,7 @@ def _get_zonal_faces_weight_at_constLat( # Iterate through all faces and their edges for face_index, face_edges in enumerate(faces_edges_cart_candidate): # Remove the Int_fill_value from the face_edges - face_edges = face_edges[all(face_edges != INT_FILL_VALUE, axis=(1, 2))] + face_edges = face_edges[np.all(face_edges != INT_FILL_VALUE, axis=(1, 2))] if is_face_GCA_list is not None: is_GCA_list = is_face_GCA_list[face_index] else: From 88a40fcea08909dbd44b655de9ab171c0b7a9c1c Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 13:20:15 -0500 Subject: [PATCH 22/24] remove local benchmark --- benchmarks/face_bounds.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/benchmarks/face_bounds.py b/benchmarks/face_bounds.py index 122ab5901..b249e7b99 100644 --- a/benchmarks/face_bounds.py +++ b/benchmarks/face_bounds.py @@ -28,16 +28,3 @@ def time_face_bounds(self, grid_path): def peakmem_face_bounds(self, grid_path): """Peak memory usage obtain ``Grid.face_bounds.""" face_bounds = self.uxgrid.bounds - -from asv_runner.benchmarks.mark import skip_benchmark_if, timeout_class_at -@timeout_class_at(1000) -class Bounds: - def setup(self): - self.uxgrid = ux.open_grid(r"C:\Users\chmie\PycharmProjects\ncar-uxarray\uxarray-hongyu\benchmarks\oQU120.grid.nc") - - def teardown(self): - del self.uxgrid - - def time_face_bounds(self): - """Time to obtain ``Grid.face_bounds``""" - self.uxgrid.bounds From ef9bc27a9f725bf5e9930903a992939b609bd9e4 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 17:53:18 -0500 Subject: [PATCH 23/24] optimize xyz_to_latlon_rad_scalar --- benchmarks/face_bounds.py | 10 +++++++ uxarray/grid/arcs.py | 19 ++++++++----- uxarray/grid/coordinates.py | 53 +++++++++++++++++++++++++++++++++-- uxarray/grid/intersections.py | 4 +-- uxarray/utils/computing.py | 11 ++++++++ 5 files changed, 86 insertions(+), 11 deletions(-) diff --git a/benchmarks/face_bounds.py b/benchmarks/face_bounds.py index b249e7b99..db3106fcb 100644 --- a/benchmarks/face_bounds.py +++ b/benchmarks/face_bounds.py @@ -28,3 +28,13 @@ def time_face_bounds(self, grid_path): def peakmem_face_bounds(self, grid_path): """Peak memory usage obtain ``Grid.face_bounds.""" face_bounds = self.uxgrid.bounds + +class Bounds: + def setup(self): + self.uxgrid = ux.open_grid(r"C:\Users\chmie\PycharmProjects\ncar-uxarray\uxarray-hongyu\benchmarks\oQU480.grid.nc") + + def teardown(self): + del self.uxgrid + + def time_bounds(self): + self.uxgrid.bounds diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index bb161bf7b..2bffa2177 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -3,8 +3,7 @@ # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place from uxarray.grid.coordinates import ( - _xyz_to_lonlat_rad, - _xyz_to_lonlat_rad_no_norm, + _xyz_to_lonlat_rad_scalar, _normalize_xyz_scalar, ) from uxarray.constants import ERROR_TOLERANCE, MACHINE_EPSILON @@ -217,12 +216,18 @@ def point_within_gca(pt, gca_cart, is_directed=False): Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons. """ # Convert the cartesian coordinates to lonlat coordinates - pt_lonlat = np.array(_xyz_to_lonlat_rad_no_norm(pt[0], pt[1], pt[2])) + pt_lonlat = np.array( + _xyz_to_lonlat_rad_scalar(pt[0], pt[1], pt[2], normalize=False) + ) GCRv0_lonlat = np.array( - _xyz_to_lonlat_rad_no_norm(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2]) + _xyz_to_lonlat_rad_scalar( + gca_cart[0][0], gca_cart[0][1], gca_cart[0][2], normalize=False + ) ) GCRv1_lonlat = np.array( - _xyz_to_lonlat_rad_no_norm(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2]) + _xyz_to_lonlat_rad_scalar( + gca_cart[1][0], gca_cart[1][1], gca_cart[1][2], normalize=False + ) ) gca_cart = np.asarray(gca_cart) @@ -363,8 +368,8 @@ def extreme_gca_latitude(gca_cart, extreme_type): else d_a_max ) # Before we make sure the grid coordinates are normalized, do not try to skip the normalization steps! - _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2]) - _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2]) + _, lat_n1 = _xyz_to_lonlat_rad_scalar(n1[0], n1[1], n1[2], normalize=True) + _, lat_n2 = _xyz_to_lonlat_rad_scalar(n2[0], n2[1], n2[2], normalize=True) if 0 < d_a_max < 1: node3 = (1 - d_a_max) * n1 + d_a_max * n2 diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index d34f0fc20..f25553321 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -68,6 +68,55 @@ def _xyz_to_lonlat_rad_no_norm( return lon, lat +@njit +def _xyz_to_lonlat_rad_scalar( + x: Union[np.ndarray, float], + y: Union[np.ndarray, float], + z: Union[np.ndarray, float], + normalize: bool = True, +): + """Converts a Cartesian x,y,z coordinates into Spherical latitude and + longitude without normalization, decorated with Numba. + + Parameters + ---------- + x : float + Cartesian x coordinate + y: float + Cartesiain y coordinate + z: float + Cartesian z coordinate + + + Returns + ------- + lon : float + Longitude in radians + lat: float + Latitude in radians + """ + + if normalize: + x, y, z = _normalize_xyz_scalar(x, y, z) + denom = np.abs(x * x + y * y + z * z) + x /= denom + y /= denom + z /= denom + + lon = math.atan2(y, x) + lat = math.asin(z) + + # set longitude range to [0, pi] + lon = np.mod(lon, 2 * np.pi) + + z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE + + lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat) + lon = np.where(z_mask, 0.0, lon) + + return lon, lat + + def _xyz_to_lonlat_rad( x: Union[np.ndarray, float], y: Union[np.ndarray, float], @@ -103,8 +152,8 @@ def _xyz_to_lonlat_rad( y /= denom z /= denom - lon = np.arctan2(y, x, dtype=np.float64) - lat = np.arcsin(z, dtype=np.float64) + lon = np.arctan2(y, x) + lat = np.arcsin(z) # set longitude range to [0, pi] lon = np.mod(lon, 2 * np.pi) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index dadf8caf3..a26136e8f 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -4,7 +4,7 @@ from uxarray.grid.arcs import point_within_gca, extreme_gca_latitude, in_between import platform import warnings -from uxarray.utils.computing import cross_fma, allclose, dot, cross +from uxarray.utils.computing import cross_fma, allclose, dot, cross, norm def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): @@ -99,7 +99,7 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True): return np.array(res) # Normalize the cross_norms - cross_norms = cross_norms / np.linalg.norm(cross_norms) + cross_norms = cross_norms / norm(cross_norms) x1 = cross_norms x2 = -x1 diff --git a/uxarray/utils/computing.py b/uxarray/utils/computing.py index f9929deb4..6baa25313 100644 --- a/uxarray/utils/computing.py +++ b/uxarray/utils/computing.py @@ -61,6 +61,17 @@ def dot(a, b): return np.dot(a, b) +@njit +def norm(x): + """Numba decorated implementation of ``np.linalg.norm()`` + + See Also + -------- + numpy.linalg.norm + """ + return np.linalg.norm(x) + + def _fmms(a, b, c, d): """ Calculate the difference of products using the FMA (fused multiply-add) operation: (a * b) - (c * d). From 6d209a432a657f7a43e67c32eea4d2f5a04055d5 Mon Sep 17 00:00:00 2001 From: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Date: Mon, 16 Sep 2024 20:16:12 -0500 Subject: [PATCH 24/24] cleanup exception handling --- uxarray/grid/intersections.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index a26136e8f..81829a68a 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -230,9 +230,7 @@ def gca_constLat_intersection( else np.vstack((res, converged_pt)) ) except RuntimeError: - print(f"Error encountered with initial guess: {p1}") - print(f"gca_cart: {gca_cart}") - raise + raise RuntimeError(f"Error encountered with initial guess: {p1}") if point_within_gca(p2, gca_cart, is_directed=is_directed): try: @@ -253,8 +251,6 @@ def gca_constLat_intersection( else np.vstack((res, converged_pt)) ) except RuntimeError: - print(f"Error encountered with initial guess: {p2}") - print(f"gca_cart: {gca_cart}") - raise + raise RuntimeError(f"Error encountered with initial guess: {p2}") return res if res is not None else np.array([])