Merge pull request #319 from Deltares/from-structured-curvilinear

From structured curvilinear

docs/api.rst

@@ -43,7 +43,7 @@ UgridDataArray
 
     UgridDataArray
     UgridDataArray.ugrid
-    UgridDataArray.from_structured
+    UgridDataArray.from_structured2d
     UgridDataArray.from_data
 
 UgridDataset
@@ -55,7 +55,7 @@ UgridDataset
     UgridDataset
     UgridDataset.ugrid
     UgridDataset.from_geodataframe
-    UgridDataset.from_structured
+    UgridDataset.from_structured2d
 
 UGRID Accessor
 --------------

docs/changelog.rst

@@ -9,6 +9,29 @@ The format is based on `Keep a Changelog`_, and this project adheres to
 Unreleased
 ----------
 
+Changed
+~~~~~~~
+
+- :meth:`xugrid.UgridDataset.from_structured` and
+  :meth:`xugrid.UgridDataArray.from_structured` are deprecated and will be
+  removed in the future; calling them will raise a FutureWarning. They have
+  been replaced by :meth:`xugrid.UgridDataset.from_structured2d` and
+  :meth:`xugrid.UgridDataArray.from_structured2d` respectively.
+
+Added
+~~~~~
+
+- :meth:`xugrid.Ugrid2d.from_structured_bounds` now accepts 3D bounds to allow
+  conversion of grids with non-monotonic x and y coordinates, such as strongly
+  curvilinear grids.
+- :meth:`xugrid.Ugrid2d.from_structured_bounds` now takes an optional
+  ``return_index`` argument to return the indices of invalid grid faces,
+  identified by one or more NaNs in its bounds.
+- This method is used in :meth:`xugrid.UgridDataArray.from_structured2d` and
+  :meth:`xugrid.UgridDataset.from_structured2d` when the optional arguments
+  ``x_bounds`` and ``y_bounds`` are provided.
+
+
 [0.12.2] 2025-01-31
 -------------------
 

tests/fixtures/fixture_regridder.py

@@ -47,7 +47,7 @@ def quads_0_25():
             "x",
         ],
     )
-    return xu.UgridDataArray.from_structured(da)
+    return xu.UgridDataArray.from_structured2d(da)
 
 
 @pytest.fixture(scope="function")
@@ -83,7 +83,7 @@ def quads_1():
             "x",
         ],
     )
-    return xu.UgridDataArray.from_structured(da)
+    return xu.UgridDataArray.from_structured2d(da)
 
 
 @pytest.fixture(scope="function")

tests/test_conversion.py

@@ -253,9 +253,9 @@ def test_infer_breaks_intervals1d_errors(structured_mesh_ascending):
         cv.infer_interval_breaks1d(up, "x")
 
 
-def test_bounds_to_vertices():
+def test_bounds1d_to_vertices():
     with pytest.raises(ValueError, match="Bounds are not monotonic"):
