Merge pull request #2391 from Parcels-code/bugfix/2383-curvilinear-periodic-point-in-cell

Adjust longitudes for antimeridian cells to account for periodicity

Author: fluidnumerics-joe

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/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)