From 62e0a9e8f01d23ba581fc2764e5d21588e9add7b Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Mon, 31 Jul 2023 21:28:47 +0200 Subject: [PATCH] Add validate_edge_node_connectivity. Fixes #129 --- docs/api.rst | 2 ++ docs/changelog.rst | 11 ++++++++++- tests/test_connectivity.py | 24 +++++++++++++++++++++++ tests/test_ugrid2d.py | 26 +++++++++--------------- xugrid/ugrid/connectivity.py | 33 ++++++++++++++++++++++++++++++- xugrid/ugrid/ugrid2d.py | 38 ++++++++++++++++++++++++++++-------- 6 files changed, 107 insertions(+), 27 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index ce4c7c73a..3d63140be 100644 --- a/docs/api.rst +++ b/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 diff --git a/docs/changelog.rst b/docs/changelog.rst index 5492b600b..b900de8f7 100644 --- a/docs/changelog.rst +++ b/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 `_. - Bug where error was thrown in the RelativeOverlapRegridder upon flipping the y coordinate. - + [0.6.1] 2023-07-07 ------------------ diff --git a/tests/test_connectivity.py b/tests/test_connectivity.py index 1e51a7dd8..f28d63eb9 100644 --- a/tests/test_connectivity.py +++ b/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( [ diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py index 6ffd82d5c..43f7ae6b2 100644 --- a/tests/test_ugrid2d.py +++ b/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( diff --git a/xugrid/ugrid/connectivity.py b/xugrid/ugrid/connectivity.py index d1e8af5c6..43ba02e06 100644 --- a/xugrid/ugrid/connectivity.py +++ b/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, diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py index e7fa8a4a9..7af3978d5 100644 --- a/xugrid/ugrid/ugrid2d.py +++ b/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],