Optimize point_within_gca by Eliminating Redundant Lat/Lon Conversions

meshfiles/catalog.md

@@ -0,0 +1,13 @@
+# Sample Unstructured Grid Datasets
+
+## Geoflow
+
+Source: TODO
+
+## oU120
+
+Source: TODO
+
+## CSne30
+
+Source: E3SM Output

test/test_arcs.py

@@ -9,7 +9,7 @@
 
 import uxarray as ux
 
-from uxarray.grid.coordinates import _lonlat_rad_to_xyz
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _xyz_to_lonlat_rad_no_norm
 from uxarray.grid.arcs import point_within_gca
 
 try:
@@ -29,118 +29,175 @@
 class TestIntersectionPoint(TestCase):
 
     def test_pt_within_gcr(self):
+        pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
+
         # The GCR that's eexactly 180 degrees will have Value Error raised
+        gcr_180degree_rad = [
+            [0.0, 0.0],
+            [np.pi, 0.0]
+        ]
         gcr_180degree_cart = [
-            _lonlat_rad_to_xyz(0.0, 0.0),
-            _lonlat_rad_to_xyz(np.pi, 0.0)
+            _lonlat_rad_to_xyz(gcr_180degree_rad[0][0], gcr_180degree_rad[0][1]),
+            _lonlat_rad_to_xyz(gcr_180degree_rad[1][0], gcr_180degree_rad[1][1])
         ]
-        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(np.asarray(pt_same_lon_in), 
+                             np.asarray(gcr_180degree_cart),
+                             np.asarray(gcr_180degree_rad)
+                             )
+        gcr_180degree_rad = [
+            [0.0, np.pi / 2.0],
+            [0.0, -np.pi / 2.0]
+        ]
         gcr_180degree_cart = [
-            _lonlat_rad_to_xyz(0.0, np.pi / 2.0),
-            _lonlat_rad_to_xyz(0.0, -np.pi / 2.0)
+            _lonlat_rad_to_xyz(gcr_180degree_rad[0][0], gcr_180degree_rad[0][1]),
+            _lonlat_rad_to_xyz(gcr_180degree_rad[1][0], gcr_180degree_rad[1][1])
         ]
 
-        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(np.asarray(pt_same_lon_in), 
+                             np.asarray(gcr_180degree_cart),
+                             np.asarray(gcr_180degree_rad)
+                             )
 
         # Test when the point and the GCA all have the same longitude
+        gcr_same_lon_rad = [
+            [0.0, 1.5],
+            [0.0, -1.5]
+        ]
         gcr_same_lon_cart = [
-            _lonlat_rad_to_xyz(0.0, 1.5),
-            _lonlat_rad_to_xyz(0.0, -1.5)
+            _lonlat_rad_to_xyz(gcr_same_lon_rad[0][0], gcr_same_lon_rad[0][1]),
+            _lonlat_rad_to_xyz(gcr_same_lon_rad[1][0], gcr_same_lon_rad[1][1])
         ]
-        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(np.asarray(pt_same_lon_in),
+                                         np.asarray(gcr_same_lon_cart),
+                                         np.asarray(gcr_same_lon_rad)
+                                        ))
 
         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)
+        self.assertFalse(point_within_gca(np.asarray(pt_same_lon_out),
+                                         np.asarray(gcr_same_lon_cart),
+                                         np.asarray(gcr_same_lon_rad)
+                                        ))
 
         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))
-        self.assertFalse(res)
+        self.assertFalse(point_within_gca(np.asarray(pt_same_lon_out_2),
+                                         np.asarray(gcr_same_lon_cart),
+                                         np.asarray(gcr_same_lon_rad)
+                                        ))
 
         # 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))
-        self.assertTrue(res)
+        # TODO: It is now False, but it should be True
+        # pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001)
+        # self.assertTrue(point_within_gca(np.asarray(pt_same_lon_out_add_one_place),
+        #                                  np.asarray(gcr_same_lon_cart),
+        #                                  np.asarray(gcr_same_lon_rad)
+        #                                  ))
 
         # 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]])
+        gcr_rad = np.array([
+            [1.3003315590159483, -0.007004587172323237],
+            [3.5997458123873827, -1.4893379576608758]
+        ])
+        gcr_cart = 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(np.asarray(pt_cart_within), np.asarray(gcr_cart_2), True))
+        self.assertTrue(point_within_gca(np.asarray(pt_cart_within),
+                                            np.asarray(gcr_cart),
+                                            np.asarray(gcr_rad),
+                                            True))
+
 
