Merge branch 'v4-dev' into interpolation_tutorial_v4

Author: reint-fischer

src/parcels/_core/field.py

@@ -21,7 +21,7 @@
 from parcels._core.utils.string import _assert_str_and_python_varname
 from parcels._core.utils.time import TimeInterval
 from parcels._core.uxgrid import UxGrid
-from parcels._core.xgrid import XGrid, _transpose_xfield_data_to_tzyx
+from parcels._core.xgrid import XGrid, _transpose_xfield_data_to_tzyx, assert_all_field_dims_have_axis
 from parcels._python import assert_same_function_signature
 from parcels._reprs import default_repr
 from parcels._typing import VectorType
@@ -103,6 +103,7 @@ def __init__(
         _assert_compatible_combination(data, grid)
 
         if isinstance(grid, XGrid):
+            assert_all_field_dims_have_axis(data, grid.xgcm_grid)
             data = _transpose_xfield_data_to_tzyx(data, grid.xgcm_grid)
 
         self.name = name

src/parcels/_core/index_search.py

@@ -102,6 +102,23 @@ def curvilinear_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray
     )
 
     px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]])
+    # Map grid and particle longitudes to [-180,180)
+    px = ((px + 180.0) % 360.0) - 180.0
+    x = ((x + 180.0) % 360.0) - 180.0
+
+    # Create a mask for antimeridian cells
+    lon_span = px.max(axis=0) - px.min(axis=0)
+    antimeridian_cell = lon_span > 180.0
+
+    if np.any(antimeridian_cell):
+        # For any antimeridian cell ...
+        #  If particle longitude is closer to 180.0, then negative cell longitudes need to be bumped by +360
+        mask = (px < 0.0) & antimeridian_cell[np.newaxis, :] & (x[np.newaxis, :] > 0.0)
+        px[mask] += 360.0
+        #  If particle longitude is closer to -180.0, then positive cell longitudes need to be bumped by -360
+        mask = (px > 0.0) & antimeridian_cell[np.newaxis, :] & (x[np.newaxis, :] < 0.0)
+        px[mask] -= 360.0
+
     py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]])
 
     a, b = np.dot(invA, px), np.dot(invA, py)
@@ -119,7 +136,6 @@ def curvilinear_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray
             ((y - py[0]) / (py[1] - py[0]) + (y - py[3]) / (py[2] - py[3])) * 0.5,
             (x - a[0] - a[2] * eta) / (a[1] + a[3] * eta),
         )
-
     is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0)
 
     return is_in_cell, np.column_stack((xsi, eta))

src/parcels/_core/xgrid.py

@@ -48,6 +48,16 @@ def _drop_field_data(ds: xr.Dataset) -> xr.Dataset:
     return ds.drop_vars(ds.data_vars)
 
 
+def assert_all_field_dims_have_axis(da: xr.DataArray, xgcm_grid: xgcm.Grid) -> None:
+    ax_dims = [(get_axis_from_dim_name(xgcm_grid.axes, dim), dim) for dim in da.dims]
+    for dim in ax_dims:
+        if dim[0] is None:
+            raise ValueError(
+                f'Dimension "{dim[1]}" has no axis attribute. '
+                f'HINT: You may want to add an {{"axis": A}} to your DataSet["{dim[1]}"], where A is one of "X", "Y", "Z" or "T"'
+            )
+
+
 def _transpose_xfield_data_to_tzyx(da: xr.DataArray, xgcm_grid: xgcm.Grid) -> xr.DataArray:
     """
     Transpose a DataArray of any shape into a 4D array of order TZYX. Uses xgcm to determine

src/parcels/interpolators.py

@@ -175,8 +175,7 @@ def CGrid_Velocity(
         py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]])
 
     if grid._mesh == "spherical":
-        px[0] = np.where(px[0] < lon - 225, px[0] + 360, px[0])
-        px[0] = np.where(px[0] > lon + 225, px[0] - 360, px[0])
+        px = ((px + 180.0) % 360.0) - 180.0
         px[1:] = np.where(px[1:] - px[0] > 180, px[1:] - 360, px[1:])
         px[1:] = np.where(-px[1:] + px[0] > 180, px[1:] + 360, px[1:])
     c1 = i_u._geodetic_distance(
@@ -291,7 +290,10 @@ def CGrid_Velocity(
 
     # check whether the grid conversion has been applied correctly
     xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3]
-    u = np.where(np.abs((xx - lon) / lon) > 1e-4, np.nan, u)
+    dlon = xx - lon
+    if grid._mesh == "spherical":
+        dlon = ((dlon + 180.0) % 360.0) - 180.0
+    u = np.where(np.abs(dlon / lon) > 1e-4, np.nan, u)
 
     if vectorfield.W:
         data = vectorfield.W.data

tests/test_advection.py

@@ -454,20 +454,17 @@ def test_nemo_curvilinear_fieldset():
     U = parcels.Field("U", ds["U"], grid, interp_method=XLinear)
     V = parcels.Field("V", ds["V"], grid, interp_method=XLinear)
     U.units = parcels.GeographicPolar()
-    V.units = parcels.GeographicPolar()  # U and V need GoegraphicPolar for  C-Grid interpolation to work correctly
+    V.units = parcels.GeographicPolar()  # U and V need GeographicPolar for C-Grid interpolation to work correctly
     UV = parcels.VectorField("UV", U, V, vector_interp_method=CGrid_Velocity)
     fieldset = parcels.FieldSet([U, V, UV])
 
     npart = 20
     lonp = 30 * np.ones(npart)
     latp = np.linspace(-70, 88, npart)
-    runtime = np.timedelta64(12, "h")  # TODO increase to 160 days
-
-    def periodicBC(particles, fieldset):  # pragma: no cover
-        particles.dlon = np.where(particles.lon > 180, particles.dlon - 360, particles.dlon)
+    runtime = np.timedelta64(160, "D")
 
     pset = parcels.ParticleSet(fieldset, lon=lonp, lat=latp)
-    pset.execute([AdvectionEE, periodicBC], runtime=runtime, dt=np.timedelta64(6, "h"))
+    pset.execute(AdvectionEE, runtime=runtime, dt=np.timedelta64(6, "h"))
     np.testing.assert_allclose(pset.lat, latp, atol=1e-1)
 
 

tests/test_xgrid.py

@@ -136,6 +136,13 @@ def test_invalid_depth():
         XGrid.from_dataset(ds)
 
 
+def test_dim_without_axis():
+    ds = xr.Dataset({"z1d": (["depth"], [0])}, coords={"depth": [0]})
+    grid = XGrid.from_dataset(ds)
+    with pytest.raises(ValueError, match='Dimension "depth" has no axis attribute*'):
+        Field("z1d", ds["z1d"], grid, XLinear)
+
+
 def test_vertical1D_field():
     nz = 11
     ds = xr.Dataset(