From 9184f3940786a2eecad0b66a2d240ca0374825cc Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 15 Aug 2023 14:45:53 +0200
Subject: [PATCH 1/6] Add perimeter property, to_aperiodic method

---
 xugrid/core/dataarray_accessor.py | 14 +++++
 xugrid/core/dataset_accessor.py   | 17 ++++++
 xugrid/ugrid/connectivity.py      | 16 ++++++
 xugrid/ugrid/ugrid2d.py           | 91 +++++++++++++++++++++++++++++++
 4 files changed, 138 insertions(+)

diff --git a/xugrid/core/dataarray_accessor.py b/xugrid/core/dataarray_accessor.py
index db64f69d1..dd7063c69 100644
--- a/xugrid/core/dataarray_accessor.py
+++ b/xugrid/core/dataarray_accessor.py
@@ -226,6 +226,20 @@ def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]) -> xr.DataArray
         )
         return self._raster(x, y, index)
 
+    def to_aperiodic(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.
+        """
+        grid, obj = self.grid.to_aperiodic(xmax=xmax, obj=self.obj)
+        return UgridDataArray(obj, grid)
+
     @property
     def crs(self):
         """
diff --git a/xugrid/core/dataset_accessor.py b/xugrid/core/dataset_accessor.py
index 9c3631014..7ea8132c4 100644
--- a/xugrid/core/dataset_accessor.py
+++ b/xugrid/core/dataset_accessor.py
@@ -281,6 +281,23 @@ 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_aperiodic(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.
+        """
+        grids = []
+        for grid in self.grids:
+            grid, obj = self.grid.to_aperiodic(xmax=xmax, obj=self.obj)
+            grids.append(grid)
+        return UgridDataset(obj, grids)
+
     def to_dataset(self, optional_attributes: bool = False):
         """
         Converts this UgridDataArray or UgridDataset into a standard
diff --git a/xugrid/ugrid/connectivity.py b/xugrid/ugrid/connectivity.py
index ca7b0b783..aaed9a668 100644
--- a/xugrid/ugrid/connectivity.py
+++ b/xugrid/ugrid/connectivity.py
@@ -515,6 +515,22 @@ def structured_connectivity(active: IntArray) -> AdjacencyMatrix:
     return AdjacencyMatrix(A.indices, A.indptr, A.nnz, n, m)
 
 
+def perimeter(
+    face_node_connectivity: IntArray,
+    fill_value: int,
+    node_x: FloatArray,
+    node_y: FloatArray,
+):
+    nodes = np.column_stack([node_x, node_y])
+    closed, _ = close_polygons(face_node_connectivity, fill_value)
+    coordinates = nodes[closed]
+    # Shift coordinates to avoid precision loss
+    coordinates[..., 0] -= node_x.mean()
+    coordinates[..., 1] -= node_y.mean()
+    dxy = np.diff(coordinates, axis=1)
+    return np.linalg.norm(dxy, axis=-1).sum(axis=1)
+
+
 def area(
     face_node_connectivity: IntArray,
     fill_value: int,
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index bbd9f557d..c28e33372 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -122,6 +122,8 @@ def __init__(
         self._meshkernel = None
         # Celltree
         self._celltree = None
+        # Perimeter
+        self._perimeter = None
         # Area
         self._area = None
         # Centroids
@@ -163,6 +165,8 @@ def _clear_geometry_properties(self):
         self._meshkernel = None
         # Celltree
         self._celltree = None
+        # Perimeter
+        self._perimeter = None
         # Area
         self._area = None
         # Centroids
@@ -518,6 +522,17 @@ def area(self) -> FloatArray:
             )
         return self._area
 
+    @property
+    def perimeter(self) -> FloatArray:
+        if self._perimeter is None:
+            self._perimeter = connectivity.perimeter(
+                self.face_node_connectivity,
+                self.fill_value,
+                self.node_x,
+                self.node_y,
+            )
+        return self._perimeter
+
     @property
     def face_bounds(self):
         """
@@ -1410,6 +1425,82 @@ def merge_partitions(grids: Sequence["Ugrid2d"]) -> "Ugrid2d":
         )
         return merged_grid, indexes
 
+    def to_aperiodic(self, xmax: float, obj=None):
+        """
+        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.
+        obj: xr.DataArray or xr.Dataset
+
+        Returns
+        -------
+        aperiodic_grid: Ugrid2d
+        aligned: xr.DataArray or xr.Dataset
+        """
+        xmin, _, xmax, _ = self.bounds
+        half_domain = 0.5 * (xmax - xmin)
+
+        x = self.face_node_coordinates[..., 0]
+        is_periodic = (np.nanmax(x, axis=1)[:, np.newaxis] - x) > half_domain
+        periodic_nodes = self.face_node_connectivity[is_periodic]
+
+        uniques, new_nodes = np.unique(periodic_nodes, return_inverse=True)
+        new_nodes += self.n_node
+        new_x = np.full(uniques.size, xmax)
+        new_y = self.node_y[uniques]
+        new_faces = self.face_node_connectivity.copy()
+        new_faces[is_periodic] = new_nodes
+
+        # edge_node_connectivity must be rederived!
+        new = Ugrid2d(
+            node_x=np.concatenate((self.node_x, new_x)),
+            node_y=np.concatenate((self.node_y, new_y)),
+            face_node_connectivity=new_faces,
+            fill_value=self.fill_value,
+            name=self.name,
+            edge_node_connectivity=None,
+            indexes=self._indexes,
+            projected=self.projected,
+            crs=self.crs,
+            attrs=self.attrs,
+        )
+
+        edge_index = None
+        if self._edge_node_connectivity is not None:
+            # Use a casting trick so we can use numpy isin.
+            edges = (
+                np.sort(self.edge_node_connectivity.astype(np.int32), axis=1)
+                .view(np.int64)
+                .ravel()
+            )
+            sorter = np.argsort(edges)
+            new_edges = (
+                new.edge_node_connectivity.astype(np.int32).view(np.int64).ravel()
+            )
+            is_old = np.isin(new_edges, edges)
+            edge_index = np.full(new.n_edge, -1)
+            edge_index[is_old] = np.searchsorted(
+                edges, new_edges[is_old], sorter=sorter
+            )
+
+        if obj is not None:
+            indexes = {
+                self.face_dimension: pd.RangeIndex(0, self.n_face),
+                self.node_dimension: pd.Index(
+                    np.concatenate((np.arange(self.n_node), new_nodes))
+                ),
+            }
+            if edge_index is not None:
+                indexes[self.edge_dimension] = pd.Index(edge_index)
+            return new, obj.isel(**indexes)
+        else:
+            return new
+
     def triangulate(self):
         """
         Triangulate this UGRID2D topology, breaks more complex polygons down

From dad7912b4d8eea974e7c70a9d89216fc1fb720a2 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 16 Aug 2023 14:53:37 +0200
Subject: [PATCH 2/6] Change name to nonperiodic, fix edge_index. First start
 with to_periodic.

---
 xugrid/core/dataarray_accessor.py |   4 +-
 xugrid/core/dataset_accessor.py   |   4 +-
 xugrid/ugrid/ugrid2d.py           | 103 ++++++++++++++++++++++++++----
 3 files changed, 95 insertions(+), 16 deletions(-)

diff --git a/xugrid/core/dataarray_accessor.py b/xugrid/core/dataarray_accessor.py
index dd7063c69..852b04f08 100644
--- a/xugrid/core/dataarray_accessor.py
+++ b/xugrid/core/dataarray_accessor.py
@@ -226,7 +226,7 @@ def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]) -> xr.DataArray
         )
         return self._raster(x, y, index)
 
