Robust implementation of _point_within_gca_cartesian using only Cartesian coordinates (#1112)

* initial commit

* initial implementation

* initial testcase setup

* required fix for ` _angle_of_2_vectors`

* update ` _angle_of_2_vectors` to return degree in 0~360

* update `_pt_on_gca`

* `test_populate_bounds_normal` failed when using njit

* update _populate_face_latlon_bound

* fix tests and add test case for case at pole

* add test case exactly at pole

* add condition to check if point is almost exactly equal to one of the end points

* add general north and south pole tests

* add tests for north and south pole

* add extra test case for north and south pole

* fix assert

* add edge case test

* Add error tolerance for values 0.0

* run pre-commit

* remove debugging code and use njit

---------

Co-authored-by: Hongyu Chen <hyvchen@ucdavis.edu>

Author: philipc2

test/test_arcs.py

@@ -10,7 +10,7 @@
 import uxarray as ux
 
 from uxarray.grid.coordinates import _lonlat_rad_to_xyz
-from uxarray.grid.arcs import point_within_gca, _point_within_gca_cartesian
+from uxarray.grid.arcs import point_within_gca
 
 try:
     import constants
@@ -31,65 +31,46 @@ class TestIntersectionPoint(TestCase):
     def test_pt_within_gcr(self):
         # The GCR that's eexactly 180 degrees will have Value Error raised
 
-        gcr_180degree_cart = [
+        gcr_180degree_cart = np.asarray([
             _lonlat_rad_to_xyz(0.0, np.pi / 2.0),
             _lonlat_rad_to_xyz(0.0, -np.pi / 2.0)
-        ]
+        ])
+        pt_same_lon_in = np.asarray(_lonlat_rad_to_xyz(0.0, 0.0))
 
-        pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
         with self.assertRaises(ValueError):
-            _point_within_gca_cartesian(pt_same_lon_in, gcr_180degree_cart)
-
+            point_within_gca(pt_same_lon_in, gcr_180degree_cart[0],gcr_180degree_cart[1] )
+        #
         # Test when the point and the GCA all have the same longitude
-        gcr_same_lon_cart = [
+        gcr_same_lon_cart = np.asarray([
             _lonlat_rad_to_xyz(0.0, 1.5),
             _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_cartesian(pt_same_lon_in, gcr_same_lon_cart))
+        ])
+        pt_same_lon_in = np.asarray(_lonlat_rad_to_xyz(0.0, 0.0))
+        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart[0], gcr_same_lon_cart[1]))
 
-        pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
-        res = _point_within_gca_cartesian(pt_same_lon_out, gcr_same_lon_cart)
+        pt_same_lon_out = np.asarray(_lonlat_rad_to_xyz(0.0, 1.5000001))
+        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart[0], gcr_same_lon_cart[1])
         self.assertFalse(res)
 
-        pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0)
-        res = _point_within_gca_cartesian(pt_same_lon_out_2, gcr_same_lon_cart)
+        pt_same_lon_out_2 = np.asarray(_lonlat_rad_to_xyz(0.1, 1.0))
+        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart[0], gcr_same_lon_cart[1])
         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_cartesian(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
-        self.assertTrue(res)
-
-        # Normal case
-        # GCR vertex0 in radian : [1.3003315590159483, -0.007004587172323237],
-        # GCR vertex1 in radian : [3.5997458123873827, -1.4893379576608758]
-        # Point in radian : [1.3005410084914981, -0.010444274637648326]
-        gcr_cart_2 = np.array([[0.267, 0.963, -0.007], [-0.073, -0.036,
-                                                        -0.997]])
-        pt_cart_within = np.array(
-            [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
-        self.assertTrue(_point_within_gca_cartesian(pt_cart_within, gcr_cart_2, True))
-
-        # Test other more complicate cases : The anti-meridian case
-
     def test_pt_within_gcr_antimeridian(self):
         # GCR vertex0 in radian : [5.163808182822441, 0.6351384888657234],
         # GCR vertex1 in radian : [0.8280410325693055, 0.42237025187091526]
         # Point in radian : [0.12574759138415173, 0.770098701904903]
         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_cartesian(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
+        self.assertTrue(
+            point_within_gca(pt_cart, gcr_cart[0], gcr_cart[1]))
+
         gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
                                                           0.593]])
-        with self.assertRaises(ValueError):
-            _point_within_gca_cartesian(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_cartesian(pt_cart, gcr_cart_flip, is_directed=False))
+            point_within_gca(pt_cart, gcr_cart_flip[0], gcr_cart_flip[1]))
 
         # 2nd anti-meridian case
         # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
