Check Grid for Partial Spherical Coverage (#899)

* Added `contains_holes` function

* o vectorize _mesh_contains_holes

* Addressed discussion points

* Updated API

* Fixed Pre-Commit

* Updated function names and docstrings

* Fixed test case

* Update uxarray/grid/grid.py

Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

* Added descriptors, benchmarks

* Update descriptors.py

* Fixed Benchmarks

* Fixed Benchmarks

---------

Co-authored-by: Rajeev Jain <rajeeja@gmail.com>
Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com>

benchmarks/mpas_ocean.py

@@ -164,3 +164,17 @@ def time_nearest_neighbor_remapping(self):
 
     def time_inverse_distance_weighted_remapping(self):
         self.uxds_480["bottomDepth"].remap.inverse_distance_weighted(self.uxds_120.uxgrid)
+
+
+class HoleEdgeIndices:
+    param_names = ['resolution']
+    params = ['480km', '120km']
+
+    def setup(self, resolution):
+        self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])
+
+    def teardown(self, resolution):
+        del self.uxds
+
+    def time_construct_hole_edge_indices(self, resolution):
+        ux.grid.geometry._construct_hole_edge_indices(self.uxds.uxgrid.edge_face_connectivity)

docs/internal_api/index.rst

@@ -141,6 +141,7 @@ Geometry
    grid.geometry._insert_pt_in_latlonbox
    grid.geometry._populate_face_latlon_bound
    grid.geometry._populate_bounds
+   grid.geometry._construct_hole_edge_indices
 
 Coordinates
 -----------

docs/user_api/index.rst

@@ -210,6 +210,7 @@ IO
    Grid.to_polycollection
    Grid.to_linecollection
    Grid.validate
+   Grid.hole_edge_indices
 
 
 Methods

test/test_grid.py

@@ -31,6 +31,7 @@
 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_holes = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc'
 
 dsfile_vortex_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_vortex.nc"
 dsfile_var2_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30_var2.nc"
@@ -48,6 +49,14 @@ def test_validate(self):
         grid_mpas = ux.open_grid(gridfile_mpas)
         assert (grid_mpas.validate())
 
+    def test_grid_with_holes(self):
+        """Test _holes_in_mesh function."""
+        grid_without_holes = ux.open_grid(gridfile_mpas)
+        grid_with_holes = ux.open_grid(gridfile_mpas_holes)
+
+        self.assertTrue(grid_with_holes.hole_edge_indices.size != 0)
+        self.assertTrue(grid_without_holes.hole_edge_indices.size == 0)
+
     def test_encode_as(self):
         """Reads a ugrid file and encodes it as `xarray.Dataset` in various
         types."""
@@ -658,8 +667,8 @@ def test_build_face_edges_connectivity_mpas(self):
 
         # construct edge nodes
         edge_nodes_output, _, _ = _build_edge_node_connectivity(mpas_grid_ux.face_node_connectivity.values,
-                                                          mpas_grid_ux.n_face,
-                                                          mpas_grid_ux.n_max_face_nodes)
+                                                                mpas_grid_ux.n_face,
+                                                                mpas_grid_ux.n_max_face_nodes)
 
         # _populate_face_edge_connectivity(mpas_grid_ux)
         # edge_nodes_output = mpas_grid_ux._ds['edge_node_connectivity'].values
@@ -932,6 +941,7 @@ def test_from_face_vertices(self):
 
 class TestLatlonBounds(TestCase):
     gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
+
     def test_populate_bounds_GCA_mix(self):
         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]]
@@ -941,11 +951,10 @@ def test_populate_bounds_GCA_mix(self):
         faces = [face_1, face_2, face_3, face_4]
 
         # Hand calculated bounds for the above faces in radians
-        expected_bounds = [[[0.17453293, 1.07370494],[0.17453293, 0.87266463]],
-                           [[0.17453293, 1.10714872],[6.10865238, 0.87266463]],
-                           [[1.04719755, 1.57079633],[3.66519143, 0.52359878]],
-                           [[1.04719755,1.57079633],[0.,         6.28318531]]]
-
+        expected_bounds = [[[0.17453293, 1.07370494], [0.17453293, 0.87266463]],
+                           [[0.17453293, 1.10714872], [6.10865238, 0.87266463]],
+                           [[1.04719755, 1.57079633], [3.66519143, 0.52359878]],
+                           [[1.04719755, 1.57079633], [0., 6.28318531]]]
 
         grid = ux.Grid.from_face_vertices(faces, latlon=True)
         bounds_xarray = grid.bounds

uxarray/conventions/descriptors.py

@@ -25,3 +25,9 @@
     "cf_role": "edge_node_distances",
     "long_name": "Distances between the nodes that make up " "each edge.",
 }
+
+HOLE_EDGE_INDICES_DIMS = ["n_edge"]
+HOLE_EDGE_INDICES_ATTRS = {
+    "cf_role": "hole_edge_indices",
+    "long_name": "Indices of edges that border a region of the grid not covered by any geometry",
+}

uxarray/grid/geometry.py

@@ -1021,3 +1021,12 @@ def _populate_bounds(
         return bounds
     else:
         grid._ds["bounds"] = bounds
+
+
+def _construct_hole_edge_indices(edge_face_connectivity):
+    """Index the missing edges on a partial grid with holes, that is a region
+    of the grid that is not covered by any geometry."""
+
+    # If an edge only has one face saddling it than the mesh has holes in it
+    edge_with_holes = np.where(edge_face_connectivity[:, 1] == INT_FILL_VALUE)[0]
+    return edge_with_holes

uxarray/grid/grid.py

@@ -54,6 +54,7 @@
     _grid_to_matplotlib_polycollection,
     _grid_to_matplotlib_linecollection,
     _populate_bounds,
+    _construct_hole_edge_indices,
 )
 
 from uxarray.grid.neighbors import (
@@ -1085,6 +1086,16 @@ def face_jacobian(self):
             _ = self.face_areas
         return self._face_jacobian
 
+    @property
+    def hole_edge_indices(self):
+        """Indices of edges that border a region of the grid not covered by any
+        geometry, know as a hole, in a partial grid."""
+        if "hole_edge_indices" not in self._ds:
+            self._ds["hole_edge_indices"] = _construct_hole_edge_indices(
+                self.edge_face_connectivity.values
+            )
+        return self._ds["hole_edge_indices"]
+
     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.

uxarray/grid/validation.py

@@ -1,6 +1,7 @@
 import numpy as np
 from warnings import warn
 
+
 from uxarray.constants import ERROR_TOLERANCE