Merge coincident nodes when constructing any Grid - #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.
There was a problem hiding this comment.
This fix looks reasonable. One more nitpick here though: your comment on the original thread #865 (comment) suggested you wanted the user to be able to specify the tolerance.
Right now, there is no way to specify the "tolerance for merging duplicate nodes" from any public-facing functions. To specify it as a larger value (e.g. 1e-4) a user would need to call _merge_coincident_grid_ds_nodes directly on grid_ds before passing to Grid.__init__, while there is presumably no way to properly specify it as a smaller value (e.g. 1e-12) because the _merge_coincident_grid_ds_nodes call inside Grid.__init__ would always rerun the process with a tolerance of ERROR_TOLERANCE=1e-8.
Should this option be added somewhere (e.g. something like a duplicates_tolerance parameter in Grid.__init__) or is "allow the user to specify the tolerance" no longer desirable?
There was a problem hiding this comment.
Can be a good edition, but not in this PR. I'll open a separate issue for a user-facing tolerance.
Adding this means threading it through open_grid and the from_* constructors..
| 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.
There was a problem hiding this comment.
I think that makes sense. Added the run-benchmarks label because this PR adds a call to _merge_coincident_grid_ds_nodes into Grid.__init__ so it might affect uxarray performance. I believe my original comment here is resolved, but I am keeping this thread open for now as a reminder to myself to consider the benchmarking results before approving this PR.
There was a problem hiding this comment.
The cache wasn't the whole story - dual mesh is still ~1.6x slower because get_dual builds a new Grid, which re-runs the coincident merge on the dual. Fixing that next, so don't go by the current benchmark numbers.
| 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}.
There was a problem hiding this comment.
The new version of the test looks much clearer and the comments do a good job explaining the steps overall.
I still don't think the number 3850 is explained clearly. In particular, why should we expect len(duplicates)==2150 for this case? I am extra suspicious because the implementation changes have changed the number from 3803 to 3850.
One way forward would be to somehow separately figure out the number you would expect here. Although, I do think that the test_duplicate_nodes_minimal_example might be a reasonable enough test for ensuring the code is identifying duplicates correctly. If you want to avoid finding a justification for the len(duplicates)==2150 number, it would be sufficient to leave a comment there to say something like "The source file just happened to have 2150 duplicates according to uxarray, but there is not necessarily any reason to suspect it must be 2150. This test does not check for correctness of which duplicates are identified, the number is hard-coded just to prevent unexpected regressions. For a correctness test of the duplicate-identifying algorithm, see test_duplicate_nodes_minimal_example."
There was a problem hiding this comment.
2150 is 2111 nodes exactly equal in lon/lat plus 39 within ERROR_TOLERANCE of an earlier node, and 10 pole nodes are held out of merging by design. The 3803 you're comparing to was the old dual.n_face assertion, not a node count; the comparable number now is 3840, and 3850 is distinct nodes.
|
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.
| 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.
There was a problem hiding this comment.
Is there any reason to avoid using the same meaning for tolerance here as in _read_structured_grid? I would be less confused by it if it meant the same thing in both places. It seems like both routines are using the input tolerance in a similar way, while building a KDTree to query for everything within a certain distance of a given point.
There was a problem hiding this comment.
Only the input unit differs: _read_structured_grid takes tol in degrees from the user so it converts, while ERROR_TOLERANCE is already a chord on the unit sphere. Both query the KDTree with the same chord radius by default.
…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
"Coincident nodes" is the domain term used by TempestRemap, MOAB, and our own validation module; "dedupe" reads like a dataframe operation. Also note in the docstring that we keep the redundant coordinates and only remap connectivity, where TempestRemap and MOAB delete and renumber.
_coincident_node_canonical_indices exempts pole nodes from merging so each face touching a pole keeps its own longitude. The mask used np.isclose(|z|, 1.0, atol=tolerance), which is wrong twice: tolerance is a chord radius everywhere else in the function, and np.isclose's default rtol=1e-5 swamped atol entirely. The carve-out spanned 1 - |z| <= 1.001e-5, a chord of 4.5e-3 or ~28 km on Earth, so any node within 28 km of a pole was silently exempt from merging. For a point at colatitude t, 1 - |z| = 1 - cos(t) = chord**2 / 2, so the correct cap is tolerance**2 / 2 with rtol pinned to zero.
construct_dual reads node_face_connectivity with no duplicate handling, so a face still referencing a coincident duplicate index yields a degenerate dual face instead of an error. Merging at construction makes that unreachable today, so the guard is a no-op safety net; check_duplicate_nodes is now deprecated and ignored, since the check always runs.
The pole/seam merge in _read_structured_grid changes the ersstv5 grid (16290 nodes to 16200) and turns pole quads into triangles, so the near-pole plots were stale. Outputs re-executed; the stray OMP stderr line is dropped.
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
The notebook's cell sources are identical to main; only the outputs differed, and those differences were an xarray coordinate-ordering repr change and a stray OMP info line. Keeping the re-run added ~1700 lines of diff with no content change.
Coincident-node merging now runs on every Grid construction, so the cost is paid by grids that have no duplicates at all. Two points within a chord of ERROR_TOLERANCE differ by at most that in x, so in x-sorted order every consecutive gap between them is also within tolerance; a point whose sorted neighbours are both further away cannot be coincident with anything. Filtering on that is exact and drops the tree build entirely for a clean grid: 1.23s -> 0.28s over 2M nodes.
rajeeja
changed the title
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have stayed the same:
Benchmarks that have got worse:
|
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
Making good progress! Still a few inline comments to resolve, see below and also I left some follow-up comments on previous inline comment threads.
Also, looking at ASV benchmark results, it looks like quite a few benchmarks have slowed down significantly, especially dual mesh construction being ~50% slower. Why would dual mesh construction slow down so much, if the search for duplicates primarily occurs when building the original grid? Tagging @cmdupuis3 to maybe comment on the benchmarks if you have additional insights. See also #1710 as a sanity check that the benchmark changes represent real changes to performance, rather than quirks of the benchmarking suite.
|
|
||
| 3---2---7 lat 1 nodes 2,3 are the shared edge | ||
| | | | nodes 7,6 are their duplicates | ||
| 0---1---6 lat 0 |
There was a problem hiding this comment.
I don't think this diagram corresponds to the words/code… did a robot draw this? Please fix! (Or, please clarify how it corresponds to the data, if I misunderstood….) Where are nodes 4 and 5?
There was a problem hiding this comment.
Also, the repeated nodes in the code below are 1,2 mapping to 6,7, not 2,3 mapping to 6,7 as claimed in the test comments above.
There was a problem hiding this comment.
Redrew the diagram and fixed the node numbering in 77b5f83.
| 0---1---6 lat 0 | ||
| """ | ||
| node_lon = np.array([0.0, 1.0, 1.0, 0.0, 2.0, 2.0, 1.0, 1.0]) | ||
| node_lat = np.array([0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0]) |
There was a problem hiding this comment.
The test suite is missing a (non-pole) test demonstrating that error tolerance is being handled properly. Adding a tiny number (1e-9?) to one of the duplicates here would be sufficient to spot check this. E.g. node 1 could keep node_lon = 1.0, but set node 6's node_lon to 1.0+1e-9.
There was a problem hiding this comment.
Added test_duplicate_nodes_tolerance in c968d01, node 6 at lon 1.0+1e-9.
There was a problem hiding this comment.
Correction: that offset was 1e-9 degrees but the cutoff is a 1e-8 chord (~5.7e-7 degrees), so the test passed far inside the boundary. Parametrized in 4d580fc over 1e-9 and 1e-7 (merge) and 1e-5 (must not).
| ) | ||
| ) | ||
| else: | ||
| return grid_ds |
There was a problem hiding this comment.
Why return grid_ds unchanged if it doesn't have x,y,z or lon,lat values? My concern is that this is a good way to silently cause incorrect results, where you expect duplicates to have been merged but can't trust that they actually were merged. Maybe crash or at least raise a warning in this case?
There was a problem hiding this comment.
Warns with a RuntimeWarning now instead of returning silently, in f1e6e5c. Kept the unchanged return so datasets with no coordinates still round-trip.
There was a problem hiding this comment.
Added test_merge_warns_when_node_locations_are_unknown in 4d580fc; the branch had no coverage.
| ) | ||
| ) | ||
| grid_ds["face_node_connectivity"] = grid_ds["face_node_connectivity"].copy( | ||
| data=_collapse_repeated_face_corners( |
There was a problem hiding this comment.
Why is something like _collapse_repeated_face_corners only relevant in the face_node_connectivity case? Maybe it is true but could you clarify with comment reasoning through the other connectivity combinations and justifying why it wouldn't happen for them?
There was a problem hiding this comment.
Added the reasoning as a comment in e61e1fd.
There was a problem hiding this comment.
Correction: node_node_connectivity can gain a repeat too ([1, 1, 3, FILL] when a row lists both a duplicate and its canonical node). Rewrote the comment in 4d580fc to cover all three groups.
Thanks, the dual mesh slowdown was the duplicate check running on every get_dual call and rescanning every node each time. It's cached on the grid now. |
get_dual constructs a new Grid from this grid's face centers, which re-ran the coincident node search on every call and made dual mesh construction ~1.6x slower. The dual's connectivity is already canonical, so the merge is skipped there via a new internal Grid flag.
…ason through corner collapsing
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.