Add validate_edge_node_connectivity. Fixes #129

docs/api.rst

@@ -254,6 +254,8 @@ UGRID2D Topology
     Ugrid2d.edge_node_connectivity
     Ugrid2d.face_edge_connectivity
     Ugrid2d.face_face_connectivity
+
+    Ugrid2d.validate_edge_node_connectivity
 
     Ugrid2d.exterior_edges
     Ugrid2d.exterior_faces

docs/changelog.rst

@@ -17,6 +17,15 @@ Added
 - Added :func:`xugrid.polygonize` to create vector polygons for all connected
   regions of a Ugrid2d topology sharing a common value. The result is a
   geopandas GeoDataFrame.
+- :meth:`xugrid.Ugrid2d.validate_edge_node_connectivity` has been added to
+  validate edge_node_connectivity by comparing with the face_node_connectivity.
+  The result can be used to define a valid subselection.
+
+Changed
+~~~~~~~
+
+- Initializing a Ugrid2d topology with an invalid edge_node_connectivity will
+  no longer raise an error.
 
 [0.6.2] 2023-07-26
 ------------------
@@ -30,7 +39,7 @@ Fixed
   See `#119 <https://github.com/Deltares/xugrid/issues/99>`_.
 - Bug where error was thrown in the RelativeOverlapRegridder upon 
   flipping the y coordinate.
-
+  
 
 [0.6.1] 2023-07-07
 ------------------

tests/test_connectivity.py

@@ -354,6 +354,30 @@ def test_edge_connectivity(mixed_mesh):
     assert np.array_equal(face_edges, expected_face_edges)
 
 
