Merge coincident pole and antimeridian nodes in structured grids

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()

uxarray/io/_structured.py

@@ -1,7 +1,7 @@
 import numpy as np
 import xarray as xr
 
-from uxarray.constants import INT_DTYPE
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE
 from uxarray.conventions import ugrid
 
 
@@ -21,7 +21,7 @@ def _read_structured_grid(lon, lat, tol=1e-10):
     lat : array_like
         1D array of latitude coordinates in degrees.
     tol : float, optional
-        Tolerance for considering nodes as identical (default is `1e-10`).
+        Tolerance in degrees for considering nodes as identical (default is `1e-10`).
 
     Returns
     -------
@@ -79,11 +79,25 @@ def _read_structured_grid(lon, lat, tol=1e-10):
     # Stack longitude and latitude for processing
     nodes = np.column_stack((node_lon, node_lat))
 
+    # Match nodes on the sphere rather than in the lon/lat plane, so that the poles
+    # (many lon values, one point) and the antimeridian seam (lon differing by 360)
+    # are recognized as coincident.
+    lon_rad = np.deg2rad(node_lon)
+    lat_rad = np.deg2rad(node_lat)
+    cos_lat = np.cos(lat_rad)
+    node_xyz = np.column_stack(
+        (cos_lat * np.cos(lon_rad), cos_lat * np.sin(lon_rad), np.sin(lat_rad))
+    )
+
     # Build KDTree
-    tree = KDTree(nodes)
+    tree = KDTree(node_xyz)
+
+    # ``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)
 
     # Find all pairs of nodes within the tolerance
-    pairs = tree.query_pairs(r=tol)
+    pairs = tree.query_pairs(r=chord_tol)
 
     n_nodes = len(nodes)
     if pairs:
@@ -138,6 +152,20 @@ def _read_structured_grid(lon, lat, tol=1e-10):
     # Stack the node indices to form face_node_connectivity
     face_node_conn = np.vstack((n1, n2, n3, n4), dtype=INT_DTYPE).T
 
+    # Merging the poles leaves their quads with a repeated corner; drop it so those
+    # faces are stored as the triangles they are, padded with the fill value.
+    keep = face_node_conn != np.roll(face_node_conn, 1, axis=1)
+    if not keep.all():
+        n_nodes_per_face = keep.sum(axis=1)
+        order = np.argsort(~keep, axis=1, kind="stable")
+        compacted = np.take_along_axis(face_node_conn, order, axis=1)
+        n_max_face_nodes = n_nodes_per_face.max()
+        compacted = compacted[:, :n_max_face_nodes]
+        compacted[np.arange(n_max_face_nodes) >= n_nodes_per_face[:, None]] = (
+            INT_FILL_VALUE
+        )
+        face_node_conn = compacted
+
     out_ds["node_lon"] = xr.DataArray(
         data=unique_node_lon, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LON_ATTRS
     )