-        cv.bounds_to_vertices(
+        cv.bounds1d_to_vertices(
             xr.DataArray(
                 data=[[0.0, 1.0], [2.0, 3.0], [1.0, 2.0]], dims=["x", "nbound"]
             )
@@ -266,15 +266,85 @@ def test_bounds_to_vertices():
     # Ascending
     x_bounds = np.column_stack((x_vertices[:-1], x_vertices[1:]))
     y_bounds = np.column_stack((y_vertices[:-1], y_vertices[1:]))
-    assert np.allclose(cv.bounds_to_vertices(x_bounds), x_vertices)
-    assert np.allclose(cv.bounds_to_vertices(y_bounds), y_vertices)
+    assert np.allclose(cv.bounds1d_to_vertices(x_bounds), x_vertices)
+    assert np.allclose(cv.bounds1d_to_vertices(y_bounds), y_vertices)
     # Descending
     xrev = x_vertices[::-1]
     yrev = y_vertices[::-1]
     x_bounds = np.column_stack((xrev[1:], xrev[:-1]))
     y_bounds = np.column_stack((yrev[1:], yrev[:-1]))
-    assert np.allclose(cv.bounds_to_vertices(x_bounds), xrev)
-    assert np.allclose(cv.bounds_to_vertices(y_bounds), yrev)
+    assert np.allclose(cv.bounds1d_to_vertices(x_bounds), xrev)
+    assert np.allclose(cv.bounds1d_to_vertices(y_bounds), yrev)
+
+
+def test_bounds2d_to_topology2d():
+    # Clockwise
+    x_bounds = np.array(
+        [[[0.0, 0.0, 1.0, 1.0], [2.0, 2.0, 3.0, 3.0], [4.0, 4.0, 5.0, 5.0]]]
+    )
+    y_bounds = np.array(
+        [[[0.0, 1.0, 1.0, 0.0], [2.0, 3.0, 3.0, 2.0], [4.0, 5.0, 5.0, 4.0]]]
+    )
+    x, y, faces, index = cv.bounds2d_to_topology2d(x_bounds, y_bounds)
+
+    assert index.all()
+    assert faces.shape == (3, 4)
+
+    # Result should be counterclockwise
+    expected_first_cell = np.array(
+        [
+            [0.0, 0.0],
+            [1.0, 0.0],
+            [1.0, 1.0],
+            [0.0, 1.0],
+        ]
+    )
+    actual = np.column_stack((x, y))[faces[0]]
+    assert np.allclose(actual, expected_first_cell)
+
+    # Test with some invalid (NaN) coordinates
+    x_bounds_with_nan = x_bounds.copy()
+    x_bounds_with_nan[0, 0, 0] = np.nan
+    x, y, faces, index = cv.bounds2d_to_topology2d(x_bounds_with_nan, y_bounds)
+
+    # Test that the first cell is marked invalid
+    assert not index[0]
+    assert index[1:].all()
+    assert faces.shape == (2, 4)
+
+    # Test bad bounds. Triangles are allowed, points and lines are not.
+    x_bounds = np.array(
+        [
+            [
+                [0.0, 0.0, 0.0, 0.0],
+                [1.0, 2.0, 2.0, 1.0],
+                [2.0, 3.0, 3.0, 2.0],
+                [2.0, 2.0, 3.0, 3.0],
+            ]
+        ]
+    )
+    y_bounds = np.array(
+        [
+            [
+                [0.0, 0.0, 0.0, 0.0],
+                [0.0, 0.0, 1.0, 1.0],
+                [0.0, 0.0, 0.0, 1.0],
+                [0.0, 0.0, 1.0, 1.0],
+            ]
+        ]
+    )
+
+    with pytest.warns(
+        UserWarning, match="A UGRID2D face requires at least three unique vertices."
+    ):
+        x, y, faces, index = cv.bounds2d_to_topology2d(x_bounds, y_bounds)
+
+    assert np.array_equal(index, [False, True, True, False])
+    expected_x = [
+        [1.0, 2.0, 2.0, 1.0],
+        [2.0, 3.0, 2.0, 3.0],  # repeated node moved to last one
+    ]
+    assert np.array_equal(x[faces], expected_x)
 
 
 def test_infer_xy_coords():

tests/test_regrid/test_regridder.py

@@ -302,14 +302,14 @@ def test_directional_dependence():
         coords={"y": [10.0, 20.0], "x": [0.0, 10.0]},
         dims=("y", "x"),
     )
-    target_uda = xu.UgridDataArray.from_structured(target_da)
+    target_uda = xu.UgridDataArray.from_structured2d(target_da)
 
     flip = slice(None, None, -1)
     flipy = da.isel(y=flip)
     flipx = da.isel(x=flip)
     flipxy = da.isel(x=flip, y=flip)
-    uda = xu.UgridDataArray.from_structured(da)
-    uda_flipxy = xu.UgridDataArray.from_structured(flipxy)
+    uda = xu.UgridDataArray.from_structured2d(da)
+    uda_flipxy = xu.UgridDataArray.from_structured2d(flipxy)
 
     # Structured target: test whether the result is the same regardless of source
     # orientation.
@@ -383,7 +383,7 @@ def test_barycentric_structured():
     regridder = xu.BarycentricInterpolator(source=da, target=target)
     out_structured = regridder.regrid(da)
 
-    target_uda = xu.UgridDataArray.from_structured(target)
+    target_uda = xu.UgridDataArray.from_structured2d(target)
     regridder = xu.BarycentricInterpolator(source=da, target=target_uda)
     out_unstructured = regridder.regrid(da)
 

tests/test_ugrid2d.py

@@ -1198,6 +1198,20 @@ def test_from_structured_bounds():
     assert isinstance(grid, xugrid.Ugrid2d)
     assert grid.n_face == 12
 
+    x_bounds = np.array(
+        [[[0.0, 0.0, 1.0, 1.0], [2.0, 2.0, 3.0, 3.0], [4.0, 4.0, 5.0, 5.0]]]
+    )
+    y_bounds = np.array(
+        [[[0.0, 1.0, 1.0, 0.0], [2.0, 3.0, 3.0, 2.0], [4.0, 5.0, 5.0, 4.0]]]
+    )
+    grid = xugrid.Ugrid2d.from_structured_bounds(x_bounds, y_bounds)
+    assert grid.n_face == 3
+    assert np.allclose(grid.area, 1.0)
+    assert grid.bounds == (0.0, 0.0, 5.0, 5.0)
+
+    with pytest.raises(ValueError, match="Bounds shapes do not match"):
+        xugrid.Ugrid2d.from_structured_bounds(x_bounds, y_bounds.T)
+
 
 def test_from_structured():
     da = xr.DataArray(

tests/test_ugrid_dataset.py

@@ -811,54 +811,60 @@ def setup(self):
             }
         )
 
