diff --git a/grid_geoflow.exo b/grid_geoflow.exo new file mode 100644 index 000000000..e49b00e4f Binary files /dev/null and b/grid_geoflow.exo differ diff --git a/test/test_helpers.py b/test/test_helpers.py index c5b923a26..079c8c3e6 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -261,6 +261,9 @@ def test_get_cartesian_face_edge_nodes_pipeline(self): node_y = grid.node_y.values node_z = grid.node_z.values + # node_lat = grid.node_lat.values + # node_lon = grid.node_lon.values + # Call the function to test face_edges_connectivity_cartesian = _get_cartesian_face_edge_nodes( face_node_conn, n_face, n_max_face_edges, node_x, node_y, node_z @@ -268,7 +271,7 @@ def test_get_cartesian_face_edge_nodes_pipeline(self): # 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', face_edges_connectivity_cartesian[0] + 'North', face_edges_connectivity_cartesian[0], ) # Assert that the result is True diff --git a/test/test_integrate.py b/test/test_integrate.py index 078f355cf..845758ca1 100644 --- a/test/test_integrate.py +++ b/test/test_integrate.py @@ -69,14 +69,22 @@ def test_get_zonal_face_interval(self): [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]], + face_edge_nodes_cart = np.array([[vertices[0], vertices[1]], [vertices[1], vertices[2]], [vertices[2], vertices[3]], [vertices[3], vertices[0]]]) + face_edge_nodes_rad = np.array([ + [vertices_lonlat[0], vertices_lonlat[1]], + [vertices_lonlat[1], vertices_lonlat[2]], + [vertices_lonlat[2], vertices_lonlat[3]], + [vertices_lonlat[3], vertices_lonlat[0]] + ]) 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, + interval_df = _get_zonal_face_interval(face_edge_nodes_cart, + face_edge_nodes_rad, + constZ, np.array([[-0.25 * np.pi, 0.25 * np.pi], [1.6 * np.pi, 0.4 * np.pi]]), is_directed=False) @@ -94,6 +102,7 @@ def test_get_zonal_face_interval(self): # Asserting almost equal arrays nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) + # TODO: ERROR 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], @@ -103,13 +112,16 @@ def test_get_zonal_face_interval_GCA_constLat(self): vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat] - face_edge_nodes = np.array([[vertices[0], vertices[1]], + face_edge_nodes_cart = np.array([[vertices[0], vertices[1]], [vertices[1], vertices[2]], [vertices[2], vertices[3]], [vertices[3], vertices[0]]]) - + face_edge_nodes_rad = np.array([[vertices_lonlat[0], vertices_lonlat[1]], + [vertices_lonlat[1], vertices_lonlat[2]], + [vertices_lonlat[2], vertices_lonlat[3]], + [vertices_lonlat[3], vertices_lonlat[0]]]) constZ = np.sin(0.20 * np.pi) - interval_df = _get_zonal_face_interval(face_edge_nodes, constZ, + interval_df = _get_zonal_face_interval(face_edge_nodes_cart, face_edge_nodes_rad, 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])) @@ -127,6 +139,7 @@ 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) + # TODO: ERROR def test_get_zonal_face_interval_equator(self): """Test that the zonal face weights are correct.""" vertices_lonlat = [[-0.4 * np.pi, 0.25 * np.pi], [-0.4 * np.pi, 0.0], @@ -134,12 +147,16 @@ def test_get_zonal_face_interval_equator(self): vertices = [_lonlat_rad_to_xyz(*v) for v in vertices_lonlat] - face_edge_nodes = np.array([[vertices[0], vertices[1]], + face_edge_nodes_cart = np.array([[vertices[0], vertices[1]], [vertices[1], vertices[2]], [vertices[2], vertices[3]], [vertices[3], vertices[0]]]) + face_edge_nodes_rad = np.array([[vertices_lonlat[0], vertices_lonlat[1]], + [vertices_lonlat[1], vertices_lonlat[2]], + [vertices_lonlat[2], vertices_lonlat[3]], + [vertices_lonlat[3], vertices_lonlat[0]]]) - interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0, + interval_df = _get_zonal_face_interval(face_edge_nodes_cart, face_edge_nodes_rad, 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])) @@ -158,7 +175,7 @@ def test_get_zonal_face_interval_equator(self): nt.assert_array_almost_equal(actual_values_sorted, expected_values_sorted, decimal=13) # Even if we change the is_GCA_list to False, the result should be the same - interval_df = _get_zonal_face_interval(face_edge_nodes, 0.0, + interval_df = _get_zonal_face_interval(face_edge_nodes_cart, 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])) diff --git a/test/test_intersections.py b/test/test_intersections.py index b134d9029..5a7184c2e 100644 --- a/test/test_intersections.py +++ b/test/test_intersections.py @@ -13,129 +13,183 @@ 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_rad = np.array([ + [np.deg2rad(170.0), np.deg2rad(89.99)], + [np.deg2rad(170.0), np.deg2rad(10.0)] + ]) GCR1_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(89.99)), - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(10.0)) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) ]) + GCR2_rad = np.array([ + [np.deg2rad(70.0), np.deg2rad(0.0)], + [np.deg2rad(179.0), np.deg2rad(0.0)] + ]) GCR2_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0), - _lonlat_rad_to_xyz(np.deg2rad(179.0), 0.0) + _lonlat_rad_to_xyz(GCR2_rad[0][0], GCR2_rad[0][1]), + _lonlat_rad_to_xyz(GCR2_rad[1][0], GCR2_rad[1][1]) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = gca_gca_intersection(GCR1_cart, GCR1_rad, GCR2_cart, GCR2_rad) # res_cart should be empty since these two GCRs are not intersecting self.assertTrue(np.array_equal(res_cart, np.array([]))) + GCR1_rad = np.array([ + [np.deg2rad(170.0), np.deg2rad(89.0)], + [np.deg2rad(170.0), np.deg2rad(-10.0)] + ]) 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)) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) + ]) + GCR2_rad = np.array([ + [np.deg2rad(70.0), 0.0], + [np.deg2rad(175.0), 0.0] ]) GCR2_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(70.0), 0.0), - _lonlat_rad_to_xyz(np.deg2rad(175.0), 0.0) + _lonlat_rad_to_xyz(GCR2_rad[0][0], GCR2_rad[0][1]), + _lonlat_rad_to_xyz(GCR2_rad[1][0], GCR2_rad[1][1]) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + res_cart = gca_gca_intersection(GCR1_cart, GCR1_rad, GCR2_cart, GCR2_rad) + # Test if the result is normalized self.assertTrue( - np.allclose(np.linalg.norm(res_cart, axis=0), - 1.0, - atol=ERROR_TOLERANCE)) + np.allclose( np.linalg.norm(res_cart, axis=0), + 1.0, + atol=ERROR_TOLERANCE + ) + ) res_lonlat_rad = _xyz_to_lonlat_rad(res_cart[0], res_cart[1], res_cart[2]) # res_cart should be [170, 0] self.assertTrue( - np.array_equal(res_lonlat_rad, - np.array([np.deg2rad(170.0), - np.deg2rad(0.0)]))) + np.array_equal( res_lonlat_rad, + np.array([np.deg2rad(170.0), np.deg2rad(0.0)]) + ) + ) def test_get_GCA_GCA_intersections_parallel(self): # Test the case where the two GCAs are parallel + GCR1_rad = np.array([ + [0.3 * np.pi, 0.0], + [0.5 * np.pi, 0.0] + ]) GCR1_cart = np.array([ - _lonlat_rad_to_xyz(0.3 * np.pi, 0.0), - _lonlat_rad_to_xyz(0.5 * np.pi, 0.0) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) + ]) + GCR2_rad = np.array([ + [0.5 * np.pi, 0.0], + [-0.5 * np.pi - 0.01, 0.0] ]) GCR2_cart = np.array([ - _lonlat_rad_to_xyz(0.5 * np.pi, 0.0), - _lonlat_rad_to_xyz(-0.5 * np.pi - 0.01, 0.0) + _lonlat_rad_to_xyz(GCR2_rad[0][0], GCR2_rad[0][1]), + _lonlat_rad_to_xyz(GCR2_rad[1][0], GCR2_rad[1][1]) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + + # Update the function call to pass both Cartesian and lat/lon in radians + res_cart = gca_gca_intersection(GCR1_cart, GCR1_rad, GCR2_cart, GCR2_rad) + self.assertTrue(np.array_equal(res_cart, np.array([]))) def test_get_GCA_GCA_intersections_perpendicular(self): # Test the case where the two GCAs are perpendicular to each other + GCR1_rad = np.array([ + [np.deg2rad(170.0), np.deg2rad(0.0)], + [np.deg2rad(170.0), np.deg2rad(10.0)] + ]) GCR1_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(0.0)), - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(10.0)) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) + ]) + GCR2_rad = np.array([ + [0.5 * np.pi, 0.0], + [-0.5 * np.pi - 0.01, 0.0] ]) GCR2_cart = np.array([ - _lonlat_rad_to_xyz(*[0.5 * np.pi, 0.0]), - _lonlat_rad_to_xyz(*[-0.5 * np.pi - 0.01, 0.0]) + _lonlat_rad_to_xyz(GCR2_rad[0][0], GCR2_rad[0][1]), + _lonlat_rad_to_xyz(GCR2_rad[1][0], GCR2_rad[1][1]) ]) - res_cart = gca_gca_intersection(GCR1_cart, GCR2_cart) + + # Update the function call to pass both Cartesian and lat/lon in radians + res_cart = gca_gca_intersection(GCR1_cart, GCR1_rad, GCR2_cart, GCR2_rad) # Test if the result is normalized self.assertTrue( - np.allclose(np.linalg.norm(res_cart, axis=0), - 1.0, - atol=ERROR_TOLERANCE)) + np.allclose(np.linalg.norm(res_cart, axis=0), + 1.0, + atol=ERROR_TOLERANCE + ) + ) res_lonlat_rad = _xyz_to_lonlat_rad(*res_cart) self.assertTrue( - np.allclose(res_lonlat_rad, - np.array([np.deg2rad(170.0), - np.deg2rad(0.0)]))) + np.allclose(res_lonlat_rad, + np.array([np.deg2rad(170.0), + np.deg2rad(0.0)]) + ) + ) class TestGCAconstLatIntersection(TestCase): def test_GCA_constLat_intersections_antimeridian(self): + GCR1_rad = np.array([ + [np.deg2rad(170.0), np.deg2rad(89.99)], + [np.deg2rad(170.0), np.deg2rad(10.0)] + ]) GCR1_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(89.99)), - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(10.0)) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) ]) - res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(60.0)), verbose=True) + # Update the function call to pass both Cartesian and lat/lon in radians + res = gca_constLat_intersection(GCR1_cart, GCR1_rad, np.sin(np.deg2rad(60.0)), verbose=True) + res_lonlat_rad = _xyz_to_lonlat_rad(*(res[0].tolist())) self.assertTrue( - np.allclose(res_lonlat_rad, - np.array([np.deg2rad(170.0), - np.deg2rad(60.0)]))) + np.allclose( + res_lonlat_rad, + np.array([np.deg2rad(170.0), np.deg2rad(60.0)]) + ) + ) def test_GCA_constLat_intersections_empty(self): + GCR1_rad = np.array([ + [np.deg2rad(170.0), np.deg2rad(89.99)], + [np.deg2rad(170.0), np.deg2rad(10.0)] + ]) GCR1_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(89.99)), - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(10.0)) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) ]) - res = gca_constLat_intersection(GCR1_cart, np.sin(np.deg2rad(-10.0)), verbose=False) + # Update the function call to pass both Cartesian and lat/lon in radians + res = gca_constLat_intersection(GCR1_cart, GCR1_rad, np.sin(np.deg2rad(-10.0)), verbose=False) + self.assertTrue(res.size == 0) + def test_GCA_constLat_intersections_two_pts(self): + GCR1_rad = np.array( + [ + [np.deg2rad(10.0), np.deg2rad(10.0)], + [np.deg2rad(170.0), np.deg2rad(10.0)] + ] + ) GCR1_cart = np.array([ - _lonlat_rad_to_xyz(np.deg2rad(10.0), - np.deg2rad(10)), - _lonlat_rad_to_xyz(np.deg2rad(170.0), - np.deg2rad(10.0)) + _lonlat_rad_to_xyz(GCR1_rad[0][0], GCR1_rad[0][1]), + _lonlat_rad_to_xyz(GCR1_rad[1][0], GCR1_rad[1][1]) ]) - max_lat = ux.grid.arcs.extreme_gca_latitude(GCR1_cart, 'max') + # Calculate the maximum latitude and query latitude + max_lat = ux.grid.arcs.extreme_gca_latitude(GCR1_cart, 'max') query_lat = (np.deg2rad(10.0) + max_lat) / 2.0 - res = gca_constLat_intersection(GCR1_cart, np.sin(query_lat), verbose=False) + # Update the function call to pass both Cartesian and lat/lon in radians + res = gca_constLat_intersection(GCR1_cart, GCR1_rad, np.sin(query_lat), verbose=False) + + # Assert the expected result self.assertTrue(res.shape[0] == 2) +