Merge pull request #154 from Deltares/aperiodic

Period and aperiodic boundaries

docs/api.rst

@@ -80,6 +80,10 @@ UgridDataArray or UgridDataset.
     UgridDatasetAccessor.to_crs
     UgridDatasetAccessor.sel
     UgridDatasetAccessor.sel_points
+    UgridDatasetAccessor.rasterize
+    UgridDatasetAccessor.rasterize_like
+    UgridDatasetAccessor.to_periodic
+    UgridDatasetAccessor.to_nonperiodic
     UgridDatasetAccessor.intersect_line
     UgridDatasetAccessor.intersect_linestring
     UgridDatasetAccessor.partition
@@ -112,6 +116,8 @@ UgridDataArray or UgridDataset.
     UgridDataArrayAccessor.partition_by_label
     UgridDataArrayAccessor.rasterize
     UgridDataArrayAccessor.rasterize_like
+    UgridDataArrayAccessor.to_periodic
+    UgridDataArrayAccessor.to_nonperiodic
     UgridDataArrayAccessor.to_geodataframe
     UgridDataArrayAccessor.binary_dilation
     UgridDataArrayAccessor.binary_erosion
@@ -254,7 +260,9 @@ UGRID2D Topology
     Ugrid2d.face_dimension
     Ugrid2d.face_coordinates
     Ugrid2d.centroids
+    Ugrid2d.circumcenters
     Ugrid2d.area
+    Ugrid2d.perimeter
     Ugrid2d.face_x
     Ugrid2d.face_y
 
@@ -298,6 +306,8 @@ UGRID2D Topology
     Ugrid2d.locate_bounding_box
     Ugrid2d.rasterize
     Ugrid2d.rasterize_like
+    Ugrid2d.to_periodic
+    Ugrid2d.to_nonperiodic
     Ugrid2d.topology_subset
     Ugrid2d.label_partitions
     Ugrid2d.partition

docs/changelog.rst

@@ -12,9 +12,24 @@ Unreleased
 Added
 ~~~~~
 
+- :meth:`xugrid.Ugrid2d.to_nonperiodic`,
+  :meth:`xugrid.UgridDataArrayAccessor.to_nonperiodic` and
+  :meth:`xugrid.UgridDatasetAccessor.to_nonperiodic` have been added to convert
+  a "periodid grid" (where the leftmost nodes are the same as the rightmost
+  nodes, e.g. a mesh for the globe) to an "ordinary" grid.
+- Conversely, :meth:`xugrid.Ugrid2d.to_periodic`,
+  :meth:`xugrid.UgridDataArrayAccessor.to_periodic` and
+  :meth:`xugrid.UgridDatasetAccessor.to_periodic` have been added to convert an
+  ordinary grid to a periodic grid.
+- :attr:`xugrid.Ugrid2d.perimeter` has been added the compute the length of the
+  face perimeters.
+
 Changed
 ~~~~~~~
 
+- UGRID 2D topologies are no longer automatically forced in counterclockwise
+  orientation during initialization.
+
 Fixed
 ~~~~~
 

tests/test_connectivity.py

@@ -537,6 +537,20 @@ def test_structured_connectivity():
     assert np.array_equal(connectivity.neighbors(A, 6), [4])
 
 
+def test_area(mixed_mesh):
+    node_x = np.array([0.0, 1.0, 1.0, 2.0, 2.0])
+    node_y = np.array([0.0, 0.0, 1.0, 0.0, 1.0])
+    actual = connectivity.area(*mixed_mesh, node_x, node_y)
+    assert np.allclose(actual, [0.5, 1.0])
+
+
+def test_perimeter(mixed_mesh):
+    node_x = np.array([0.0, 1.0, 1.0, 2.0, 2.0])
+    node_y = np.array([0.0, 0.0, 1.0, 0.0, 1.0])
+    actual = connectivity.perimeter(*mixed_mesh, node_x, node_y)
+    assert np.allclose(actual, [2.0 + np.sqrt(2), 4.0])
+
+
 def test_triangulate(mixed_mesh):
     faces, fill_value = mixed_mesh
     actual_triangles, actual_faces = connectivity.triangulate_dense(faces, fill_value)

