Merge duplicate nodes when constructing the dual mesh - #1692
Conversation
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.
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.
There was a problem hiding this comment.
Fixed — matching goes through _coincident_node_canonical_indices, which is a KDTree query within ERROR_TOLERANCE on the sphere, not exact equality.
| 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.)
There was a problem hiding this comment.
The dict only holds duplicates, not every node, so it stays small. The real cost was _check_duplicate_nodes_indices looping over every face in Python; that is now one np.isin over the connectivity array.
| """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)
There was a problem hiding this comment.
Fixed — _remap_node_connectivity builds a lookup array and indexes into it once, no per-duplicate array scan.
| 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.
There was a problem hiding this comment.
Added assert grid.n_node == 6000 and len(duplicates) == 2150, and derived 3840 in the test: 3850 distinct nodes remain after the merge, ten touched by one face only, so no dual cell. Also added test_duplicate_nodes_minimal_example — two quads, eight nodes at six locations, asserting the map is {6: 1, 7: 2}.
|
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?) |
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.
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.
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
Suggestion: please change the title of this PR to reflect the full scope. Duplicate nodes are now handled directly whenever constructing a Grid, not just when constructing the dual mesh. I was confused about why _check_duplicate_nodes_indices checks were removed despite construct_dual() being unchanged. I think the reason is because all Grid objects are now guaranteed to not have duplicate nodes.
I tried leaving a full review but I kept getting the feeling that something weird was happening, a nagging feeling like "hey I think I've looked at this code before and left an inline comment, why am I reviewing it again?" Then I realized many of the changes here are also in #1690 which hasn't been merged yet. That led to duplicating review work and will probably lead to needing to apply fixes multiple times. I'm not sure the cleanest way forwards at this point, but I might suggest waiting to merge this until after #1690 gets merged. At least, I will want to do another close review of this PR after that PR merges, because I lost track of which things I already looked at closely and which I need to consider again.
| def test_global_structured_grid_merges_poles_and_seam(): | ||
| """Nodes coincident on the sphere must be merged, even though their | ||
| (lon, lat) pairs differ. Regression test for issue #1689.""" | ||
| import numpy as np |
There was a problem hiding this comment.
Please move numpy imports to top of file; numpy imports should always be at top of file throughout all parts of uxarray.
| return canonical | ||
|
|
||
| tree = KDTree(points_xyz[mergeable_indices]) | ||
| pairs = tree.query_pairs(r=tolerance, output_type="ndarray") |
There was a problem hiding this comment.
Why does this tolerance use ERROR_TOLERANCE for a radius, but other tolerances convert to chord_tol?
There was a problem hiding this comment.
ERROR_TOLERANCE is already a chord distance on the unit sphere, so it is the radius directly. Conversion is only needed where the tolerance arrives in degrees, as in _read_structured_grid.
|
|
||
| # ``tol`` is an angle in degrees; on the unit sphere the matching radius is the | ||
| # chord subtended by that angle, so the threshold keeps its documented meaning. | ||
| chord_tol = 2.0 * np.sin(np.deg2rad(tol) / 2.0) |
There was a problem hiding this comment.
This computation is redundant with the _DEFAULT_STRUCTURED_TOL_DEG computation above; the result will be incorrect because the formula is being applied twice.
…des' into rajeeja/coincident-nodes # Conflicts: # test/io/test_structured.py # uxarray/core/dataarray.py # uxarray/core/dataset.py # uxarray/grid/grid.py # uxarray/io/_structured.py
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.