+def test_validate_edge_connectivity(mixed_mesh):
+    faces, fill_value = mixed_mesh
+    edges = np.array([[0, 1]])
+    with pytest.raises(ValueError, match="face_node_connectivity defines 6 edges"):
+        connectivity.validate_edge_node_connectivity(faces, fill_value, edges)
+
+    edges = np.array(
+        [
+            [0, 1],  # T
+            [0, 1],  # F
+            [1, 0],  # F
+            [0, 2],  # T
+            [1, 2],  # T
+            [1, 3],  # T
+            [2, 4],  # T
+            [3, 4],  # T
+            [0, 4],  # F
+        ]
+    )
+    actual = connectivity.validate_edge_node_connectivity(faces, fill_value, edges)
+    expected = np.array([True, False, False, True, True, True, True, True, False])
+    assert np.array_equal(actual, expected)
+
+
 def test_node_node_connectivity():
     edge_nodes = np.array(
         [

tests/test_ugrid2d.py

@@ -1,4 +1,3 @@
-import re
 from typing import NamedTuple
 
 import geopandas as gpd
@@ -108,22 +107,6 @@ def test_ugrid2d_init():
     assert grid._face_edge_connectivity is None
 
 
-def test_ugrid2d_init__invalid_edge_nodes():
-    edge_nodes = np.vstack((EDGE_NODES, [6, 7]))
-    msg = re.escape(
-        "edge_node_connectivity is invalid, the following edges "
-        "are not associated with any face: [10]"
-    )
-    with pytest.raises(ValueError, match=msg):
-        xugrid.Ugrid2d(
-            node_x=VERTICES[:, 0],
-            node_y=VERTICES[:, 1],
-            fill_value=-1,
-            face_node_connectivity=FACES,
-            edge_node_connectivity=edge_nodes,
-        )
-
-
 def test_ugrid2d_alternative_init():
     custom_attrs = {
         "node_dimension": "nNetNode",
@@ -172,6 +155,15 @@ def test_ugrid2d_properties():
     assert isinstance(grid.attrs, dict)
 
 
+def test_validate_edge_node_connectivity():
+    # Full test at test_connectivity
+    grid = grid2d()
+    valid = grid.validate_edge_node_connectivity()
+    assert isinstance(valid, np.ndarray)
+    assert valid.size == grid.n_edge
+    assert valid.all()
+
+
 def test_ugrid2d_edge_bounds():
     grid = grid2d()
     expected = np.array(

xugrid/ugrid/connectivity.py

@@ -2,6 +2,7 @@
 
 import numba as nb
 import numpy as np
+import pandas as pd
 from scipy import sparse
 
 from xugrid.constants import BoolArray, FloatArray, IntArray, IntDType, SparseMatrix
@@ -403,7 +404,9 @@ def edge_connectivity(
         unique, index = np.unique(np.sort(prior, axis=1), axis=0, return_index=True)
         # Check whether everything looks okay:
         if not np.array_equal(unique, edge_node_connectivity):
-            raise ValueError("Invalid edge_node_connectivity")
+            raise ValueError(
+                "Invalid edge_node_connectivity. Run .validate_edge_node_connectivity()."
+            )
         inverse_indices = index[inverse_indices]
         edge_node_connectivity = prior
 
@@ -414,6 +417,34 @@ def edge_connectivity(
     return edge_node_connectivity, face_edge_connectivity
 
 
+def validate_edge_node_connectivity(
+    face_node_connectivity: IntArray,
+    fill_value: int,
+    edge_node_connectivity: IntArray,
+) -> BoolArray:
+    """
+    * Is the edge defined by the face_node_connectivity?
+    * Are the edges unique?
+    """
+    new, _ = edge_connectivity(face_node_connectivity, fill_value)
+    old = np.sort(edge_node_connectivity, axis=1)
+
+    new1d = new.astype(np.int32).view(np.int64).ravel()
+    old1d = old.astype(np.int32).view(np.int64).ravel()
+
+    s = pd.Series(old1d)
+    n_edge = len(new1d)
+    n_old = s.nunique()
+    if n_old < n_edge:
+        raise ValueError(
+            f"face_node_connectivity defines {n_edge} edges, but "
+            f"edge_node_connectivity defines only {n_old} edges."
+        )
+
+    valid = np.isin(old1d, new1d) & (~s.duplicated().to_numpy())
+    return valid
+
+
 def face_face_connectivity(
     edge_face_connectivity: IntArray,
     fill_value: int,

xugrid/ugrid/ugrid2d.py

@@ -442,14 +442,6 @@ def edge_node_connectivity(self) -> IntArray:
 
     @edge_node_connectivity.setter
     def edge_node_connectivity(self, value):
-        if value is not None:
-            associated = np.isin(value, self.face_node_connectivity)
-            associated = associated[:, 0] & associated[:, 1]
-            if not associated.all():
-                raise ValueError(
-                    "edge_node_connectivity is invalid, the following edges are not "
-                    f"associated with any face: {np.flatnonzero(~associated)}"
-                )
         self._edge_node_connectivity = value
 
     @property
@@ -782,6 +774,36 @@ def celltree(self):
             )
         return self._celltree
 
+    def validate_edge_node_connectivity(self):
+        """
+        Mark valid edges, by comparing face_node_connectivity and
+        edge_node_connectivity. Edges that are not part of a face, as well as
+        duplicate edges are marked ``False``.
+
+        An error is raised if the face_node_connectivity defines more unique
+        edges than the edge_node_connectivity.
+
+        Returns
+        -------
+        valid: np.ndarray of bool
+            Marks for every edge whether it is valid.
+
+        Examples
+        --------
+
+        To purge invalid edges and associated data from a dataset that contains
+        un-associated or duplicate edges:
+
+        >>> uds = xugrid.open_dataset("example.nc")
+        >>> valid = uds.ugrid.grid.validate_edge_node_connectivity()
+        >>> purged = uds.isel({grid.edge_dimension: valid})
+        """
+        return connectivity.validate_edge_node_connectivity(
+            self.face_node_connectivity,
+            self.fill_value,
+            self.edge_node_connectivity,
+        )
+
     def assign_face_coords(
         self,
         obj: Union[xr.DataArray, xr.Dataset],