Check Grid for Partial Spherical Coverage

[aaronzedwick]

Closes #898

Overview

Allows for a grid to be checked for holes.

Expected Usage

import uxarray as ux

grid_path = "/path/to/grid.nc"

uxds = ux.open_grid(grid_path)

uxds.contains_holes()

PR Checklist

General

• An issue is linked created and linked
• Add appropriate labels
• Filled out Overview and Expected Usage (if applicable) sections

Testing

• Adequate tests are created if there is new functionality
• Tests cover all possible logical paths in your function
• Tests are not too basic (such as simply calling a function and nothing else)

Documentation

• Docstrings have been added to all new functions
• Docstrings have updated with any function changes
• Internal functions have a preceding underscore (_) and have been added to docs/internal_api/index.rst
• User functions have been added to docs/user_api/index.rst

[rajeeja]

looks good adds closes #141

[hongyuchen1030]

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]

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]

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]

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!

[hongyuchen1030]

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]

@aaronzedwick

If all we are interested in knowing is whether there are empty patches with no geometry, such as continents, we should be able to check the total face area of the unit sphere.

If the area is less than 4pi, it means that not the entire surface is covered with geometries.

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

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

[philipc2]

I don't see an issue with this check, as the reasoning that @aaronzedwick explained is sound. If there exists any edge that does not have two neighboring faces, then there must exist some "empty" area adjacent to that edge. Illustration below for reference: the edges marked in red each neighbor an empty region, which means that we can't construct an edge between the centers of the two faces indicated in the edge_face_connectivity

This causes issues with the Dual Mesh construction, since node-centered data variables become face-centered. However, there will be no faces constructed centered around these values. Given our interpretation of the UGRID conventions, we have no way to construct a face for this data, meaning that the data is essentially "lost" from the dual mesh construction process.

Also, I will suggest caching the result once we've run it, so that we don't need to run the repeated searches over the edge_face_connectivity. Something as simple as a Grid._partial_coverage attribute would suffice.

[hongyuchen1030]

I don't see an issue with this check, as the reasoning that @aaronzedwick explained is sound. If there exists any edge that does not have two neighboring faces, then there must exist some "empty" area adjacent to that edge. Illustration below for reference: the edges marked in red each neighbor an empty region, which means that we can't construct an edge between the centers of the two faces indicated in the edge_face_connectivity

We have already gone through this part in the previous discussion above. I understand what you mean by "missing a face" here. My current concern is: after we have examined the entire grid to check if there's a partial coverage, what's the next step to deal with the missing faces? Do we need to rerun a similar check to identify which faces don't have neighbors? So basically two O(N) operations?

Also, I will suggest caching the result once we've run it, so that we don't need to run the repeated searches over the edge_face_connectivity. Something as simple as a Grid._partial_coverage attribute would suffice.

This is a good point. As we agreed in the long conversation from the zonal_mean PR, we need to provide the sudo codes for the algorithm. So that we can better adapt it to the UXarray API and the python vectorization styles.

[philipc2]

My current concern is: after we have examined the entire grid to check if there's a partial coverage, what's the next step to deal with the missing faces? Do we need to rerun a similar check to identify which faces don't have neighbors? So basically two O(N) operations?

This is a good point that I overlooked. We should be keeping track of which indices are considered to border an empty area, which can be achieved by removing the np.any call

edge_face_connectivity[:, 1] == INT_FILL_VALUE

Once we have this, we need to identify which nodes (indicated in black) will need to be removed from the new grid that is constructure from the Dual Grid.

The number of faces in the new Grid should equal the original number of nodes (n_node) minus the number of nodes that share one or more of the "invalid" edges shown above.

[hongyuchen1030]

edge_face_connectivity[:, 1] == INT_FILL_VALUE

This sounds good. If we need to check if there's a partial coverage, don't just return a boolean value, but instead directly return where is the partial coverage. And if the return is None, then we know is all covered.

The number of faces in the new Grid should equal the original number of nodes (n_node) minus the number of nodes that share one or more of the "invalid" edges shown above.

It sounds good to me though I haven't thought deep into this yet. But I will keep you updated if anything other things comes into my head.

[aaronzedwick]