@@ -100,22 +81,9 @@ def test_pt_within_gcr_antimeridian(self):
         pt_cart_within = np.array(
             [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
         self.assertFalse(
-            _point_within_gca_cartesian(pt_cart_within, gcr_cart_1, is_directed=True))
-        self.assertFalse(
-            _point_within_gca_cartesian(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]
-        v2_rad = [2 * np.pi - 0.1, 0.0]
-        v1_cart = _lonlat_rad_to_xyz(v1_rad[0], v1_rad[1])
-        v2_cart = _lonlat_rad_to_xyz(v2_rad[0], v1_rad[1])
-        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_cartesian(pt_cart, gcr_cart, is_directed=True)
-        gcr_car_flipped = np.array([v2_cart, v1_cart])
-        self.assertTrue(
-            _point_within_gca_cartesian(pt_cart, gcr_car_flipped, is_directed=True))
+            point_within_gca(pt_cart_within, gcr_cart_1[0], gcr_cart_1[1]))
+
+
 
     def test_pt_within_gcr_cross_pole(self):
         gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]])
@@ -125,15 +93,4 @@ def test_pt_within_gcr_cross_pole(self):
         # Normalize the point abd the GCA
         pt_cart = pt_cart / np.linalg.norm(pt_cart)
         gcr_cart = np.array([x / np.linalg.norm(x) for x in gcr_cart])
-        self.assertTrue(_point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=False))
-
-        gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, -0.6]])
-        pt_cart = np.array(
-            [0.10, 0.0, 0.8])
-
-        # When the point is not within the GCA
-        pt_cart = pt_cart / np.linalg.norm(pt_cart)
-        gcr_cart = np.array([x / np.linalg.norm(x) for x in gcr_cart])
-        self.assertFalse(_point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=False))
-        with self.assertRaises(ValueError):
-            _point_within_gca_cartesian(pt_cart, gcr_cart, is_directed=True)
+        self.assertTrue(point_within_gca(pt_cart, gcr_cart[0], gcr_cart[1]))

test/test_helpers.py

@@ -238,6 +238,31 @@ def test_angle_of_2_vectors(self):
         v2 = np.array([1.0, 0.0, 0.0])
         self.assertAlmostEqual(_angle_of_2_vectors(v1, v2), 0.0)
 
+    def test_angle_of_2_vectors_180_degree(self):
+        GCR1_cart = np.array([
+            _lonlat_rad_to_xyz(np.deg2rad(0.0),
+                                    np.deg2rad(0.0)),
+            _lonlat_rad_to_xyz(np.deg2rad(181.0),
+                                    np.deg2rad(0.0))
+        ])
+
+        res = _angle_of_2_vectors( GCR1_cart[0],  GCR1_cart[1])
+
+        # The angle between the two vectors should be 181 degree
+        self.assertAlmostEqual(res, np.deg2rad(181.0), places=8)
+
+        GCR1_cart = np.array([
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(89.0)),
+            _lonlat_rad_to_xyz(np.deg2rad(170.0),
+                                    np.deg2rad(-10.0))
+        ])
+
+        res = _angle_of_2_vectors( GCR1_cart[0],  GCR1_cart[1])
+
+        # The angle between the two vectors should be 181 degree
+        self.assertAlmostEqual(res, np.deg2rad(89.0+10.0), places=8)
+
 
 class TestFaceEdgeConnectivityHelper(TestCase):
 

test/test_integrate.py

@@ -81,10 +81,9 @@ def test_get_faces_constLat_intersection_info_one_intersection(self):
         ])
 
         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)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 1
         self.assertEqual(len(unique_intersections), 1)
 
@@ -111,10 +110,10 @@ def test_get_faces_constLat_intersection_info_encompass_pole(self):
         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)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length should be no greater than 2* n_edges
         self.assertLessEqual(len(unique_intersections), 2*len(face_edges_cart))
 
@@ -133,10 +132,9 @@ def test_get_faces_constLat_intersection_info_on_pole(self):
              [-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)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 2)
 
@@ -156,10 +154,9 @@ def test_get_faces_constLat_intersection_info_near_pole(self):
         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)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 1)
 
