Merge duplicate nodes when constructing the dual mesh - #1692
Conversation
Canonicalize duplicate node indices in the face-node connectivity before building the dual, so grids with repeated nodes produce a correct dual instead of being rejected. Also vectorize the duplicate lookup and deprecate the now redundant check_duplicate_nodes argument.
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
Mostly looks like a good fix! I think it still needs a little bit of extra work, and I left some inline comments accordingly.
| """Map duplicate node indices to the first index with the same coordinates.""" | ||
| node_coordinates = np.column_stack((grid.node_lon.values, grid.node_lat.values)) | ||
| _, first_indices, inverse_indices = np.unique( | ||
| node_coordinates, axis=0, return_index=True, return_inverse=True |
There was a problem hiding this comment.
Is "exact equality" the correct way to go here? My intuition originally was that there should probably be some sort of tolerance here, e.g. if values agree to within 1e-12, they are probably the same node, right?
There was a problem hiding this comment.
In fact, looking back at the original issue thread, it looks like you were the person who originally suggested there be a tolerance in the first place! So, I would actually now assert that exact inequality is not the desired implementation, there should be a tolerance, as clarified in thread of #865.
| np.arange(grid.n_node, dtype=INT_DTYPE) != first_indices[inverse_indices] | ||
| ) | ||
| return { | ||
| INT_DTYPE(index): INT_DTYPE(first_indices[inverse_indices[index]]) |
There was a problem hiding this comment.
I would guess that creating a dict here is extremely inefficient… maybe that doesn't matter if there are only ever a tiny number of duplicate nodes. Not necessarily blocking, but have you looked into how long this takes to run for any larger grids containing duplicate nodes? (Or, do you expect a very limited number of duplicate nodes in most cases? I would guess it is probably not worth worrying about if there are less than ~1000 duplicates or so.)
| """Return a copy of connectivity with duplicate node indices canonicalized.""" | ||
| remapped_connectivity = connectivity.copy() | ||
| for duplicate_index, source_index in duplicate_node_indices.items(): | ||
| remapped_connectivity[remapped_connectivity == duplicate_index] = source_index |
There was a problem hiding this comment.
(There is also probably some cleverer way to do this using numpy indexing without a loop through a dictionary… but maybe not necessary, see my comment related to efficiency in validation.py)
| def test_dual_duplicate(gridpath): | ||
| """Test dual mesh creation with duplicate grids.""" | ||
| dataset = ux.open_dataset(gridpath("ugrid", "geoflow-small", "grid.nc"), gridpath("ugrid", "geoflow-small", "grid.nc")) | ||
| """Test dual mesh creation with duplicate node indices.""" |
There was a problem hiding this comment.
Could you include something in this test to assert there are actually duplicate nodes in the original grid? That would help to prove this test is actually testing what it claims to be testing. Right now I just have to trust that geoflow-small grid happens to contain duplicates, but I can't see from these lines if that's really true, or how many duplicates there are.
Can you also clarify with a comment where the number 3803 comes from?
Extra helpful, but not necessarily required, would be if you are able to construct a tiny example inline here which clearly has some duplicate nodes, something small enough to directly reason through how they should be handled.
| @@ -129,7 +129,14 @@ def test_dual_mesh_mpas(gridpath): | |||
|
|
|||
|
|
|||
| def test_dual_duplicate(gridpath): | |||
There was a problem hiding this comment.
Please add tests which also check the other grids mentioned in the original issue report: geos-cs/c1/test-c12.native.nc4 and esmf/ne30/ne30pg3.grid.nc. It would be fine to add them as part of this function, unless you prefer to make new test functions for them.
|
Actually, apologies for not including this during the original review, comment but one more thought: does this actually fully close the original issue? The issue writeup makes it sound to me like duplicate nodes should be handled immediately upon constructing the grid, not just during one functionality (get_dual). Is there a reason that duplicate nodes should be handled only during get_dual, instead of immediately? (Do all other current/planned functions work properly regardless of whether there are duplicate nodes?) |
Match nodes in Cartesian space rather than the lon/lat plane so pole and antimeridian nodes are recognized as the same point, and store the resulting polar faces as triangles instead of quads with a repeated corner.
float32 input (e.g. real climate datasets) silently ran the whole xyz/tolerance pipeline at float32 precision, causing pole/antimeridian merges to fail or merge only partially.
Extend #865's fix beyond the dual mesh: canonicalize duplicate/coincident node indices in connectivity for every Grid construction path, not just construct_dual. Detection is now tolerance-based (unit-sphere chordal distance) instead of exact lon/lat match, so pole-degenerate duplicates are also caught. Node coordinate arrays are left untouched by design; only connectivity is remapped to canonical indices, with any resulting repeated face corners collapsed.
construct_dual no longer needs its own per-call duplicate detection and remap, and get_dual() no longer needs to hard-gate on duplicate node indices, since Grid construction now canonicalizes them structurally before any of this code runs.
Since duplicate node coordinates are intentionally left unreferenced by connectivity, a node KDTree/BallTree built over the raw coordinate array could select an index no face actually points to, silently returning empty or wrong nearest-neighbor results. Build the "nodes" tree only over live (referenced) indices and translate query results back to original index space.
polars' unique() with maintain_order unset does not guarantee row order across runs, so the node index assigned to a given corner coordinate could vary between reads of the same file. This is normally harmless, but it made canonical-node selection for coincident duplicates (e.g. pole points with differing longitude) flaky from run to run.
test_dual_duplicate: validate() now succeeds since connectivity is fully canonicalized (duplicate coordinates remain by design, but nothing references a dead index anymore). test_grid_nn_subset: max valid k for a node search is now bounded by the live node count, not raw node count. test_to_geodataframe_preserves_antimeridian_faces: pole-coincident corners with differing longitude are now also merged, shifting the antimeridian face count.
Merging pole-adjacent duplicate nodes was collapsing each face's own locally-meaningful longitude at the pole into one arbitrary canonical value, which corrupted lat/lon bounds and broke zonal weight computation for cube-sphere grids near the poles.
2 participants
Closes #865
geoflow-smallpreviously could not produce a dual at all and now yields 3803 faces; the test asserts this and fails on main._find_duplicate_nodesis vectorized withnp.uniqueinstead of a per-node dict.Grid.get_dual(check_duplicate_nodes=...)is now ignored and deprecated rather than removed, so existing callers keep working.UxDataArray.get_dualandUxDataset.get_dualstill raiseGridInvalidError, since node-centered data cannot be remapped onto a merged node set.(lon, lat), such as poles and the antimeridian; that isGrid.from_structureddoes not merge coincident pole and antimeridian nodes #1689 / Merge coincident pole and antimeridian nodes in structured grids #1690.