Repository: UXARRAY/uxarray

Issue:
Function to Check For Meshes With Holes
### Proposed new feature or change:

Function that returns whether or not a mesh has holes in it. Something like this:

`Grid.contains_holes()`

Returns True if any holes exists and False if not.

Hints from issue comments:
@rajeeja: Why can one function not work for both of them?

On Thu, Aug 22, 2024, 9:42 AM Aaron Zedwick ***@***.***>
wrote:

> Proposed new feature or change:
>
> Function that returns whether or not a mesh has holes in it. Something
> like this:
>
> Grid.validate.holes_in_mesh()
>
> Returns True if any holes exists and False if not.
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/UXARRAY/uxarray/issues/898>, or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AALF6AR4XZLIY5REUY2JHE3ZSX2GVAVCNFSM6AAAAABM6J4QE6VHI2DSMVQWIX3LMV43ASLTON2WKOZSGQ4DAOJXGI2DMNY>
> .
> You are receiving this because you are subscribed to this thread.Message
> ID: ***@***.***>
>
@aaronzedwick: > Why can one function not work for both of them?
> […](#)
> On Thu, Aug 22, 2024, 9:42 AM Aaron Zedwick ***@***.***> wrote: Proposed new feature or change: Function that returns whether or not a mesh has holes in it. Something like this: Grid.validate.holes_in_mesh() Returns True if any holes exists and False if not. — Reply to this email directly, view it on GitHub <#898>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AALF6AR4XZLIY5REUY2JHE3ZSX2GVAVCNFSM6AAAAABM6J4QE6VHI2DSMVQWIX3LMV43ASLTON2WKOZSGQ4DAOJXGI2DMNY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

A grid with holes in it is a valid grid, so putting under the normal validation function doesn't make sense. I just need a function to return whether a grid has holes or not for the dual mesh, because if a grid has holes in it, it should be approached differently especially when it comes to DataArrays. @rajeeja: > 
> A grid with holes in it is a valid grid, so putting under the normal validation function doesn't make sense. I just need a function to return whether a grid has holes or not for the dual mesh, because if a grid has holes in it, it should be approached differently especially when it comes to DataArrays.

`validation` if you see the file `grid/validation.py` there are functions such as `_check_connectivity _check_duplicates` etc., you can add something to detect holes `_detect_holes` a new function is not required.

The problem is the dual mesh must be constructed in a way that respects the topology of the original mesh, which can be non-trivial. What you really need is the idea of "geometry" and "mesh", it might be doable without the computational geometry here as we limit ourselves to spherical and earth type topology/meshes.  
@aaronzedwick: > > A grid with holes in it is a valid grid, so putting under the normal validation function doesn't make sense. I just need a function to return whether a grid has holes or not for the dual mesh, because if a grid has holes in it, it should be approached differently especially when it comes to DataArrays.
> 
> `validation` if you see the file `grid/validation.py` there are functions such as `_check_connectivity _check_duplicates` etc., you can add something to detect holes `_detect_holes` a new function is not required.
> 
> The problem is the dual mesh must be constructed in a way that respects the topology of the original mesh, which can be non-trivial. What you really need is the idea of "geometry" and "mesh", it might be doable without the computational geometry here as we limit ourselves to spherical and earth type topology/meshes.

I opened a PR earlier with my implementation idea, could you give it a review so I can better understand how you want it to be organized? I added a similar thing as what you said above, but I don’t know that it’s exactly what you are saying. @aaronzedwick: #899 is the PR number, I changed the name of the issue as I think it might have been poorly worded and possibly confusing. I don't want to add a new validation function, I simply just want to add a check for holes function inside of `validation.py`@philipc2: This is a much needed check that we should include. 

Looping @hongyuchen1030 in for any inputs on this. We can chat about the algorithm in next weeks UXarray meeting. @hongyuchen1030: > This is a much needed check that we should include.
> 
> Looping @hongyuchen1030 in for any inputs on this. We can chat about the algorithm in next weeks UXarray meeting.

All face in the Uxarray grid should be **convex**. And this mean there will not be any "hole" inside
And before answering the  "is Hole inside" algorithm, please think about the following question first:
1. How to detect if a point is inside a face on the spherical surface. (I know the immediate answer is ray-casting, but we are on the spherical surface, so the ray will go back!)
2. How to know an interval is "inside" a face。
3. Why we need this "is Hole inside" function?@hongyuchen1030: @aaronzedwick And I also add a comment for the current implementation in the PR https://github.com/UXARRAY/uxarray/pull/899#issuecomment-2307786474

Such implementation idea is not checking the hole, but to see if a face has **FEWER** nodes than other faces. @rajeeja: > @aaronzedwick And I also add a comment for the current implementation in the PR [#899 (comment)](https://github.com/UXARRAY/uxarray/pull/899#issuecomment-2307786474)
> 
> Such implementation idea is not checking the hole, but to see if a face has **FEWER** nodes than other faces.

Oh, in that case it'd come out as True for a mesh with faces with different numbers of edges (mixed mesh, eg. tri/quad/hex etc.)@hongyuchen1030: > @aaronzedwick And I also add a comment for the current implementation in the PR [#899 (comment)](https://github.com/UXARRAY/uxarray/pull/899#issuecomment-2307786474)
> 
> Such implementation idea is not checking the hole, but to see if a face has **FEWER** nodes than other faces.

Sorry for the misunderstanding about the purpose of such a function. The updated comment listed here. https://github.com/UXARRAY/uxarray/pull/899#issuecomment-2307975194

Basically my question is why this function and what's the following steps for such function.

Solution patch:
diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py
index d2baaf86b..51b0c99d2 100644
--- a/benchmarks/mpas_ocean.py
+++ b/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)
diff --git a/docs/internal_api/index.rst b/docs/internal_api/index.rst
index b263f4ec8..6dd40b2c0 100644
--- a/docs/internal_api/index.rst
+++ b/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
 -----------
diff --git a/docs/user_api/index.rst b/docs/user_api/index.rst
index ec3278b6d..da904fe0d 100644
--- a/docs/user_api/index.rst
+++ b/docs/user_api/index.rst
@@ -210,6 +210,7 @@ IO
    Grid.to_polycollection
    Grid.to_linecollection
    Grid.validate
+   Grid.hole_edge_indices
 
 
 Methods
diff --git a/uxarray/conventions/descriptors.py b/uxarray/conventions/descriptors.py
index 2cd0eb2a7..1d020b913 100644
--- a/uxarray/conventions/descriptors.py
+++ b/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",
+}
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index 08229d1a0..afa31909d 100644
--- a/uxarray/grid/geometry.py
+++ b/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
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 97b9f5d1f..b1ff3e43c 100644
--- a/uxarray/grid/grid.py
+++ b/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.
diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py
index 1a4852711..4f74249a5 100644
--- a/uxarray/grid/validation.py
+++ b/uxarray/grid/validation.py
@@ -1,6 +1,7 @@
 import numpy as np
 from warnings import warn
 
+
 from uxarray.constants import ERROR_TOLERANCE
 
 


Test patch:
diff --git a/test/test_grid.py b/test/test_grid.py
index 57eb3f3b8..6b80aa87a 100644
--- a/test/test_grid.py
+++ b/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