Merge pull request #2461 from Parcels-code/add_XLinearVelocity

Adding XLinear_Velocity interpolator

Author: erikvansebille

docs/user_guide/examples/explanation_kernelloop.md

@@ -71,7 +71,12 @@ ds_fields["VWind"] = xr.DataArray(
 fieldset = parcels.FieldSet.from_copernicusmarine(ds_fields)
 
 # Create a vecorfield for the wind
-windvector = parcels.VectorField("Wind", fieldset.UWind, fieldset.VWind)
+windvector = parcels.VectorField(
+    "Wind",
+    fieldset.UWind,
+    fieldset.VWind,
+    vector_interp_method=parcels.interpolators.XLinear_Velocity
+)
 fieldset.add_field(windvector)
 ```
 

docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb

@@ -147,11 +147,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "u, v = fieldset.UV.eval(\n",
-    "    np.array([0]), np.array([0]), np.array([20]), np.array([50]), apply_conversion=False\n",
-    ")  # do not convert m/s to deg/s\n",
+    "u, v = fieldset.UV.eval(np.array([0]), np.array([0]), np.array([60]), np.array([50]))\n",
+    "u *= 1852 * 60 * np.cos(np.deg2rad(60))  # convert from degrees s^-1 to m s^-1\n",
+    "v *= 1852 * 60  # convert from degrees s^-1 to m s\n",
     "print(f\"(u, v) = ({u[0]:.3f}, {v[0]:.3f})\")\n",
-    "assert np.isclose(u, 1.0, atol=1e-3)"
+    "assert np.isclose(u, 1.0, atol=1e-3)\n",
+    "assert np.isclose(v, 0.0, atol=1e-3)"
    ]
   },
   {

docs/user_guide/examples/tutorial_unitconverters.ipynb

@@ -15,7 +15,9 @@
    "source": [
     "In most applications, Parcels works with `spherical` meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use `flat` meshes, where longitude and latitude are given in meters (note that the dimensions are then still called `longitude` and `latitude` for consistency reasons).\n",
     "\n",
-    "In all cases, velocities are given in m/s. Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this guide explains how this works.\n"
+    "In all cases, velocities are given in m s<sup>-1</sup>. But for advection (such as `particles.dlat = v * particles.dt`), the velocity needs to be in degrees s<sup>-1</sup>. \n",
+    "\n",
+    "Parcels seamlessly converts between meters s<sup>-1</sup> and degrees s<sup>-1</sup>, under the hood. For transparency, this guide explains how this works.\n"
    ]
   },
   {
@@ -130,7 +132,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Indeed, if we multiply the value of the U field with 1852 \\* 60 \\* cos(lat) (the number of meters in 1 degree of longitude) and the value of the V field with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m s<sup>-1</sup> for `u` and `v`.\n"
+    "Indeed, if we multiply the value of `u` with 1852 \\* 60 \\* cos(lat) (the number of meters in 1 degree of longitude) and the value of `v` with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m s<sup>-1</sup> for `u` and `v`.\n"
    ]
   },
   {
@@ -149,31 +151,6 @@
     "assert np.isclose(v, 1.0)"
    ]
   },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "You can also interpolate the Field to these values directly by using the `eval()` method and setting `applyConversion=False`, as below\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "print(\n",
-    "    fieldset.UV.eval(\n",
-    "        time,\n",
-    "        z,\n",
-    "        lat,\n",
-    "        lon,\n",
-    "        apply_conversion=False,\n",
-    "    )\n",
-    ")"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -185,7 +162,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Sampling `U` and `V` separately will _not_ convert to degrees s<sup>-1</sup>, so these velocities cannot be used directly for advection on spherical coordinates. This is one of the main reasons to always use the `UV` VectorField for velocity sampling in Parcels.\n"
+    "Sampling `U` and `V` separately will _not_ convert to degrees s<sup>-1</sup>, so these velocities _should not_ be used directly for advection on spherical coordinates. This is one of the main reasons to always use the `UV` VectorField for velocity sampling in Parcels.\n"
    ]
   },
   {
@@ -204,7 +181,7 @@
    "source": [
     "## Unit conversion for other fields such as diffusivity\n",
     "\n",
-    "For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a `spherical` mesh. For example, if we define a diffusivity field with value 10 m<sup>2</sup> s<sup>-1</sup>, Parcels will not convert this to degrees<sup>2</sup> s<sup>-1</sup> under the hood. \n",
+    "For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a `spherical` mesh. For example, if we define a diffusivity field in m<sup>2</sup> s<sup>-1</sup>, Parcels will _not_ convert this to degrees<sup>2</sup> s<sup>-1</sup> under the hood. \n",
     "\n",
     "If you want to work with diffusivity in degrees<sup>2</sup> s<sup>-1</sup> (for example to move particles using a random walk), you will have to convert this yourself in your kernel. \n",
     "\n",

docs/user_guide/v4-migration.md

@@ -40,7 +40,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
 - `Field.eval()` returns an array of floats instead of a single float (related to the vectorization)
 - `Field.eval()` does not throw OutOfBounds or other errors
 - The `NestedField` class has been removed. See the Nested Grids how-to guide for how to set up Nested Grids in v4.
-- `applyConversion` has been renamed to `apply_conversion` and only works for VectorFields. Conversion of units should be handled in Kernels.
+- `applyConversion` has been removed. Interpolation on VectorFields automatically converts from m/s to degrees/s for spherical meshes. Other conversion of units should be handled in Interpolators or Kernels.
 
 ## GridSet
 

src/parcels/_core/field.py

@@ -247,8 +247,10 @@ def __init__(
         W: Field | None = None,  # noqa: N803
         vector_interp_method: Callable | None = None,
     ):
-        _assert_str_and_python_varname(name)
+        if vector_interp_method is None:
+            raise ValueError("vector_interp_method must be provided for VectorField initialization.")
 
+        _assert_str_and_python_varname(name)
         self.name = name
         self.U = U
         self.V = V
@@ -268,12 +270,8 @@ def __init__(
         else:
             self.vector_type = "2D"
 
-        # Setting the interpolation method dynamically
-        if vector_interp_method is None:
-            self._vector_interp_method = None
-        else:
-            assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector, context="Interpolation")
-            self._vector_interp_method = vector_interp_method
+        assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector, context="Interpolation")
+        self._vector_interp_method = vector_interp_method
 
     def __repr__(self):
         return vectorfield_repr(self)
@@ -287,7 +285,7 @@ def vector_interp_method(self, method: Callable):
         assert_same_function_signature(method, ref=ZeroInterpolator_Vector, context="Interpolation")
         self._vector_interp_method = method
 
-    def eval(self, time: datetime, z, y, x, particles=None, apply_conversion=True):
+    def eval(self, time: datetime, z, y, x, particles=None):
         """Interpolate vectorfield values in space and time.
 
         Parameters
@@ -300,10 +298,6 @@ def eval(self, time: datetime, z, y, x, particles=None, apply_conversion=True):
         particles : ParticleSet, optional
             If provided, used to associate results with particle indices and to
             update particle state and element indices. Defaults to None.
-        apply_conversion : bool, default True
-            If True and the underlying grid is spherical, the horizontal velocity
-            components (U, V) are converted from m s^-1 to degrees s^-1.
-            If False, raw interpolated values are returned.
 
         Returns
         -------
@@ -327,26 +321,7 @@ def eval(self, time: datetime, z, y, x, particles=None, apply_conversion=True):
 
         particle_positions, grid_positions = _get_positions(self.U, time, z, y, x, particles, _ei)
 
-        if self._vector_interp_method is None:
-            u = self.U._interp_method(particle_positions, grid_positions, self.U)
-            v = self.V._interp_method(particle_positions, grid_positions, self.V)
-
-            if apply_conversion and self.U.grid._mesh == "spherical":
-                u /= 1852 * 60.0 * np.cos(np.deg2rad(y))
-                v /= 1852 * 60.0
-
-            if "3D" in self.vector_type:
-                w = self.W._interp_method(particle_positions, grid_positions, self.W)
-            else:
-                w = 0.0
-        else:
-            (u, v, w) = self._vector_interp_method(particle_positions, grid_positions, self)
-
-            if apply_conversion:
-                if self.U.grid._mesh == "spherical":
-                    meshJac = 1852 * 60.0 * np.cos(np.deg2rad(y))
-                    u = u / meshJac
-                    v = v / meshJac
+        (u, v, w) = self._vector_interp_method(particle_positions, grid_positions, self)
 
         for vel in (u, v, w):
             _update_particle_states_interp_value(particles, vel)

src/parcels/_core/fieldset.py

@@ -20,7 +20,14 @@
 from parcels._logger import logger
 from parcels._reprs import fieldset_repr
 from parcels._typing import Mesh
-from parcels.interpolators import UxPiecewiseConstantFace, UxPiecewiseLinearNode, XConstantField, XLinear
+from parcels.interpolators import (
+    Ux_Velocity,
+    UxPiecewiseConstantFace,
+    UxPiecewiseLinearNode,
+    XConstantField,
+    XLinear,
+    XLinear_Velocity,
+)
 
 if TYPE_CHECKING:
     from parcels._core.basegrid import BaseGrid
@@ -266,9 +273,11 @@ def from_fesom2(cls, ds: ux.UxDataset):
 
             if "W" in ds.data_vars:
                 fields["W"] = Field("W", ds["W"], grid, _select_uxinterpolator(ds["U"]))
-                fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"])
+                fields["UVW"] = VectorField(
+                    "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=Ux_Velocity
+                )
             else:
-                fields["UV"] = VectorField("UV", fields["U"], fields["V"])
+                fields["UV"] = VectorField("UV", fields["U"], fields["V"], vector_interp_method=Ux_Velocity)
 
         for varname in set(ds.data_vars) - set(fields.keys()):
             fields[varname] = Field(varname, ds[varname], grid, _select_uxinterpolator(ds[varname]))
@@ -350,11 +359,13 @@ def from_sgrid_conventions(
         if "U" in ds.data_vars and "V" in ds.data_vars:
             fields["U"] = Field("U", ds["U"], grid, XLinear)
             fields["V"] = Field("V", ds["V"], grid, XLinear)
-            fields["UV"] = VectorField("UV", fields["U"], fields["V"])
+            fields["UV"] = VectorField("UV", fields["U"], fields["V"], vector_interp_method=XLinear_Velocity)
 
             if "W" in ds.data_vars:
                 fields["W"] = Field("W", ds["W"], grid, XLinear)
-                fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"])
+                fields["UVW"] = VectorField(
+                    "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=XLinear_Velocity
+                )
 
         for varname in set(ds.data_vars) - set(fields.keys()) - skip_vars:
             fields[varname] = Field(varname, ds[varname], grid, XLinear)

src/parcels/_datasets/structured/generated.py

@@ -15,6 +15,7 @@
 
 def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh="spherical"):
     max_lon = 180.0 if mesh == "spherical" else 1e6
+    max_lat = 90.0 if mesh == "spherical" else 1e6
 
     return xr.Dataset(
         {"U": (["time", "depth", "YG", "XG"], np.zeros(dims)), "V": (["time", "depth", "YG", "XG"], np.zeros(dims))},
@@ -25,7 +26,7 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh="spherical"):
             "YG": (["YG"], np.arange(dims[2]), {"axis": "Y", "c_grid_axis_shift": -0.5}),
             "XC": (["XC"], np.arange(dims[3]) + 0.5, {"axis": "X"}),
             "XG": (["XG"], np.arange(dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}),
-            "lat": (["YG"], np.linspace(-90, 90, dims[2]), {"axis": "Y", "c_grid_axis_shift": -0.5}),
+            "lat": (["YG"], np.linspace(-max_lat, max_lat, dims[2]), {"axis": "Y", "c_grid_axis_shift": -0.5}),
             "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}),
         },
     ).pipe(

src/parcels/interpolators.py

@@ -20,10 +20,12 @@
     "CGrid_Velocity",
     "UxPiecewiseConstantFace",
     "UxPiecewiseLinearNode",
+    "Ux_Velocity",
     "XConstantField",
     "XFreeslip",
     "XLinear",
     "XLinearInvdistLandTracer",
+    "XLinear_Velocity",
     "XNearest",
     "XPartialslip",
     "ZeroInterpolator",
@@ -146,6 +148,25 @@ def XConstantField(
     return field.data[0, 0, 0, 0].values
 
 
+def XLinear_Velocity(
+    particle_positions: dict[str, float | np.ndarray],
+    grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]],
+    vectorfield: VectorField,
+):
+    """Trilinear interpolation on a regular grid for VectorFields of velocity."""
+    u = XLinear(particle_positions, grid_positions, vectorfield.U)
+    v = XLinear(particle_positions, grid_positions, vectorfield.V)
+    if vectorfield.grid._mesh == "spherical":
+        u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"]))
+        v /= 1852 * 60
+
+    if vectorfield.W:
+        w = XLinear(particle_positions, grid_positions, vectorfield.W)
+    else:
+        w = 0.0
+    return u, v, w
+
+
 def CGrid_Velocity(
     particle_positions: dict[str, float | np.ndarray],
     grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]],
@@ -268,9 +289,10 @@ def CGrid_Velocity(
     U = (1 - xsi) * U0 + xsi * U1
     V = (1 - eta) * V0 + eta * V1
 
-    meshJac = 1852 * 60.0 if grid._mesh == "spherical" else 1
-
-    jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
+    if grid._mesh == "spherical":
+        jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * 1852 * 60.0
+    else:
+        jac = i_u._compute_jacobian_determinant(py, px, eta, xsi)
 
     u = (
         (-(1 - eta) * U - (1 - xsi) * V) * px[0]
@@ -288,6 +310,11 @@ def CGrid_Velocity(
         u = u.compute()
         v = v.compute()
 
+    if grid._mesh == "spherical":
+        conversion = 1852 * 60.0 * np.cos(np.deg2rad(particle_positions["lat"]))
+        u /= conversion
+        v /= conversion
+
     # 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]
     dlon = xx - lon
@@ -682,3 +709,22 @@ def UxPiecewiseLinearNode(
     zk = field.grid.z.values[zi]
     zkp1 = field.grid.z.values[zi + 1]
     return (fzk * (zkp1 - z) + fzkp1 * (z - zk)) / (zkp1 - zk)  # Linear interpolation in the vertical direction
+
+
+def Ux_Velocity(
+    particle_positions: dict[str, float | np.ndarray],
+    grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]],
+    vectorfield: VectorField,
+):
+    """Interpolation kernel for Vectorfields of velocity on a UxGrid."""
+    u = vectorfield.U._interp_method(particle_positions, grid_positions, vectorfield.U)
+    v = vectorfield.V._interp_method(particle_positions, grid_positions, vectorfield.V)
+    if vectorfield.grid._mesh == "spherical":
+        u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"]))
+        v /= 1852 * 60
+
+    if "3D" in vectorfield.vector_type:
+        w = vectorfield.W._interp_method(particle_positions, grid_positions, vectorfield.W)
+    else:
+        w = 0.0
+    return u, v, w

tests-v3/test_fieldset_sampling.py

@@ -29,10 +29,6 @@ def SampleUV(particle, fieldset, time):  # pragma: no cover
     (particle.u, particle.v) = fieldset.UV[time, particle.depth, particle.lat, particle.lon]
 
 
-def SampleUVNoConvert(particle, fieldset, time):  # pragma: no cover
-    (particle.u, particle.v) = fieldset.UV.eval(time, particle.depth, particle.lat, particle.lon, applyConversion=False)
-
-
 def SampleP(particle, fieldset, time):  # pragma: no cover
     particle.p = fieldset.P[particle]
 
@@ -439,24 +435,6 @@ def test_fieldset_sample_geographic(fieldset_geometric):
     assert np.allclose(pset.u, lat, rtol=1e-6)
 
 
-@pytest.mark.v4alpha
-@pytest.mark.xfail(reason="GH1946")
-def test_fieldset_sample_geographic_noconvert(fieldset_geometric):
-    """Sample a fieldset without conversion to geographic units."""
-    npart = 120
-    fieldset = fieldset_geometric
-    lon = np.linspace(-170, 170, npart)
-    lat = np.linspace(-80, 80, npart)
-
-    pset = ParticleSet(fieldset, pclass=pclass(), lon=lon, lat=np.zeros(npart) + 70.0)
-    pset.execute(pset.Kernel(SampleUVNoConvert), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.v, lon * 1000 * 1.852 * 60, rtol=1e-6)
-
-    pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0)
-    pset.execute(pset.Kernel(SampleUVNoConvert), endtime=1.0, dt=1.0)
-    assert np.allclose(pset.u, lat * 1000 * 1.852 * 60, rtol=1e-6)
-
-
 @pytest.mark.v4alpha
 @pytest.mark.xfail(reason="GH1946")
 def test_fieldset_sample_geographic_polar(fieldset_geometric_polar):