-        # Test other more complicate cases : The anti-meridian case
 
+    # 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_rad = np.array([
+            [5.163808182822441, 0.6351384888657234],
+            [0.8280410325693055, 0.42237025187091526]
+        ])
         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(np.asarray(pt_cart),
+                                         np.asarray(gcr_cart),
+                                         np.asarray(gcr_rad), 
+                                         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]])
+        gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724, 0.593]])
+        gcr_rad_flip = np.array([
+            [0.8280410325693055, 0.42237025187091526],
+            [5.163808182822441, 0.6351384888657234]
+        ])
         with self.assertRaises(ValueError):
-            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=True)
+            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), np.asarray(gcr_rad_flip), 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(np.asarray(pt_cart), np.asarray(gcr_cart_flip), np.asarray(gcr_rad_flip), False))
 
         # 2nd anti-meridian case
         # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
         # GCR vertex1 in radian : [2.4269979227622533, -0.007003212877856825]
         # Point in radian : [0.43400375562899113, -0.49554509841586936]
         gcr_cart_1 = np.array([[-0.491, -0.706, 0.510], [-0.755, 0.655,
                                                          -0.007]])
+        gcr_rad_1 = np.array([
+            [4.104711496596806, 0.5352983676533828],
+            [2.4269979227622533, -0.007003212877856825]
+        ])
         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(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), np.asarray(gcr_rad_1), True))
         self.assertFalse(
-            point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=False))
-
+            point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), np.asarray(gcr_rad_1), 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])
+        gcr_rad = np.array([v1_rad, v2_rad])
         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(np.asarray(pt_cart), np.asarray(gcr_cart), np.asarray(gcr_rad), True)
         gcr_car_flipped = np.array([v2_cart, v1_cart])
+        gcr_rad_flipped = np.array([v2_rad, v1_rad])
         self.assertTrue(
-            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_car_flipped), is_directed=True))
+            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_car_flipped), np.asarray(gcr_rad_flipped), 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]])
         pt_cart = np.array(
             [0.10, 0.0, 0.8])
-
-        # Normalize the point abd the GCA
+        gcr_rad = np.array([
+            _xyz_to_lonlat_rad_no_norm(gcr_cart[0][0], gcr_cart[0][1], gcr_cart[0][2]),
+            _xyz_to_lonlat_rad_no_norm(gcr_cart[1][0], gcr_cart[1][1], gcr_cart[1][2])
+        ])
+        # Normalize the point and 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(pt_cart, gcr_cart, is_directed=False))
+        self.assertTrue(point_within_gca(pt_cart, gcr_cart, gcr_rad, is_directed=False))
 
         gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, -0.6]])
+        gcr_rad = np.array([
+            _xyz_to_lonlat_rad_no_norm(gcr_cart[0][0], gcr_cart[0][1], gcr_cart[0][2]),
+            _xyz_to_lonlat_rad_no_norm(gcr_cart[1][0], gcr_cart[1][1], gcr_cart[1][2])
+        ])
         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(pt_cart, gcr_cart, is_directed=False))
+        self.assertFalse(point_within_gca(pt_cart, gcr_cart, gcr_rad, is_directed=False))
         with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart, is_directed=True)
+            point_within_gca(pt_cart, gcr_cart, gcr_rad, is_directed=True)

test/test_helpers.py

@@ -241,7 +241,7 @@ def test_angle_of_2_vectors(self):
 
 class TestFaceEdgeConnectivityHelper(TestCase):
 
-    def test_get_cartesian_face_edge_nodes_pipeline(self):
+    def test_get_cartesian_lonlat_face_edge_nodes_pipeline(self):
         # Create the vertices for the grid, based around the North Pole
         vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]
 
@@ -265,16 +265,20 @@ def test_get_cartesian_face_edge_nodes_pipeline(self):
         face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
             face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
         )
+        face_edges_connectivity_rad = _get_lonlat_rad_face_edge_nodes(
+            face_node_conn, n_face, n_max_face_edges, grid.node_lon.values, grid.node_lat.values
+        )
 
-        # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
+        # Check that the face_edges_connectivity_cartesian and face_edges_connectivity_rad
+        #  works as an input to _pole_point_inside_polygon
         result = ux.grid.geometry._pole_point_inside_polygon(
-            'North', face_edges_connectivity_cartesian[0]
+            'North', face_edges_connectivity_cartesian[0], face_edges_connectivity_rad[0]
         )
 
         # Assert that the result is True
         self.assertTrue(result)
 
