From 4027f8077be5457546642dfc902906e93dbdbcc4 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 20 May 2026 16:16:44 +0200 Subject: [PATCH 01/10] Rename _datasets.structured.generic::datasets to datasets_comodo --- src/parcels/_datasets/structured/generic.py | 8 ++-- tests/datasets/test_structured.py | 6 +-- tests/datasets/test_utils.py | 4 +- tests/test_field.py | 2 +- tests/test_fieldset.py | 2 +- tests/test_index_search.py | 4 +- tests/test_kernel.py | 2 +- tests/test_particlefile.py | 4 +- tests/test_particleset.py | 2 +- tests/test_particleset_execute.py | 2 +- tests/test_particlesetview.py | 2 +- tests/test_spatialhash.py | 8 ++-- tests/test_xgrid.py | 50 ++++++++++----------- 13 files changed, 48 insertions(+), 48 deletions(-) diff --git a/src/parcels/_datasets/structured/generic.py b/src/parcels/_datasets/structured/generic.py index ff729c8a13..88ab7609bc 100644 --- a/src/parcels/_datasets/structured/generic.py +++ b/src/parcels/_datasets/structured/generic.py @@ -5,7 +5,7 @@ from . import T, X, Y, Z -__all__ = ["T", "X", "Y", "Z", "datasets"] +__all__ = ["T", "X", "Y", "Z", "datasets_comodo"] TIME = xr.date_range("2000", "2001", T) @@ -139,7 +139,7 @@ def _unrolled_cone_curvilinear_grid(): ) -datasets = { +datasets_comodo = { "2d_left_rotated": _rotated_curvilinear_grid(), "ds_2d_left": xr.Dataset( # MITgcm indexing style { @@ -238,7 +238,7 @@ def _unrolled_cone_curvilinear_grid(): } datasets_sgrid = { "ds_2d_padded_high": ( - datasets["ds_2d_left"] + datasets_comodo["ds_2d_left"] .pipe( sgrid._attach_sgrid_metadata, sgrid.SGrid2DMetadata( @@ -258,7 +258,7 @@ def _unrolled_cone_curvilinear_grid(): ) ), "ds_2d_padded_low": ( - datasets["ds_2d_right"] + datasets_comodo["ds_2d_right"] .pipe( sgrid._attach_sgrid_metadata, sgrid.SGrid2DMetadata( diff --git a/tests/datasets/test_structured.py b/tests/datasets/test_structured.py index 2750583f58..9ca6195f10 100644 --- a/tests/datasets/test_structured.py +++ b/tests/datasets/test_structured.py @@ -1,12 +1,12 @@ import xgcm from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS -from parcels._datasets.structured.generic import datasets +from parcels._datasets.structured.generic import datasets_comodo def test_left_indexed_dataset(): """Checks that 'ds_2d_left' is right indexed on all variables.""" - ds = datasets["ds_2d_left"] + ds = datasets_comodo["ds_2d_left"] grid = xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS) for _axis_name, axis in grid.axes.items(): @@ -16,7 +16,7 @@ def test_left_indexed_dataset(): def test_right_indexed_dataset(): """Checks that 'ds_2d_right' is right indexed on all variables.""" - ds = datasets["ds_2d_right"] + ds = datasets_comodo["ds_2d_right"] grid = xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS) for _axis_name, axis in grid.axes.items(): for pos, _dim_name in axis.coords.items(): diff --git a/tests/datasets/test_utils.py b/tests/datasets/test_utils.py index b376283e30..914e4055bf 100644 --- a/tests/datasets/test_utils.py +++ b/tests/datasets/test_utils.py @@ -3,7 +3,7 @@ import xarray as xr from parcels._datasets import utils -from parcels._datasets.structured.generic import datasets +from parcels._datasets.structured.generic import datasets_comodo @pytest.fixture @@ -28,7 +28,7 @@ def nonzero_ds(): ) -@pytest.mark.parametrize("ds", [pytest.param(v, id=k) for k, v in datasets.items()]) +@pytest.mark.parametrize("ds", [pytest.param(v, id=k) for k, v in datasets_comodo.items()]) @pytest.mark.parametrize("except_for", [None, "coords"]) def test_replace_arrays_with_zeros(ds, except_for): # make sure doesn't error with range of datasets diff --git a/tests/test_field.py b/tests/test_field.py index cc264beae3..70fc6a7530 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -7,7 +7,7 @@ from parcels import Field, UxGrid, VectorField, XGrid from parcels._datasets.structured.generic import T as T_structured -from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.structured.generic import datasets_comodo as datasets_structured from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels.interpolators import ( UxConstantFaceConstantZC, diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 6eeef20a6d..198e1e54b6 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -10,7 +10,7 @@ from parcels import Field, ParticleFile, ParticleSet, VectorField, XGrid, convert from parcels._core.fieldset import CalendarError, FieldSet, _datetime_to_msg from parcels._datasets.structured.generic import T as T_structured -from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.structured.generic import datasets_comodo 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, XLinear_Velocity diff --git a/tests/test_index_search.py b/tests/test_index_search.py index b85e6e3b63..153eafafe1 100644 --- a/tests/test_index_search.py +++ b/tests/test_index_search.py @@ -5,13 +5,13 @@ import parcels.tutorial from parcels import Field, XGrid from parcels._core.index_search import _latlon_rad_to_xyz, _search_indices_curvilinear_2d -from parcels._datasets.structured.generic import datasets +from parcels._datasets.structured.generic import datasets_comodo from parcels.interpolators import XLinear @pytest.fixture def field_cone(): - ds = datasets["2d_left_unrolled_cone"] + ds = datasets_comodo["2d_left_unrolled_cone"] grid = XGrid.from_dataset(ds, mesh="flat") field = Field( name="test_field", diff --git a/tests/test_kernel.py b/tests/test_kernel.py index b4a3de9225..691e889384 100644 --- a/tests/test_kernel.py +++ b/tests/test_kernel.py @@ -8,7 +8,7 @@ XGrid, ) from parcels._core.kernel import Kernel -from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.structured.generic import datasets_comodo as datasets_structured from parcels.interpolators import XLinear from parcels.kernels import AdvectionRK4, AdvectionRK45 from tests.common_kernels import MoveEast, MoveNorth diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index b0938db64e..2f658c55c8 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -26,7 +26,7 @@ from parcels._core.particlefile import _get_schema 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._datasets.structured.generic import datasets_comodo from parcels.convert import copernicusmarine_to_sgrid from parcels.interpolators import XLinear, XLinear_Velocity from parcels.kernels import AdvectionRK4 @@ -36,7 +36,7 @@ @pytest.fixture def fieldset() -> FieldSet: # TODO v4: Move into a `conftest.py` file and remove duplicates """Fixture to create a FieldSet object for testing.""" - ds = datasets["ds_2d_left"] + ds = datasets_comodo["ds_2d_left"] 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) diff --git a/tests/test_particleset.py b/tests/test_particleset.py index 5c66341b5e..565655039c 100644 --- a/tests/test_particleset.py +++ b/tests/test_particleset.py @@ -15,7 +15,7 @@ Variable, XGrid, ) -from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.structured.generic import datasets_comodo as datasets_structured from parcels.interpolators import XLinear from tests.common_kernels import DoNothing from tests.utils import round_and_hash_float_array diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index 87a6fabd13..012d14f2ce 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -21,7 +21,7 @@ ) from parcels._core.utils.time import timedelta_to_float from parcels._datasets.structured.generated import simple_UV_dataset -from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.structured.generic import datasets_comodo as datasets_structured from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels.interpolators import ( Ux_Velocity, diff --git a/tests/test_particlesetview.py b/tests/test_particlesetview.py index 2f0532543c..fb28f1e05d 100644 --- a/tests/test_particlesetview.py +++ b/tests/test_particlesetview.py @@ -3,7 +3,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._datasets.structured.generic import datasets_comodo as datasets_structured from parcels.interpolators import XLinear, XLinear_Velocity diff --git a/tests/test_spatialhash.py b/tests/test_spatialhash.py index 6beedfbf6f..4d87882cbc 100644 --- a/tests/test_spatialhash.py +++ b/tests/test_spatialhash.py @@ -1,18 +1,18 @@ import numpy as np from parcels import XGrid -from parcels._datasets.structured.generic import datasets +from parcels._datasets.structured.generic import datasets_comodo def test_spatialhash_init(): - ds = datasets["2d_left_rotated"] + ds = datasets_comodo["2d_left_rotated"] grid = XGrid.from_dataset(ds, mesh="flat") spatialhash = grid.get_spatial_hash() assert spatialhash is not None def test_invalid_positions(): - ds = datasets["2d_left_rotated"] + ds = datasets_comodo["2d_left_rotated"] grid = XGrid.from_dataset(ds, mesh="flat") j, i, _ = grid.get_spatial_hash().query([np.nan, np.inf], [np.nan, np.inf]) @@ -21,7 +21,7 @@ def test_invalid_positions(): def test_mixed_positions(): - ds = datasets["2d_left_rotated"] + ds = datasets_comodo["2d_left_rotated"] grid = XGrid.from_dataset(ds, mesh="flat") lat = grid.lat.mean() lon = grid.lon.mean() diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index 53b0124e8a..c9e3881b6d 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -17,24 +17,24 @@ XGrid, _transpose_xfield_data_to_tzyx, ) -from parcels._datasets.structured.generic import X, Y, Z, datasets, datasets_sgrid +from parcels._datasets.structured.generic import X, Y, Z, datasets_comodo, datasets_sgrid from parcels.interpolators import XLinear from tests import utils GridTestCase = namedtuple("GridTestCase", ["ds", "attr", "expected"]) test_cases = [ - GridTestCase(datasets["ds_2d_left"], "lon", datasets["ds_2d_left"].XG.values), - GridTestCase(datasets["ds_2d_left"], "lat", datasets["ds_2d_left"].YG.values), - GridTestCase(datasets["ds_2d_left"], "depth", datasets["ds_2d_left"].ZG.values), + GridTestCase(datasets_comodo["ds_2d_left"], "lon", datasets_comodo["ds_2d_left"].XG.values), + GridTestCase(datasets_comodo["ds_2d_left"], "lat", datasets_comodo["ds_2d_left"].YG.values), + GridTestCase(datasets_comodo["ds_2d_left"], "depth", datasets_comodo["ds_2d_left"].ZG.values), GridTestCase( - datasets["ds_2d_left"], + datasets_comodo["ds_2d_left"], "time", - datasets["ds_2d_left"].time.values.astype(np.float64) / 1e9, + datasets_comodo["ds_2d_left"].time.values.astype(np.float64) / 1e9, ), - GridTestCase(datasets["ds_2d_left"], "xdim", X - 1), - GridTestCase(datasets["ds_2d_left"], "ydim", Y - 1), - GridTestCase(datasets["ds_2d_left"], "zdim", Z - 1), + GridTestCase(datasets_comodo["ds_2d_left"], "xdim", X - 1), + GridTestCase(datasets_comodo["ds_2d_left"], "ydim", Y - 1), + GridTestCase(datasets_comodo["ds_2d_left"], "zdim", Z - 1), ] @@ -48,7 +48,7 @@ def assert_equal(actual, expected): assert_allclose(actual, expected) -@pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) +@pytest.mark.parametrize("ds", [datasets_comodo["ds_2d_left"]]) def test_grid_init_param_types(ds): with pytest.raises(ValueError, match="Invalid value 'invalid'. Valid options are.*"): XGrid.from_dataset(ds, mesh="invalid") @@ -61,25 +61,25 @@ def test_xgrid_properties_ground_truth(ds, attr, expected): assert_equal(actual, expected) -@pytest.mark.parametrize("ds", [pytest.param(ds, id=key) for key, ds in datasets.items()]) +@pytest.mark.parametrize("ds", [pytest.param(ds, id=key) for key, ds in datasets_comodo.items()]) def test_xgrid_from_dataset_on_generic_datasets(ds): XGrid.from_dataset(ds, mesh="flat") -@pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) +@pytest.mark.parametrize("ds", [datasets_comodo["ds_2d_left"]]) def test_xgrid_axes(ds): grid = XGrid.from_dataset(ds, mesh="flat") assert grid.axes == ["Z", "Y", "X"] -@pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) +@pytest.mark.parametrize("ds", [datasets_comodo["ds_2d_left"]]) @pytest.mark.parametrize("mesh", ["flat", "spherical"]) def test_uxgrid_mesh(ds, mesh): grid = XGrid.from_dataset(ds, mesh=mesh) assert grid._mesh == mesh -@pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) +@pytest.mark.parametrize("ds", [datasets_comodo["ds_2d_left"]]) def test_transpose_xfield_data_to_tzyx(ds): da = ds["data_g"] grid = XGrid.from_dataset(ds, mesh="flat") @@ -93,7 +93,7 @@ def test_transpose_xfield_data_to_tzyx(ds): utils.assert_valid_field_data(da_test, grid) -@pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) +@pytest.mark.parametrize("ds", [datasets_comodo["ds_2d_left"]]) def test_xgrid_get_axis_dim(ds): grid = XGrid.from_dataset(ds, mesh="flat") assert grid.get_axis_dim("Z") == Z - 1 @@ -108,7 +108,7 @@ def test_invalid_xgrid_field_array(): def test_invalid_lon_lat(): """Stress test the grid initialiser by creating incompatible datasets that test the edge cases""" - ds = datasets["ds_2d_left"].copy() + ds = datasets_comodo["ds_2d_left"].copy() ds["lon"], ds["lat"] = xr.broadcast(ds["YC"], ds["XC"]) with pytest.raises( @@ -117,7 +117,7 @@ def test_invalid_lon_lat(): ): XGrid.from_dataset(ds, mesh="flat") - ds = datasets["ds_2d_left"].copy() + ds = datasets_comodo["ds_2d_left"].copy() ds["lon"], _ = xr.broadcast(ds["YG"], ds["XG"]) with pytest.raises( ValueError, @@ -125,7 +125,7 @@ def test_invalid_lon_lat(): ): XGrid.from_dataset(ds, mesh="flat") - ds = datasets["ds_2d_left"].copy() + ds = datasets_comodo["ds_2d_left"].copy() ds["lon"], ds["lat"] = xr.broadcast(ds["YG"], ds["XG"]) ds["lon"], ds["lat"] = ds["lon"].transpose(), ds["lat"].transpose() @@ -137,7 +137,7 @@ def test_invalid_lon_lat(): def test_invalid_depth(): - ds = datasets["ds_2d_left"].copy() + ds = datasets_comodo["ds_2d_left"].copy() ds = ds.reindex({"ZG": ds.ZG[::-1]}) with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"): @@ -202,8 +202,8 @@ def test_time1D_field(): @pytest.mark.parametrize( "ds", [ - pytest.param(datasets["ds_2d_left"], id="1D lon/lat"), - pytest.param(datasets["2d_left_rotated"], id="2D lon/lat"), + pytest.param(datasets_comodo["ds_2d_left"], id="1D lon/lat"), + pytest.param(datasets_comodo["2d_left_rotated"], id="2D lon/lat"), ], ) # for key, ds in datasets.items()]) def test_xgrid_search_cpoints(ds): @@ -282,7 +282,7 @@ def test_search_1d_array_some_out_of_bounds(array, x, expected_xi): "ds, da_name, expected", [ pytest.param( - datasets["ds_2d_left"], + datasets_comodo["ds_2d_left"], "U_C_grid", { "XG": (np.int64(0), np.float64(0.0)), @@ -292,7 +292,7 @@ def test_search_1d_array_some_out_of_bounds(array, x, expected_xi): id="MITgcm indexing style U_C_grid", ), pytest.param( - datasets["ds_2d_left"], + datasets_comodo["ds_2d_left"], "V_C_grid", { "XC": (np.int64(-1), np.float64(0.5)), @@ -302,7 +302,7 @@ def test_search_1d_array_some_out_of_bounds(array, x, expected_xi): id="MITgcm indexing style V_C_grid", ), pytest.param( - datasets["ds_2d_right"], + datasets_comodo["ds_2d_right"], "U_C_grid", { "XG": (np.int64(0), np.float64(0.0)), @@ -312,7 +312,7 @@ def test_search_1d_array_some_out_of_bounds(array, x, expected_xi): id="NEMO indexing style U_C_grid", ), pytest.param( - datasets["ds_2d_right"], + datasets_comodo["ds_2d_right"], "V_C_grid", { "XC": (np.int64(0), np.float64(0.5)), From 66236a9553c26bedf3e6b1f066aa0eca493eac1b Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 20 May 2026 16:18:30 +0200 Subject: [PATCH 02/10] Update __all__ --- src/parcels/_datasets/structured/generic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_datasets/structured/generic.py b/src/parcels/_datasets/structured/generic.py index 88ab7609bc..5d2da3d63a 100644 --- a/src/parcels/_datasets/structured/generic.py +++ b/src/parcels/_datasets/structured/generic.py @@ -5,7 +5,7 @@ from . import T, X, Y, Z -__all__ = ["T", "X", "Y", "Z", "datasets_comodo"] +__all__ = ["T", "X", "Y", "Z", "datasets_comodo", "datasets_sgrid"] TIME = xr.date_range("2000", "2001", T) From c62007611580b40c1e71b050d1c2fdf49581a3dd Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 21 May 2026 14:10:06 +0200 Subject: [PATCH 03/10] Cleanup test_length1dimensions --- tests/test_advection.py | 84 ++++++++++++++++++----------------------- 1 file changed, 36 insertions(+), 48 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index e5c000132e..f137e48223 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -26,6 +26,7 @@ simple_UV_dataset, stommel_gyre_dataset, ) +from parcels._datasets.structured.generic import datasets_sgrid from parcels.interpolators import CGrid_Velocity, XLinear, XLinear_Velocity from parcels.kernels import ( AdvectionDiffusionEM, @@ -188,61 +189,48 @@ def SubmergeParticle(particles, fieldset): # pragma: no cover assert len(pset) == 0 -@pytest.mark.parametrize("u", [-0.3, np.array(0.2)]) -@pytest.mark.parametrize("v", [0.2, np.array(1)]) -@pytest.mark.parametrize("w", [None, -0.2, np.array(0.7)]) -def test_length1dimensions(u, v, w): # TODO: Refactor this test to be more readable (and isolate test setup) - (lon, xdim) = (np.linspace(-10, 10, 21), 21) if isinstance(u, np.ndarray) else (np.array([0]), 1) - (lat, ydim) = (np.linspace(-15, 15, 31), 31) if isinstance(v, np.ndarray) else (np.array([-4]), 1) - (depth, zdim) = ( - (np.linspace(-5, 5, 11), 11) if (isinstance(w, np.ndarray) and w is not None) else (np.array([3]), 1) - ) - - tdim = 2 # TODO make this also work for length-1 time dimensions - dims = (tdim, zdim, ydim, xdim) - U = u * np.ones(dims, dtype=np.float32) - V = v * np.ones(dims, dtype=np.float32) - if w is not None: - W = w * np.ones(dims, dtype=np.float32) - - ds = xr.Dataset( - { - "U": (["time", "depth", "YG", "XG"], U), - "V": (["time", "depth", "YG", "XG"], V), - }, - coords={ - "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "s")], {"axis": "T"}), - "depth": (["depth"], depth, {"axis": "Z"}), - "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), - "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), - "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), - "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), - "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), - "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), - }, - ) - if w: - ds["W"] = (["time", "depth", "YG", "XG"], W) +@pytest.mark.parametrize( + "u_value, x_slice", [(-0.3, slice(0, 1)), (0.2, slice(None))], ids=["single_u_layer", "full_u"] +) +@pytest.mark.parametrize("v_value, y_slice", [(0.2, slice(0, 1)), (1.0, slice(None))], ids=["single_v_layer", "full_v"]) +@pytest.mark.parametrize( + "w_value, z_slice", + [(None, None), (-0.2, slice(0, 1)), (0.7, slice(None))], + ids=["no_vertical", "single_w_layer", "full_w"], +) +def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice): + ds = datasets_sgrid["ds_2d_padded_high"].copy()[["U_A_grid", "grid"]] + ds = ds.rename({"U_A_grid": "U"}) + ds["U"] = xr.full_like(ds["U"], u_value) + ds["V"] = xr.full_like(ds["U"], v_value) + if w_value is not None: + ds["W"] = xr.full_like(ds["U"], w_value) + + # adjust spatial extent + ds["depth"] -= 5 + ds["lon"] *= 4 + ds["lon"] -= 15 + ds["lat"] *= 4 + ds["lat"] -= 15 + + # Slice dataset + indexers = {"node_dimension1": x_slice, "node_dimension2": y_slice} + if w_value: + indexers.update({"vertical_dimensions_dim1": z_slice}) + ds = ds.sgrid.isel(indexers) - 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, vector_interp_method=XLinear_Velocity)] - if w: - W = Field("W", ds["W"], grid, interp_method=XLinear) - fields.append(VectorField("UVW", U, V, W, vector_interp_method=XLinear_Velocity)) - fieldset = FieldSet(fields) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") x0, y0, z0 = 2, 8, -4 pset = ParticleSet(fieldset, lon=x0, lat=y0, z=z0) - kernel = AdvectionRK4 if w is None else AdvectionRK4_3D + kernel = AdvectionRK4 if w_value is None else AdvectionRK4_3D pset.execute(kernel, runtime=np.timedelta64(4, "s"), dt=np.timedelta64(1, "s")) assert len(pset.lon) == len([p.lon for p in pset]) - np.testing.assert_allclose(np.array([p.lon - x0 for p in pset]), 4 * u, atol=1e-6) - np.testing.assert_allclose(np.array([p.lat - y0 for p in pset]), 4 * v, atol=1e-6) - if w: - np.testing.assert_allclose(np.array([p.z - z0 for p in pset]), 4 * w, atol=1e-6) + np.testing.assert_allclose(np.array([p.lon - x0 for p in pset]), 4 * u_value, atol=1e-6) + np.testing.assert_allclose(np.array([p.lat - y0 for p in pset]), 4 * v_value, atol=1e-6) + if w_value: + np.testing.assert_allclose(np.array([p.z - z0 for p in pset]), 4 * w_value, atol=1e-6) def test_radialrotation(npart=10): From fba5867e8f3c08d21d9b8f74d042423fcf575e9d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 21 May 2026 14:54:00 +0200 Subject: [PATCH 04/10] Fix spatial extent for test Needed to relax the tolerance here. Is that acceptable? Might need to revisit the original datasets defined in generic.py to assure the spatial extents are acceptable --- tests/test_advection.py | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index f137e48223..e00591e029 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -190,12 +190,14 @@ def SubmergeParticle(particles, fieldset): # pragma: no cover @pytest.mark.parametrize( - "u_value, x_slice", [(-0.3, slice(0, 1)), (0.2, slice(None))], ids=["single_u_layer", "full_u"] + "u_value, x_slice", [(-0.03, slice(0, 1)), (0.02, slice(None))], ids=["single_u_layer", "full_u"] +) +@pytest.mark.parametrize( + "v_value, y_slice", [(0.02, slice(0, 1)), (0.1, slice(None))], ids=["single_v_layer", "full_v"] ) -@pytest.mark.parametrize("v_value, y_slice", [(0.2, slice(0, 1)), (1.0, slice(None))], ids=["single_v_layer", "full_v"]) @pytest.mark.parametrize( "w_value, z_slice", - [(None, None), (-0.2, slice(0, 1)), (0.7, slice(None))], + [(None, None), (-0.02, slice(0, 1)), (0.07, slice(None))], ids=["no_vertical", "single_w_layer", "full_w"], ) def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice): @@ -206,13 +208,6 @@ def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice) if w_value is not None: ds["W"] = xr.full_like(ds["U"], w_value) - # adjust spatial extent - ds["depth"] -= 5 - ds["lon"] *= 4 - ds["lon"] -= 15 - ds["lat"] *= 4 - ds["lat"] -= 15 - # Slice dataset indexers = {"node_dimension1": x_slice, "node_dimension2": y_slice} if w_value: @@ -221,16 +216,16 @@ def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") - x0, y0, z0 = 2, 8, -4 + x0, y0, z0 = 3, 3, 20 pset = ParticleSet(fieldset, lon=x0, lat=y0, z=z0) kernel = AdvectionRK4 if w_value is None else AdvectionRK4_3D pset.execute(kernel, runtime=np.timedelta64(4, "s"), dt=np.timedelta64(1, "s")) assert len(pset.lon) == len([p.lon for p in pset]) - np.testing.assert_allclose(np.array([p.lon - x0 for p in pset]), 4 * u_value, atol=1e-6) - np.testing.assert_allclose(np.array([p.lat - y0 for p in pset]), 4 * v_value, atol=1e-6) + np.testing.assert_allclose(np.array([p.lon - x0 for p in pset]), 4 * u_value, atol=1e-5) + np.testing.assert_allclose(np.array([p.lat - y0 for p in pset]), 4 * v_value, atol=1e-5) if w_value: - np.testing.assert_allclose(np.array([p.z - z0 for p in pset]), 4 * w_value, atol=1e-6) + np.testing.assert_allclose(np.array([p.z - z0 for p in pset]), 4 * w_value, atol=1e-5) def test_radialrotation(npart=10): From e0d333c7ce4848d0f4b1f3c69a29841cdbe68189 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 21 May 2026 15:02:27 +0200 Subject: [PATCH 05/10] Update structured/generated.py datasets with sgrid metadata --- src/parcels/_datasets/structured/generated.py | 63 +++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/src/parcels/_datasets/structured/generated.py b/src/parcels/_datasets/structured/generated.py index 665bcaf0f4..a86efab900 100644 --- a/src/parcels/_datasets/structured/generated.py +++ b/src/parcels/_datasets/structured/generated.py @@ -75,6 +75,19 @@ def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square field "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, + ).pipe( + sgrid._attach_sgrid_metadata, + sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + sgrid.FaceNodePadding("XC", "XG", sgrid.Padding.LOW), + sgrid.FaceNodePadding("YC", "YG", sgrid.Padding.HIGH), + ), + vertical_dimensions=(sgrid.FaceNodePadding("ZC", "depth", sgrid.Padding.BOTH),), + ), ) @@ -111,6 +124,19 @@ def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with x "u_g": u_g, "f": f, }, + ).pipe( + sgrid._attach_sgrid_metadata, + sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + sgrid.FaceNodePadding("XC", "XG", sgrid.Padding.LOW), + sgrid.FaceNodePadding("YC", "YG", sgrid.Padding.HIGH), + ), + vertical_dimensions=(sgrid.FaceNodePadding("ZC", "depth", sgrid.Padding.BOTH),), + ), ) @@ -161,6 +187,19 @@ def decaying_moving_eddy_dataset(xdim=2, ydim=2): "gamma": gamma, "gamma_g": gamma_g, }, + ).pipe( + sgrid._attach_sgrid_metadata, + sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + sgrid.FaceNodePadding("XC", "XG", sgrid.Padding.LOW), + sgrid.FaceNodePadding("YC", "YG", sgrid.Padding.HIGH), + ), + vertical_dimensions=(sgrid.FaceNodePadding("ZC", "depth", sgrid.Padding.BOTH),), + ), ) @@ -244,6 +283,18 @@ def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, + ).pipe( + sgrid._attach_sgrid_metadata, + sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + sgrid.FaceNodePadding("XC", "XG", sgrid.Padding.LOW), + sgrid.FaceNodePadding("YC", "YG", sgrid.Padding.HIGH), + ), + ), ) @@ -300,4 +351,16 @@ def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, + ).pipe( + sgrid._attach_sgrid_metadata, + sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + sgrid.FaceNodePadding("XC", "XG", sgrid.Padding.LOW), + sgrid.FaceNodePadding("YC", "YG", sgrid.Padding.HIGH), + ), + ), ) From 02ebfd259e7ea4f65cf6aee41929de3129441767 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 21 May 2026 15:17:00 +0200 Subject: [PATCH 06/10] Update test_advection.py tests to use SGRID ingestion --- tests/test_advection.py | 52 ++++++++++------------------------------- 1 file changed, 12 insertions(+), 40 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index e00591e029..8c27edf44f 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -6,15 +6,12 @@ import parcels import parcels.tutorial from parcels import ( - Field, FieldSet, Particle, ParticleFile, ParticleSet, StatusCode, Variable, - VectorField, - XGrid, convert, ) from parcels._core.utils.time import timedelta_to_float @@ -27,7 +24,6 @@ stommel_gyre_dataset, ) from parcels._datasets.structured.generic import datasets_sgrid -from parcels.interpolators import CGrid_Velocity, XLinear, XLinear_Velocity from parcels.kernels import ( AdvectionDiffusionEM, AdvectionDiffusionM1, @@ -230,11 +226,7 @@ def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice) def test_radialrotation(npart=10): ds = radial_rotation_dataset() - 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, vector_interp_method=XLinear_Velocity) - fieldset = parcels.FieldSet([U, V, UV]) + fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh="flat") dt = np.timedelta64(30, "s") lon = np.linspace(32, 50, npart) @@ -268,22 +260,18 @@ def test_radialrotation(npart=10): ) def test_moving_eddy(kernel, rtol): ds = moving_eddy_dataset() - 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) 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, vector_interp_method=XLinear_Velocity) - fieldset = FieldSet([U, V, W, UVW]) - else: - UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) - fieldset = FieldSet([U, V, UV]) + ds["W"] = ds["V"] + if kernel in [AdvectionDiffusionEM, AdvectionDiffusionM1]: # Add zero diffusivity field for diffusion kernels - ds["Kh"] = (["time", "depth", "YG", "XG"], np.full(ds["U"].shape, 0)) - fieldset.add_field(Field("Kh", ds["Kh"], grid, interp_method=XLinear), "Kh_zonal") - fieldset.add_field(Field("Kh", ds["Kh"], grid, interp_method=XLinear), "Kh_meridional") + ds["Kh_zonal"] = (["time", "depth", "YG", "XG"], np.full(ds["U"].shape, 0)) + ds["Kh_meridional"] = ds["Kh_zonal"] + + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + + if kernel in [AdvectionDiffusionEM, AdvectionDiffusionM1]: fieldset.add_constant("dres", 0.1) start_lon, start_lat, start_z = 12000, 12500, 12500 @@ -322,11 +310,7 @@ def truth_moving(x_0, y_0, t): ) def test_decaying_moving_eddy(kernel, rtol): ds = decaying_moving_eddy_dataset() - 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, vector_interp_method=XLinear_Velocity) - fieldset = FieldSet([U, V, UV]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") start_lon, start_lat = 10000, 10000 dt = np.timedelta64(60, "m") @@ -371,14 +355,7 @@ def truth_moving(x_0, y_0, t): @pytest.mark.parametrize("grid_type", ["A", "C"]) 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 = 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) - UV = VectorField("UV", U, V, vector_interp_method=vector_interp_method) - fieldset = FieldSet([U, V, P, UV]) + fieldset = FieldSet.from_sgrid_conventions(stommel_gyre_dataset(grid_type=grid_type), mesh="flat") dt = np.timedelta64(30, "m") runtime = np.timedelta64(1, "D") @@ -413,12 +390,7 @@ def UpdateP(particles, fieldset): # pragma: no cover def test_peninsula_fieldset(kernel, rtol, grid_type): npart = 2 ds = peninsula_dataset(grid_type=grid_type) - 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) - P = Field("P", ds["P"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) - fieldset = FieldSet([U, V, P, UV]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") dt = np.timedelta64(30, "m") runtime = np.timedelta64(23, "h") From 9a70166f4627c3f1cb49e33dcd24d1ead6124636 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 21 May 2026 17:54:41 +0200 Subject: [PATCH 07/10] Update test_write_fieldset_without_time --- tests/test_particlefile.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index 2f658c55c8..2335455595 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -77,8 +77,7 @@ def test_compression(fieldset, tmp_parquet, compression): def test_write_fieldset_without_time(tmp_parquet): ds = peninsula_dataset() # DataSet without time assert "time" not in ds.dims - grid = XGrid.from_dataset(ds, mesh="flat") - fieldset = FieldSet([Field("U", ds["U"], grid, XLinear)]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") pset = ParticleSet(fieldset, pclass=Particle, lon=0, lat=0) From 9c60e452471472648bc62631a9be6c34c1987a28 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 22 May 2026 11:19:49 +0200 Subject: [PATCH 08/10] Add comment back --- tests/test_advection.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_advection.py b/tests/test_advection.py index 8c27edf44f..9bb383ad7a 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -198,6 +198,7 @@ def SubmergeParticle(particles, fieldset): # pragma: no cover ) def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice): ds = datasets_sgrid["ds_2d_padded_high"].copy()[["U_A_grid", "grid"]] + ds = ds.isel(time=slice(2)) # TODO make this also work for length-1 time dimensions ds = ds.rename({"U_A_grid": "U"}) ds["U"] = xr.full_like(ds["U"], u_value) ds["V"] = xr.full_like(ds["U"], v_value) From 9c772c0e81f0dc05985bf2442aa3104f34c3cf3d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 22 May 2026 09:22:28 +0000 Subject: [PATCH 09/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_advection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index 9bb383ad7a..17511e82e2 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -198,7 +198,7 @@ def SubmergeParticle(particles, fieldset): # pragma: no cover ) def test_length1dimensions(u_value, x_slice, v_value, y_slice, w_value, z_slice): ds = datasets_sgrid["ds_2d_padded_high"].copy()[["U_A_grid", "grid"]] - ds = ds.isel(time=slice(2)) # TODO make this also work for length-1 time dimensions + ds = ds.isel(time=slice(2)) # TODO make this also work for length-1 time dimensions ds = ds.rename({"U_A_grid": "U"}) ds["U"] = xr.full_like(ds["U"], u_value) ds["V"] = xr.full_like(ds["U"], v_value) From 862875ef883b0320d593eb1ed211ce72c2b0715f Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 22 May 2026 14:14:42 +0200 Subject: [PATCH 10/10] Fix docs --- docs/development/posting-issues.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/development/posting-issues.md b/docs/development/posting-issues.md index 0d47aea8ae..ad29ed1efb 100644 --- a/docs/development/posting-issues.md +++ b/docs/development/posting-issues.md @@ -46,8 +46,8 @@ As a user with access to your dataset, you would do: # Generate an example dataset to zip. The user would use their own. import xarray as xr -from parcels._datasets.structured.generic import datasets -datasets['ds_2d_left'].to_netcdf("my_dataset.nc") +from parcels._datasets.structured.generic import datasets_comodo +datasets_comodo['ds_2d_left'].to_netcdf("my_dataset.nc") ``` ```{code-cell}