+    def test_future_warnings(self):
+        with pytest.warns(FutureWarning):
+            xugrid.UgridDataArray.from_structured(self.da2d)
+        with pytest.warns(FutureWarning):
+            xugrid.UgridDataset.from_structured(self.ds)
+
     def test_error_1d(self):
         with pytest.raises(
             ValueError, match="DataArray must have at least two spatial dimensions"
         ):
-            xugrid.UgridDataArray.from_structured(self.da1d)
+            xugrid.UgridDataArray.from_structured2d(self.da1d)
 
     def test_error_x_xor_y(self):
         with pytest.raises(ValueError, match="Provide both x and y, or neither."):
-            xugrid.UgridDataArray.from_structured(self.da2d, x="this")
+            xugrid.UgridDataArray.from_structured2d(self.da2d, x="this")
 
     def test_missing_xy(self):
         with pytest.raises(ValueError, match="Coordinates xc and yc are not present."):
-            xugrid.UgridDataArray.from_structured(self.da2d, x="xc", y="yc")
+            xugrid.UgridDataArray.from_structured2d(self.da2d, x="xc", y="yc")
 
     def test_from_dataarray(self):
-        uda = xugrid.UgridDataArray.from_structured(self.da2d)
+        uda = xugrid.UgridDataArray.from_structured2d(self.da2d)
         assert isinstance(uda, xugrid.UgridDataArray)
         assert uda.name == "grid"
         assert uda.dims == ("layer", "mesh2d_nFaces")
         assert uda.shape == (2, 12)
         assert np.allclose(uda.ugrid.sel(x=2.0, y=5.0), [[0], [12]])
         # Check whether flipping the y-axis doesn't cause any problems
         flipped = self.da2d.isel(y=slice(None, None, -1))
-        uda = xugrid.UgridDataArray.from_structured(flipped)
+        uda = xugrid.UgridDataArray.from_structured2d(flipped)
         assert np.allclose(uda.ugrid.sel(x=2.0, y=5.0), [[0], [12]])
 
         # And test transposed.
         daT = self.da2d.transpose()
-        uda = xugrid.UgridDataArray.from_structured(daT)
+        uda = xugrid.UgridDataArray.from_structured2d(daT)
         assert isinstance(uda, xugrid.UgridDataArray)
         assert uda.name == "grid"
         assert uda.dims == ("layer", "mesh2d_nFaces")
         assert uda.shape == (2, 12)
         assert np.allclose(uda.ugrid.sel(x=2.0, y=5.0), [[0], [12]])
 
     def test_from_multicoord(self):
-        uda = xugrid.UgridDataArray.from_structured(self.da_coords2d)
+        uda = xugrid.UgridDataArray.from_structured2d(self.da_coords2d)
         assert isinstance(uda, xugrid.UgridDataArray)
         assert np.array_equal(np.unique(uda.ugrid.grid.node_x), [-0.5, 0.5, 1.5])
         assert np.array_equal(uda.data, [0, 1, 2, 3])
 
