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.

Discussion:
hongyuchen1030: ```python
def _mesh_contains_holes(edge_face_connectivity):
    """Check if a mesh has holes in it."""

    # If an edge only has one face saddling it than the mesh has holes in it
    return np.any(edge_face_connectivity[:, 1] == INT_FILL_VALUE)
```

This function is basically checking if the `edge_face_connectivity` has a **FILL VALUE**, which means it has fewer nodes than other faces in the mesh files, **it doesn't mean it has a hole.**
aaronzedwick: > ```python
> def _mesh_contains_holes(edge_face_connectivity):
>     """Check if a mesh has holes in it."""
> 
>     # If an edge only has one face saddling it than the mesh has holes in it
>     return np.any(edge_face_connectivity[:, 1] == INT_FILL_VALUE)
> ```
> 
> This function is basically checking if the `edge_face_connectivity` has a **FILL VALUE**, which means it has fewer nodes than other faces in the mesh files, **it doesn't mean it has a hole.**

Doesn’t this connective hold for each edge what face saddles it? If it is the case that an edge has only one face and the other side of the edge has a fill value, that means it lacks a face there and therefore there is a hole in the mesh? For each edge it can only ever have 2 faces that saddle it, and if there is only one face saddling it, the logic goes that there must be a hole. Is that assumption incorrect?
hongyuchen1030: > > ```python
> > def _mesh_contains_holes(edge_face_connectivity):
> >     """Check if a mesh has holes in it."""
> > 
> >     # If an edge only has one face saddling it than the mesh has holes in it
> >     return np.any(edge_face_connectivity[:, 1] == INT_FILL_VALUE)
> > ```
> > 
> > 
> >     
> >       
> >     
> > 
> >       
> >     
> > 
> >     
> >   
> > This function is basically checking if the `edge_face_connectivity` has a **FILL VALUE**, which means it has fewer nodes than other faces in the mesh files, **it doesn't mean it has a hole.**
> 
> Doesn’t this connective hold for each edge what face saddles it? If it is the case that an edge has only one face and the other side of the edge has a fill value, that means it lacks a face there and therefore there is a hole in the mesh? For each edge it can only ever have 2 faces that saddle it, and if there is only one face saddling it, the logic goes that there must be a hole. Is that assumption incorrect?

Oh, so you want to see if a given **mesh file** contains holes, instead of a single face has holes inside it? 
So yeah, your implementation is propably correct for this purpose then. 

However, why we need to run this check over the entire grid?

So basically what happens is: you check the ENTIRE grid to see if it has a hole or not, which will still at least needs an O(log_2 N_Edges) time complexity. And then what will happen if we detect there's at least one hole in the face? What's the following procedure? Are you planning to do another query about where the hole is? 
aaronzedwick: > > > ```python
> > > def _mesh_contains_holes(edge_face_connectivity):
> > >     """Check if a mesh has holes in it."""
> > > 
> > >     # If an edge only has one face saddling it than the mesh has holes in it
> > >     return np.any(edge_face_connectivity[:, 1] == INT_FILL_VALUE)
> > > ```
> > > 
> > > 
> > >     
> > >       
> > >     
> > > 
> > >       
> > >     
> > > 
> > >     
> > >   
> > > This function is basically checking if the `edge_face_connectivity` has a **FILL VALUE**, which means it has fewer nodes than other faces in the mesh files, **it doesn't mean it has a hole.**
> > 
> > 
> > Doesn’t this connective hold for each edge what face saddles it? If it is the case that an edge has only one face and the other side of the edge has a fill value, that means it lacks a face there and therefore there is a hole in the mesh? For each edge it can only ever have 2 faces that saddle it, and if there is only one face saddling it, the logic goes that there must be a hole. Is that assumption incorrect?
> 
> Oh, so you want to see if a given **mesh file** contains holes, instead of a single face has holes inside it? So yeah, your implementation is propably correct for this purpose then.
> 
> However, why we need to run this check over the entire grid?
> 
> So basically what happens is: you check the ENTIRE grid to see if it has a hole or not, which will still at least needs an O(log_2 N_Edges) time complexity. And then what will happen if we detect there's at least one hole in the face? What's the following procedure? Are you planning to do another query about where the hole is?

