Merge Duplicate Node Indices

[philipc2]

Closes #865

Overview

• Adds the merge_duplicate_node_indices() function that replaces the index of nodes that have the exact latitude and longitude to only reference one node.

Expected Usage

import uxarray as ux

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

uxgrid = ux.open_grid(grid_path)

# inplace
uxgrid.merge_duplicate_node_indices()

# return a new grid
new_grid = uxgrid.merge_duplicate_node_indices(inplace=False)

[aaronzedwick]

Do you know what the problem was with the previous implementation? If I remember it was still having problems?

[philipc2]

Do you know what the problem was with the previous implementation? If I remember it was still having problems?

I haven't been able to pinpoint the problem, but everything appears to be working correctly here according to the tests.

[philipc2]

@aaronzedwick

I'm going to merge these changes into the branch used in #859 to see the compatibility.

[philipc2]

@aaronzedwick

May you take a look at the case using the geos dataset? Looks like something is breaking in the construct_faces function call.

[rajeeja]

what tests are failing? we may want to issue a warning on how and what is supported.
Supporting all formats may be harder and might break with corner cases.

[philipc2]

what tests are failing? we may want to issue a warning on how and what is supported. Supporting all formats may be harder and might break with corner cases.

Something with the Dual Mesh Construction on the Cube Sphere grid is breaking. It's complaining about a dimension miss match when assigned the faces.

I don't believe it's an issue with how we are merging ndoes.

[aaronzedwick]

@aaronzedwick

May you take a look at the case using the geos dataset? Looks like something is breaking in the construct_faces function call.

Oh yeah, I know the problem. Let me fix it for you.

[aaronzedwick]

@philipc2 There you go!

[philipc2]

@philipc2 There you go!

Thanks for fixing this so quickly!

[philipc2]

@aaronzedwick

There still appear to be some underlying issues, at least on the Cube Sphere grid.

Running the following executes, however when I try to plot it I get the following error:

uxgrid.merge_duplicate_node_indices(inplace=True)
dual = uxgrid.compute_dual()

---

IndexError Traceback (most recent call last)
Cell In[18], line 1
----> 1 dual.plot()

File ~/anaconda3/envs/dev-uxarray/lib/python3.11/site-packages/uxarray/plot/accessor.py:30, in GridPlotAccessor.call(self, **kwargs)
29 def call(self, **kwargs) -> Any:
---> 30 return grid_plot.plot(self._uxgrid, **kwargs)

File ~/anaconda3/envs/dev-uxarray/lib/python3.11/site-packages/uxarray/plot/grid_plot.py:16, in plot(grid, **kwargs)
14 def plot(grid: Grid, **kwargs):
15 """Default Plotting Method for Grid."""
---> 16 return mesh(grid, **kwargs)

File ~/anaconda3/envs/dev-uxarray/lib/python3.11/site-packages/uxarray/plot/grid_plot.py:43, in mesh(uxgrid, backend, exclude_antimeridian, width, height, xlabel, ylabel, **kwargs)
19 def mesh(
20 uxgrid: Grid,
21 backend: Optional[str] = "bokeh",
(...)
27 **kwargs,
28 ):
29 """Vector Line Plot of the edges that make up each face.
30
31 Parameters
(...)
40 Plot Width for Bokeh Backend
41 """
---> 43 gdf = uxgrid.to_geodataframe(exclude_antimeridian=exclude_antimeridian)
45 hv_paths = hv.Path(gdf)
47 uxarray.plot.utils.backend.assign(backend=backend)

File ~/anaconda3/envs/dev-uxarray/lib/python3.11/site-packages/uxarray/grid/grid.py:1304, in Grid.to_geodataframe(self, override, cache, exclude_antimeridian)
1301 return self._gdf
1303 # construct a geodataframe with the faces stored as polygons as the geometry
-> 1304 gdf = _grid_to_polygon_geodataframe(
1305 self, exclude_antimeridian=exclude_antimeridian
1306 )
1308 # cache computed geodataframe
1309 if cache:

File ~/anaconda3/envs/dev-uxarray/lib/python3.11/site-packages/uxarray/grid/geometry.py:95, in _grid_to_polygon_geodataframe(grid, exclude_antimeridian)
91 def _grid_to_polygon_geodataframe(grid, exclude_antimeridian):
92 """Converts the faces of a Grid into a spatialpandas.GeoDataFrame
93 with a geometry column of polygons."""
---> 95 polygon_shells = _build_polygon_shells(
96 grid.node_lon.values,
97 grid.node_lat.values,
98 grid.face_node_connectivity.values,
99 grid.n_face,
100 grid.n_max_face_nodes,
101 grid.n_nodes_per_face.values,
102 )
104 antimeridian_face_indices = _build_antimeridian_face_indices(
105 polygon_shells[:, :, 0]
106 )
108 if grid.n_face > GDF_POLYGON_THRESHOLD:

File ~/anaconda3/envs/dev-uxarray/lib/python3.11/site-packages/uxarray/grid/geometry.py:82, in _build_polygon_shells(node_lon, node_lat, face_node_connectivity, n_face, n_max_face_nodes, n_nodes_per_face, projection)
77 node_lon = lonlat_proj[:, 0]
78 node_lat = lonlat_proj[:, 1]
80 polygon_shells = (
81 np.array(
---> 82 [node_lon[closed_face_nodes], node_lat[closed_face_nodes]], dtype=np.float32
83 )
84 .swapaxes(0, 1)
85 .swapaxes(1, 2)
86 )
88 return polygon_shells

IndexError: index -9223372036854775808 is out of bounds for axis 0 with size 864

Some of the entries in the face_node_connectivity are now incorrect. There are some entries that are an entire row of fill values.

The number of entirely zero entries corresponds to how many duplicates were found.

Instead of having rows of full fill values, we should ensure that each duplicate has the same entry in the connectivity.

[aaronzedwick]

@philipc2 I am not exactly sure what's happening, but I don't believe its a problem with my code. The only reason you are getting that error is because the connectivity that you are passing in is broken I think. Comment out the line in construct_faces that says if n_edges[i] < 3: continue and then compute dual will run, and the plot is still incorrect. All those fill value entries look to only have one node for some reason around which it is trying to make a dual face. Something with the node_face_connectivity I suppose?

[philipc2]

@philipc2 I am not exactly sure what's happening, but I don't believe its a problem with my code. The only reason you are getting that error is because the connectivity that you are passing in is broken I think. Comment out the line in construct_faces that says if n_edges[i] < 3: continue and then compute dual will run, and the plot is still incorrect. All those fill value entries look to only have one node for some reason around which it is trying to make a dual face. Something with the node_face_connectivity I suppose?

Yeah that seems to be the case, I'll look into it!

[philipc2]

@philipc2 I am not exactly sure what's happening, but I don't believe its a problem with my code. The only reason you are getting that error is because the connectivity that you are passing in is broken I think. Comment out the line in construct_faces that says if n_edges[i] < 3: continue and then compute dual will run, and the plot is still incorrect. All those fill value entries look to only have one node for some reason around which it is trying to make a dual face. Something with the node_face_connectivity I suppose?

Something still appears to be off, I'm getting the same error now (even though the connectivity appears to have been corrected)

[aaronzedwick]

Can you push the changes that fixed the connectivity?

[rajeeja]

just one minor test and this is good go.

[rajeeja]

add an assert to check that compute_dual fails or produces incorrect stuff without this operation.