Merge duplicate nodes when constructing the dual mesh

test/grid/grid/test_core.py

@@ -129,7 +129,53 @@ 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."""
+    from uxarray.grid.validation import _check_duplicate_nodes_indices, _find_duplicate_nodes
+
+    grid_path = gridpath("ugrid", "geoflow-small", "grid.nc")
+    grid = ux.open_grid(grid_path)
+
+    # source file has duplicate node coordinates; connectivity should already
+    # be canonicalized to a single index per coincident group
+    assert len(_find_duplicate_nodes(grid)) > 0
+    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
+    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)."""
+    from uxarray.grid.validation import _check_duplicate_nodes_indices, _find_duplicate_nodes
+
+    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_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."""
+    from uxarray.grid.validation import _find_duplicate_nodes
+
+    grid_path = gridpath("esmf", "ne30", "ne30pg3.grid.nc")
+    grid = ux.open_grid(grid_path)
+
+    assert len(_find_duplicate_nodes(grid)) == 0

test/io/test_structured.py

@@ -74,3 +74,52 @@ def test_from_xarray_with_grid_from_latlon(ds_name):
     subset = uxds["air"].isel(time=0).subset.bounding_circle((-100.0, 40.0), 5)
     assert "n_face" in subset.dims
     assert subset.sizes["n_face"] > 0
+
+
+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
+
+    n_lon, n_lat = 36, 18
+    d_lat = 180.0 / n_lat
+    lon = np.linspace(-180, 180, n_lon, endpoint=False)
+    lat = np.linspace(-90 + d_lat / 2, 90 - d_lat / 2, n_lat)
+
+    uxgrid = ux.Grid.from_structured(lon=lon, lat=lat)
+
+    # Every duplicated pole node and antimeridian node must be gone.
+    assert uxgrid.n_node < (n_lon + 1) * (n_lat + 1)
+    assert np.isclose(uxgrid.node_lat.values, 90.0).sum() == 1
+    assert np.isclose(uxgrid.node_lat.values, -90.0).sum() == 1
+
+    # A closed sphere: V - E + F == 2.
+    assert uxgrid.n_node - uxgrid.n_edge + uxgrid.n_face == 2
+
+    # The pole is now a real singularity touching every longitude column, and
+    # its faces are triangles rather than quads with a repeated corner.
+    face_nodes = uxgrid.face_node_connectivity.values
+    n_nodes_per_face = uxgrid.n_nodes_per_face.values
+    assert (n_nodes_per_face == 3).sum() == 2 * n_lon
+
+    for face, n_nodes in zip(face_nodes, n_nodes_per_face):
+        nodes = face.tolist()[:n_nodes]
+        assert len(set(nodes)) == n_nodes
+
+    pole = int(np.flatnonzero(np.isclose(uxgrid.node_lat.values, 90.0))[0])
+    assert (face_nodes == pole).any(axis=1).sum() == n_lon
+
+
+def test_regional_structured_grid_is_unchanged():
+    """A grid that touches neither pole nor the antimeridian must keep every
+    node and stay entirely quadrilateral."""
+    import numpy as np
+
+    lon = np.linspace(-50, -10, 20)
+    lat = np.linspace(10, 40, 15)
+
+    uxgrid = ux.Grid.from_structured(lon=lon, lat=lat)
+
+    assert uxgrid.n_node == 21 * 16
+    assert uxgrid.n_face == 20 * 15
+    assert (uxgrid.n_nodes_per_face.values == 4).all()

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

@@ -34,7 +34,6 @@
 from uxarray.formatting_html import array_repr
 from uxarray.grid import Grid
 from uxarray.grid.dual import construct_dual
-from uxarray.grid.validation import _check_duplicate_nodes_indices
 from uxarray.io._healpix import get_zoom_from_cells
 from uxarray.plot.accessor import UxDataArrayPlotAccessor
 from uxarray.remap.accessor import RemapAccessor
@@ -2177,9 +2176,6 @@ def get_dual(self):
             Dual Mesh `uxda` constructed
         """
 
-        if _check_duplicate_nodes_indices(self.uxgrid):
-            raise GridInvalidError("Duplicate nodes found, cannot construct dual")
-
         if self.uxgrid.partial_sphere_coverage:
             warn(
                 "This mesh is partial, which could cause inconsistent results and data will be lost",

uxarray/core/dataset.py

@@ -18,7 +18,6 @@
 from uxarray.formatting_html import dataset_repr
 from uxarray.grid import Grid
 from uxarray.grid.dual import construct_dual
-from uxarray.grid.validation import _check_duplicate_nodes_indices
 from uxarray.io._healpix import get_zoom_from_cells
 from uxarray.plot.accessor import UxDatasetPlotAccessor
 from uxarray.remap.accessor import RemapAccessor
@@ -706,9 +705,6 @@ def get_dual(self):
             Dual Mesh `uxds` constructed
         """
 
-        if _check_duplicate_nodes_indices(self.uxgrid):
-            raise GridInvalidError("Duplicate nodes found, cannot construct dual")
-
         if self.uxgrid.partial_sphere_coverage:
             warn(
                 "This mesh is partial, which could cause inconsistent results and data will be lost",

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,158 @@ 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 _dedupe_grid_ds_nodes(grid_ds, tolerance=ERROR_TOLERANCE):
+    """Canonicalize duplicate (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 duplicate nodes are remapped to a single canonical
+    (lowest-indexed) node.
+    """
+    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