-    def to_aperiodic(self, xmax: float):
+    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
@@ -237,7 +237,7 @@ def to_aperiodic(self, xmax: float):
         xmax: float
             The x-value of the newly created rightmost boundary nodes.
         """
-        grid, obj = self.grid.to_aperiodic(xmax=xmax, obj=self.obj)
+        grid, obj = self.grid.to_nonperiodic(xmax=xmax, obj=self.obj)
         return UgridDataArray(obj, grid)
 
     @property
diff --git a/xugrid/core/dataset_accessor.py b/xugrid/core/dataset_accessor.py
index 7ea8132c4..bc5a2bb92 100644
--- a/xugrid/core/dataset_accessor.py
+++ b/xugrid/core/dataset_accessor.py
@@ -281,7 +281,7 @@ 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_aperiodic(self, xmax: float):
+    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
@@ -294,7 +294,7 @@ def to_aperiodic(self, xmax: float):
         """
         grids = []
         for grid in self.grids:
-            grid, obj = self.grid.to_aperiodic(xmax=xmax, obj=self.obj)
+            grid, obj = self.grid.to_nonperiodic(xmax=xmax, obj=self.obj)
             grids.append(grid)
         return UgridDataset(obj, grids)
 
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index c28e33372..d7372626b 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -1425,7 +1425,65 @@ def merge_partitions(grids: Sequence["Ugrid2d"]) -> "Ugrid2d":
         )
         return merged_grid, indexes
 