@@ -182,10 +179,9 @@ def test_get_faces_constLat_intersection_info_2(self):
                                      [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)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 2)
 
@@ -208,10 +204,9 @@ def test_get_faces_constLat_intersection_info_2(self):
                                      [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)
+        unique_intersections, pt_lon_min, pt_lon_max = _get_faces_constLat_intersection_info(face_edges_cart, latitude_cart, is_GCA_list, is_latlonface)
         # The expected unique_intersections length is 2
         self.assertEqual(len(unique_intersections), 2)
 
@@ -240,8 +235,7 @@ def test_get_zonal_face_interval(self):
         # 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)
+                                                                                         0.4 * np.pi]]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -285,7 +279,7 @@ def test_get_zonal_face_interval_empty_interval(self):
             [3.14159265, 3.2321175]
         ])
 
-        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds, is_directed=False)
+        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds)
         expected_res = pd.DataFrame({"start": [0.0], "end": [0.0]})
         pd.testing.assert_frame_equal(res, expected_res)
 
@@ -322,7 +316,7 @@ def test_get_zonal_face_interval_encompass_pole(self):
         })
 
         # Call the function to get the result
-        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds, is_directed=False)
+        res = _get_zonal_face_interval(face_edges_cart, latitude_cart, face_latlon_bounds)
 
         # Assert the result matches the expected DataFrame
         pd.testing.assert_frame_equal(res,  expected_df)
@@ -349,8 +343,7 @@ def test_get_zonal_face_interval_FILL_VALUE(self):
         # 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)
+                                                                                         0.4 * np.pi]]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -383,7 +376,7 @@ def test_get_zonal_face_interval_GCA_constLat(self):
         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, is_GCA_list=np.array([True, False, True, False]))
+                                               is_GCA_list=np.array([True, False, True, False]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -415,7 +408,7 @@ def test_get_zonal_face_interval_equator(self):
         interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
                                                                                          0.4 * np.pi]]),
-                                               is_directed=False, is_GCA_list=np.array([True, True, True, True]))
+                                               is_GCA_list=np.array([True, True, True, True]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -434,7 +427,7 @@ def test_get_zonal_face_interval_equator(self):
         interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0,
                                                np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi,
                                                                                          0.4 * np.pi]]),
-                                               is_directed=False, is_GCA_list=np.array([True, False, True, False]))
+                                               is_GCA_list=np.array([True, False, True, False]))
         expected_interval_df = pd.DataFrame({
             'start': [1.6 * np.pi, 0.0],
             'end': [2.0 * np.pi, 00.4 * np.pi]
@@ -605,7 +598,7 @@ 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)
 
         nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
@@ -666,7 +659,7 @@ def test_get_zonal_faces_weight_at_constLat_regular(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
-        ]), np.sin(0.1 * np.pi), latlon_bounds, is_directed=False)
+        ]), np.sin(0.1 * np.pi), latlon_bounds)
 
         nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
 
@@ -693,7 +686,7 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_one_face(self):
         ])
         constLat_cart = -1
 
-        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds, is_directed=False)
+        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
         # Define the expected DataFrame
         expected_weight_df = pd.DataFrame({"face_index": [0], "weight": [1.0]})
 
@@ -740,7 +733,7 @@ def test_get_zonal_faces_weight_at_constLat_on_pole_faces(self):
 
         constLat_cart = 1.0
 
-        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds, is_directed=False)
+        weight_df = _get_zonal_faces_weight_at_constLat(face_edges_cart, constLat_cart, face_bounds)
         # Define the expected DataFrame
         expected_weight_df = pd.DataFrame({
             'face_index': [0, 1, 2, 3],
@@ -774,7 +767,7 @@ def test_get_zonal_face_interval_pole(self):
         ])
         constLat_cart = -0.9986295347545738
 
-        weight_df = _get_zonal_face_interval(face_edges_cart, constLat_cart, face_bounds, is_directed=False)
+        weight_df = _get_zonal_face_interval(face_edges_cart, constLat_cart, face_bounds)
         # No Nan values should be present in the weight_df
         self.assertFalse(weight_df.isnull().values.any())
 
@@ -827,7 +820,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_latlonface=True)
 
 
         nt.assert_array_almost_equal(weight_df, expected_weight_df, decimal=3)
@@ -839,4 +832,4 @@ def test_get_zonal_faces_weight_at_constLat_latlonface(self):
         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)
+            ]), np.deg2rad(20), latlon_bounds)