-        uda = xugrid.UgridDataArray.from_structured(self.da_coords2d, x="xc", y="yc")
+        uda = xugrid.UgridDataArray.from_structured2d(self.da_coords2d, x="xc", y="yc")
         assert isinstance(uda, xugrid.UgridDataArray)
         assert np.array_equal(np.unique(uda.ugrid.grid.node_x), [9.0, 11.0, 13.0])
         assert np.array_equal(uda.data, [0, 1, 2, 3])
 
     def test_from_dataset(self):
-        uds = xugrid.UgridDataset.from_structured(self.ds)
+        uds = xugrid.UgridDataset.from_structured2d(self.ds)
         assert isinstance(uds, xugrid.UgridDataset)
         assert set(uds.data_vars) == {"a", "b", "c"}
         assert uds["a"].dims == ("layer", "mesh2d_nFaces")
@@ -873,13 +879,13 @@ def test_from_multicoord_dataset(self):
         da = self.da_coords2d.rename({"x": "x1", "y": "y1"})
         ds["d"] = da
         # Unspecified: it'll only infer x and y.
-        uds = xugrid.UgridDataset.from_structured(ds)
+        uds = xugrid.UgridDataset.from_structured2d(ds)
         assert isinstance(uds, xugrid.UgridDataset)
         assert uds["a"].dims == ("layer", "mesh2d_nFaces")
         assert uds["d"].dims == ("y1", "x1")
         assert len(uds.ugrid.grids) == 1
         # Now specify separate topologies.
-        uds = xugrid.UgridDataset.from_structured(
+        uds = xugrid.UgridDataset.from_structured2d(
             ds, {"mesh2d_0": ("x", "y"), "mesh2d_1": ("xc", "yc")}
         )
         assert isinstance(uds, xugrid.UgridDataset)
@@ -889,6 +895,63 @@ def test_from_multicoord_dataset(self):
         assert uds["d"].dims == ("mesh2d_1_nFaces",)
         assert len(uds.ugrid.grids) == 2
 
+    def test_from_bounds(self):
+        # Construct bounds from ugrid form, it's easier...
+        uda = xugrid.UgridDataArray.from_structured2d(self.da2d)
+        grid = uda.ugrid.grid
+        bounds_x = xr.DataArray(
+            grid.face_node_coordinates[..., 0].reshape(3, 4, 4),
+            dims=("y", "x", "bound"),
+        )
+        bounds_y = xr.DataArray(
+            grid.face_node_coordinates[..., 1].reshape(3, 4, 4),
+            dims=("y", "x", "bound"),
+        )
+        uda2 = xugrid.UgridDataArray.from_structured2d(
+            self.da2d, "x", "y", bounds_x, bounds_y
+        )
+        assert uda.identical(uda2)
+        # Try inconsistent dim order on bounds:
+        uda3 = xugrid.UgridDataArray.from_structured2d(
+            self.da2d, "x", "y", bounds_x.transpose(), bounds_y
+        )
+        assert uda.identical(uda3)
+
+        with pytest.raises(ValueError, match="x and y must be provided for bounds"):
+            xugrid.UgridDataArray.from_structured2d(
+                self.da2d, x_bounds=bounds_x, y_bounds=bounds_y
+            )
+
+        # Now dataset
+        ds = self.ds.copy()
+        ds["grid_x"] = bounds_x
+        ds["grid_y"] = bounds_y
+        uds = xugrid.UgridDataset.from_structured2d(
+            ds,
+            topology={
+                "mesh2d": {
+                    "x": "x",
+                    "y": "y",
+                    "x_bounds": "grid_x",
+                    "y_bounds": "grid_y",
+                }
+            },
+        )
+        assert set(uds.data_vars) == {"a", "b", "c", "grid_x", "grid_y"}
+        assert uds["a"].dims == ("layer", "mesh2d_nFaces")
+        assert uds["b"].dims == ("x",)
+        assert uds["c"].dims == ()
+
+        with pytest.raises(TypeError):
+            xugrid.UgridDataset.from_structured2d(
+                ds, topology={"x_bounds": "grid_x", "y_bounds": "grid_y"}
+            )
+
+        with pytest.raises(ValueError, match="x and y must be provided for bounds"):
+            xugrid.UgridDataset.from_structured2d(
+                ds, topology={"mesh2d": {"x_bounds": "grid_x", "y_bounds": "grid_y"}}
+            )
+
 
 def test_multiple_coordinates():
     grid = GRID()