-    def to_aperiodic(self, xmax: float, obj=None):
+    def to_periodic(self, obj=None):
+        """
+        Convert this grid to a periodic grid, where the rightmost nodes are
+        equal to the leftmost nodes. Note: for this to work, the y-coordinates
+        on the left boundary must match those on the right boundary exactly.
+
+        Returns
+        -------
+        periodic_grid: Ugrid2d
+        aligned: xr.DataArray or xr.Dataset
+        """
+        xmin, _, xmax, _ = self.bounds
+        coordinates = self.node_coordinates
+        is_right = np.isclose(coordinates[:, 0], xmax)
+        is_left = np.isclose(coordinates[:, 0], xmin)
+
+        node_y = coordinates[:, 1]
+        if not np.allclose(np.sort(node_y[is_left]), np.sort(node_y[is_right])):
+            raise ValueError(
+                "y-coordinates of the left and right boundaries do not match"
+            )
+
+        coordinates[is_right, 0] = xmin
+        new_xy, node_index, inverse = np.unique(
+            coordinates, return_index=True, return_inverse=True
+        )
+        new_faces = inverse[self.face_node_connectivity]
+        new_edges = None
+        if self._edge_node_connectivity is not None:
+            new_edges = inverse[self.edge_node_connectivity]
+            new_edges.sort(axis=1)
+            new_edges, edge_index = np.unique(new_edges, axis=0, return_index=True)
+
+        new = Ugrid2d(
+            node_x=new_xy[:, 0],
+            node_y=new_xy[:, 1],
+            face_node_connectivity=new_faces,
+            fill_value=self.fill_value,
+            name=self.name,
+            edge_node_connectivity=None,
+            indexes=self._indexes,
+            projected=self.projected,
+            crs=self.crs,
+            attrs=self.attrs,
+        )
+
+        if obj is not None:
+            indexes = {
+                self.face_dimension: pd.RangeIndex(0, self.n_face),
+                self.node_dimension: pd.Index(node_index),
+            }
+            if edge_index is not None:
+                indexes[self.edge_dimension] = pd.Index(edge_index)
+            indexes = {k: v for k, v in indexes.items() if k in obj.dims}
+            return new, obj.isel(**indexes)
+        else:
+            return new
+
+    def to_nonperiodic(self, xmax: float, obj=None):
         """
         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
@@ -1439,7 +1497,7 @@ def to_aperiodic(self, xmax: float, obj=None):
 
         Returns
         -------
-        aperiodic_grid: Ugrid2d
+        nonperiodic_grid: Ugrid2d
         aligned: xr.DataArray or xr.Dataset
         """
         xmin, _, xmax, _ = self.bounds