tests/test_ugrid2d.py

@@ -131,7 +131,6 @@ def test_ugrid2d_alternative_init():
 
 
 def test_ugrid2d_properties():
-    # These are defined in the base class
     grid = grid2d()
     assert grid.edge_dimension == f"{NAME}_nEdges"
     assert grid.node_dimension == f"{NAME}_nNodes"
@@ -149,6 +148,8 @@ def test_ugrid2d_properties():
     face_node_coords = grid.face_node_coordinates
     assert edge_node_coords.shape == (10, 2, 2)
     assert face_node_coords.shape == (4, 4, 2)
+    assert grid.area.shape == (grid.n_face,)
+    assert grid.perimeter.shape == (grid.n_face,)
     are_nan = np.isnan(face_node_coords)
     assert are_nan[2:, -1:, :].all()
     assert not are_nan[:, :-1, :].any()
@@ -1064,3 +1065,145 @@ def test_ugrid2d_rename_with_dataset():
         "two",
     ]
     assert sorted(dataset.coords) == ["__renamed_node_x", "__renamed_node_y"]
+
+
+class TestPeriodicGridConversion:
+    @pytest.fixture(autouse=True)
+    def setup(self):
+        self.vertices = np.array(
+            [
+                [0.0, 0.0],  # 0
+                [1.0, 0.0],  # 1
+                [2.0, 0.0],  # 2
+                [3.0, 0.0],  # 3
+                [0.0, 1.0],  # 4
+                [1.0, 1.0],  # 5
+                [2.0, 1.0],  # 6
+                [3.0, 1.0],  # 7
+                [0.0, 2.0],  # 8
+                [1.0, 2.0],  # 9
+                [2.0, 2.0],  # 10
+                [3.0, 2.0],  # 11
+            ]
+        )
+        self.faces = np.array(
+            [
+                [0, 1, 5, 4],
+                [1, 2, 6, 5],
+                [2, 3, 7, 6],
+                [4, 5, 9, 8],
+                [5, 6, 10, 9],
+                [6, 7, 11, 10],
+            ]
+        )
+        grid = xugrid.Ugrid2d(*self.vertices.T, -1, self.faces)
+        ds = xr.Dataset()
+        ds["a"] = xr.DataArray(np.arange(grid.n_node), dims=(grid.node_dimension,))
+        ds["b"] = xr.DataArray(np.arange(grid.n_edge), dims=(grid.edge_dimension,))
+        ds["c"] = xr.DataArray(np.arange(grid.n_face), dims=(grid.face_dimension,))
+        self.ds = ds
+        self.grid = grid
+
+    def test_to_periodic(self):
+        grid = self.grid.copy()
+
+        # Trigger edge node connectivity
+        _ = grid.edge_node_connectivity
+        # Convert
+        new, new_ds = grid.to_periodic(obj=self.ds)
+
+        # Absent vertices: 3, 7, 11
+        expected_vertices = self.vertices[[0, 1, 2, 4, 5, 6, 8, 9, 10]]
+        expected_faces = np.array(
+            [
+                [0, 1, 4, 3],
+                [1, 2, 5, 4],
+                [2, 0, 3, 5],
+                [3, 4, 7, 6],
+                [4, 5, 8, 7],
+                [5, 3, 6, 8],
+            ]
+        )
+        expected_edges = np.array(
+            [
+                [0, 1],
+                [0, 3],
+                [1, 2],
+                [1, 4],
+                [0, 2],
+                [2, 5],
+                [3, 4],
+                [3, 6],
+                [4, 5],
+                [4, 7],
+                [3, 5],
+                [5, 8],
+                [6, 7],
+                [7, 8],
+                [6, 8],
+            ]
+        )
+        assert np.array_equal(
+            new.face_node_connectivity,
+            expected_faces,
+        )
+        assert np.allclose(
+            new.node_coordinates,
+            expected_vertices,
+        )
+        assert np.array_equal(
+            new.edge_node_connectivity,
+            expected_edges,
+        )
+        # Remove nodes (3 & 7 & 11) and edges (6 & 13)
+        expected_a = np.arange(grid.n_node).tolist()
+        expected_a.remove(3)
+        expected_a.remove(7)
+        expected_a.remove(11)
+        expected_b = np.arange(grid.n_edge).tolist()
+        expected_b.remove(6)
+        expected_b.remove(13)
+        assert np.array_equal(new_ds["a"], expected_a)
+        assert np.array_equal(new_ds["b"], expected_b)
+        assert np.array_equal(new_ds["c"], [0, 1, 2, 3, 4, 5])
+
+        # Test whether it also works without an object provided.
+        new = grid.to_periodic()
+        assert np.array_equal(
+            new.face_node_connectivity,
+            expected_faces,
+        )
+        assert np.allclose(new.node_coordinates, expected_vertices)
+        assert np.array_equal(new.edge_node_connectivity, expected_edges)
+
+    def test_to_nonperiodic(self):
+        grid = self.grid.copy()
+        _ = grid.edge_node_connectivity  # trigger generation of edge nodes
+        periodic_grid, new_ds = grid.to_periodic(obj=self.ds)
+        back = periodic_grid.to_nonperiodic(xmax=3.0)
+
+        expected_vertices = self.vertices[[0, 1, 2, 4, 5, 6, 8, 9, 10, 3, 7, 11]]
+        expected_faces = np.array(
+            [
+                [0, 1, 4, 3],
+                [1, 2, 5, 4],
+                [2, 9, 10, 5],
+                [3, 4, 7, 6],
+                [4, 5, 8, 7],
+                [5, 10, 11, 8],
+            ]
+        )
+        back, back_ds = periodic_grid.to_nonperiodic(xmax=3.0, obj=new_ds)
+        assert np.allclose(back.node_coordinates, expected_vertices)
+        assert np.array_equal(back.face_node_connectivity, expected_faces)
+        assert back.edge_node_connectivity.shape == (17, 2)
+        assert np.array_equal(back_ds["a"], [0, 1, 2, 4, 5, 6, 8, 9, 10, 0, 4, 8])
+        assert np.array_equal(
+            back_ds["b"], [0, 1, 2, 3, 3, 4, 5, 7, 8, 9, 10, 10, 11, 12, 14, 15, 16]
+        )
+        assert np.array_equal(back_ds["c"], [0, 1, 2, 3, 4, 5])
+
+        back = periodic_grid.to_nonperiodic(xmax=3.0)
+        assert np.allclose(back.node_coordinates, expected_vertices)
+        assert np.array_equal(back.face_node_connectivity, expected_faces)
+        assert back.edge_node_connectivity.shape == (17, 2)