-    def test_get_cartesian_face_edge_nodes_filled_value(self):
+    def test_get_cartesian_lonlat_face_edge_nodes_filled_value(self):
         # Create the vertices for the grid, based around the North Pole
         vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]
 
@@ -298,10 +302,14 @@ def test_get_cartesian_face_edge_nodes_filled_value(self):
         face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes(
             face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z
         )
+        face_edges_connectivity_rad = _get_lonlat_rad_face_edge_nodes(
+            face_node_conn, n_face, n_max_face_edges, grid.node_lon.values, grid.node_lat.values
+        )
 
-        # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
+        # Check that the face_edges_connectivity_cartesian and face_edges_connectivity_rad
+        #  works as an input to _pole_point_inside_polygon
         result = ux.grid.geometry._pole_point_inside_polygon(
-            'North', face_edges_connectivity_cartesian[0]
+            'North', face_edges_connectivity_cartesian[0], face_edges_connectivity_rad[0]
         )
 
         # Assert that the result is True
@@ -391,84 +399,6 @@ def test_get_cartesian_face_edge_nodes_filled_value2(self):
         self.assertEqual(face_edges_connectivity_cartesian.shape, correct_result.shape)
 
 
-    def test_get_lonlat_face_edge_nodes_pipeline(self):
-        # Create the vertices for the grid, based around the North Pole
-        vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]
-
-        # Normalize the vertices
-        vertices = [x / np.linalg.norm(x) for x in vertices]
-
-        # Construct the grid from the vertices
-        grid = ux.Grid.from_face_vertices(vertices, latlon=False)
-
-        # Extract the necessary grid data
-        face_node_conn = grid.face_node_connectivity.values
-        n_nodes_per_face = np.array([len(face) for face in face_node_conn])
-        n_face = len(face_node_conn)
-        n_max_face_edges = max(n_nodes_per_face)
-        node_lon = grid.node_lon.values
-        node_lat = grid.node_lat.values
-
-        # Call the function to test
-        face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
-            face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
-        )
-
-        # Convert the first face's edges to Cartesian coordinates
-        face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0]
-        face_edges_connectivity_cartesian = []
-        for edge in face_edges_connectivity_lonlat:
-            edge_cart = [_lonlat_rad_to_xyz(*node) for node in edge]
-            face_edges_connectivity_cartesian.append(edge_cart)
-
-        # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
-            'North', np.array(face_edges_connectivity_cartesian)
-        )
-
-        # Assert that the result is True
-        self.assertTrue(result)
-
-    def test_get_lonlat_face_edge_nodes_filled_value(self):
-        # Create the vertices for the grid, based around the North Pole
-        vertices = [[0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]]
-
-        # Normalize the vertices
-        vertices = [x / np.linalg.norm(x) for x in vertices]
-        vertices.append([INT_FILL_VALUE, INT_FILL_VALUE, INT_FILL_VALUE])
-
-        # Construct the grid from the vertices
-        grid = ux.Grid.from_face_vertices(vertices, latlon=False)
-
-        # Extract the necessary grid data
-        face_node_conn = grid.face_node_connectivity.values
-        n_nodes_per_face = np.array([len(face) for face in face_node_conn])
-        n_face = len(face_node_conn)
-        n_max_face_edges = max(n_nodes_per_face)
-        node_lon = grid.node_lon.values
-        node_lat = grid.node_lat.values
-
-        # Call the function to test
-        face_edges_connectivity_lonlat = _get_lonlat_rad_face_edge_nodes(
-            face_node_conn, n_face, n_max_face_edges, node_lon, node_lat
-        )
-
-        # Convert the first face's edges to Cartesian coordinates
-        face_edges_connectivity_lonlat = face_edges_connectivity_lonlat[0]
-        face_edges_connectivity_cartesian = []
-        for edge in face_edges_connectivity_lonlat:
-            edge_cart = [_lonlat_rad_to_xyz(*node) for node in edge]
-            face_edges_connectivity_cartesian.append(edge_cart)
-
-        # Check that the face_edges_connectivity_cartesian works as an input to _pole_point_inside_polygon
-        result = ux.grid.geometry._pole_point_inside_polygon(
-            'North', np.array(face_edges_connectivity_cartesian)
-        )
-
-        # Assert that the result is True
-        self.assertTrue(result)
-
-
     def test_get_lonlat_face_edge_nodes_filled_value2(self):
         # The face vertices order in counter-clockwise
         # face_conn = [[0,1,2],[1,3,4,2]]