@@ -1456,7 +1514,8 @@ def to_aperiodic(self, xmax: float, obj=None):
         new_faces = self.face_node_connectivity.copy()
         new_faces[is_periodic] = new_nodes
 
-        # edge_node_connectivity must be rederived!
+        # edge_node_connectivity must be rederived, since we've added a number
+        # of new edges and new nodes.
         new = Ugrid2d(
             node_x=np.concatenate((self.node_x, new_x)),
             node_y=np.concatenate((self.node_y, new_y)),
@@ -1472,21 +1531,40 @@ def to_aperiodic(self, xmax: float, obj=None):
 
         edge_index = None
         if self._edge_node_connectivity is not None:
-            # Use a casting trick so we can use numpy isin.
+            # If there is edge associated data, we need to duplicate the data
+            # of the edges. It is impossible(?) to do this on the edges
+            # directly, due to the possible presence of "symmetric" edges:
+            #     2
+            #    /|\
+            #   / | \
+            #  0__1__0
+            #
+            # (0, 1) and (1, 0) are topologically distinct, but only in the
+            # face definition. In the new grid, the 0 on the right will have
+            # become node 3, creating distinct edges.
+            #
+            # Note that any data with the edge is only stored once, which is
+            # incorrect(!), but a given for these grids and would be a problem
+            # for the simulation code producing these results.
+            #
+            # We use a casting trick to collapse two integers into one so we
+            # can use searchsorted easily.
             edges = (
-                np.sort(self.edge_node_connectivity.astype(np.int32), axis=1)
+                np.sort(self.edge_node_connectivity, axis=1)
+                .astype(np.int32)
                 .view(np.int64)
                 .ravel()
             )
-            sorter = np.argsort(edges)
+            # Create a mapping of the new nodes created above, to the original nodes.
+            # Then, find the new edges in the old using searchsorted.
+            mapping = np.concatenate((np.arange(self.n_node), uniques))
             new_edges = (
-                new.edge_node_connectivity.astype(np.int32).view(np.int64).ravel()
-            )
-            is_old = np.isin(new_edges, edges)
-            edge_index = np.full(new.n_edge, -1)
-            edge_index[is_old] = np.searchsorted(
-                edges, new_edges[is_old], sorter=sorter
+                np.sort(mapping[new.edge_node_connectivity], axis=1)
+                .astype(np.int32)
+                .view(np.int64)
+                .ravel()
             )
+            edge_index = np.searchsorted(edges, new_edges, sorter=np.argsort(edges))
 
         if obj is not None:
             indexes = {
@@ -1497,6 +1575,7 @@ def to_aperiodic(self, xmax: float, obj=None):
             }
             if edge_index is not None:
                 indexes[self.edge_dimension] = pd.Index(edge_index)
+            indexes = {k: v for k, v in indexes.items() if k in obj.dims}
             return new, obj.isel(**indexes)
         else:
             return new

From 4c6a66b20a6afabf3537a1d3ca762157986470a1 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 4 Sep 2023 13:34:29 +0200
Subject: [PATCH 3/6] Fix node indexing, iterate over grids

---
 xugrid/core/dataset_accessor.py | 4 ++--
 xugrid/ugrid/ugrid2d.py         | 5 ++---
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/xugrid/core/dataset_accessor.py b/xugrid/core/dataset_accessor.py
index bc5a2bb92..32fa4e4d3 100644
--- a/xugrid/core/dataset_accessor.py
+++ b/xugrid/core/dataset_accessor.py
@@ -294,8 +294,8 @@ def to_nonperiodic(self, xmax: float):
         """
         grids = []
         for grid in self.grids:
-            grid, obj = self.grid.to_nonperiodic(xmax=xmax, obj=self.obj)
-            grids.append(grid)
+            new_grid, obj = grid.to_nonperiodic(xmax=xmax, obj=self.obj)
+            grids.append(new_grid)
         return UgridDataset(obj, grids)
 
     def to_dataset(self, optional_attributes: bool = False):
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index d7372626b..70e63f6f9 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -1508,11 +1508,10 @@ def to_nonperiodic(self, xmax: float, obj=None):
         periodic_nodes = self.face_node_connectivity[is_periodic]
 
         uniques, new_nodes = np.unique(periodic_nodes, return_inverse=True)
-        new_nodes += self.n_node
         new_x = np.full(uniques.size, xmax)
         new_y = self.node_y[uniques]
         new_faces = self.face_node_connectivity.copy()
-        new_faces[is_periodic] = new_nodes
+        new_faces[is_periodic] = new_nodes + self.n_node
 
         # edge_node_connectivity must be rederived, since we've added a number
         # of new edges and new nodes.
@@ -1570,7 +1569,7 @@ def to_nonperiodic(self, xmax: float, obj=None):
             indexes = {
                 self.face_dimension: pd.RangeIndex(0, self.n_face),
                 self.node_dimension: pd.Index(
-                    np.concatenate((np.arange(self.n_node), new_nodes))
+                    np.concatenate((np.arange(self.n_node), uniques))
                 ),
             }
             if edge_index is not None:

From 3c94af377cb054847a764f7f92670f4921243525 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 5 Sep 2023 13:43:37 +0200
Subject: [PATCH 4/6] Do no overwrite xmax variable

---
 xugrid/ugrid/ugrid2d.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index 70e63f6f9..88b5b08c1 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -1500,9 +1500,12 @@ def to_nonperiodic(self, xmax: float, obj=None):
         nonperiodic_grid: Ugrid2d
         aligned: xr.DataArray or xr.Dataset
         """
-        xmin, _, xmax, _ = self.bounds
-        half_domain = 0.5 * (xmax - xmin)
+        xleft, _, xright, _ = self.bounds
+        half_domain = 0.5 * (xright - xleft)
 
+        # Extract all x coordinates for every face. Then identify the nodes
+        # which have a value of e.g. -180, while the max x value for the face
+        # is 180.0. These nodes should be duplicated.
         x = self.face_node_coordinates[..., 0]
         is_periodic = (np.nanmax(x, axis=1)[:, np.newaxis] - x) > half_domain
         periodic_nodes = self.face_node_connectivity[is_periodic]

From 517980bd0db2e1a1285c7539924921a4d277eaa6 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 5 Oct 2023 17:57:06 +0200
Subject: [PATCH 5/6] Add tests for to_periodic, to_nonperiodic Add tests for
 perimeter, area Do not coercing to counterclockwise during initialization:
 the current criterium fails on a "periodic" grid; added TO DO's.

Arguably, no coercing should happen during initialization at all. It should be available as a utility and in validation.
---
 tests/test_connectivity.py   |  14 ++++
 tests/test_ugrid2d.py        | 145 ++++++++++++++++++++++++++++++++++-
 xugrid/ugrid/connectivity.py |   5 +-
 xugrid/ugrid/ugrid2d.py      |  37 +++++++--
 4 files changed, 190 insertions(+), 11 deletions(-)

diff --git a/tests/test_connectivity.py b/tests/test_connectivity.py
index b9ef1e651..5065ab7dc 100644
--- a/tests/test_connectivity.py
+++ b/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)
diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py
index 7dac6544a..44ce41c7c 100644
--- a/tests/test_ugrid2d.py
+++ b/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)
diff --git a/xugrid/ugrid/connectivity.py b/xugrid/ugrid/connectivity.py
index 0100b7540..2504a49df 100644
--- a/xugrid/ugrid/connectivity.py
+++ b/xugrid/ugrid/connectivity.py
@@ -354,6 +354,7 @@ def reverse_orientation(face_node_connectivity: IntArray, fill_value: int):
 def counterclockwise(
     face_node_connectivity: IntArray, fill_value: int, nodes: FloatArray
 ) -> IntArray:
