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/grid/grid.py

@@ -89,7 +89,10 @@
 from uxarray.io._icon import _read_icon
 from uxarray.io._mpas import _read_mpas
 from uxarray.io._scrip import _encode_scrip, _read_scrip
-from uxarray.io._structured import _read_structured_grid
+from uxarray.io._structured import (
+    _DEFAULT_STRUCTURED_TOL_DEG,
+    _read_structured_grid,
+)
 from uxarray.io._topology import _read_topology
 from uxarray.io._ugrid import (
     _encode_ugrid,
@@ -501,7 +504,11 @@ def from_topology(
 
     @classmethod
     def from_structured(
-        cls, ds: xr.Dataset = None, lon=None, lat=None, tol: float | None = 1e-10
+        cls,
+        ds: xr.Dataset = None,
+        lon=None,
+        lat=None,
+        tol: float | None = _DEFAULT_STRUCTURED_TOL_DEG,
     ):
         """
         Converts a structured ``xarray.Dataset`` or longitude and latitude coordinates into an unstructured ``uxarray.Grid``.
@@ -525,8 +532,9 @@ def from_structured(
             Should be a one-dimensional or two-dimensional array following CF conventions.
 
         tol : float, optional
-            Tolerance for considering nodes as identical when constructing the grid from longitude and latitude.
-            Default is `1e-10`.
+            Tolerance in degrees for considering nodes as identical when constructing the grid from
+            longitude and latitude. Defaults to the angle whose chord length on the unit sphere equals
+            ``uxarray.constants.ERROR_TOLERANCE``.
 
         Returns
         -------

uxarray/io/_structured.py

@@ -1,11 +1,17 @@
 import numpy as np
 import xarray as xr
 
-from uxarray.constants import INT_DTYPE
+from uxarray.constants import ERROR_TOLERANCE, INT_DTYPE, INT_FILL_VALUE
 from uxarray.conventions import ugrid
+from uxarray.grid.coordinates import _lonlat_rad_to_xyz
 
+# ``ERROR_TOLERANCE`` is defined as a Cartesian distance on the unit sphere; convert
+# it to the angular tolerance `tol` is documented in so the default tracks the same
+# precision assumption used everywhere else in the codebase.
+_DEFAULT_STRUCTURED_TOL_DEG = np.rad2deg(2.0 * np.arcsin(ERROR_TOLERANCE / 2.0))
 
-def _read_structured_grid(lon, lat, tol=1e-10):
+
+def _read_structured_grid(lon, lat, tol=_DEFAULT_STRUCTURED_TOL_DEG):
     """
     Constructs an unstructured grid dataset from structured longitude and latitude coordinates.
 
@@ -21,7 +27,9 @@ 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. Defaults to the angle
+        whose chord length on the unit sphere equals ``uxarray.constants.ERROR_TOLERANCE``,
+        matching the precision assumption used elsewhere in the codebase.
 
     Returns
     -------
@@ -42,6 +50,13 @@ def _read_structured_grid(lon, lat, tol=1e-10):
 
     out_ds = xr.Dataset()
 
+    # Coincidence detection below relies on float64 precision (~1e-16); real-world
+    # datasets often store lon/lat as float32 (~1e-7), which silently propagates
+    # through this pipeline and causes pole/antimeridian merges to fail or merge
+    # only partially, regardless of ``tol``.
+    lon = np.asarray(lon, dtype=np.float64)
+    lat = np.asarray(lat, dtype=np.float64)
+
     sorted_indices = np.argsort(lon)
     lon = lon[sorted_indices]
 
@@ -79,11 +94,22 @@ 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)
+    node_xyz = np.column_stack(_lonlat_rad_to_xyz(lon_rad, 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 +164,34 @@ 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
 
+    # No new faces are created here -- this only shrinks the width of existing rows
+    # in face_node_conn for faces that became degenerate after the pole merge above.
+    #
+    # A face touching the pole is built from 2 distinct edge-longitudes at the pole
+    # latitude, e.g. corners (n1, n2, n3, n4) = (A, P, P, B), where P is the single
+    # merged pole node that both pole-row corners now point to (n2 == n3). That is a
+    # triangle A-P-B stored as a 4-column quad with one corner repeated, so:
+    #   1. `keep` marks, per face, which corners differ from their cyclic predecessor
+    #      (n2 == n3 above means the P at position 2 is dropped from that row).
+    #   2. The kept corners are pushed to the front of each row (`order`), giving
+    #      (A, P, B, B) instead of (A, P, P, B) -- still 4 columns, but the last
+    #      column is now the padding slot for a 3-node face.
+    #   3. `n_max_face_nodes` is the largest node count any face still needs (3 here,
+    #      unless some other face in the grid still has 4 distinct corners, in which
+    #      case nothing is trimmed and this is a no-op). Columns beyond each face's
+    #      own count are set to `INT_FILL_VALUE`, giving (A, P, B, FILL).
+    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
     )