Sounds perfect to me, @philipc2 @hongyuchen1030 thanks for the suggestions on this one! Let me know if you have anymore thoughts and I will get the changes mentioned here implemented!

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

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.

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.

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

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.

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

[philipc2]

May we include a basic benchmark in the mpas_ocean file for computing Grid.hole_edge_indices

[aaronzedwick]

May we include a basic benchmark in the mpas_ocean file for computing Grid.hole_edge_indices

Yeah, that's a good idea!

[github-actions]

ASV Benchmarking

Benchmark Comparison Results

Benchmarks that have improved:

Change	Before [f71e11b]	After [641d1f0]	Ratio	Benchmark (Parameter)
	failed	214±8μs	n/a	mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
	failed	108±2μs	n/a	mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
-	406M	354M	0.87	mpas_ocean.Integrate.peakmem_integrate('480km')

Benchmarks that have stayed the same:

Change	Before [f71e11b]	After [641d1f0]	Ratio	Benchmark (Parameter)
	379M	380M	1	face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
	380M	381M	1	face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
	383M	384M	1	face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
	381M	382M	1	face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
	1.33±0s	1.35±0s	1.02	face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
	175±2ms	178±1ms	1.02	face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
	1.60±0.01s	1.60±0s	1	face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
	7.80±0.07ms	7.77±0.1ms	1	face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
	1.69±0.01s	1.70±0.01s	1.01	import.Imports.timeraw_import_uxarray
	629±5ms	615±4ms	0.98	mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
	40.3±0.5ms	39.0±0.2ms	0.97	mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
	1.74±0.07ms	1.65±0.3ms	0.95	mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
	484±20μs	482±10μs	1	mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
	296±2ns	302±7ns	1.02	mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
	819±10ns	819±5ns	1	mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
	279±2ns	285±3ns	1.02	mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
	399M	399M	1	mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('120km', False)
	387M	387M	1	mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('120km', True)
	362M	362M	1	mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('480km', False)
	361M	361M	1	mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('480km', True)
	1.03±0.01s	1.02±0.01s	0.99	mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
	55.9±0.6ms	57.1±1ms	1.02	mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
	75.8±0.6ms	76.4±0.6ms	1.01	mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
	5.12±0.1ms	5.17±0.2ms	1.01	mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
	263M	263M	1	mpas_ocean.Gradient.peakmem_gradient('120km')
	243M	240M	0.99	mpas_ocean.Gradient.peakmem_gradient('480km')
	2.66±0.01ms	2.67±0.01ms	1	mpas_ocean.Gradient.time_gradient('120km')
	291±3μs	282±0.7μs	0.97	mpas_ocean.Gradient.time_gradient('480km')
	370M	370M	1	mpas_ocean.Integrate.peakmem_integrate('120km')
	175±1ms	175±1ms	1	mpas_ocean.Integrate.time_integrate('120km')
	11.8±0.09ms	11.6±0.06ms	0.99	mpas_ocean.Integrate.time_integrate('480km')
	340±0.7ms	344±4ms	1.01	mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
	342±2ms	348±2ms	1.02	mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
	343±0.8ms	344±3ms	1.01	mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
	21.6±0.6ms	21.7±0.3ms	1	mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
	21.2±0.3ms	21.8±0.3ms	1.03	mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
	21.8±0.2ms	21.8±0.2ms	1	mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
	54.8±0.07ms	55.0±0.3ms	1	mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
	44.5±0.2ms	44.2±0.08ms	0.99	mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
	360±1ms	362±0.6ms	1	mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
	262±0.6ms	263±0.6ms	1	mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
	236M	236M	1	quad_hexagon.QuadHexagon.peakmem_open_dataset
	236M	236M	1	quad_hexagon.QuadHexagon.peakmem_open_grid
	6.41±0.02ms	6.45±0.02ms	1.01	quad_hexagon.QuadHexagon.time_open_dataset
	5.53±0.05ms	5.58±0.1ms	1.01	quad_hexagon.QuadHexagon.time_open_grid

Benchmarks that have got worse:

Change	Before [f71e11b]	After [641d1f0]	Ratio	Benchmark (Parameter)
+	1.19±0μs	1.34±0μs	1.12	mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')

[philipc2]

The implementation here looks great and the execution time is quick! Thanks for putting this together.