Merge duplicate nodes when constructing the dual mesh

test/grid/grid/test_core.py

@@ -5,6 +5,11 @@
 
 import uxarray as ux
 from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE
+from uxarray.grid.validation import (
+    _check_duplicate_nodes_indices,
+    _find_duplicate_nodes,
+)
+from uxarray.errors import GridInvalidError
 
 
 def test_grid_with_holes(gridpath):
@@ -129,7 +134,142 @@ def test_dual_mesh_mpas(gridpath):
 
 
 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"))
-    with pytest.raises(ux.errors.GridInvalidError):
-        dataset.get_dual()
+    """Test dual mesh creation on a grid whose source file has duplicate
+    (coincident) node indices, merged at construction time."""
+    grid_path = gridpath("ugrid", "geoflow-small", "grid.nc")
+    grid = ux.open_grid(grid_path)
+
+    # The source file really does contain duplicates: 6000 node coordinates for
+    # 3850 distinct locations, so 2150 indices are coincident with an earlier one.
+    duplicates = _find_duplicate_nodes(grid)
+    assert grid.n_node == 6000
+    assert len(duplicates) == 2150
+
+    # Connectivity is canonicalized to a single index per coincident group, so no
+    # face references any of those 2150 duplicate indices.
+    assert not _check_duplicate_nodes_indices(grid)
+    # duplicate coordinates are left in place by design, but connectivity is
+    # fully canonicalized, so validation passes
+    assert grid.validate()
+
+    dual = grid.get_dual()
+
+    assert dual.n_node == grid.n_face
+
+    # One dual face per node that is a corner of at least three faces. After the
+    # merge, 3850 distinct nodes remain, ten of which are touched by a single face
+    # only and so produce no dual cell, leaving 3840.
+    face_nodes = grid.face_node_connectivity.values
+    faces_per_node = np.bincount(
+        face_nodes[face_nodes != INT_FILL_VALUE], minlength=grid.n_node
+    )
+    assert grid.n_node - len(duplicates) == 3850
+    assert (faces_per_node >= 3).sum() == 3840
+    assert dual.n_face == 3840
+
+    dataset = ux.open_dataset(grid_path, grid_path)
+    dual_ds = dataset.get_dual()
+    assert dual_ds.uxgrid.n_face == dual.n_face
+
+
+def test_dual_duplicate_geos_cs(gridpath):
+    """Test dual mesh creation on a cube-sphere grid with duplicate node
+    indices (issue #865)."""
+    grid_path = gridpath("geos-cs", "c12", "test-c12.native.nc4")
+    grid = ux.open_grid(grid_path)
+
+    assert len(_find_duplicate_nodes(grid)) > 0
+    assert not _check_duplicate_nodes_indices(grid)
+
+    dual = grid.get_dual()
+    assert dual.n_node == grid.n_face
+    assert dual.n_face > 0
+
+
+def test_duplicate_nodes_minimal_example():
+    """Two quads that share an edge, but whose shared corners are stored twice.
+
+    Nodes 2 and 3 are repeated as nodes 6 and 7, so the file describes 8 nodes at
+    6 distinct locations. Node 6 must canonicalize to node 2 and node 7 to node 3,
+    leaving the second face pointing at the first face's corners.
+
+        3---2---7      lat 1   nodes 2,3 are the shared edge
+        |   |   |              nodes 7,6 are their duplicates
+        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])
+    #                    left quad        right quad, via the duplicates
+    face_node_connectivity = np.array([[0, 1, 2, 3], [6, 4, 5, 7]])
+
+    grid = ux.Grid.from_topology(node_lon, node_lat, face_node_connectivity)
+
+    duplicates = _find_duplicate_nodes(grid)
+    assert duplicates == {6: 1, 7: 2}
+
+    # No face may still reference a duplicate index.
+    assert not _check_duplicate_nodes_indices(grid)
+    nt.assert_equal(
+        grid.face_node_connectivity.values, np.array([[0, 1, 2, 3], [1, 4, 5, 2]])
+    )
+
+
+def test_get_dual_rejects_faces_referencing_duplicate_nodes():
+    """``construct_dual`` reads ``node_face_connectivity`` with no duplicate
+    handling, so a face still pointing at a dead duplicate index would yield a
+    degenerate dual face instead of an error. Merging at construction makes this
+    unreachable today; the guard keeps it that way."""
+    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])
+    unmerged = np.array([[0, 1, 2, 3], [6, 4, 5, 7]])
+
+    grid = ux.Grid.from_topology(node_lon, node_lat, unmerged)
+    # Construction canonicalized the connectivity; put the duplicates back.
+    grid.face_node_connectivity = xr.DataArray(
+        unmerged, dims=grid.face_node_connectivity.dims
+    )
+
+    assert _check_duplicate_nodes_indices(grid)
+    with pytest.raises(GridInvalidError):
+        grid.get_dual()
+
+
+def test_pole_exception_uses_a_chord_tolerance():
+    """The pole carve-out must be a chord radius, not a raw ``|z|`` deviation.
+
+    ``np.isclose(|z|, 1.0, atol=tolerance)`` also carries numpy's default
+    ``rtol=1e-5``, so the carve-out spanned ``1 - |z| <= 1.001e-5`` -- a chord of
+    4.5e-3, or ~28 km on Earth. Every node within that cap was exempted from
+    merging. Only nodes at the pole itself may be exempt.
+    """
+    from uxarray.grid.validation import _coincident_node_canonical_indices
+
+    # Colatitude chosen so 1 - z = 1e-6: well inside the old carve-out, and far
+    # outside a chord of ERROR_TOLERANCE (whose cap is 1 - z <= 5e-17).
+    z = 1.0 - 1e-6
+    x = np.sqrt(1.0 - z * z)
+
+    points_xyz = np.array(
+        [
+            [0.0, 0.0, 1.0],  # north pole, kept distinct from the next node
+            [0.0, 0.0, 1.0],  # same location, its own face-specific longitude
+            [x, 0.0, z],  # near the pole, genuinely coincident with the next
+            [x, 0.0, z],
+        ]
+    )
+
+    canonical = _coincident_node_canonical_indices(points_xyz)
+
+    # Nodes at a pole are still never merged with one another.
+    nt.assert_equal(canonical[:2], np.array([0, 1]))
+    # Near-pole coincident nodes now merge; before the fix they were exempt.
+    nt.assert_equal(canonical[2:], np.array([2, 2]))
+
+
+def test_no_duplicate_nodes_ne30pg3(gridpath):
+    """``esmf/ne30/ne30pg3.grid.nc`` no longer reproduces issue #865's
+    duplicate-node bug; this only checks the general fix is a safe no-op."""
+    grid_path = gridpath("esmf", "ne30", "ne30pg3.grid.nc")
+    grid = ux.open_grid(grid_path)
+
+    assert len(_find_duplicate_nodes(grid)) == 0

test/test_subset.py

@@ -78,8 +78,12 @@ def test_grid_nn_subset(gridpath):
     for grid_path in GRID_PATHS:
         grid = ux.open_grid(grid_path)
 
-        # corner-nodes
-        ks = [1, 2, grid.n_node - 1]
+        # corner-nodes -- k is bounded by the number of *live* (non-duplicate)
+        # nodes, since the node search tree excludes dead coincident indices
+        from uxarray.grid.validation import _live_node_indices
+
+        n_live_nodes = len(_live_node_indices(grid))
+        ks = [1, 2, n_live_nodes - 1]
         for coord in coord_locs:
             for k in ks:
                 grid_subset = grid.subset.nearest_neighbor(coord,

uxarray/core/dataarray.py

@@ -2172,8 +2172,12 @@ def get_dual(self):
         --------
         dual : uxda
             Dual Mesh `uxda` constructed
-        """
 