tests/test_ugrid_dataset.py

@@ -1117,3 +1117,55 @@ def test_fm_facenodeconnectivity_fillvalue():
 
     # xugrid 0.6.0 has -2 values in the array
     assert (uds.grid.face_node_connectivity != -2).all()
+
+
+def test_periodic_conversion():
+    vertices = np.array(
+        [
+            [0.0, 0.0],  # 0
+            [1.0, 0.0],  # 1
+            [2.0, 0.0],  # 2
+            [3.0, 0.0],  # 3
+            [0.0, 1.0],  # 4
+            [1.0, 1.0],  # 5
+            [2.0, 1.0],  # 6
+            [3.0, 1.0],  # 7
+            [0.0, 2.0],  # 8
+            [1.0, 2.0],  # 9
+            [2.0, 2.0],  # 10
+            [3.0, 2.0],  # 11
+        ]
+    )
+    faces = np.array(
+        [
+            [0, 1, 5, 4],
+            [1, 2, 6, 5],
+            [2, 3, 7, 6],
+            [4, 5, 9, 8],
+            [5, 6, 10, 9],
+            [6, 7, 11, 10],
+        ]
+    )
+    grid = xugrid.Ugrid2d(*vertices.T, -1, faces)
+    da = xr.DataArray([0, 1, 2, 3, 4, 5], dims=(grid.face_dimension,))
+    uda = xugrid.UgridDataArray(da, grid)
+    periodic = uda.ugrid.to_periodic()
+    back = periodic.ugrid.to_nonperiodic(xmax=3.0)
+    assert isinstance(periodic, xugrid.UgridDataArray)
+    assert isinstance(back, xugrid.UgridDataArray)
+    back_grid = back.ugrid.grid
+    assert back_grid.n_face == grid.n_face
+    assert back_grid.n_edge == grid.n_edge
+    assert back_grid.n_node == grid.n_node
+
+    # Also test a multi-topology dataset The 1D grid should be skipped: it
+    # doesn't implement anything for these conversions, but should simply be
+    # added as-is to the result.
+    uds = ugrid1d_ds()
+    uds["a2d"] = uda
+    periodic_ds = uds.ugrid.to_periodic()
+    back_ds = periodic_ds.ugrid.to_nonperiodic(xmax=3.0)
+    assert isinstance(periodic_ds, xugrid.UgridDataset)
+    assert isinstance(back_ds, xugrid.UgridDataset)
+    assert "a1d" in back_ds
+    assert "a2d" in back_ds

