From 1b6224e07ae485663563f07d0b37ef1a5a903c90 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Wed, 19 Aug 2026 16:51:58 -0500 Subject: [PATCH 1/4] Merge coincident nodes in structured grids on the sphere 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. --- test/io/test_structured.py | 49 ++++++++++++++++++++++++++++++++++++++ uxarray/io/_structured.py | 28 ++++++++++++++++++++-- 2 files changed, 75 insertions(+), 2 deletions(-) diff --git a/test/io/test_structured.py b/test/io/test_structured.py index 54e735e4b..166387a18 100644 --- a/test/io/test_structured.py +++ b/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() diff --git a/uxarray/io/_structured.py b/uxarray/io/_structured.py index e3b8cab1e..f88d1b4d2 100644 --- a/uxarray/io/_structured.py +++ b/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 @@ -79,8 +79,18 @@ 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) # Find all pairs of nodes within the tolerance pairs = tree.query_pairs(r=tol) @@ -138,6 +148,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 ) From 3b1016efd05137b6e70183756960bbc05d4de8fe Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Wed, 19 Aug 2026 17:02:14 -0500 Subject: [PATCH 2/4] Keep tol in degrees when matching nodes on the sphere --- uxarray/io/_structured.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/uxarray/io/_structured.py b/uxarray/io/_structured.py index f88d1b4d2..7f7c47f24 100644 --- a/uxarray/io/_structured.py +++ b/uxarray/io/_structured.py @@ -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 ------- @@ -92,8 +92,12 @@ def _read_structured_grid(lon, lat, tol=1e-10): # Build KDTree 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: From 70cade5acb9265e31fd788c2358b5a31e30fa84e Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 21 Aug 2026 11:00:33 -0500 Subject: [PATCH 3/4] Address review: derive tol from ERROR_TOLERANCE, reuse xyz helper, clarify comment --- uxarray/grid/grid.py | 16 ++++++++++++---- uxarray/io/_structured.py | 37 ++++++++++++++++++++++++++++--------- 2 files changed, 40 insertions(+), 13 deletions(-) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index 711b30389..894718cd6 100644 --- a/uxarray/grid/grid.py +++ b/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 ------- diff --git a/uxarray/io/_structured.py b/uxarray/io/_structured.py index 7f7c47f24..39a45a928 100644 --- a/uxarray/io/_structured.py +++ b/uxarray/io/_structured.py @@ -1,11 +1,17 @@ import numpy as np import xarray as xr -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 +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 in degrees 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 ------- @@ -84,10 +92,7 @@ def _read_structured_grid(lon, lat, tol=1e-10): # 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)) - ) + node_xyz = np.column_stack(_lonlat_rad_to_xyz(lon_rad, lat_rad)) # Build KDTree tree = KDTree(node_xyz) @@ -152,8 +157,22 @@ 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. + # 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) From dde8a1811808a47d3ad03be39096e8f06dac9ee8 Mon Sep 17 00:00:00 2001 From: Rajeev Jain Date: Fri, 21 Aug 2026 11:01:01 -0500 Subject: [PATCH 4/4] Upcast lon/lat to float64 before coincidence detection 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. --- uxarray/io/_structured.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/uxarray/io/_structured.py b/uxarray/io/_structured.py index 39a45a928..cced26e20 100644 --- a/uxarray/io/_structured.py +++ b/uxarray/io/_structured.py @@ -50,6 +50,13 @@ def _read_structured_grid(lon, lat, tol=_DEFAULT_STRUCTURED_TOL_DEG): 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]