Adding XLinearInvdistLandTracer interpolator

Note that the interpolator now requires nested-for-loops. So either has to be vectorized or written with numba, for speed

Author: erikvansebille

src/parcels/interpolators.py

@@ -21,6 +21,7 @@
     "UXPiecewiseLinearNode",
     "XFreeslip",
     "XLinear",
+    "XLinearInvdistLandTracer",
     "XNearest",
     "XPartialslip",
     "ZeroInterpolator",
@@ -570,6 +571,52 @@ def XNearest(
     return value.compute() if is_dask_collection(value) else value
 
 
+def XLinearInvdistLandTracer(
+    particle_positions: dict[str, float | np.ndarray],
+    grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]],
+    field: Field,
+):
+    """Linear spatial interpolation on a regular grid, where points on land are not used."""
+    values = XLinear(particle_positions, grid_positions, field)
+
+    on_land = np.argwhere(np.isnan(values))
+
+    xi, xsi = grid_positions["X"]["index"], grid_positions["X"]["bcoord"]
+    yi, eta = grid_positions["Y"]["index"], grid_positions["Y"]["bcoord"]
+    zi, zeta = grid_positions["Z"]["index"], grid_positions["Z"]["bcoord"]
+    ti, tau = grid_positions["T"]["index"], grid_positions["T"]["bcoord"]
+
+    axis_dim = field.grid.get_axis_dim_mapping(field.data.dims)
+    lenT = 2 if np.any(tau > 0) else 1
+    lenZ = 2 if np.any(zeta > 0) else 1
+
+    corner_data = _get_corner_data_Agrid(field.data, ti, zi, yi, xi, lenT, lenZ, len(xsi), axis_dim)
+
+    def is_land(p: int):
+        value = corner_data[:, :, :, :, p]
+        return np.where(np.isnan(value), True, False)
+
+    for p in on_land:
+        land = is_land(p)
+        nb_land = np.sum(land)
+        if nb_land == 4 * lenZ * lenT:
+            values[p] = 0.0
+        else:
+            val = 0
+            w_sum = 0
+            for t in range(lenT):
+                for k in range(lenZ):
+                    for j in range(2):
+                        for i in range(2):
+                            if land[t][k][j][i] == 0:
+                                distance = pow((eta[p] - j), 2) + pow((xsi[p] - i), 2)
+                                val += corner_data[t, k, j, i, p] / distance
+                                w_sum += 1 / distance
+            values[p] = val / w_sum
+
+    return values.compute() if is_dask_collection(values) else values
+
+
 def UXPiecewiseConstantFace(
     particle_positions: dict[str, float | np.ndarray],
     grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]],

tests/test_interpolation.py

@@ -21,6 +21,7 @@
     UXPiecewiseLinearNode,
     XFreeslip,
     XLinear,
+    XLinearInvdistLandTracer,
     XNearest,
     XPartialslip,
     ZeroInterpolator,
@@ -80,6 +81,7 @@ def field():
             [1.49, 6.49, 13.99],
             id="Linear-3",
         ),
+        pytest.param(XLinearInvdistLandTracer, np.timedelta64(1, "s"), 2.5, 0.49, 0.51, 13.99, id="LinearInvDistLand"),
         pytest.param(
             XNearest,
             [np.timedelta64(0, "s"), np.timedelta64(3, "s")],
@@ -131,6 +133,38 @@ def test_spatial_slip_interpolation(field, func, t, z, y, x, expected):
     np.testing.assert_array_almost_equal(velocities, expected)
 
 
+@pytest.mark.parametrize(
+    "func, t, z, y, x, expected",
+    [
+        (XLinearInvdistLandTracer, np.timedelta64(1, "s"), 0, 0.5, 0.5, 1.0),
+        (XLinearInvdistLandTracer, np.timedelta64(1, "s"), 0, 1.5, 1.5, 0.0),
+        (
+            XLinearInvdistLandTracer,
+            [np.timedelta64(0, "s"), np.timedelta64(1, "s")],
+            [0, 2],
+            [0.5, 0.5],
+            [0.5, 0.5],
+            1.0,
+        ),
+        (
+            XLinearInvdistLandTracer,
+            [np.timedelta64(0, "s"), np.timedelta64(1, "s")],
+            [0, 2],
+            [0.5, 1.5],
+            [0.5, 1.5],
+            [1.0, 0.0],
+        ),
+    ],
+)
+def test_invdistland_interpolation(field, func, t, z, y, x, expected):
+    field.data[:] = 1.0
+    field.data[:, :, 1:3, 1:3] = np.nan  # Set NaN land value to test inv_dist
+    field.interp_method = func
+
+    value = field[t, z, y, x]
+    np.testing.assert_array_almost_equal(value, expected)
+
+
 @pytest.mark.parametrize("mesh", ["spherical", "flat"])
 def test_interpolation_mesh_type(mesh, npart=10):
     ds = simple_UV_dataset(mesh=mesh)