Add test for crossing antimeridian and fix cartesian barycentric coordinate calc

Author: fluidnumericsJoe

parcels/_datasets/unstructured/generic.py

@@ -208,7 +208,109 @@ def _fesom2_square_delaunay_uniform_z_coordinate():
     return ux.UxDataset({"U": u, "V": v, "W": w, "p": p}, uxgrid=uxgrid)
 
 
+def _fesom2_square_delaunay_antimeridian():
+    """
+    Delaunay grid that crosses the antimeridian with uniform z-coordinate, mimicking a FESOM2 dataset.
+    This dataset consists of a square domain with closed boundaries, where the grid is generated using Delaunay triangulation.
+    The bottom topography is flat and uniform, and the vertical grid spacing is constant with 10 layers spanning [0,1000.0]
+    The lateral velocity field components are non-zero constant, and the vertical velocity component is zero.
+    The pressure field is constant.
+    All fields are placed on location consistent with FESOM2 variable placement conventions
+    """
+    lon, lat = np.meshgrid(
+        np.linspace(-210.0, -150.0, Nx, dtype=np.float32), np.linspace(0, 60.0, Nx, dtype=np.float32)
+    )
+    # wrap longitude from [-180,180]
+    lon = np.where(lon < -180, lon + 360, lon)
+    lon_flat = lon.ravel()
+    lat_flat = lat.ravel()
+    zf = np.linspace(0.0, 1000.0, 10, endpoint=True, dtype=np.float32)  # Vertical element faces
+    zc = 0.5 * (zf[:-1] + zf[1:])  # Vertical element centers
+    nz = zf.size
+    nz1 = zc.size
+
+    # mask any point on one of the boundaries
+    mask = (
+        np.isclose(lon_flat, 0.0) | np.isclose(lon_flat, 60.0) | np.isclose(lat_flat, 0.0) | np.isclose(lat_flat, 60.0)
+    )
+
+    boundary_points = np.flatnonzero(mask)
+
+    uxgrid = ux.Grid.from_points(
+        (lon_flat, lat_flat),
+        method="regional_delaunay",
+        boundary_points=boundary_points,
+    )
+    uxgrid.attrs["Conventions"] = "UGRID-1.0"
+
+    # Define arrays U (zonal), V (meridional) and P (sea surface height)
+    U = np.ones(
+        (T, nz1, uxgrid.n_face), dtype=np.float64
+    )  # Lateral velocity is on the element centers and face centers
+    V = np.ones(
+        (T, nz1, uxgrid.n_face), dtype=np.float64
+    )  # Lateral velocity is on the element centers and face centers
+    W = np.zeros(
+        (T, nz, uxgrid.n_node), dtype=np.float64
+    )  # Vertical velocity is on the element faces and face vertices
+    P = np.ones((T, nz1, uxgrid.n_node), dtype=np.float64)  # Pressure is on the element centers and face vertices
+
+    u = ux.UxDataArray(
+        data=U,
+        name="U",
+        uxgrid=uxgrid,
+        dims=["time", "nz1", "n_face"],
+        coords=dict(
+            time=(["time"], TIME),
+            nz1=(["nz1"], zc),
+        ),
+        attrs=dict(
+            description="zonal velocity", units="m/s", location="face", mesh="delaunay", Conventions="UGRID-1.0"
+        ),
+    )
+    v = ux.UxDataArray(
+        data=V,
+        name="V",
+        uxgrid=uxgrid,
+        dims=["time", "nz1", "n_face"],
+        coords=dict(
+            time=(["time"], TIME),
+            nz1=(["nz1"], zc),
+        ),
+        attrs=dict(
+            description="meridional velocity", units="m/s", location="face", mesh="delaunay", Conventions="UGRID-1.0"
+        ),
+    )
+    w = ux.UxDataArray(
+        data=W,
+        name="w",
+        uxgrid=uxgrid,
+        dims=["time", "nz", "n_node"],
+        coords=dict(
+            time=(["time"], TIME),
+            nz=(["nz"], zf),
+        ),
+        attrs=dict(
+            description="vertical velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0"
+        ),
+    )
+    p = ux.UxDataArray(
+        data=P,
+        name="p",
+        uxgrid=uxgrid,
+        dims=["time", "nz1", "n_node"],
+        coords=dict(
+            time=(["time"], TIME),
+            nz1=(["nz1"], zc),
+        ),
+        attrs=dict(description="pressure", units="N/m^2", location="node", mesh="delaunay", Conventions="UGRID-1.0"),
+    )
+
+    return ux.UxDataset({"U": u, "V": v, "W": w, "p": p}, uxgrid=uxgrid)
+
+
 datasets = {
     "stommel_gyre_delaunay": _stommel_gyre_delaunay(),
     "fesom2_square_delaunay_uniform_z_coordinate": _fesom2_square_delaunay_uniform_z_coordinate(),
+    "fesom2_square_delaunay_antimeridian": _fesom2_square_delaunay_antimeridian(),
 }