+        Raises
+        ------
+        GridInvalidError
+            If any face still references a coincident duplicate node.
+        """
         if _check_duplicate_nodes_indices(self.uxgrid):
             raise GridInvalidError("Duplicate nodes found, cannot construct dual")
 

uxarray/core/dataset.py

@@ -755,8 +755,12 @@ def get_dual(self):
         --------
         dual : uxds
             Dual Mesh `uxds` constructed
-        """
 
+        Raises
+        ------
+        GridInvalidError
+            If any face still references a coincident duplicate node.
+        """
         if _check_duplicate_nodes_indices(self.uxgrid):
             raise GridInvalidError("Duplicate nodes found, cannot construct dual")
 

uxarray/grid/connectivity.py

@@ -2,7 +2,7 @@
 import xarray as xr
 from numba import njit
 
-from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
+from uxarray.constants import ERROR_TOLERANCE, INT_DTYPE, INT_FILL_VALUE
 from uxarray.conventions import ugrid
 
 
@@ -299,6 +299,161 @@ def _build_face_edge_connectivity(inverse_indices, n_face, n_max_face_nodes):
     return inverse_indices
 
 
+def _remap_node_connectivity(connectivity, duplicate_node_map, n_node):
+    """Return a copy of connectivity with duplicate node indices canonicalized."""
+    if not duplicate_node_map:
+        return connectivity
+
+    lookup = np.arange(n_node, dtype=INT_DTYPE)
+    keys = np.fromiter(
+        duplicate_node_map.keys(), dtype=INT_DTYPE, count=len(duplicate_node_map)
+    )
+    vals = np.fromiter(
+        duplicate_node_map.values(), dtype=INT_DTYPE, count=len(duplicate_node_map)
+    )
+    lookup[keys] = vals
+
+    remapped_connectivity = connectivity.copy()
+    valid = connectivity != INT_FILL_VALUE
+    remapped_connectivity[valid] = lookup[connectivity[valid]]
+    return remapped_connectivity
+
+
+def _collapse_repeated_face_corners(face_node_connectivity, canonical_values):
+    """Collapse consecutive (cyclically) repeated corners in each face row.
+
+    Remapping two originally-distinct, now-coincident corners of the same
+    face to a single canonical node can leave that node referenced twice in
+    a row, e.g. a quad (A, P, P, B) at a merged pole -- a triangle stored as
+    a 4-column row with one corner repeated. This pads it back down to
+    (A, P, B, FILL) so it is treated as the triangle it actually is.
+
+    Only rows containing a value in ``canonical_values`` (nodes that
+    absorbed at least one duplicate) are inspected, since no other row can
+    have gained a repeat from the remap.
+    """
+    if len(canonical_values) == 0:
+        return face_node_connectivity
+
+    affected_rows = np.flatnonzero(
+        np.isin(face_node_connectivity, canonical_values).any(axis=1)
+    )
+    if len(affected_rows) == 0:
+        return face_node_connectivity
+
+    face_node_connectivity = face_node_connectivity.copy()
+    for row_index in affected_rows:
+        row = face_node_connectivity[row_index]
+        valid = row != INT_FILL_VALUE
+        n_valid = int(valid.sum())
+        if n_valid <= 1:
+            continue
+
+        corners = row[:n_valid]
+        keep = corners != np.roll(corners, 1)
+        if keep.all():
+            continue
+
+        compacted = corners[keep]
+        new_row = np.full_like(row, INT_FILL_VALUE)
+        new_row[: len(compacted)] = compacted
+        face_node_connectivity[row_index] = new_row
+
+    return face_node_connectivity
+
+
+# node-index-valued connectivity: safe to remap element-wise in place
+_NODE_INDEX_CONNECTIVITY_TO_REMAP = ("face_node_connectivity", "node_node_connectivity")
+
+# connectivity derived from (and referencing) node indices, but whose rows must stay
+# unique (e.g. edge_node_connectivity) -- dropped rather than remapped, so the
+# existing lazy `@property` getters rebuild them cleanly from the corrected
+# face_node_connectivity instead of leaving phantom duplicate rows behind.
+_DERIVED_CONNECTIVITY_TO_INVALIDATE = (
+    "edge_node_connectivity",
+    "face_edge_connectivity",
+    "edge_face_connectivity",
+    "face_face_connectivity",
+    "node_edge_connectivity",
+    "node_face_connectivity",
+)
+
+
+def _merge_coincident_grid_ds_nodes(grid_ds, tolerance=ERROR_TOLERANCE):
+    """Canonicalize coincident (within ``tolerance``) node indices in a raw grid
+    dataset's connectivity, before it is wrapped in a ``Grid``.
+
+    Per issue #865, node coordinate/data arrays are left untouched -- only
+    connectivity references to coincident nodes are remapped to a single canonical
+    (lowest-indexed) node. Note this differs from TempestRemap and MOAB, which
+    delete the redundant nodes and renumber; keeping them preserves round-trip
+    fidelity and leaves node-centered data index-aligned, at the cost of leaving
+    unreferenced coordinates behind (see ``_live_node_indices``).
+    """
+    from uxarray.grid.coordinates import _lonlat_rad_to_xyz
+    from uxarray.grid.validation import _coincident_node_canonical_indices
+
+    if "face_node_connectivity" not in grid_ds:
+        return grid_ds
+
+    if {"node_x", "node_y", "node_z"} <= set(grid_ds.variables):
+        points_xyz = np.column_stack(
+            (
+                grid_ds["node_x"].values,
+                grid_ds["node_y"].values,
+                grid_ds["node_z"].values,
+            )
+        )
+    elif "node_lon" in grid_ds and "node_lat" in grid_ds:
+        points_xyz = np.column_stack(
+            _lonlat_rad_to_xyz(
+                np.deg2rad(grid_ds["node_lon"].values),
+                np.deg2rad(grid_ds["node_lat"].values),
+            )
+        )
+    else:
+        return grid_ds
+
+    n_node = points_xyz.shape[0]
+    canonical = _coincident_node_canonical_indices(points_xyz, tolerance)
+    duplicate_node_map = {
+        INT_DTYPE(index): INT_DTYPE(canonical[index])
+        for index in np.flatnonzero(canonical != np.arange(n_node, dtype=INT_DTYPE))
+    }
+    if not duplicate_node_map:
+        return grid_ds
+
+    grid_ds = grid_ds.copy()
+
+    for name in _NODE_INDEX_CONNECTIVITY_TO_REMAP:
+        if name in grid_ds:
+            grid_ds[name] = grid_ds[name].copy(
+                data=_remap_node_connectivity(
+                    grid_ds[name].values, duplicate_node_map, n_node
+                )
+            )
+
+    if "face_node_connectivity" in grid_ds:
+        canonical_values = np.unique(
+            np.fromiter(
+                duplicate_node_map.values(),
+                dtype=INT_DTYPE,
+                count=len(duplicate_node_map),
+            )
+        )
+        grid_ds["face_node_connectivity"] = grid_ds["face_node_connectivity"].copy(
+            data=_collapse_repeated_face_corners(
+                grid_ds["face_node_connectivity"].values, canonical_values
+            )
+        )
+
+    for name in _DERIVED_CONNECTIVITY_TO_INVALIDATE:
+        if name in grid_ds:
+            grid_ds = grid_ds.drop_vars(name)
+
+    return grid_ds
+
+
 def _populate_node_face_connectivity(grid):
     """Constructs the UGRID connectivity variable (``node_face_connectivity``)
     and stores it within the internal (``Grid._ds``) and through the attribute