From 1e17a46f0148ec2f8c51262105b34b268c0a6b7f Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 9 Jan 2026 16:15:01 +0100 Subject: [PATCH 01/11] Adding XLinear_Velocity interpolator --- src/parcels/interpolators.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 8e23d49521..d2deaa9315 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -146,6 +146,20 @@ 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, +): + u = XLinear(particle_positions, grid_positions, vectorfield.U) + v = XLinear(particle_positions, grid_positions, vectorfield.V) + 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]], From 49a0a823deebe549d6e1ab5f851829ae44d22f6e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Sun, 11 Jan 2026 20:40:33 +0100 Subject: [PATCH 02/11] Making unit tests of zonal and meridional flow stricter Following the discussion around zonal and meridional spherical conversion in #2455 --- src/parcels/_datasets/structured/generated.py | 3 +- tests/test_advection.py | 40 +++++++++++++++---- 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/src/parcels/_datasets/structured/generated.py b/src/parcels/_datasets/structured/generated.py index 20440717c2..1272607bec 100644 --- a/src/parcels/_datasets/structured/generated.py +++ b/src/parcels/_datasets/structured/generated.py @@ -14,6 +14,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))}, @@ -24,7 +25,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( diff --git a/tests/test_advection.py b/tests/test_advection.py index 0a0b7484aa..7afecea177 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -30,22 +30,27 @@ @pytest.mark.parametrize("mesh", ["spherical", "flat"]) def test_advection_zonal(mesh, npart=10): - """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`.""" + """Particles at high latitude move geographically faster due to the pole correction.""" ds = simple_UV_dataset(mesh=mesh) ds["U"].data[:] = 1.0 fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) - pset = ParticleSet(fieldset, lon=np.zeros(npart) + 20.0, lat=np.linspace(0, 80, npart)) - pset.execute(AdvectionRK4, runtime=np.timedelta64(2, "h"), dt=np.timedelta64(15, "m")) + runtime = 7200 + startlat = np.linspace(0, 80, npart) + startlon = 20.0 + np.zeros(npart) + pset = ParticleSet(fieldset, lon=startlon, lat=startlat) + pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) + expected_dlon = runtime if mesh == "spherical": - assert (np.diff(pset.lon) > 1.0e-4).all() - else: - assert (np.diff(pset.lon) < 1.0e-4).all() + expected_dlon /= 1852 * 60 * np.cos(np.deg2rad(pset.lat)) + + np.testing.assert_allclose(pset.lon - startlon, expected_dlon, atol=1e-5) + np.testing.assert_allclose(pset.lat, startlat, atol=1e-5) def test_advection_zonal_with_particlefile(tmp_store): - """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`.""" + """Particles at high latitude move geographically faster due to the pole correction.""" npart = 10 ds = simple_UV_dataset(mesh="flat") ds["U"].data[:] = 1.0 @@ -88,6 +93,27 @@ def test_advection_zonal_periodic(): np.testing.assert_allclose(pset.lat, 0.5, atol=1e-5) +@pytest.mark.parametrize("mesh", ["spherical", "flat"]) +def test_advection_meridional(mesh, npart=10): + """All particles move the same in meridional direction, regardless of latitude.""" + ds = simple_UV_dataset(mesh=mesh) + ds["V"].data[:] = 1.0 + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) + + runtime = 7200 + startlat = np.linspace(0, 80, npart) + startlon = 20.0 + np.zeros(npart) + pset = ParticleSet(fieldset, lon=startlon, lat=startlat) + pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) + + expected_dlat = runtime + if mesh == "spherical": + expected_dlat /= 1852 * 60 + + np.testing.assert_allclose(pset.lon, startlon, atol=1e-5) + np.testing.assert_allclose(pset.lat - startlat, expected_dlat, atol=1e-4) + + def test_horizontal_advection_in_3D_flow(npart=10): """Flat 2D zonal flow that increases linearly with z from 0 m/s to 1 m/s.""" ds = simple_UV_dataset(mesh="flat") From 34570ebc5cefcdd0166054d351f8dba846fe0425 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Sun, 11 Jan 2026 20:47:58 +0100 Subject: [PATCH 03/11] Also testing horizontal advection in 3D for spherical mesh To confirm unit conversion works --- tests/test_advection.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index 7afecea177..0a81c04de9 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -114,17 +114,20 @@ def test_advection_meridional(mesh, npart=10): np.testing.assert_allclose(pset.lat - startlat, expected_dlat, atol=1e-4) -def test_horizontal_advection_in_3D_flow(npart=10): - """Flat 2D zonal flow that increases linearly with z from 0 m/s to 1 m/s.""" - ds = simple_UV_dataset(mesh="flat") +@pytest.mark.parametrize("mesh", ["spherical", "flat"]) +def test_horizontal_advection_in_3D_flow(mesh, npart=10): + """2D zonal flow that increases linearly with z from 0 m/s to 1 m/s.""" + ds = simple_UV_dataset(mesh=mesh) ds["U"].data[:] = 1.0 ds["U"].data[:, 0, :, :] = 0.0 # Set U to 0 at the surface - fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.zeros(npart), z=np.linspace(0.1, 0.9, npart)) pset.execute(AdvectionRK4, runtime=np.timedelta64(2, "h"), dt=np.timedelta64(15, "m")) expected_lon = pset.z * pset.time + if mesh == "spherical": + expected_lon /= 1852 * 60 * np.cos(np.deg2rad(pset.lat)) np.testing.assert_allclose(pset.lon, expected_lon, atol=1.0e-1) From b6079cbb91b704d0553b7152d59a4e976cb1cf6e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 12 Jan 2026 18:42:37 +0100 Subject: [PATCH 04/11] Adding XLinear_Velocity as default vector_interp_method --- src/parcels/_core/fieldset.py | 14 +++++++++++--- src/parcels/interpolators.py | 2 ++ tests/test_advection.py | 18 +++++++++--------- 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 6b6c14682f..6a64e60fc6 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -20,7 +20,13 @@ 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 ( + UxPiecewiseConstantFace, + UxPiecewiseLinearNode, + XConstantField, + XLinear, + XLinear_Velocity, +) if TYPE_CHECKING: from parcels._core.basegrid import BaseGrid @@ -358,11 +364,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) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index e59d7df779..83e82406a0 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -24,6 +24,7 @@ "XFreeslip", "XLinear", "XLinearInvdistLandTracer", + "XLinear_Velocity", "XNearest", "XPartialslip", "ZeroInterpolator", @@ -151,6 +152,7 @@ def XLinear_Velocity( 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.W: diff --git a/tests/test_advection.py b/tests/test_advection.py index 1ce0d5180c..c73b82cab6 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -14,7 +14,7 @@ simple_UV_dataset, stommel_gyre_dataset, ) -from parcels.interpolators import CGrid_Velocity, XLinear +from parcels.interpolators import CGrid_Velocity, XLinear, XLinear_Velocity from parcels.kernels import ( AdvectionDiffusionEM, AdvectionDiffusionM1, @@ -213,10 +213,10 @@ def test_length1dimensions(u, v, w): # TODO: Refactor this test to be more read grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) - fields = [U, V, VectorField("UV", U, V)] + fields = [U, V, VectorField("UV", U, V, vector_interp_method=XLinear_Velocity)] if w: W = Field("W", ds["W"], grid, interp_method=XLinear) - fields.append(VectorField("UVW", U, V, W)) + fields.append(VectorField("UVW", U, V, W, vector_interp_method=XLinear_Velocity)) fieldset = FieldSet(fields) x0, y0, z0 = 2, 8, -4 @@ -236,7 +236,7 @@ def test_radialrotation(npart=10): grid = XGrid.from_dataset(ds, mesh="flat") U = parcels.Field("U", ds["U"], grid, interp_method=XLinear) V = parcels.Field("V", ds["V"], grid, interp_method=XLinear) - UV = parcels.VectorField("UV", U, V) + UV = parcels.VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = parcels.FieldSet([U, V, UV]) dt = np.timedelta64(30, "s") @@ -277,10 +277,10 @@ def test_moving_eddy(kernel, rtol): if kernel in [AdvectionRK2_3D, AdvectionRK4_3D]: # Using W to test 3D advection (assuming same velocity as V) W = Field("W", ds["V"], grid, interp_method=XLinear) - UVW = VectorField("UVW", U, V, W) + UVW = VectorField("UVW", U, V, W, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, W, UVW]) else: - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, UV]) if kernel in [AdvectionDiffusionEM, AdvectionDiffusionM1]: # Add zero diffusivity field for diffusion kernels @@ -328,7 +328,7 @@ def test_decaying_moving_eddy(kernel, rtol): grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, UV]) start_lon, start_lat = 10000, 10000 @@ -376,7 +376,7 @@ def test_stommelgyre_fieldset(kernel, rtol, grid_type): npart = 2 ds = stommel_gyre_dataset(grid_type=grid_type) grid = XGrid.from_dataset(ds, mesh="flat") - vector_interp_method = None if grid_type == "A" else CGrid_Velocity + vector_interp_method = XLinear_Velocity if grid_type == "A" else CGrid_Velocity U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) P = Field("P", ds["P"], grid, interp_method=XLinear) @@ -420,7 +420,7 @@ def test_peninsula_fieldset(kernel, rtol, grid_type): U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) P = Field("P", ds["P"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, P, UV]) dt = np.timedelta64(30, "m") From ac21aa2b3d8afe77c68a97533f39edf91b87cb63 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 12 Jan 2026 18:57:24 +0100 Subject: [PATCH 05/11] forcing vector_interp_method to not be None --- src/parcels/_core/field.py | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index ffbf7731cb..3bf371590e 100644 --- a/src/parcels/_core/field.py +++ b/src/parcels/_core/field.py @@ -272,12 +272,10 @@ 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 + raise ValueError("vector_interp_method must be provided for VectorField initialization.") + 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) @@ -308,15 +306,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 "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) + (u, v, w) = self._vector_interp_method(particle_positions, grid_positions, self) if apply_conversion: u = self.U.units.to_target(u, z, y, x) From 842b543d0905a23816f70000a413b3e53d2e4638 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 12 Jan 2026 18:58:09 +0100 Subject: [PATCH 06/11] Adding Ux_Velocity interpolator --- src/parcels/_core/fieldset.py | 7 +++++-- src/parcels/interpolators.py | 16 ++++++++++++++++ tests/test_fieldset.py | 7 +++---- 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 6a64e60fc6..dfd3a285cb 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -21,6 +21,7 @@ from parcels._reprs import fieldset_repr from parcels._typing import Mesh from parcels.interpolators import ( + Ux_Velocity, UxPiecewiseConstantFace, UxPiecewiseLinearNode, XConstantField, @@ -280,9 +281,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])) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 83e82406a0..a7324e65bf 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -20,6 +20,7 @@ "CGrid_Velocity", "UxPiecewiseConstantFace", "UxPiecewiseLinearNode", + "Ux_Velocity", "XConstantField", "XFreeslip", "XLinear", @@ -698,3 +699,18 @@ 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 "3D" in vectorfield.vector_type: + w = vectorfield.W._interp_method(particle_positions, grid_positions, vectorfield.W) + else: + w = 0.0 + return u, v, w diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 131826647f..b2443a51c0 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -13,7 +13,7 @@ from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.structured.generic import datasets_sgrid from parcels._datasets.unstructured.generic import datasets as datasets_unstructured -from parcels.interpolators import XLinear +from parcels.interpolators import XLinear, XLinear_Velocity from tests import utils ds = datasets_structured["ds_2d_left"] @@ -25,7 +25,7 @@ def fieldset() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet( [U, V, UV], @@ -165,7 +165,6 @@ def test_fieldset_init_incompatible_calendars(): grid = XGrid.from_dataset(ds1, mesh="flat") U = Field("U", ds1["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds1["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) ds2 = ds.copy() ds2["time"] = ( @@ -177,7 +176,7 @@ def test_fieldset_init_incompatible_calendars(): incompatible_calendar = Field("test", ds2["data_g"], grid2, interp_method=XLinear) with pytest.raises(CalendarError, match="Expected field '.*' to have calendar compatible with datetime object"): - FieldSet([U, V, UV, incompatible_calendar]) + FieldSet([U, V, incompatible_calendar]) def test_fieldset_add_field_incompatible_calendars(fieldset): From 4c50a6aef17e821797c2f02f103899d034d7c88f Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 13 Jan 2026 13:00:15 +0100 Subject: [PATCH 07/11] Updating apply_conversion descriptions --- .../examples/tutorial_nemo_curvilinear.ipynb | 9 ++--- .../examples/tutorial_unitconverters.ipynb | 35 ++++--------------- docs/user_guide/v4-migration.md | 2 +- src/parcels/_core/field.py | 6 +--- 4 files changed, 13 insertions(+), 39 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index ffe6f6a787..795743b863 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/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)" ] }, { diff --git a/docs/user_guide/examples/tutorial_unitconverters.ipynb b/docs/user_guide/examples/tutorial_unitconverters.ipynb index 4bcf19ede6..e75a686425 100644 --- a/docs/user_guide/examples/tutorial_unitconverters.ipynb +++ b/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-1. But for advection (such as `particles.dlat = v * particles.dt`), the velocity needs to be in degrees s-1. \n", + "\n", + "Parcels seamlessly converts between meters s-1 and degrees s-1, 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-1 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-1 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-1, 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-1, 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 m2 s-1, Parcels will not convert this to degrees2 s-1 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 m2 s-1, Parcels will _not_ convert this to degrees2 s-1 under the hood. \n", "\n", "If you want to work with diffusivity in degrees2 s-1 (for example to move particles using a random walk), you will have to convert this yourself in your kernel. \n", "\n", diff --git a/docs/user_guide/v4-migration.md b/docs/user_guide/v4-migration.md index af0056d2d9..6e2f4743df 100644 --- a/docs/user_guide/v4-migration.md +++ b/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 diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index 6a8f73f062..e404ad3627 100644 --- a/src/parcels/_core/field.py +++ b/src/parcels/_core/field.py @@ -285,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 @@ -298,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 ------- From 53b847ea1ad78ea1ef10fe149672e14af419ba46 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 13 Jan 2026 13:01:35 +0100 Subject: [PATCH 08/11] Converting velocities inside Interpolators --- src/parcels/interpolators.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index a7324e65bf..61b58650b0 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -156,6 +156,10 @@ def XLinear_Velocity( """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: @@ -285,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] @@ -305,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 @@ -709,6 +719,10 @@ def Ux_Velocity( """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: From cd4dcee4a2fee3c6810d9df8d5d79c32723c9793 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 13 Jan 2026 13:13:42 +0100 Subject: [PATCH 09/11] Adding vector_interp_method to all tests --- tests-v3/test_fieldset_sampling.py | 22 ---------------------- tests/test_particlefile.py | 4 ++-- tests/test_particleset_execute.py | 12 ++++++------ tests/test_particlesetview.py | 4 ++-- tests/test_uxarray_fieldset.py | 6 +++++- 5 files changed, 15 insertions(+), 33 deletions(-) diff --git a/tests-v3/test_fieldset_sampling.py b/tests-v3/test_fieldset_sampling.py index 1a3c68623e..291c27b885 100644 --- a/tests-v3/test_fieldset_sampling.py +++ b/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): diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index 2f6df98a5f..5b433c3e0a 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -25,7 +25,7 @@ from parcels._core.utils.time import TimeInterval, timedelta_to_float from parcels._datasets.structured.generated import peninsula_dataset from parcels._datasets.structured.generic import datasets -from parcels.interpolators import XLinear +from parcels.interpolators import XLinear, XLinear_Velocity from parcels.kernels import AdvectionRK4 from tests.common_kernels import DoNothing @@ -37,7 +37,7 @@ def fieldset() -> FieldSet: # TODO v4: Move into a `conftest.py` file and remov grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, XLinear) V = Field("V", ds["V_A_grid"], grid, XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet( [U, V, UV], diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index b989356dcd..5cab91ee33 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -23,7 +23,7 @@ from parcels._datasets.structured.generated import simple_UV_dataset from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.unstructured.generic import datasets as datasets_unstructured -from parcels.interpolators import UxPiecewiseConstantFace, UxPiecewiseLinearNode, XLinear +from parcels.interpolators import Ux_Velocity, UxPiecewiseConstantFace, UxPiecewiseLinearNode, XLinear, XLinear_Velocity from parcels.kernels import AdvectionEE, AdvectionRK2, AdvectionRK4, AdvectionRK4_3D, AdvectionRK45 from tests.common_kernels import DoNothing from tests.utils import DEFAULT_PARTICLES @@ -35,7 +35,7 @@ def fieldset() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet([U, V, UV]) @@ -47,7 +47,7 @@ def fieldset_no_time_interval() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet([U, V, UV]) @@ -494,7 +494,7 @@ def test_uxstommelgyre_pset_execute(): grid=grid, interp_method=UxPiecewiseConstantFace, ) - UV = VectorField(name="UV", U=U, V=V) + UV = VectorField(name="UV", U=U, V=V, vector_interp_method=Ux_Velocity) fieldset = FieldSet([UV, UV.U, UV.V, P]) pset = ParticleSet( fieldset, @@ -540,7 +540,7 @@ def test_uxstommelgyre_multiparticle_pset_execute(): grid=grid, interp_method=UxPiecewiseConstantFace, ) - UVW = VectorField(name="UVW", U=U, V=V, W=W) + UVW = VectorField(name="UVW", U=U, V=V, W=W, vector_interp_method=Ux_Velocity) fieldset = FieldSet([UVW, UVW.U, UVW.V, UVW.W, P]) pset = ParticleSet( fieldset, @@ -579,7 +579,7 @@ def test_uxstommelgyre_pset_execute_output(): grid=grid, interp_method=UxPiecewiseConstantFace, ) - UV = VectorField(name="UV", U=U, V=V) + UV = VectorField(name="UV", U=U, V=V, vector_interp_method=Ux_Velocity) fieldset = FieldSet([UV, UV.U, UV.V, P]) pset = ParticleSet( fieldset, diff --git a/tests/test_particlesetview.py b/tests/test_particlesetview.py index 847878f54e..d5d60161a2 100644 --- a/tests/test_particlesetview.py +++ b/tests/test_particlesetview.py @@ -4,7 +4,7 @@ from parcels import Field, FieldSet, Particle, ParticleSet, Variable, VectorField, XGrid from parcels._core.statuscodes import StatusCode from parcels._datasets.structured.generic import datasets as datasets_structured -from parcels.interpolators import XLinear +from parcels.interpolators import XLinear, XLinear_Velocity @pytest.fixture @@ -13,7 +13,7 @@ def fieldset() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet([U, V, UV]) diff --git a/tests/test_uxarray_fieldset.py b/tests/test_uxarray_fieldset.py index cf3c56dd92..b6895b52dc 100644 --- a/tests/test_uxarray_fieldset.py +++ b/tests/test_uxarray_fieldset.py @@ -13,6 +13,7 @@ ) from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels.interpolators import ( + Ux_Velocity, UxPiecewiseConstantFace, UxPiecewiseLinearNode, ) @@ -47,6 +48,7 @@ def uv_fesom_channel(ds_fesom_channel) -> VectorField: grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"], mesh="flat"), interp_method=UxPiecewiseConstantFace, ), + vector_interp_method=Ux_Velocity, ) return UV @@ -73,6 +75,7 @@ def uvw_fesom_channel(ds_fesom_channel) -> VectorField: grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"], mesh="flat"), interp_method=UxPiecewiseLinearNode, ), + vector_interp_method=Ux_Velocity, ) return UVW @@ -118,11 +121,12 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): U=Field(name="U", data=ds.U, grid=grid, interp_method=UxPiecewiseConstantFace), V=Field(name="V", data=ds.V, grid=grid, interp_method=UxPiecewiseConstantFace), W=Field(name="W", data=ds.W, grid=grid, interp_method=UxPiecewiseLinearNode), + vector_interp_method=Ux_Velocity, ) P = Field(name="p", data=ds.p, grid=grid, interp_method=UxPiecewiseLinearNode) fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W]) - (u, v, w) = fieldset.UVW.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], apply_conversion=False) + (u, v, w) = fieldset.UVW.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0]) assert np.allclose([u.item(), v.item(), w.item()], [1.0, 1.0, 0.0], rtol=1e-3, atol=1e-6) assert np.isclose( From 4248b3fa0492283f7ed2edec2321f5dc386b090d Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 13 Jan 2026 14:29:31 +0100 Subject: [PATCH 10/11] Adding vector_interp_method to explanation_kernelloop --- docs/user_guide/examples/explanation_kernelloop.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/user_guide/examples/explanation_kernelloop.md b/docs/user_guide/examples/explanation_kernelloop.md index 143c4aab81..025305ce32 100644 --- a/docs/user_guide/examples/explanation_kernelloop.md +++ b/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) ``` From 30e1245e6d8392383ea624d17a8c8cf737fe79e6 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 13 Jan 2026 15:06:48 +0100 Subject: [PATCH 11/11] Moving ValueError when no vector_interp_method to top of initi --- src/parcels/_core/field.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index e404ad3627..76e2c20cdc 100644 --- a/src/parcels/_core/field.py +++ b/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,8 +270,6 @@ def __init__( else: self.vector_type = "2D" - if vector_interp_method is None: - raise ValueError("vector_interp_method must be provided for VectorField initialization.") assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector, context="Interpolation") self._vector_interp_method = vector_interp_method