+    # TODO: in case of "periodic grids", we need to wrap around the grid.
     closed, _ = close_polygons(face_node_connectivity, fill_value)
     p = nodes[closed]
     dxy = np.diff(p, axis=1)
@@ -525,8 +526,8 @@ def perimeter(
     closed, _ = close_polygons(face_node_connectivity, fill_value)
     coordinates = nodes[closed]
     # Shift coordinates to avoid precision loss
-    coordinates[..., 0] -= node_x.mean()
-    coordinates[..., 1] -= node_y.mean()
+    xy0 = coordinates[:, 0]
+    coordinates -= xy0[:, np.newaxis]
     dxy = np.diff(coordinates, axis=1)
     return np.linalg.norm(dxy, axis=-1).sum(axis=1)
 
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index d848d4913..16e6eaba0 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -109,9 +109,13 @@ def __init__(
                 "face_node_connectivity should be an array of integers or a sparse matrix"
             )
 
-        self.face_node_connectivity = connectivity.counterclockwise(
-            face_node_connectivity, self.fill_value, self.node_coordinates
-        )
+        self.face_node_connectivity = face_node_connectivity
+
+        # TODO: do this in validation instead. While UGRID conventions demand it,
+        # where does it go wrong?
+        # self.face_node_connectivity = connectivity.counterclockwise(
+        #    face_node_connectivity, self.fill_value, self.node_coordinates
+        # )
 
         self._initialize_indexes_attrs(name, dataset, indexes, attrs)
         self._dataset = dataset
@@ -1529,16 +1533,28 @@ def to_periodic(self, obj=None):
                 "y-coordinates of the left and right boundaries do not match"
             )
 