Exactly. The specific use case for this is within the dual mesh. If the grid has holes (like say it doesn’t include landmasses) then the dual mesh of a data array would give incorrect results, as some face centers and nodes won’t be used in the dual mesh (imagine if there was a line of squares in a row, no dual mesh could be made from just that and therefore all data stored there would be lost). Also, for bilinear remapping, if the mesh has holes in it, the construction of the dual mesh needs to be handled differently, in accordance to what Paul does in his implementation. So inside the bilinear remapping a call to this function must be made to determine if a simply dual mesh will work, or if a different approach to constructing the dual must be taken. I hope that helps clear things up!
hongyuchen1030: > > > > ```python
> > > > def _mesh_contains_holes(edge_face_connectivity):
> > > >     """Check if a mesh has holes in it."""
> > > > 
> > > >     # If an edge only has one face saddling it than the mesh has holes in it
> > > >     return np.any(edge_face_connectivity[:, 1] == INT_FILL_VALUE)
> > > > ```
> > > > 
> > > > 
> > > >     
> > > >       
> > > >     
> > > > 
> > > >       
> > > >     
> > > > 
> > > >     
> > > >   
> > > > This function is basically checking if the `edge_face_connectivity` has a **FILL VALUE**, which means it has fewer nodes than other faces in the mesh files, **it doesn't mean it has a hole.**
> > > 
> > > 
> > > Doesn’t this connective hold for each edge what face saddles it? If it is the case that an edge has only one face and the other side of the edge has a fill value, that means it lacks a face there and therefore there is a hole in the mesh? For each edge it can only ever have 2 faces that saddle it, and if there is only one face saddling it, the logic goes that there must be a hole. Is that assumption incorrect?
> > 
> > 
> > Oh, so you want to see if a given **mesh file** contains holes, instead of a single face has holes inside it? So yeah, your implementation is propably correct for this purpose then.
> > However, why we need to run this check over the entire grid?
> > So basically what happens is: you check the ENTIRE grid to see if it has a hole or not, which will still at least needs an O(log_2 N_Edges) time complexity. And then what will happen if we detect there's at least one hole in the face? What's the following procedure? Are you planning to do another query about where the hole is?
> 
> Exactly. The specific use case for this is within the dual mesh. If the grid has holes (like say it doesn’t include landmasses) then the dual mesh of a data array would give incorrect results, as some face centers and nodes won’t be used in the dual mesh (imagine if there was a line of squares in a row, no dual mesh could be made from just that and therefore all data stored there would be lost). Also, for bilinear remapping, if the mesh has holes in it, the construction of the dual mesh needs to be handled differently, in accordance to what Paul does in his implementation. So inside the bilinear remapping a call to this function must be made to determine if a global dual mesh will work, or if a different local approach to constructing the dual must be taken. I hope that helps clear things up!

So yes, my concern was true. After we call this boolean function to check the entire grid and return False, we need to check the entire grid again to find where the holes are missing. So two same traverse for the entire grid, which is very computationally expensive.

Why don't we just return the clean-up meshes or identify the missing faces during the re-mapping process?

During the `zonal_mean` implementation, we encountered a similar problem and it can be solved using the sweep-line algorithm without checking if such grid contains a missing face
philipc2: I'd also avoid referring to this behavior as "holes", as it can be confused with what we've discussed. We are essentially just checking if the entire sphere is covered with our grid. A grid with partial coverage, such as the grids used in ocean models, only cover a fraction of the sphere. 
hongyuchen1030: > I'd also avoid referring to this behavior as "holes", as it can be confused with what we've discussed. We are essentially just checking if the entire sphere is covered with our grid. A grid with partial coverage, such as the grids used in ocean models, only cover a fraction of the sphere.

And partial grid should not impact the re-gridding behavior as well.
aaronzedwick: > > I'd also avoid referring to this behavior as "holes", as it can be confused with what we've discussed. We are essentially just checking if the entire sphere is covered with our grid. A grid with partial coverage, such as the grids used in ocean models, only cover a fraction of the sphere.
> 
> And partial grid should not impact the re-gridding behavior as well.

https://github.com/ClimateGlobalChange/tempestremap/blob/6f1b7839a919482606caee6ca5ffb114f2241e29/src/LinearRemapFV.cpp#L2996

Paul does a check for it in his code, unless I am misunderstanding what he does there.
philipc2: > > > I'd also avoid referring to this behavior as "holes", as it can be confused with what we've discussed. We are essentially just checking if the entire sphere is covered with our grid. A grid with partial coverage, such as the grids used in ocean models, only cover a fraction of the sphere.
> > 
> > And partial grid should not impact the re-gridding behavior as well.
> 
> https://github.com/ClimateGlobalChange/tempestremap/blob/6f1b7839a919482606caee6ca5ffb114f2241e29/src/LinearRemapFV.cpp#L2996
> 
> Paul does a check for it in his code, unless I am misunderstanding what he does there. Unless he is checking for holes inside a face but looking at his code for `DoesMeshContain` holes it seems he does a similar thing as me. Also, @philipc2 the face area idea would work but seems to me that it would take much longer as you need to calculate the area for each face.

