diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index aadbd45a1..7f7c78a03 100644 --- a/benchmarks/mpas_ocean.py +++ b/benchmarks/mpas_ocean.py @@ -3,6 +3,7 @@ from pathlib import Path import numpy as np +from polars.testing.parametric import dtypes import uxarray as ux @@ -194,6 +195,31 @@ def time_const_lat(self, resolution, lat_step): self.uxgrid.cross_section.constant_latitude(lat) +class PointInPolygon: + param_names = ['resolution'] + params = ['480km', '120km'] + + def setup(self, resolution): + self.uxgrid = ux.open_grid(file_path_dict[resolution][0]) + self.uxgrid.normalize_cartesian_coordinates() + + # Construct variables needed to ensure that the benchmark doesn't measure construction time + _ = self.uxgrid.face_edge_connectivity + _ = self.uxgrid.face_x.values + _ = self.uxgrid.face_lon.values + + point = np.array([0.0, 0.0, 1.0]) + res = self.uxgrid.get_faces_containing_point(point) + + def teardown(self, resolution): + del self.uxgrid + + def time_face_search(self, resolution): + point_xyz = np.array([self.uxgrid.face_x[0].values, self.uxgrid.face_y[0].values, self.uxgrid.face_z[0].values], dtype=np.float64) + point_lonlat = np.array([self.uxgrid.face_lon[0].values, self.uxgrid.face_lat.values[0]], dtype=np.float64) + self.uxgrid.get_faces_containing_point(point_xyz=point_xyz, point_lonlat=point_lonlat) + + class ZonalAverage(DatasetBenchmark): def setup(self, resolution, *args, **kwargs): self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1]) diff --git a/docs/api.rst b/docs/api.rst index e5ca434aa..2c9170447 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -139,6 +139,7 @@ Descriptors Grid.partial_sphere_coverage Grid.global_sphere_coverage Grid.triangular + Grid.max_face_radius Attributes ~~~~~~~~~~ @@ -159,6 +160,7 @@ Methods Grid.calculate_total_face_area Grid.normalize_cartesian_coordinates Grid.construct_face_centers + Grid.get_faces_containing_point Inheritance of Xarray Functionality ----------------------------------- diff --git a/test/test_geometry.py b/test/test_geometry.py index 7de043479..c066f612a 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -9,12 +9,15 @@ import uxarray as ux from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE import uxarray.utils.computing as ac_utils -from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad +from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad, \ + _xyz_to_lonlat_deg, _xyz_to_lonlat_rad_scalar from uxarray.grid.arcs import extreme_gca_latitude, extreme_gca_z from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes -from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, stereographic_projection, inverse_stereographic_projection +from uxarray.grid.geometry import _populate_face_latlon_bound, _populate_bounds, _pole_point_inside_polygon_cartesian, \ + stereographic_projection, inverse_stereographic_projection, point_in_face, haversine_distance +from sklearn.metrics.pairwise import haversine_distances current_path = Path(os.path.dirname(os.path.realpath(__file__))) @@ -31,9 +34,12 @@ grid_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc" grid_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc" +grid_mpas_2 = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" + # List of grid files to test grid_files_latlonBound = [grid_quad_hex, grid_geoflow, gridfile_CSne8, grid_mpas] + def test_antimeridian_crossing(): verts = [[[-170, 40], [180, 30], [165, 25], [-170, 20]]] @@ -44,6 +50,7 @@ def test_antimeridian_crossing(): assert len(uxgrid.antimeridian_face_indices) == 1 assert len(gdf['geometry']) == 1 + def test_antimeridian_point_on(): verts = [[[-170, 40], [180, 30], [-170, 20]]] @@ -51,10 +58,12 @@ def test_antimeridian_point_on(): assert len(uxgrid.antimeridian_face_indices) == 1 + def test_linecollection_execution(): uxgrid = ux.open_grid(gridfile_CSne8) lines = uxgrid.to_linecollection() + def test_pole_point_inside_polygon_from_vertice_north(): 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]] @@ -74,6 +83,7 @@ def test_pole_point_inside_polygon_from_vertice_north(): result = _pole_point_inside_polygon_cartesian('South', face_edge_cart) assert not result, "South pole should not be inside the polygon" + def test_pole_point_inside_polygon_from_vertice_south(): vertices = [[0.5, 0.5, -0.5], [-0.5, 0.5, -0.5], [0.0, 0.0, -1.0]] @@ -91,6 +101,7 @@ def test_pole_point_inside_polygon_from_vertice_south(): result = _pole_point_inside_polygon_cartesian('South', face_edge_cart) assert result, "South pole should be inside the polygon" + def test_pole_point_inside_polygon_from_vertice_pole(): vertices = [[0, 0, 1], [-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5]] @@ -110,6 +121,7 @@ def test_pole_point_inside_polygon_from_vertice_pole(): result = _pole_point_inside_polygon_cartesian('South', face_edge_cart) assert not result, "South pole should not be inside the polygon" + def test_pole_point_inside_polygon_from_vertice_cross(): vertices = [[0.6, -0.3, 0.5], [0.2, 0.2, -0.2], [-0.5, 0.1, -0.2], [-0.1, -0.2, 0.2]] @@ -198,7 +210,7 @@ def _max_latitude_rad_iterative(gca_cart): # Update maximum latitude and section if needed max_lat = max(max_lat, w1_lonlat[1], w2_lonlat[1]) if np.abs(w2_lonlat[1] - - w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[ + w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[ 1] == max_lat == w2_lonlat[1]: max_section = [w1_new, w2_new] break @@ -215,6 +227,7 @@ def _max_latitude_rad_iterative(gca_cart): return np.average([b_lonlat[1], c_lonlat[1]]) + def _min_latitude_rad_iterative(gca_cart): """Calculate the minimum latitude of a great circle arc defined by two points. @@ -285,7 +298,7 @@ def _min_latitude_rad_iterative(gca_cart): # Update minimum latitude and section if needed min_lat = min(min_lat, w1_lonlat[1], w2_lonlat[1]) if np.abs(w2_lonlat[1] - - w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[ + w1_lonlat[1]) <= ERROR_TOLERANCE or w1_lonlat[ 1] == min_lat == w2_lonlat[1]: min_section = [w1_new, w2_new] break @@ -302,6 +315,7 @@ def _min_latitude_rad_iterative(gca_cart): return np.average([b_lonlat[1], c_lonlat[1]]) + def test_extreme_gca_latitude_max(): gca_cart = np.array([ _normalize_xyz(*[0.5, 0.5, 0.5]), @@ -317,6 +331,7 @@ def test_extreme_gca_latitude_max(): expected_max_latitude = 1.0 assert np.isclose(max_latitude, expected_max_latitude, atol=ERROR_TOLERANCE) + def test_extreme_gca_latitude_max_short(): # 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]]) @@ -346,6 +361,7 @@ def test_extreme_gca_latitude_min(): expected_min_latitude = -np.pi / 2 assert np.isclose(min_latitude, expected_min_latitude, atol=ERROR_TOLERANCE) + def test_get_latlonbox_width(): gca_latlon = np.array([[0.0, 0.0], [0.0, 3.0]]) width = ux.grid.geometry._get_latlonbox_width(gca_latlon) @@ -355,6 +371,7 @@ def test_get_latlonbox_width(): width = ux.grid.geometry._get_latlonbox_width(gca_latlon) assert width == 2.0 + def test_insert_pt_in_latlonbox_non_periodic(): old_box = np.array([[0.1, 0.2], [0.3, 0.4]]) # Radians new_pt = np.array([0.15, 0.35]) @@ -362,6 +379,7 @@ def test_insert_pt_in_latlonbox_non_periodic(): result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt, False) np.testing.assert_array_equal(result, expected) + def test_insert_pt_in_latlonbox_periodic(): old_box = np.array([[0.1, 0.2], [6.0, 0.1]]) # Radians, periodic new_pt = np.array([0.15, 6.2]) @@ -369,6 +387,7 @@ def test_insert_pt_in_latlonbox_periodic(): result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt, True) np.testing.assert_array_equal(result, expected) + def test_insert_pt_in_latlonbox_pole(): old_box = np.array([[0.1, 0.2], [0.3, 0.4]]) new_pt = np.array([np.pi / 2, np.nan]) # Pole point @@ -376,6 +395,7 @@ def test_insert_pt_in_latlonbox_pole(): result = ux.grid.geometry.insert_pt_in_latlonbox(old_box, new_pt) np.testing.assert_array_equal(result, expected) + def test_insert_pt_in_empty_state(): old_box = np.array([[np.nan, np.nan], [np.nan, np.nan]]) # Empty state @@ -385,7 +405,8 @@ def test_insert_pt_in_empty_state(): np.testing.assert_array_equal(result, expected) -def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z): +def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid, + node_x, node_y, node_z): """Construct an array to hold the edge Cartesian coordinates connectivity for a face in a grid.""" mask = face_edges_ind != INT_FILL_VALUE valid_edges = face_edges_ind[mask] @@ -406,7 +427,9 @@ def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ return cartesian_coordinates -def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat): + +def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes_ind, face_edges_ind, edge_nodes_grid, + node_lon, node_lat): """Construct an array to hold the edge lat lon in radian connectivity for a face in a grid.""" mask = face_edges_ind != INT_FILL_VALUE valid_edges = face_edges_ind[mask] @@ -433,6 +456,7 @@ def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca(face_nodes return lonlat_coordinates + def test_populate_bounds_normal_latlon_bounds_gca(): vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -460,6 +484,7 @@ def test_populate_bounds_normal_latlon_bounds_gca(): bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_antimeridian_latlon_bounds_gca(): vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -487,6 +512,7 @@ def test_populate_bounds_antimeridian_latlon_bounds_gca(): bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_equator_latlon_bounds_gca(): face_edges_cart = np.array([ [[0.99726469, -0.05226443, -0.05226443], [0.99862953, 0.0, -0.05233596]], @@ -501,6 +527,7 @@ def test_populate_bounds_equator_latlon_bounds_gca(): expected_bounds = np.array([[-0.05235988, 0], [6.23082543, 0]]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_south_sphere_latlon_bounds_gca(): face_edges_cart = np.array([ [[-1.04386773e-01, -5.20500333e-02, -9.93173799e-01], [-1.04528463e-01, -1.28010448e-17, -9.94521895e-01]], @@ -516,6 +543,7 @@ def test_populate_bounds_south_sphere_latlon_bounds_gca(): expected_bounds = np.array([[-1.51843645, -1.45388627], [3.14159265, 3.92699082]]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_near_pole_latlon_bounds_gca(): face_edges_cart = np.array([ [[3.58367950e-01, 0.00000000e+00, -9.33580426e-01], [3.57939780e-01, 4.88684203e-02, -9.32465008e-01]], @@ -531,6 +559,7 @@ def test_populate_bounds_near_pole_latlon_bounds_gca(): expected_bounds = np.array([[-1.20427718, -1.14935491], [0, 0.13568803]]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_near_pole2_latlon_bounds_gca(): face_edges_cart = np.array([ [[3.57939780e-01, -4.88684203e-02, -9.32465008e-01], [3.58367950e-01, 0.00000000e+00, -9.33580426e-01]], @@ -546,6 +575,7 @@ def test_populate_bounds_near_pole2_latlon_bounds_gca(): expected_bounds = np.array([[-1.20427718, -1.14935491], [6.147497, 4.960524e-16]]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_long_face_latlon_bounds_gca(): face_edges_cart = np.array([ [[9.9999946355819702e-01, -6.7040475551038980e-04, 8.0396590055897832e-04], @@ -570,6 +600,7 @@ def test_populate_bounds_long_face_latlon_bounds_gca(): # The expected bounds should not contain the south pole [0,-0.5*np.pi] assert bounds[1][0] != 0.0 + def test_populate_bounds_node_on_pole_latlon_bounds_gca(): vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -596,6 +627,7 @@ def test_populate_bounds_node_on_pole_latlon_bounds_gca(): bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_edge_over_pole_latlon_bounds_gca(): vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -622,6 +654,7 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_gca(): bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_pole_inside_latlon_bounds_gca(): vertices_lonlat = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -649,7 +682,8 @@ def test_populate_bounds_pole_inside_latlon_bounds_gca(): np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) -def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z): +def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind, + edge_nodes_grid, node_x, node_y, node_z): """Construct an array to hold the edge Cartesian coordinates connectivity for a face in a grid.""" mask = face_edges_ind != INT_FILL_VALUE valid_edges = face_edges_ind[mask] @@ -670,7 +704,9 @@ def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face return cartesian_coordinates -def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat): + +def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(face_nodes_ind, face_edges_ind, + edge_nodes_grid, node_lon, node_lat): """Construct an array to hold the edge lat lon in radian connectivity for a face in a grid.""" mask = face_edges_ind != INT_FILL_VALUE valid_edges = face_edges_ind[mask] @@ -697,6 +733,7 @@ def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_latlonface(fac return lonlat_coordinates + def test_populate_bounds_normal_latlon_bounds_latlonface(): vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -718,9 +755,11 @@ def test_populate_bounds_normal_latlon_bounds_latlonface(): grid.edge_node_connectivity.values, grid.node_lon.values, grid.node_lat.values) expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]]) - bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True) + bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, + is_latlonface=True) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_antimeridian_latlon_bounds_latlonface(): vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -742,9 +781,11 @@ def test_populate_bounds_antimeridian_latlon_bounds_latlonface(): grid.edge_node_connectivity.values, grid.node_lon.values, grid.node_lat.values) expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]]) - bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True) + bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, + is_latlonface=True) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_node_on_pole_latlon_bounds_latlonface(): vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -766,9 +807,11 @@ def test_populate_bounds_node_on_pole_latlon_bounds_latlonface(): grid.edge_node_connectivity.values, grid.node_lon.values, grid.node_lat.values) expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]]) - bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True) + bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, + is_latlonface=True) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_edge_over_pole_latlon_bounds_latlonface(): vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -790,9 +833,11 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_latlonface(): grid.edge_node_connectivity.values, grid.node_lon.values, grid.node_lat.values) expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]]) - bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True) + bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, + is_latlonface=True) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_pole_inside_latlon_bounds_latlonface(): vertices_lonlat = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -814,11 +859,13 @@ def test_populate_bounds_pole_inside_latlon_bounds_latlonface(): grid.edge_node_connectivity.values, grid.node_lon.values, grid.node_lat.values) expected_bounds = np.array([[lat_min, lat_max], [lon_min, lon_max]]) - bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, is_latlonface=True) + bounds = _populate_face_latlon_bound(face_edges_connectivity_cartesian, face_edges_connectivity_lonlat, + is_latlonface=True) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) -def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_x, node_y, node_z): +def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind, + edge_nodes_grid, node_x, node_y, node_z): """Construct an array to hold the edge Cartesian coordinates connectivity for a face in a grid.""" mask = face_edges_ind != INT_FILL_VALUE valid_edges = face_edges_ind[mask] @@ -839,7 +886,9 @@ def _get_cartesian_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_n return cartesian_coordinates -def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind, edge_nodes_grid, node_lon, node_lat): + +def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_nodes_ind, face_edges_ind, + edge_nodes_grid, node_lon, node_lat): """Construct an array to hold the edge lat lon in radian connectivity for a face in a grid.""" mask = face_edges_ind != INT_FILL_VALUE valid_edges = face_edges_ind[mask] @@ -866,6 +915,7 @@ def _get_lonlat_rad_face_edge_nodes_testcase_helper_latlon_bounds_gca_list(face_ return lonlat_coordinates + def test_populate_bounds_normal_latlon_bounds_gca_list(): vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -891,6 +941,7 @@ def test_populate_bounds_normal_latlon_bounds_gca_list(): is_GCA_list=[True, False, True, False]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_antimeridian_latlon_bounds_gca_list(): vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -916,6 +967,7 @@ def test_populate_bounds_antimeridian_latlon_bounds_gca_list(): is_GCA_list=[True, False, True, False]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_node_on_pole_latlon_bounds_gca_list(): vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -941,6 +993,7 @@ def test_populate_bounds_node_on_pole_latlon_bounds_gca_list(): is_GCA_list=[True, False, True, False]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_edge_over_pole_latlon_bounds_gca_list(): vertices_lonlat = [[210.0, 80.0], [350.0, 60.0], [10.0, 60.0], [30.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -966,6 +1019,7 @@ def test_populate_bounds_edge_over_pole_latlon_bounds_gca_list(): is_GCA_list=[True, False, True, False]) np.testing.assert_allclose(bounds, expected_bounds, atol=ERROR_TOLERANCE) + def test_populate_bounds_pole_inside_latlon_bounds_gca_list(): vertices_lonlat = [[200.0, 80.0], [350.0, 60.0], [10.0, 60.0], [40.0, 80.0]] vertices_lonlat = np.array(vertices_lonlat) @@ -1010,6 +1064,7 @@ def test_populate_bounds_GCA_mix_latlon_bounds_mix(): for i in range(len(faces)): np.testing.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) + def test_populate_bounds_LatlonFace_mix_latlon_bounds_mix(): face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] face_2 = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] @@ -1029,6 +1084,7 @@ def test_populate_bounds_LatlonFace_mix_latlon_bounds_mix(): for i in range(len(faces)): np.testing.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) + def test_populate_bounds_GCAList_mix_latlon_bounds_mix(): face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] face_2 = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] @@ -1051,6 +1107,7 @@ def test_populate_bounds_GCAList_mix_latlon_bounds_mix(): for i in range(len(faces)): np.testing.assert_allclose(face_bounds[i], expected_bounds[i], atol=ERROR_TOLERANCE) + def test_face_bounds_latlon_bounds_files(): """Test to ensure ``Grid.face_bounds`` works correctly for all grid files.""" for grid_path in grid_files_latlonBound: @@ -1071,16 +1128,19 @@ def test_face_bounds_latlon_bounds_files(): # Clean up the grid object del uxgrid + def test_engine_geodataframe(): uxgrid = ux.open_grid(gridfile_geoflow) for engine in ['geopandas', 'spatialpandas']: gdf = uxgrid.to_geodataframe(engine=engine) + def test_periodic_elements_geodataframe(): uxgrid = ux.open_grid(gridfile_geoflow) for periodic_elements in ['ignore', 'exclude', 'split']: gdf = uxgrid.to_geodataframe(periodic_elements=periodic_elements) + def test_to_gdf_geodataframe(): uxgrid = ux.open_grid(gridfile_geoflow) @@ -1088,6 +1148,7 @@ def test_to_gdf_geodataframe(): gdf_without_am = uxgrid.to_geodataframe(exclude_antimeridian=True) + def test_cache_and_override_geodataframe(): """Tests the cache and override functionality for GeoDataFrame conversion.""" uxgrid = ux.open_grid(gridfile_geoflow) @@ -1114,6 +1175,200 @@ def test_cache_and_override_geodataframe(): assert gdf_f is not gdf_e + +# Test point_in_face function +def test_point_inside(): + """Test the function `point_in_face`, where the points are all inside the face""" + + # Open grid + grid = ux.open_grid(grid_mpas_2) + + # Get the face edges of all faces in the grid + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + # Loop through each face + for i in range(grid.n_face): + # Set the point as the face center of the polygon + point_xyz = np.array([grid.face_x[i].values, grid.face_y[i].values, grid.face_z[i].values]) + + # Assert that the point is in the polygon + assert point_in_face(faces_edges_cartesian[i], point_xyz, inclusive=True) + + +def test_point_outside(): + """Test the function `point_in_face`, where the point is outside the face""" + + # Open grid + grid = ux.open_grid(grid_mpas_2) + + # Get the face edges of all faces in the grid + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + # Set the point as the face center of a different face than the face tested + point_xyz = np.array([grid.face_x[1].values, grid.face_y[1].values, grid.face_z[1].values]) + + # Assert that the point is not in the face tested + assert not point_in_face(faces_edges_cartesian[0], point_xyz, inclusive=True) + + +def test_point_on_node(): + """Test the function `point_in_face`, when the point is on the node of the polygon""" + + # Open grid + grid = ux.open_grid(grid_mpas_2) + + # Get the face edges of all faces in the grid + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + # Set the point as a node + point_xyz = np.array([*faces_edges_cartesian[0][0][0]]) + + # Assert that the point is in the face when inclusive is true + assert point_in_face(faces_edges_cartesian[0], point_xyz, inclusive=True) + + # Assert that the point is not in the face when inclusive is false + assert not point_in_face(faces_edges_cartesian[0], point_xyz, inclusive=False) + + +def test_point_inside_close(): + """Test the function `point_in_face`, where the point is inside the face, but very close to the edge""" + + # Create a square + vertices_lonlat = [[-10.0, 10.0], [-10.0, -10.0], [10.0, -10.0], [10.0, 10.0]] + vertices_lonlat = np.array(vertices_lonlat) + + # Choose a point just inside the square + point = np.array(_lonlat_rad_to_xyz(np.deg2rad(0.0), np.deg2rad(-9.8))) + + # Create the grid and face edges + grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + # Use point in face to determine if the point is inside or out of the face + assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False) + + +def test_point_outside_close(): + """Test the function `point_in_face`, where the point is outside the face, but very close to the edge""" + + # Create a square + vertices_lonlat = [[-10.0, 10.0], [-10.0, -10.0], [10.0, -10.0], [10.0, 10.0]] + vertices_lonlat = np.array(vertices_lonlat) + + # Choose a point just outside the square + point = np.array(_lonlat_rad_to_xyz(np.deg2rad(0.0), np.deg2rad(-10.2))) + + # Create the grid and face edges + grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + # Use point in face to determine if the point is inside or out of the face + assert not point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=False) + + +def test_face_at_pole(): + """Test the function `point_in_face`, when the face is at the North Pole""" + + # Generate a face that is at a pole + vertices_lonlat = [[10.0, 90.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] + vertices_lonlat = np.array(vertices_lonlat) + + point = np.array(_lonlat_rad_to_xyz(np.deg2rad(25), np.deg2rad(30))) + + # Create the grid and face edges + grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True) + + +def test_face_at_antimeridian(): + """Test the function `point_in_face`, where the face crosses the antimeridian""" + + # Generate a face crossing the antimeridian + vertices_lonlat = [[350, 60.0], [350, 10.0], [50.0, 10.0], [50.0, 60.0]] + vertices_lonlat = np.array(vertices_lonlat) + point = np.array(_lonlat_rad_to_xyz(np.deg2rad(25), np.deg2rad(30))) + + # Create the grid and face edges + grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True) + + +def test_face_normal_face(): + """Test the function `point_in_face`, where the face is a normal face, not crossing the antimeridian or the + poles""" + + # Generate a normal face that is not crossing the antimeridian or the poles + vertices_lonlat = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] + vertices_lonlat = np.array(vertices_lonlat) + point = np.array(_lonlat_rad_to_xyz(np.deg2rad(25), np.deg2rad(30))) + + # Create the grid and face edges + grid = ux.Grid.from_face_vertices(vertices_lonlat, latlon=True) + faces_edges_cartesian = _get_cartesian_face_edge_nodes( + grid.face_node_connectivity.values, + grid.n_face, + grid.n_max_face_edges, + grid.node_x.values, + grid.node_y.values, + grid.node_z.values, + ) + + assert point_in_face(faces_edges_cartesian[0], point_xyz=point, inclusive=True) + + def test_stereographic_projection_stereographic_projection(): lon = np.array(0) lat = np.array(0) @@ -1128,3 +1383,17 @@ def test_stereographic_projection_stereographic_projection(): assert np.array_equal(lon, new_lon) assert np.array_equal(lat, new_lat) assert np.array_equal(x, y) and x == 0 + + +def test_haversine_distance_creation(): + """Tests the use of `haversine_distance`""" + + # Create two points + point_a = [np.deg2rad(-34.8), np.deg2rad(-58.5)] + point_b = [np.deg2rad(49.0), np.deg2rad(2.6)] + + result = haversine_distances([point_a, point_b]) + + distance = haversine_distance(point_a[1], point_a[0], point_b[1], point_b[0]) + + assert np.isclose(result[0][1], distance, atol=1e-6) diff --git a/test/test_grid.py b/test/test_grid.py index 4c8d98c76..df9b9dddd 100644 --- a/test/test_grid.py +++ b/test/test_grid.py @@ -13,13 +13,14 @@ from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, \ _build_edge_node_connectivity, _build_face_face_connectivity, _populate_face_face_connectivity -from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz +from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz, _xyz_to_lonlat_rad_scalar from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE from uxarray.grid.arcs import extreme_gca_latitude from uxarray.grid.validation import _find_duplicate_nodes +from .test_gradient import quad_hex_grid_path try: import constants @@ -34,8 +35,6 @@ gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" gridfile_fesom = current_path / "meshfiles" / "ugrid" / "fesom" / "fesom.mesh.diag.nc" gridfile_geoflow = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc" -gridfile_mpas = current_path / 'meshfiles' / "mpas" / "QU" / 'mesh.QU.1920km.151026.nc' -gridfile_mpas_two = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc' gridfile_geos = current_path / 'meshfiles' / "geos-cs" / "c12" / 'test-c12.native.nc4' gridfile_mpas_holes = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc' @@ -44,13 +43,10 @@ shp_filename = current_path / "meshfiles" / "shp" / "grid_fire.shp" - - grid_CSne30 = ux.open_grid(gridfile_CSne30) grid_RLL1deg = ux.open_grid(gridfile_RLL1deg) grid_RLL10deg_CSne4 = ux.open_grid(gridfile_RLL10deg_CSne4) - mpas_filepath = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" exodus_filepath = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g" ugrid_filepath_01 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" @@ -80,6 +76,7 @@ gridfile_ugrid = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc" gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" +gridfile_mpas_two = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc' gridfile_exodus = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g" gridfile_scrip = current_path / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc" @@ -89,6 +86,7 @@ def test_grid_validate(): grid_mpas = ux.open_grid(gridfile_mpas) assert grid_mpas.validate() + def test_grid_with_holes(): """Test _holes_in_mesh function.""" grid_without_holes = ux.open_grid(gridfile_mpas) @@ -97,6 +95,7 @@ def test_grid_with_holes(): assert grid_with_holes.partial_sphere_coverage assert grid_without_holes.global_sphere_coverage + def test_grid_encode_as(): """Reads a ugrid file and encodes it as `xarray.Dataset` in various types.""" grid_CSne30.encode_as("UGRID") @@ -107,6 +106,7 @@ def test_grid_encode_as(): grid_RLL1deg.encode_as("Exodus") grid_RLL10deg_CSne4.encode_as("Exodus") + def test_grid_init_verts(): """Create a uxarray grid from multiple face vertices with duplicate nodes and saves a ugrid file.""" cart_x = [ @@ -131,7 +131,7 @@ def test_grid_init_verts(): [5, 4, 7, 6], # back face [4, 0, 3, 7], # left face [3, 2, 6, 7], # top face - [4, 5, 1, 0] # bottom face + [4, 5, 1, 0] # bottom face ] faces_coords = [] @@ -163,6 +163,7 @@ def test_grid_init_verts(): assert vgrid.n_node == 6 vgrid.encode_as("UGRID") + def test_grid_init_verts_different_input_datatype(): """Create a uxarray grid from multiple face vertices with different datatypes (ndarray, list, tuple) and saves a ugrid file.""" faces_verts_ndarray = np.array([ @@ -176,8 +177,8 @@ def test_grid_init_verts_different_input_datatype(): vgrid.encode_as("UGRID") faces_verts_list = [[[150, 10], [160, 20], [150, 30], [135, 30], [125, 20], [135, 10]], - [[125, 20], [135, 30], [125, 60], [110, 60], [100, 30], [105, 20]], - [[95, 10], [105, 20], [100, 30], [85, 30], [75, 20], [85, 10]]] + [[125, 20], [135, 30], [125, 60], [110, 60], [100, 30], [105, 20]], + [[95, 10], [105, 20], [100, 30], [85, 30], [75, 20], [85, 10]]] vgrid = ux.open_grid(faces_verts_list, latlon=True) assert vgrid.n_face == 3 assert vgrid.n_node == 14 @@ -195,6 +196,7 @@ def test_grid_init_verts_different_input_datatype(): assert vgrid.validate() vgrid.encode_as("UGRID") + def test_grid_init_verts_fill_values(): faces_verts_filled_values = [[[150, 10], [160, 20], [150, 30], [135, 30], [125, 20], [135, 10]], @@ -208,6 +210,7 @@ def test_grid_init_verts_fill_values(): assert vgrid.n_face == 3 assert vgrid.n_node == 12 + def test_grid_properties(): """Tests to see if accessing variables through set properties is equal to using the dict.""" xr.testing.assert_equal(grid_CSne30.node_lon, grid_CSne30._ds["node_lon"]) @@ -234,27 +237,32 @@ def test_grid_properties(): assert n_faces == grid_geoflow.n_face assert n_face_nodes == grid_geoflow.n_max_face_nodes + def test_read_shpfile(): """Reads a shape file and write ugrid file.""" with pytest.raises(ValueError): grid_shp = ux.open_grid(shp_filename) + def test_read_scrip(): """Reads a scrip file.""" grid_CSne8 = ux.open_grid(gridfile_CSne8) # tests from scrip + def test_operators_eq(): """Test Equals ('==') operator.""" grid_CSne30_01 = ux.open_grid(gridfile_CSne30) grid_CSne30_02 = ux.open_grid(gridfile_CSne30) assert grid_CSne30_01 == grid_CSne30_02 + def test_operators_ne(): """Test Not Equals ('!=') operator.""" grid_CSne30_01 = ux.open_grid(gridfile_CSne30) grid_RLL1deg = ux.open_grid(gridfile_RLL1deg) assert grid_CSne30_01 != grid_RLL1deg + def test_face_areas_calculate_total_face_area_triangle(): """Create a uxarray grid from vertices and saves an exodus file.""" verts = [[[0.57735027, -5.77350269e-01, -0.57735027], @@ -275,11 +283,13 @@ def test_face_areas_calculate_total_face_area_triangle(): quadrature_rule="triangular", order=4) nt.assert_almost_equal(area_triangular, constants.TRI_AREA, decimal=1) + def test_face_areas_calculate_total_face_area_file(): """Create a uxarray grid from vertices and saves an exodus file.""" area = ux.open_grid(gridfile_CSne30).calculate_total_face_area() nt.assert_almost_equal(area, constants.MESH30_AREA, decimal=3) + def test_face_areas_calculate_total_face_area_sphere(): """Computes the total face area of an MPAS mesh that lies on a unit sphere, with an expected total face area of 4pi.""" mpas_grid_path = current_path / 'meshfiles' / "mpas" / "QU" / 'mesh.QU.1920km.151026.nc' @@ -293,11 +303,13 @@ def test_face_areas_calculate_total_face_area_sphere(): nt.assert_almost_equal(primal_face_area, constants.UNIT_SPHERE_AREA, decimal=3) nt.assert_almost_equal(dual_face_area, constants.UNIT_SPHERE_AREA, decimal=3) + def test_face_areas_compute_face_areas_geoflow_small(): """Checks if the GeoFlow Small can generate a face areas output.""" grid_geoflow = ux.open_grid(gridfile_geoflow) grid_geoflow.compute_face_areas() + def test_face_areas_verts_calc_area(): faces_verts_ndarray = np.array([ np.array([[150, 10, 0], [160, 20, 0], [150, 30, 0], [135, 30, 0], @@ -311,11 +323,12 @@ def test_face_areas_verts_calc_area(): face_verts_areas = verts_grid.face_areas nt.assert_almost_equal(face_verts_areas.sum(), constants.FACE_VERTS_AREA, decimal=3) + def test_populate_coordinates_populate_cartesian_xyz_coord(): # The following testcases are generated through the matlab cart2sph/sph2cart functions lon_deg = [ 45.0001052295749, 45.0001052295749, 360 - 45.0001052295749, - 360 - 45.0001052295749 + 360 - 45.0001052295749 ] lat_deg = [ 35.2655522903022, -35.2655522903022, 35.2655522903022, @@ -342,10 +355,11 @@ def test_populate_coordinates_populate_cartesian_xyz_coord(): nt.assert_almost_equal(vgrid.node_y.values[i], cart_y[i], decimal=12) nt.assert_almost_equal(vgrid.node_z.values[i], cart_z[i], decimal=12) + def test_populate_coordinates_populate_lonlat_coord(): lon_deg = [ 45.0001052295749, 45.0001052295749, 360 - 45.0001052295749, - 360 - 45.0001052295749 + 360 - 45.0001052295749 ] lat_deg = [ 35.2655522903022, -35.2655522903022, 35.2655522903022, @@ -374,8 +388,8 @@ def test_populate_coordinates_populate_lonlat_coord(): def _revert_edges_conn_to_face_nodes_conn(edge_nodes_connectivity: np.ndarray, - face_edges_connectivity: np.ndarray, - original_face_nodes_connectivity: np.ndarray): + face_edges_connectivity: np.ndarray, + original_face_nodes_connectivity: np.ndarray): """Utilize the edge_nodes_connectivity and face_edges_connectivity to generate the res_face_nodes_connectivity in the counter-clockwise order. The counter-clockwise order will be enforced by the passed in @@ -440,6 +454,7 @@ def _revert_edges_conn_to_face_nodes_conn(edge_nodes_connectivity: np.ndarray, return np.array(res_face_nodes_connectivity) + def test_connectivity_build_n_nodes_per_face(): """Tests the construction of the ``n_nodes_per_face`` variable.""" grids = [grid_mpas, grid_exodus, grid_ugrid] @@ -457,6 +472,7 @@ def test_connectivity_build_n_nodes_per_face(): expected_nodes_per_face = np.array([6, 3, 4, 6, 6, 4, 4], dtype=int) nt.assert_equal(grid_from_verts.n_nodes_per_face.values, expected_nodes_per_face) + def test_connectivity_edge_nodes_euler(): """Verifies that (``n_edge``) follows euler's formula.""" grid_paths = [exodus_filepath, ugrid_filepath_01, ugrid_filepath_02, ugrid_filepath_03] @@ -470,6 +486,7 @@ def test_connectivity_edge_nodes_euler(): assert (n_face == n_edge - n_node + 2) + def test_connectivity_build_face_edges_connectivity_mpas(): """Tests the construction of (``Mesh2_edge_nodes``) on an MPAS grid with known edge nodes.""" from uxarray.grid.connectivity import _build_edge_node_connectivity @@ -492,6 +509,7 @@ def test_connectivity_build_face_edges_connectivity_mpas(): assert (n_face == n_edge - n_node + 2) + def test_connectivity_build_face_edges_connectivity(): """Generates Grid.Mesh2_edge_nodes from Grid.face_node_connectivity.""" ug_filename_list = [ugrid_filepath_01, ugrid_filepath_02, ugrid_filepath_03] @@ -522,6 +540,7 @@ def test_connectivity_build_face_edges_connectivity(): for i in range(len(reverted_mesh2_edge_nodes)): assert np.array_equal(reverted_mesh2_edge_nodes[i], original_face_nodes_connectivity[i]) + def test_connectivity_build_face_edges_connectivity_fillvalues(): verts = [f0_deg, f1_deg, f2_deg, f3_deg, f4_deg, f5_deg, f6_deg] uds = ux.open_grid(verts) @@ -543,6 +562,7 @@ def test_connectivity_build_face_edges_connectivity_fillvalues(): assert np.array_equal(res_face_nodes_connectivity, uds._ds["face_node_connectivity"].values) + def test_connectivity_node_face_connectivity_from_verts(): """Test generating Grid.Mesh2_node_faces from array input.""" face_nodes_conn_lonlat_degree = [[162., 30], [216., 30], [70., 30], @@ -573,6 +593,7 @@ def test_connectivity_node_face_connectivity_from_verts(): assert np.array_equal(vgrid.node_face_connectivity.values, expected) + def test_connectivity_node_face_connectivity_from_files(): """Test generating Grid.Mesh2_node_faces from file input.""" grid_paths = [exodus_filepath, ugrid_filepath_01, ugrid_filepath_02, ugrid_filepath_03] @@ -593,12 +614,14 @@ def test_connectivity_node_face_connectivity_from_files(): for i in range(grid_ux.n_node): face_index_from_sparse_matrix = grid_ux.node_face_connectivity.values[i] - valid_face_index_from_sparse_matrix = face_index_from_sparse_matrix[face_index_from_sparse_matrix != grid_ux.node_face_connectivity.attrs["_FillValue"]] + valid_face_index_from_sparse_matrix = face_index_from_sparse_matrix[ + face_index_from_sparse_matrix != grid_ux.node_face_connectivity.attrs["_FillValue"]] valid_face_index_from_sparse_matrix.sort() face_index_from_dict = node_face_connectivity[i] face_index_from_dict.sort() assert np.array_equal(valid_face_index_from_sparse_matrix, face_index_from_dict) + def test_connectivity_edge_face_connectivity_mpas(): """Tests the construction of ``Mesh2_face_edges`` to the expected results of an MPAS grid.""" uxgrid = ux.open_grid(mpas_filepath) @@ -611,6 +634,7 @@ def test_connectivity_edge_face_connectivity_mpas(): nt.assert_array_equal(edge_faces_output, edge_faces_gold) + def test_connectivity_edge_face_connectivity_sample(): """Tests the construction of ``Mesh2_face_edges`` on an example with one shared edge, and the remaining edges only being part of one face.""" verts = [[(0.0, -90.0), (180, 0.0), (0.0, 90)], @@ -633,6 +657,7 @@ def test_connectivity_edge_face_connectivity_sample(): assert n_solo == uxgrid.n_edge - n_shared assert n_invalid == 0 + def test_connectivity_face_face_connectivity_construction(): """Tests the construction of face-face connectivity.""" grid = ux.open_grid(mpas_filepath) @@ -645,6 +670,7 @@ def test_connectivity_face_face_connectivity_construction(): nt.assert_array_equal(face_face_conn_new_sorted, face_face_conn_old_sorted) + def test_class_methods_from_dataset(): # UGRID xrds = xr.open_dataset(gridfile_ugrid) @@ -663,6 +689,7 @@ def test_class_methods_from_dataset(): xrds = xr.open_dataset(gridfile_scrip) uxgrid = ux.Grid.from_dataset(xrds) + def test_class_methods_from_face_vertices(): single_face_latlon = [(0.0, 90.0), (-180, 0.0), (0.0, -90)] uxgrid = ux.Grid.from_face_vertices(single_face_latlon, latlon=True) @@ -673,6 +700,7 @@ def test_class_methods_from_face_vertices(): single_face_cart = [(0.0,)] + def test_latlon_bounds_populate_bounds_GCA_mix(): gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc" face_1 = [[10.0, 60.0], [10.0, 10.0], [50.0, 10.0], [50.0, 60.0]] @@ -691,11 +719,13 @@ def test_latlon_bounds_populate_bounds_GCA_mix(): bounds_xarray = grid.bounds nt.assert_allclose(bounds_xarray.values, expected_bounds, atol=ERROR_TOLERANCE) + def test_latlon_bounds_populate_bounds_MPAS(): gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc" uxgrid = ux.open_grid(gridfile_mpas) bounds_xarray = uxgrid.bounds + def test_dual_mesh_mpas(): grid = ux.open_grid(gridfile_mpas, use_dual=False) mpas_dual = ux.open_grid(gridfile_mpas, use_dual=True) @@ -708,11 +738,13 @@ def test_dual_mesh_mpas(): nt.assert_equal(dual.face_node_connectivity.values, mpas_dual.face_node_connectivity.values) + def test_dual_duplicate(): dataset = ux.open_dataset(gridfile_geoflow, gridfile_geoflow) with pytest.raises(RuntimeError): dataset.get_dual() + def test_normalize_existing_coordinates_non_norm_initial(): gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc" from uxarray.grid.validation import _check_normalization @@ -726,9 +758,82 @@ def test_normalize_existing_coordinates_non_norm_initial(): uxgrid.normalize_cartesian_coordinates() assert _check_normalization(uxgrid) + def test_normalize_existing_coordinates_norm_initial(): gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug" from uxarray.grid.validation import _check_normalization uxgrid = ux.open_grid(gridfile_CSne30) assert _check_normalization(uxgrid) + + +def test_number_of_faces_found(): + """Test function for `self.get_face_containing_point`, + to ensure the correct number of faces is found, depending on where the point is.""" + grid = ux.open_grid(gridfile_mpas) + partial_grid = ux.open_grid(quad_hex_grid_path) + + # For a face center only one face should be found + point_xyz = np.array([grid.face_x[100].values, grid.face_y[100].values, grid.face_z[100].values], dtype=np.float64) + + assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 1 + + # For an edge two faces should be found + point_xyz = np.array([grid.edge_x[100].values, grid.edge_y[100].values, grid.edge_z[100].values], dtype=np.float64) + + assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 2 + + # For a node three faces should be found + point_xyz = np.array([grid.node_x[100].values, grid.node_y[100].values, grid.node_z[100].values], dtype=np.float64) + + assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 3 + + partial_grid.normalize_cartesian_coordinates() + + # Test for a node on the edge where only 2 faces should be found + point_xyz = np.array([partial_grid.node_x[1].values, partial_grid.node_y[1].values, partial_grid.node_z[1].values], dtype=np.float64) + + assert len(partial_grid.get_faces_containing_point(point_xyz=point_xyz)) == 2 + + +def test_whole_grid(): + """Tests `self.get_faces_containing_point`on an entire grid, + checking that for each face center, one face is found to contain it""" + + grid = ux.open_grid(gridfile_mpas_two) + grid.normalize_cartesian_coordinates() + + # Ensure a face is found on the grid for every face center + for i in range(len(grid.face_x.values)): + point_xyz = np.array([grid.face_x[i].values, grid.face_y[i].values, grid.face_z[i].values], dtype=np.float64) + + assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 1 + +def test_point_types(): + """Tests that `self.get_faces_containing_point` works with cartesian and lonlat""" + + # Open the grid + grid = ux.open_grid(gridfile_mpas) + + # Assign a cartesian point and a lon/lat point + point_xyz = np.array([grid.node_x[100].values, grid.node_y[100].values, grid.node_z[100].values], dtype=np.float64) + point_lonlat = np.array([grid.node_lon[100].values, grid.node_lat[100].values]) + + # Test both points find faces + assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) != 0 + assert len(grid.get_faces_containing_point(point_lonlat=point_lonlat)) !=0 + +def test_point_along_arc(): + node_lon = np.array([-40, -40, 40, 40]) + node_lat = np.array([-20, 20, 20, -20]) + face_node_connectivity = np.array([[0, 1, 2, 3]], dtype=np.int64) + + uxgrid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity) + + # point at exactly 20 degrees latitude + out1 = uxgrid.get_faces_containing_point(point_lonlat=np.array([0, 20], dtype=np.float64)) + + # point at 25.41 degrees latitude (max along the great circle arc) + out2 = uxgrid.get_faces_containing_point(point_lonlat=np.array([0, 25.41], dtype=np.float64)) + + nt.assert_array_equal(out1, out2) diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py index 2d78b978a..714de0fce 100644 --- a/uxarray/grid/coordinates.py +++ b/uxarray/grid/coordinates.py @@ -100,8 +100,8 @@ def _xyz_to_lonlat_rad( z: Union[np.ndarray, float], normalize: bool = True, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Converts Cartesian x, y, z coordinates in Spherical latitude and - longitude coordinates in degrees. + """Converts Cartesian x, y, z coordinates in Spherical longitude and + latitude coordinates in radians. Parameters ---------- @@ -279,7 +279,7 @@ def _populate_face_centroids(grid, repopulate=False): # Convert to xyz if there are latlon centroids already stored centroid_lon, centroid_lat = grid.face_lon.values, grid.face_lat.values centroid_x, centroid_y, centroid_z = _lonlat_rad_to_xyz( - centroid_lon, centroid_lat + np.deg2rad(centroid_lon), np.deg2rad(centroid_lat) ) # Populate the centroids @@ -648,7 +648,7 @@ def _populate_edge_centroids(grid, repopulate=False): # Convert to xyz if there are latlon centroids already stored centroid_lon, centroid_lat = grid.edge_lon.values, grid.edge_lat.values centroid_x, centroid_y, centroid_z = _lonlat_rad_to_xyz( - centroid_lon, centroid_lat + np.deg2rad(centroid_lon), np.deg2rad(centroid_lat) ) # Populate the centroids diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py index f39120fed..1c5660da6 100644 --- a/uxarray/grid/geometry.py +++ b/uxarray/grid/geometry.py @@ -18,9 +18,16 @@ INT_DTYPE, INT_FILL_VALUE, ) -from uxarray.grid.arcs import extreme_gca_latitude, point_within_gca +from uxarray.grid.arcs import ( + extreme_gca_latitude, + point_within_gca, +) + from uxarray.grid.coordinates import _xyz_to_lonlat_rad -from uxarray.grid.intersections import gca_gca_intersection + +from uxarray.grid.intersections import ( + gca_gca_intersection, +) from uxarray.grid.utils import ( _get_cartesian_face_edge_nodes, _get_lonlat_rad_face_edge_nodes, @@ -31,6 +38,10 @@ "North": np.array([0.0, 0.0, 1.0]), "South": np.array([0.0, 0.0, -1.0]), } + +REF_POINT_NORTH_XYZ = np.array([0.01745241, 0.0, 0.9998477], dtype=np.float64) +REF_POINT_SOUTH_XYZ = np.array([0.01745241, 0.0, -0.9998477], dtype=np.float64) + POLE_POINTS_LONLAT = { "North": np.array([0.0, np.pi / 2]), "South": np.array([0.0, -np.pi / 2]), @@ -1561,3 +1572,234 @@ def inverse_stereographic_projection(x, y, central_lon, central_lat): ) return lon, lat + + +@njit(cache=True) +def point_in_face( + edges_xyz, + point_xyz, + inclusive=True, +): + """Determines if a point lies inside a face. + + Parameters + ---------- + edges_xyz : numpy.ndarray + Cartesian coordinates of each point in the face + point_xyz : numpy.ndarray + Cartesian coordinate of the point + inclusive : bool + Flag to determine whether to include points on the nodes and edges of the face + + Returns + ------- + bool + True if point is inside face, False otherwise + """ + + # Validate the inputs + if len(edges_xyz[0][0]) != 3: + raise ValueError("`edges_xyz` vertices must be in Cartesian coordinates.") + + if len(point_xyz) != 3: + raise ValueError("`point_xyz` must be a single [3] Cartesian coordinate.") + + # Initialize the intersection count + intersection_count = 0 + + # Set to hold unique intersections + unique_intersections = set() + + location = _classify_polygon_location(edges_xyz) + + if location == 1: + ref_point_xyz = REF_POINT_SOUTH_XYZ + elif location == -1: + ref_point_xyz = REF_POINT_NORTH_XYZ + else: + ref_point_xyz = REF_POINT_SOUTH_XYZ + + # Initialize the points arc between the point and the reference point + gca_cart = np.empty((2, 3), dtype=np.float64) + gca_cart[0] = point_xyz + gca_cart[1] = ref_point_xyz + + # Loop through the face's edges, checking each one for intersection + for ind in range(len(edges_xyz)): + # If the point lies on an edge, return True if inclusive + if point_within_gca( + point_xyz, + edges_xyz[ind][0], + edges_xyz[ind][1], + ): + if inclusive: + return True + else: + return False + + # Get the number of intersections between the edge and the point arc + intersections = gca_gca_intersection(edges_xyz[ind], gca_cart) + + # Add any unique intersections to the intersection_count + for intersection in intersections: + intersection_tuple = ( + intersection[0], + intersection[1], + intersection[2], + ) + if intersection_tuple not in unique_intersections: + unique_intersections.add(intersection_tuple) + intersection_count += 1 + + # Return True if the number of intersections is odd, False otherwise + return intersection_count % 2 == 1 + + +@njit(cache=True) +def _find_faces(face_edge_cartesian, point_xyz, inverse_indices): + """Finds the faces that contain a given point, inside a subset `face_edge_cartesian` + Parameters + ---------- + face_edge_cartesian : numpy.ndarray + Cartesian coordinates of all the faces according to their edges + point_xyz : numpy.ndarray + Cartesian coordinate of the point + inverse_indices : numpy.ndarray + The original indices of the subsetted grid + + Returns + ------- + index : array + The index of the face that contains the point + """ + + index = [] + + # Run for each face in the subset + for i, face in enumerate(inverse_indices): + # Check to see if the face contains the point + contains_point = point_in_face( + face_edge_cartesian[i], + point_xyz, + inclusive=True, + ) + + # If the point is found, add it to the index array + if contains_point: + index.append(face) + + # Return the index array + return index + + +def _populate_max_face_radius(self): + """Populates `max_face_radius` + + Returns + ------- + max_distance : np.float64 + The max distance from a node to a face center + """ + + # Parse all variables needed for `njit` functions + face_node_connectivity = self.face_node_connectivity.values + node_lats_rad = np.deg2rad(self.node_lat.values) + node_lons_rad = np.deg2rad(self.node_lon.values) + face_lats_rad = np.deg2rad(self.face_lat.values) + face_lons_rad = np.deg2rad(self.face_lon.values) + + # Get the max distance + max_distance = calculate_max_face_radius( + face_node_connectivity, + node_lats_rad, + node_lons_rad, + face_lats_rad, + face_lons_rad, + ) + + # Return the max distance, which is the `max_face_radius` + return np.rad2deg(max_distance) + + +@njit(cache=True) +def calculate_max_face_radius( + face_node_connectivity, node_lats_rad, node_lons_rad, face_lats_rad, face_lons_rad +): + """Finds the max face radius in the mesh. + Parameters + ---------- + face_node_connectivity : numpy.ndarray + Cartesian coordinates of all the faces according to their edges + node_lats_rad : numpy.ndarray + The `Grid.node_lat` array in radians + node_lons_rad : numpy.ndarray + The `Grid.node_lon` array in radians + face_lats_rad : numpy.ndarray + The `Grid.face_lat` array in radians + face_lons_rad : numpy.ndarray + The `Grid.face_lon` array in radians + + Returns + ------- + The max distance from a node to a face center + """ + + # Array to store all distances of each face to it's furthest node. + end_distances = np.zeros(len(face_node_connectivity)) + + # Loop over each face and its nodes + for ind, face in enumerate(face_node_connectivity): + # Filter out INT_FILL_VALUE + valid_nodes = face[face != INT_FILL_VALUE] + + # Get the face lat/lon of this face + face_lat = face_lats_rad[ind] + face_lon = face_lons_rad[ind] + + # Get the node lat/lon of this face + node_lat_rads = node_lats_rad[valid_nodes] + node_lon_rads = node_lons_rad[valid_nodes] + + # Calculate Haversine distances for all nodes in this face + distances = haversine_distance(node_lon_rads, node_lat_rads, face_lon, face_lat) + + # Store the max distance for this face + end_distances[ind] = np.max(distances) + + # Return the maximum distance found across all faces + return np.max(end_distances) + + +@njit(cache=True) +def haversine_distance(lon_a, lat_a, lon_b, lat_b): + """Calculates the haversine distance between two points. + + Parameters + ---------- + lon_a : np.float64 + The longitude of the first point + lat_a : np.float64 + The latitude of the first point + lon_b : np.float64 + The longitude of the second point + lat_b : np.float64 + The latitude of the second point + + Returns + ------- + distance : np.float64 + The distance between the two points + """ + + # Differences in latitudes and longitudes + dlat = lat_b - lat_a + dlon = lon_b - lon_a + + # Haversine formula + equation_in_sqrt = (np.sin(dlat / 2) ** 2) + np.cos(lat_a) * np.cos(lat_b) * ( + np.sin(dlon / 2) ** 2 + ) + distance = 2 * np.arcsin(np.sqrt(equation_in_sqrt)) + + # Return the gotten distance + return distance diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index f927dc4f9..44bdc37fa 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -14,6 +14,8 @@ Tuple, ) +from uxarray.grid.utils import _get_cartesian_face_edge_nodes + # reader and writer imports from uxarray.io._exodus import _read_exodus, _encode_exodus from uxarray.io._mpas import _read_mpas @@ -51,6 +53,8 @@ _populate_node_xyz, _normalize_xyz, prepare_points, + _lonlat_rad_to_xyz, + _xyz_to_lonlat_deg, ) from uxarray.grid.connectivity import ( _populate_edge_node_connectivity, @@ -69,6 +73,8 @@ _populate_bounds, _construct_boundary_edge_indices, compute_temp_latlon_array, + _find_faces, + _populate_max_face_radius, ) from uxarray.grid.neighbors import ( @@ -118,7 +124,7 @@ import copy -from uxarray.constants import INT_FILL_VALUE +from uxarray.constants import INT_FILL_VALUE, ERROR_TOLERANCE from uxarray.grid.dual import construct_dual @@ -1582,6 +1588,13 @@ def is_subset(self): """Returns `True` if the Grid is a subset, 'False' otherwise.""" return self._is_subset + @property + def max_face_radius(self): + """Maximum face radius of the grid in degrees""" + if "max_face_radius" not in self._ds: + self._ds["max_face_radius"] = _populate_max_face_radius(self) + return self._ds["max_face_radius"] + def chunk(self, n_node="auto", n_edge="auto", n_face="auto"): """Converts all arrays to dask arrays with given chunks across grid dimensions in-place. @@ -2485,3 +2498,118 @@ def get_faces_between_latitudes(self, lats: Tuple[float, float]): """ return faces_within_lat_bounds(lats, self.face_bounds_lat.values) + + def get_faces_containing_point( + self, point_xyz=None, point_lonlat=None, tolerance=ERROR_TOLERANCE + ): + """Identifies the indices of faces that contain a given point. + + Parameters + ---------- + point_xyz : numpy.ndarray + A point in cartesian coordinates. Best performance if + point_lonlat : numpy.ndarray + A point in spherical coordinates. + tolerance : numpy.ndarray + An optional error tolerance for points that lie on the nodes of a face. + + Returns + ------- + index : numpy.ndarray + Array of the face indices containing point. Empty if no face is found. This function will typically return + a single face, unless the point falls directly on a corner or edge, where there will be multiple values. + + Examples + -------- + Open a grid from a file path: + + >>> import uxarray as ux + >>> uxgrid = ux.open_grid("grid_filename.nc") + + Define a spherical point: + + >>> import numpy as np + >>> point_lonlat = np.array([45.2, 32.6], dtype=np.float64) + + Define a cartesian point: + + >>> point_xyz = np.array([0.0, 0.0, 1.0], dtype=np.float64) + + Find the indices of the faces that contain the given point: + + >>> lonlat_point_face_indices = uxgrid.get_faces_containing_point( + ... point_lonlat=point_lonlat + ... ) + >>> xyz_point_face_indices = uxgrid.get_faces_containing_point( + ... point_xyz=point_xyz + ... ) + + """ + if point_xyz is None and point_lonlat is None: + raise ValueError("Either `point_xyz` or `point_lonlat` must be passed in.") + + # Depending on the provided point coordinates, convert to get all needed coordinate systems + if point_xyz is None: + point_lonlat = np.asarray(point_lonlat, dtype=np.float64) + point_xyz = np.array( + _lonlat_rad_to_xyz(*np.deg2rad(point_lonlat)), dtype=np.float64 + ) + elif point_lonlat is None: + point_xyz = np.asarray(point_xyz, dtype=np.float64) + point_lonlat = np.array(_xyz_to_lonlat_deg(*point_xyz), dtype=np.float64) + + # Get the maximum face radius of the grid, plus a small adjustment for if the point is this exact radius away + max_face_radius = self.max_face_radius.values + 0.0001 + + # Try to find a subset in which the point resides + try: + subset = self.subset.bounding_circle( + r=max_face_radius, + center_coord=point_lonlat, + element="face centers", + inverse_indices=True, + ) + # If no subset is found, warn the user + except ValueError: + # If the grid is partial, let the user know the point likely lies outside the grid region + if self.partial_sphere_coverage: + warn( + "No faces found. The grid has partial spherical coverage, and the point may be outside the defined region of the grid." + ) + else: + warn("No faces found. Try adjusting the tolerance.") + return np.empty(0, dtype=np.int64) + + # Get the faces in terms of their edges + face_edge_nodes_xyz = _get_cartesian_face_edge_nodes( + subset.face_node_connectivity.values, + subset.n_face, + subset.n_max_face_nodes, + subset.node_x.values, + subset.node_y.values, + subset.node_z.values, + ) + + # Get the original face indices from the subset + inverse_indices = subset.inverse_indices.face.values + + # Check to see if the point is on the nodes of any face + lies_on_node = np.isclose( + face_edge_nodes_xyz, + point_xyz[None, None, :], # Expands dimensions for broadcasting + rtol=tolerance, + atol=tolerance, + ) + + edge_matches = np.all(lies_on_node, axis=-1) + face_matches = np.any(edge_matches, axis=1) + face_indices = inverse_indices[np.any(face_matches, axis=1)] + + # If a face is in face_indices, return that as the point was found to lie on a node + if len(face_indices) != 0: + return face_indices + else: + # Check if any of the faces in the subset contain the point + face_indices = _find_faces(face_edge_nodes_xyz, point_xyz, inverse_indices) + + return face_indices