Merge pull request #313 from Deltares/replace-numpy-cross

cross2d, xarray dataset dims

tests/test_ugrid2d.py

@@ -1158,7 +1158,9 @@ def test_from_structured_intervals1d(xflip: bool, yflip: bool):
 
     # Make sure the orientation is still ccw if x is decreasing, y is increasing.
     dxy = np.diff(grid.face_node_coordinates, axis=1)
-    clockwise = (np.cross(dxy[:, :-1], dxy[:, 1:])).sum(axis=1) < 0
+    clockwise = (
+        xugrid.ugrid.connectivity.cross2d(dxy[:, :-1], dxy[:, 1:]).sum(axis=1) < 0
+    )
     assert not clockwise.any()
 
 

tests/test_ugrid_dataset.py

@@ -235,7 +235,7 @@ def test_intersect_line(self):
         actual = self.uda.ugrid.intersect_line(start=p0, end=p1)
         sqrt2 = np.sqrt(2.0)
         assert isinstance(actual, xr.DataArray)
-        assert actual.dims == ("mesh2d_nFaces",)
+        assert set(actual.dims) == {"mesh2d_nFaces"}
         assert np.allclose(actual["mesh2d_x"], [0.5, 1.25])
         assert np.allclose(actual["mesh2d_y"], [0.5, 1.25])
         assert np.allclose(actual["mesh2d_s"], [0.5 * sqrt2, 1.25 * sqrt2])
@@ -250,7 +250,7 @@ def test_intersect_linestring(self):
         )
         actual = self.uda.ugrid.intersect_linestring(linestring)
         assert isinstance(actual, xr.DataArray)
-        assert actual.dims == ("mesh2d_nFaces",)
+        assert set(actual.dims) == {"mesh2d_nFaces"}
         assert np.allclose(actual["mesh2d_x"], [0.75, 1.25, 1.5, 1.5])
         assert np.allclose(actual["mesh2d_y"], [0.5, 0.5, 0.75, 1.25])
         assert np.allclose(actual["mesh2d_s"], [0.25, 0.75, 1.25, 1.75])

xugrid/ugrid/connectivity.py

@@ -17,6 +17,11 @@
 )
 
 
+def cross2d(arr1: np.ndarray, arr2: np.ndarray) -> np.ndarray:
+    """Cross product for arrays with 2D coordinates."""
+    return arr1[..., 0] * arr2[..., 1] - arr1[..., 1] * arr2[..., 0]
+
+
 def argsort_rows(array: np.ndarray) -> IntArray:
     if array.ndim != 2:
         raise ValueError(f"Array is not 2D, but has shape: {array.shape}")
@@ -388,7 +393,7 @@ def counterclockwise(face_node_connectivity: IntArray, nodes: FloatArray) -> Int
     closed, _ = close_polygons(face_node_connectivity)
     p = nodes[closed]
     dxy = np.diff(p, axis=1)
-    reverse = (np.cross(dxy[:, :-1], dxy[:, 1:])).sum(axis=1) < 0
+    reverse = cross2d(dxy[:, :-1], dxy[:, 1:]).sum(axis=1) < 0
     ccw = face_node_connectivity.copy()
     if reverse.any():
         ccw[reverse] = reverse_orientation(face_node_connectivity[reverse])
@@ -566,7 +571,7 @@ def area_from_coordinates(
     xy0 = coordinates[:, 0]
     a = coordinates[:, :-1] - xy0[:, np.newaxis]
     b = coordinates[:, 1:] - xy0[:, np.newaxis]
-    determinant = np.cross(a, b)
+    determinant = cross2d(a, b)
     return 0.5 * abs(determinant.sum(axis=1))
 
 
@@ -604,7 +609,7 @@ def centroids(
         a = coordinates[:, :-1] - xy0[:, np.newaxis]
         b = coordinates[:, 1:] - xy0[:, np.newaxis]
         c = a + b
-        determinant = np.cross(a, b)
+        determinant = cross2d(a, b)
         area_weight = 1.0 / (3.0 * determinant.sum(axis=1))
         centroid_coordinates[:, 0] = area_weight * (c[..., 0] * determinant).sum(axis=1)
         centroid_coordinates[:, 1] = area_weight * (c[..., 1] * determinant).sum(axis=1)