That's a fair point. However if we already have the face area, it would be an O(n) operation.
hongyuchen1030: > > > I'd also avoid referring to this behavior as "holes", as it can be confused with what we've discussed. We are essentially just checking if the entire sphere is covered with our grid. A grid with partial coverage, such as the grids used in ocean models, only cover a fraction of the sphere.
> > 
> > 
> > And partial grid should not impact the re-gridding behavior as well.
> 
> https://github.com/ClimateGlobalChange/tempestremap/blob/6f1b7839a919482606caee6ca5ffb114f2241e29/src/LinearRemapFV.cpp#L2996
> 
> Paul does a check for it in his code, unless I am misunderstanding what he does there. Unless he is checking for holes inside a face but looking at his code for `DoesMeshContain` holes it seems he does a similar thing as me. Also, @philipc2 the face area idea would work but seems to me that it would take much longer as you need to calculate the area for each face.

I understand that you are trying to follow what the tempest extreme does here. My point is: 
1. The tempest extreme is written in C/C++ while Uxarray is written in python. Also the API and data structure are different as well.
2. Can you write down the sudo codes for the entire re-mapping procedure so we can look into what's the ideal way to get it done in the Uxarray.

These two points are very important and introduce a recent important discussion for our developing process.
hongyuchen1030: And just some information from Computational Geometry, we normally assign the uncovered area a face index on the spherical surface as well. And when we said is neighboring with this face index, we know it is next by the uncovered area. Hope this point will help.
philipc2: > And just some information from Computational Geometry, we normally assign the uncovered area a face index on the spherical surface as well. And when we said is neighboring with this face index, we know it is next by the uncovered area. Hope this point will help.

This does make sense, however the meshes that we work with that have a partial spherical coverage (i.e. ocean meshes from climate models) don't assign an index to the uncovered area. 

