Point In Face

benchmarks/mpas_ocean.py

@@ -183,3 +183,10 @@ def teardown(self, resolution, lat_step):
     def time_const_lat(self, resolution, lat_step):
         for lat in self.lats:
             self.uxgrid.cross_section.constant_latitude(lat)
+
+
+class PointInPolygon(GridBenchmark):
+    def time_whole_grid(self, resolution):
+        point_xyz = np.array([self.uxgrid.face_x[0].values, self.uxgrid.face_y[0].values, self.uxgrid.face_z[0].values])
+
+        self.uxgrid.get_faces_containing_point(point_xyz=point_xyz)

test/test_grid.py

@@ -13,7 +13,7 @@
 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
 
@@ -44,13 +44,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"
@@ -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,39 @@ 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():
+    grid = ux.open_grid(gridfile_mpas)
+
+    # 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])
+
+    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])
+
+    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])
+
+    assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 3
+
+
+def test_whole_grid():
+    grid = ux.open_grid(gridfile_mpas)
+
+    # 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])
+
+        assert len(grid.get_faces_containing_point(point_xyz=point_xyz)) == 1

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
     ----------