parcels/uxgrid.py

@@ -74,8 +74,8 @@ def try_face(fid):
             if (bcoords >= 0).all() and (bcoords <= 1).all() and err < tol:
                 return bcoords
             else:
-                bcoords, err = self._get_barycentric_coordinates_cartesian(y, x, fid)
-                if (bcoords >= 0).all() and (bcoords <= 1).all() and err < tol:
+                bcoords = self._get_barycentric_coordinates_cartesian(y, x, fid)
+                if (bcoords >= 0).all() and (bcoords <= 1).all():
                     return bcoords
 
             return None
@@ -134,27 +134,16 @@ def _get_barycentric_coordinates_cartesian(self, y, x, fi):
         # Second attempt to find barycentric coordinates using cartesian coordinates
         nodes = np.stack(
             (
-                self._source_grid.node_x[node_ids].values(),
-                self._source_grid.node_y[node_ids].values(),
-                self._source_grid.node_z[node_ids].values(),
+                self.uxgrid.node_x[node_ids].values,
+                self.uxgrid.node_y[node_ids].values,
+                self.uxgrid.node_z[node_ids].values,
             ),
             axis=-1,
         )
 
         bcoord = np.asarray(_barycentric_coordinates_cartesian(nodes, cart_coord))
-        proj_uv = np.dot(bcoord, nodes)
-        err = np.linalg.norm(proj_uv - coord)
-        face_center = np.stack(
-            (
-                self._source_grid.face_x[fi].values(),
-                self._source_grid.face_y[fi].values(),
-                self._source_grid.face_z[fi].values(),
-            ),
-            axis=-1,
-        )
-        # Compute and remove the local projection error
-        err -= np.abs(_local_projection_error(nodes, face_center))
-        return bcoord, err
+
+        return bcoord
 
 
 def _barycentric_coordinates_cartesian(nodes, point, min_area=1e-8):
@@ -195,23 +184,6 @@ def _barycentric_coordinates_cartesian(nodes, point, min_area=1e-8):
     return barycentric_coords
 
 
-def _local_projection_error(nodes, point):
-    """
-    Computes the size of the local projection error that arises from
-    assuming planar faces. Effectively, a planar face on a spherical
-    manifold is local linearization of the spherical coordinate
-    transformation. Since query points and nodes are converted to
-    cartesian coordinates using the full spherical coordinate transformation,
-    the local projection error will likely be non-zero but related to the discretiztaion.
-    """
-    a = nodes[1] - nodes[0]
-    b = nodes[2] - nodes[0]
-    normal = np.cross(a, b)
-    normal /= np.linalg.norm(normal)
-    d = point - nodes[0]
-    return abs(np.dot(d, normal))
-
-
 def _triangle_area_cartesian(A, B, C):
     """Compute the area of a triangle given by three points."""
     d1 = B - A

tests/v4/test_uxarray_fieldset.py

@@ -115,10 +115,26 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval():
         V=Field(name="V", data=ds.V, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseConstantFace),
         W=Field(name="W", data=ds.W, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseLinearNode),
     )
-    P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseConstantFace)
+    P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseLinearNode)
     fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W])
 
     assert fieldset.U.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0
     assert fieldset.V.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0
     assert fieldset.W.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 0.0
     assert fieldset.p.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0
+
+
+def test_fesom2_square_delaunay_antimeridian_eval():
+    """
+    Test the evaluation of a fieldset with a FESOM2 square Delaunay grid that crosses the antimeridian.
+    Ensures that the fieldset can be created and evaluated correctly.
+    Since the underlying data is constant, we can check that the values are as expected.
+    """
+    ds = datasets_unstructured["fesom2_square_delaunay_antimeridian"]
+    P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseLinearNode)
+    fieldset = FieldSet([P])
+
+    assert fieldset.p.eval(time=ds.time[0].values, z=1.0, y=30.0, x=-170.0, applyConversion=False) == 1.0
+    assert fieldset.p.eval(time=ds.time[0].values, z=1.0, y=30.0, x=-180.0, applyConversion=False) == 1.0
+    assert fieldset.p.eval(time=ds.time[0].values, z=1.0, y=30.0, x=180.0, applyConversion=False) == 1.0
+    assert fieldset.p.eval(time=ds.time[0].values, z=1.0, y=30.0, x=170.0, applyConversion=False) == 1.0