![image](https://github.com/user-attachments/assets/c4e20d0f-a231-423e-ae9c-3c79f3520013)

> And when we said is neighboring with this face index, we know it is next by the uncovered area.

This is roughly what we are trying to achieve with the checks added by @aaronzedwick here. Since we don't have any explicit definitions of the geometry of the uncovered areas, the only way we are able to determine if a cell neighbors an empty region is by performing a search over the `edge_face_connectivity`. An index of `INT_FILL_VALUE` in this connectivity array indicates that a given edge is saddled between a valid face, and an empty region. 
aaronzedwick: > > And just some information from Computational Geometry, we normally assign the uncovered area a face index on the spherical surface as well. And when we said is neighboring with this face index, we know it is next by the uncovered area. Hope this point will help.
> 
> This does make sense, however the meshes that we work with that have a partial spherical coverage (i.e. ocean meshes from climate models) don't assign an index to the uncovered area.
> 
> ![image](https://private-user-images.githubusercontent.com/67855069/361475177-c4e20d0f-a231-423e-ae9c-3c79f3520013.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjQ2OTk1NDYsIm5iZiI6MTcyNDY5OTI0NiwicGF0aCI6Ii82Nzg1NTA2OS8zNjE0NzUxNzctYzRlMjBkMGYtYTIzMS00MjNlLWFlOWMtM2M3OWYzNTIwMDEzLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MjYlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODI2VDE5MDcyNlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTE1NWFlOTdmMGFhODI4NTNkMTU2MTY1OGNlZDIyMTBlYjViY2Y0MDc2NTZjMDVkMGU2MzFmM2FlMzQxYzJmZDgmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.yVXst8Kp_6CvK3FeuNZRRQ7zsPBLUXfQNJ5Y6-Cqs9Y)
> 
> > And when we said is neighboring with this face index, we know it is next by the uncovered area.
> 
> This is roughly what we are trying to achieve with the checks added by @aaronzedwick here. Since we don't have any explicit definitions of the geometry of the uncovered areas, the only way we are able to determine if a cell neighbors an empty region is by performing a search over the `edge_face_connectivity`. An index of `INT_FILL_VALUE` in this connectivity array indicates that a given edge is saddled between a valid face, and an empty region.

Oddly enough, even without the current checks for the data, the dual mesh construction does work with a partially covered data array. I am not sure exactly how the data array is constructed inside UXarray, but I was confused by this.

![image](https://github.com/user-attachments/assets/d9bfb4b4-ae69-4a14-9f53-2a7733b23b0f)

hongyuchen1030: > > > And just some information from Computational Geometry, we normally assign the uncovered area a face index on the spherical surface as well. And when we said is neighboring with this face index, we know it is next by the uncovered area. Hope this point will help.
> > 
> > 
> > This does make sense, however the meshes that we work with that have a partial spherical coverage (i.e. ocean meshes from climate models) don't assign an index to the uncovered area.
> > ![image](https://private-user-images.githubusercontent.com/67855069/361475177-c4e20d0f-a231-423e-ae9c-3c79f3520013.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjQ2OTk1NDYsIm5iZiI6MTcyNDY5OTI0NiwicGF0aCI6Ii82Nzg1NTA2OS8zNjE0NzUxNzctYzRlMjBkMGYtYTIzMS00MjNlLWFlOWMtM2M3OWYzNTIwMDEzLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MjYlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODI2VDE5MDcyNlomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTE1NWFlOTdmMGFhODI4NTNkMTU2MTY1OGNlZDIyMTBlYjViY2Y0MDc2NTZjMDVkMGU2MzFmM2FlMzQxYzJmZDgmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.yVXst8Kp_6CvK3FeuNZRRQ7zsPBLUXfQNJ5Y6-Cqs9Y)
> > > And when we said is neighboring with this face index, we know it is next by the uncovered area.
> > 
> > 
> > This is roughly what we are trying to achieve with the checks added by @aaronzedwick here. Since we don't have any explicit definitions of the geometry of the uncovered areas, the only way we are able to determine if a cell neighbors an empty region is by performing a search over the `edge_face_connectivity`. An index of `INT_FILL_VALUE` in this connectivity array indicates that a given edge is saddled between a valid face, and an empty region.
> 
> Oddly enough, even without the current checks for the data, the dual mesh construction does work with a partially covered data array. I am not sure exactly how the data array is constructed inside UXarray, but I was confused by this.
> 
> ![image](https://private-user-images.githubusercontent.com/95507181/361547281-d9bfb4b4-ae69-4a14-9f53-2a7733b23b0f.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjQ3MDAxMDcsIm5iZiI6MTcyNDY5OTgwNywicGF0aCI6Ii85NTUwNzE4MS8zNjE1NDcyODEtZDliZmI0YjQtYWU2OS00YTE0LTlmNTMtMmE3NzMzYjIzYjBmLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MjYlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODI2VDE5MTY0N1omWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTZiNWYzZDIxNGIzMmVmMTU1OTY0YjZmNDc4N2Y1NmFkMzIzMjUyZjJkYmEyYTg1ZWJiNzlkNjg3MTc0Mzk1MTgmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.x1XixuHPI-52vjNbJChiA_s3q_7Rm63ckOYbRj-HrWQ)

Is it possible that it uses the `FILL_VALUE` as the face center value? Have you looked into the output value?
aaronzedwick: @hongyuchen1030 It doesn't seem to, no, the values in the data array are all actually data points without any fill values, and the number of values matches the number of faces.
philipc2: > @hongyuchen1030 It doesn't seem to, no, the values in the data array are all actual data points without any fill values, and the number of values matches the number of faces.

Seconding this. `INT_FILL_VALUE` is `int64`, while the type of the data variable isn't constant (however it is typically a `float`)


aaronzedwick: > May we include a basic benchmark in the `mpas_ocean` file for computing `Grid.hole_edge_indices`

Yeah, that's a good idea!

Patch:
diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py
index d2baaf86..51b0c99d 100644
--- a/benchmarks/mpas_ocean.py
+++ b/benchmarks/mpas_ocean.py
@@ -164,3 +164,17 @@ class RemapUpsample:
 
     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 b263f4ec..6dd40b2c 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 ec3278b6..da904fe0 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 2cd0eb2a..1d020b91 100644
--- a/uxarray/conventions/descriptors.py
+++ b/uxarray/conventions/descriptors.py
@@ -25,3 +25,9 @@ EDGE_NODE_DISTANCES_ATTRS = {
     "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 08229d1a..afa31909 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 97b9f5d1..b1ff3e43 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -54,6 +54,7 @@ from uxarray.grid.geometry import (
     _grid_to_matplotlib_polycollection,
     _grid_to_matplotlib_linecollection,
     _populate_bounds,
+    _construct_hole_edge_indices,
 )
 
 from uxarray.grid.neighbors import (
@@ -1085,6 +1086,16 @@ class Grid:
             _ = 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 1a485271..4f74249a 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 57eb3f3b..6b80aa87 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -31,6 +31,7 @@ gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne3
 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 @@ class TestGrid(TestCase):
         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 @@ class TestConnectivity(TestCase):
 
         # 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 @@ class TestClassMethods(TestCase):
 
 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 @@ class TestLatlonBounds(TestCase):
         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