+        # Discard the rightmost nodes. Preserve the order in the faces, and the
+        # order of the nodes.
         coordinates[is_right, 0] = xmin
-        new_xy, node_index, inverse = np.unique(
-            coordinates, return_index=True, return_inverse=True
+        _, node_index, inverse = np.unique(
+            coordinates, return_index=True, return_inverse=True, axis=0
         )
-        new_faces = inverse[self.face_node_connectivity]
+        # Create a mapping of the inverse index to the new node index.
+        new_index = connectivity.renumber(node_index)
+        new_faces = new_index[inverse[self.face_node_connectivity]]
+        # Get the selection of nodes, and keep the order.
+        node_index.sort()
+        new_xy = self.node_coordinates[node_index]
+
+        # Preserve the order of the edge_node_connectivity if it is present.
         new_edges = None
+        edge_index = None
         if self._edge_node_connectivity is not None:
             new_edges = inverse[self.edge_node_connectivity]
             new_edges.sort(axis=1)
-            new_edges, edge_index = np.unique(new_edges, axis=0, return_index=True)
+            _, edge_index = np.unique(new_edges, axis=0, return_index=True)
+            edge_index.sort()
+            new_edges = new_index[new_edges][edge_index]
 
         new = Ugrid2d(
             node_x=new_xy[:, 0],
@@ -1546,7 +1562,7 @@ def to_periodic(self, obj=None):
             face_node_connectivity=new_faces,
             fill_value=self.fill_value,
             name=self.name,
-            edge_node_connectivity=None,
+            edge_node_connectivity=new_edges,
             indexes=self._indexes,
             projected=self.projected,
             crs=self.crs,
@@ -1649,6 +1665,11 @@ def to_nonperiodic(self, xmax: float, obj=None):
                 .ravel()
             )
             edge_index = np.searchsorted(edges, new_edges, sorter=np.argsort(edges))
+            # Reshuffle to keep the original order as intact as possible; how
+            # much benefit does this actually give?
+            sorter = np.argsort(edge_index)
+            new._edge_node_connectivity = new._edge_node_connectivity[sorter]
+            edge_index = edge_index[sorter]
 
         if obj is not None:
             indexes = {

From 0d848553465783a3c78a30124235b4b45fc9558f Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 6 Oct 2023 16:33:00 +0200
Subject: [PATCH 6/6] Add tests for to_periodic, to_nonperiodic.

Also fix some docstrings, api.

Update changelog.
---
 docs/api.rst                      | 10 ++++++
 docs/changelog.rst                | 15 +++++++++
 tests/test_ugrid_dataset.py       | 52 +++++++++++++++++++++++++++++++
 xugrid/core/dataarray_accessor.py | 16 ++++++++++
 xugrid/core/dataset_accessor.py   | 27 ++++++++++++++--
 xugrid/ugrid/ugrid1d.py           |  6 ++++
 xugrid/ugrid/ugrid2d.py           | 10 ++++++
 7 files changed, 133 insertions(+), 3 deletions(-)

diff --git a/docs/api.rst b/docs/api.rst
index 6ce18b15d..a49b0ed95 100644
--- a/docs/api.rst
+++ b/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
diff --git a/docs/changelog.rst b/docs/changelog.rst
index f546332c4..60e1969f2 100644
--- a/docs/changelog.rst
+++ b/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
 ~~~~~
 
diff --git a/tests/test_ugrid_dataset.py b/tests/test_ugrid_dataset.py
index 08273975e..18fbdb17f 100644
--- a/tests/test_ugrid_dataset.py
+++ b/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
diff --git a/xugrid/core/dataarray_accessor.py b/xugrid/core/dataarray_accessor.py
index 3d60f0ac2..3427f468c 100644
--- a/xugrid/core/dataarray_accessor.py
+++ b/xugrid/core/dataarray_accessor.py
@@ -230,6 +230,18 @@ 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
@@ -240,6 +252,10 @@ def to_nonperiodic(self, xmax: float):
         ----------
         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)
diff --git a/xugrid/core/dataset_accessor.py b/xugrid/core/dataset_accessor.py
index ab5eaefc2..791e2e290 100644
--- a/xugrid/core/dataset_accessor.py
+++ b/xugrid/core/dataset_accessor.py
@@ -292,9 +292,25 @@ 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 this grid from a periodic grid (where the rightmost boundary shares its
+        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.
 
@@ -302,12 +318,17 @@ def to_nonperiodic(self, xmax: float):
         ----------
         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, obj = grid.to_nonperiodic(xmax=xmax, obj=self.obj)
+            new_grid, result = grid.to_nonperiodic(xmax=xmax, obj=result)
             grids.append(new_grid)
-        return UgridDataset(obj, grids)
+        return UgridDataset(result, grids)
 
     def intersect_line(
         self, start: Sequence[float], end: Sequence[float]
diff --git a/xugrid/ugrid/ugrid1d.py b/xugrid/ugrid/ugrid1d.py
index ac9c6427f..c56561fb8 100644
--- a/xugrid/ugrid/ugrid1d.py
+++ b/xugrid/ugrid/ugrid1d.py
@@ -529,6 +529,12 @@ def intersect_line(self, obj, start, stop):
     def intersect_linestring(self, obj, linestring):
         return obj
 
+    def to_periodic(self, obj):
+        return self, obj
+
+    def to_nonperiodic(self, xmax, obj):
+        return self, obj
+
     def topological_sort_by_dfs(self) -> IntArray:
         """
         Returns an array of vertices in topological order.
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index 16e6eaba0..26b2f2ee1 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -506,6 +506,10 @@ def centroids(self) -> FloatArray:
 
     @property
     def circumcenters(self):
+        """
+        Circumenter (x, y) of every face; only works for fully triangular
+        grids.
+        """
         if self._circumcenters is None:
             self._circumcenters = connectivity.circumcenters(
                 self.face_node_connectivity,
@@ -517,6 +521,9 @@ def circumcenters(self):
 
     @property
     def area(self) -> FloatArray:
+        """
+        Area of every face.
+        """
         if self._area is None:
             self._area = connectivity.area(
                 self.face_node_connectivity,
@@ -528,6 +535,9 @@ def area(self) -> FloatArray:
 
     @property
     def perimeter(self) -> FloatArray:
+        """
+        Perimeter length of every face.
+        """
         if self._perimeter is None:
             self._perimeter = connectivity.perimeter(
                 self.face_node_connectivity,