xugrid/core/dataarray_accessor.py

@@ -230,6 +230,36 @@ def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]) -> xr.DataArray
         )
         return self._raster(x, y, index)
 
+    def to_periodic(self):
+        """
+        Convert this grid to a periodic grid, where the rightmost boundary
+        shares its nodes with the leftmost boundary.
+
+        Returns
+        -------
+        periodic: UgridDataArray
+        """
+        grid, obj = self.grid.to_periodic(obj=self.obj)
+        return UgridDataArray(obj, grid)
+
+    def to_nonperiodic(self, xmax: float):
+        """
+        Convert this grid from a periodic grid (where the rightmost boundary shares its
+        nodes with the leftmost boundary) to an aperiodic grid, where the leftmost nodes
+        are separate from the rightmost nodes.
+
+        Parameters
+        ----------
+        xmax: float
+            The x-value of the newly created rightmost boundary nodes.
+
+        Returns
+        -------
+        nonperiodic: UgridDataArray
+        """
+        grid, obj = self.grid.to_nonperiodic(xmax=xmax, obj=self.obj)
+        return UgridDataArray(obj, grid)
+
     def intersect_line(
         self, start: Sequence[float], end: Sequence[float]
     ) -> xr.DataArray:

xugrid/core/dataset_accessor.py

@@ -292,6 +292,44 @@ def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]) -> xr.Dataset:
             datasets.append(self._raster(xx, yy, index))
         return xr.merge(datasets)
 
+    def to_periodic(self):
+        """
+        Converts every grid to a periodic grid, where the rightmost boundary
+        shares its nodes with the leftmost boundary.
+
+        Returns
+        -------
+        periodic: UgridDataset
+        """
+        grids = []
+        result = self.obj
+        for grid in self.grids:
+            new_grid, result = grid.to_periodic(obj=result)
+            grids.append(new_grid)
+        return UgridDataset(result, grids)
+
+    def to_nonperiodic(self, xmax: float):
+        """
+        Convert the grid from a periodic grid (where the rightmost boundary shares its
+        nodes with the leftmost boundary) to an aperiodic grid, where the leftmost nodes
+        are separate from the rightmost nodes.
+
+        Parameters
+        ----------
+        xmax: float
+            The x-value of the newly created rightmost boundary nodes.
+
+        Returns
+        -------
+        nonperiodic: UgridDataset
+        """
+        grids = []
+        result = self.obj
+        for grid in self.grids:
+            new_grid, result = grid.to_nonperiodic(xmax=xmax, obj=result)
+            grids.append(new_grid)
+        return UgridDataset(result, grids)
+
     def intersect_line(
         self, start: Sequence[float], end: Sequence[float]
     ) -> xr.Dataset: