diff --git a/docs/user_guide/examples/tutorial_diffusion.ipynb b/docs/user_guide/examples/tutorial_diffusion.ipynb index faae27bac..f9a60fef1 100644 --- a/docs/user_guide/examples/tutorial_diffusion.ipynb +++ b/docs/user_guide/examples/tutorial_diffusion.ipynb @@ -192,7 +192,7 @@ "source": [ "from parcels._datasets.structured.generated import simple_UV_dataset\n", "\n", - "ds = simple_UV_dataset(dims=(1, 1, Ny, 1), mesh=\"flat\").isel(time=0, depth=0)\n", + "ds = simple_UV_dataset(dims=(1, 1, Ny, 1), mesh=\"flat\")\n", "ds[\"lat\"][:] = np.linspace(-0.01, 1.01, Ny)\n", "ds[\"lon\"][:] = np.ones(len(ds.XG))\n", "ds[\"Kh_meridional\"] = ([\"YG\", \"XG\"], Kh_meridional[:, None])\n", @@ -205,20 +205,8 @@ "metadata": {}, "outputs": [], "source": [ - "grid = parcels.XGrid.from_dataset(ds, mesh=\"flat\")\n", - "U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "UV = parcels.VectorField(\"UV\", U, V)\n", - "\n", - "Kh_meridional_field = parcels.Field(\n", - " \"Kh_meridional\",\n", - " ds[\"Kh_meridional\"],\n", - " grid,\n", - " interp_method=parcels.interpolators.XLinear,\n", - ")\n", - "fieldset = parcels.FieldSet([U, V, UV, Kh_meridional_field])\n", + "fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh=\"flat\")\n", "fieldset.add_constant_field(\"Kh_zonal\", 1, mesh=\"flat\")\n", - "\n", "fieldset.add_constant(\"dres\", 0.00005)" ] }, diff --git a/docs/user_guide/examples/tutorial_interpolation.ipynb b/docs/user_guide/examples/tutorial_interpolation.ipynb index b64ed72c7..afced17c9 100644 --- a/docs/user_guide/examples/tutorial_interpolation.ipynb +++ b/docs/user_guide/examples/tutorial_interpolation.ipynb @@ -46,12 +46,12 @@ "source": [ "from parcels._datasets.structured.generated import simple_UV_dataset\n", "\n", - "ds = simple_UV_dataset(dims=(1, 1, 5, 4), mesh=\"flat\").isel(time=0, depth=0)\n", + "ds = simple_UV_dataset(dims=(1, 1, 5, 4), mesh=\"flat\")\n", "ds[\"lat\"][:] = np.linspace(0.0, 1.0, len(ds.YG))\n", "ds[\"lon\"][:] = np.linspace(0.0, 1.0, len(ds.XG))\n", "dx, dy = 1.0 / len(ds.XG), 1.0 / len(ds.YG)\n", "ds[\"P\"] = ds[\"U\"] + np.random.rand(5, 4) + 0.1\n", - "ds[\"P\"][1, 1] = 0\n", + "ds[\"P\"][:, :, 1, 1] = 0\n", "ds" ] }, @@ -59,7 +59,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "From this dataset we create a {py:obj}`parcels.FieldSet`. Parcels requires an interpolation method to be set for each {py:obj}`parcels.Field`, which we will later adapt to see the effects of the different interpolators. A common interpolator for fields on structured grids is (tri)linear, implemented in {py:obj}`parcels.interpolators.XLinear`." + "From this dataset we create a {py:obj}`parcels.FieldSet` using the {py:meth}`parcels.FieldSet.from_sgrid_conventions` constructor, which automatically sets up the grid and fields according to Parcels' s-grid conventions." ] }, { @@ -68,12 +68,7 @@ "metadata": {}, "outputs": [], "source": [ - "grid = parcels.XGrid.from_dataset(ds, mesh=\"flat\")\n", - "U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "UV = parcels.VectorField(\"UV\", U, V)\n", - "P = parcels.Field(\"P\", ds[\"P\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "fieldset = parcels.FieldSet([U, V, UV, P])" + "fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh=\"flat\")" ] }, { diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index 7cc465100..9ce434822 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -440,12 +440,7 @@ "source": [ "fields = [GridID]\n", "for i, ds in enumerate(ds_in):\n", - " # TODO : use FieldSet.from_sgrid_convetion here once #2437 is merged\n", - " grid = parcels.XGrid.from_dataset(ds, mesh=\"spherical\")\n", - " U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", - " V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", - " UV = parcels.VectorField(\"UV\", U, V)\n", - " fset = parcels.FieldSet([U, V, UV])\n", + " fset = parcels.FieldSet.from_sgrid_conventions(ds, mesh=\"spherical\")\n", "\n", " for fld in fset.fields.values():\n", " fld.name = f\"{fld.name}{i}\"\n", @@ -469,32 +464,19 @@ "outputs": [], "source": [ "def AdvectEE_NestedGrids(particles, fieldset):\n", - " particles.gridID = fieldset.GridID[particles]\n", - "\n", - " # TODO because of KernelParticle bug (GH #2143), we need to copy lon/lat/time to local variables\n", - " time = particles.time\n", - " z = particles.z\n", - " lat = particles.lat\n", - " lon = particles.lon\n", " u = np.zeros_like(particles.lon)\n", " v = np.zeros_like(particles.lat)\n", "\n", + " particles.gridID = fieldset.GridID[particles]\n", " unique_ids = np.unique(particles.gridID)\n", " for gid in unique_ids:\n", " mask = particles.gridID == gid\n", " UVField = getattr(fieldset, f\"UV{gid}\")\n", - " (u[mask], v[mask]) = UVField[time[mask], z[mask], lat[mask], lon[mask]]\n", + " (u[mask], v[mask]) = UVField[particles[mask]]\n", "\n", " particles.dlon += u * particles.dt\n", " particles.dlat += v * particles.dt\n", "\n", - " # TODO particle states have to be updated manually because UVField is not called with `particles` argument (becaise of GH #2143)\n", - " particles.state = np.where(\n", - " np.isnan(u) | np.isnan(v),\n", - " parcels.StatusCode.ErrorInterpolation,\n", - " particles.state,\n", - " )\n", - "\n", "\n", "lat = np.linspace(-17, 35, 10)\n", "lon = np.full(len(lat), -5)\n", diff --git a/docs/user_guide/examples/tutorial_statuscodes.md b/docs/user_guide/examples/tutorial_statuscodes.md index b887f6ea6..557cb8063 100644 --- a/docs/user_guide/examples/tutorial_statuscodes.md +++ b/docs/user_guide/examples/tutorial_statuscodes.md @@ -63,19 +63,14 @@ Let's add the `KeepInOcean` Kernel to an particle simulation where particles mov import numpy as np from parcels._datasets.structured.generated import simple_UV_dataset -ds = simple_UV_dataset(dims=(1, 2, 5, 4), mesh="flat").isel(time=0) +ds = simple_UV_dataset(dims=(1, 2, 5, 4), mesh="flat") dx, dy = 1.0 / len(ds.XG), 1.0 / len(ds.YG) # Add W velocity that pushes through surface ds["W"] = ds["U"] - 0.1 # 0.1 m/s towards the surface -grid = parcels.XGrid.from_dataset(ds, mesh="flat") -U = parcels.Field("U", ds["U"], grid, interp_method=parcels.interpolators.XLinear) -V = parcels.Field("V", ds["V"], grid, interp_method=parcels.interpolators.XLinear) -W = parcels.Field("W", ds["W"], grid, interp_method=parcels.interpolators.XLinear) -UVW = parcels.VectorField("UVW", U, V, W) -fieldset = parcels.FieldSet([U, V, W, UVW]) +fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh="flat") ``` If we advect particles with the `AdvectionRK2_3D` kernel, Parcels will raise a `FieldOutOfBoundSurfaceError`: diff --git a/docs/user_guide/examples/tutorial_unitconverters.ipynb b/docs/user_guide/examples/tutorial_unitconverters.ipynb index a844833ea..773c339b8 100644 --- a/docs/user_guide/examples/tutorial_unitconverters.ipynb +++ b/docs/user_guide/examples/tutorial_unitconverters.ipynb @@ -49,7 +49,7 @@ "\n", "nlat = 10\n", "nlon = 18\n", - "ds = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"spherical\").isel(time=0, depth=0)\n", + "ds = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"spherical\")\n", "ds[\"temperature\"] = ds[\"U\"] + 20 # add temperature field of 20 deg\n", "ds[\"U\"].data[:] = 1.0 # set U to 1 m/s\n", "ds[\"V\"].data[:] = 1.0 # set V to 1 m/s\n", @@ -61,7 +61,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To create a `parcels.FieldSet` object, we define the `parcels.Field`s and the structured grid (`parcels.XGrid`) the fields are defined on. We add the argument `mesh='spherical'` to the `parcels.XGrid` to signal that all longitudes and latitudes are in degrees.\n", + "To create a `parcels.FieldSet` object, we use the `parcels.FieldSet.from_sgrid_conventions` constructor. We add the argument `mesh='spherical'` to signal that all longitudes and latitudes are in degrees.\n", "\n", "```{note}\n", "When using a `FieldSet` method for a specific dataset, such as `from_copernicusmarine()`, the grid information is known and parsed by Parcels, so we do not have to add the `mesh` argument.\n", @@ -76,14 +76,7 @@ "metadata": {}, "outputs": [], "source": [ - "grid = parcels.XGrid.from_dataset(ds, mesh=\"spherical\")\n", - "U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "UV = parcels.VectorField(\"UV\", U, V)\n", - "temperature = parcels.Field(\n", - " \"temperature\", ds[\"temperature\"], grid, interp_method=parcels.interpolators.XLinear\n", - ")\n", - "fieldset = parcels.FieldSet([U, V, UV, temperature])\n", + "fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh=\"spherical\")\n", "\n", "plt.pcolormesh(\n", " fieldset.U.grid.lon,\n", @@ -199,7 +192,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the XGrid object.\n" + "If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the `FieldSet` object.\n" ] }, { @@ -208,21 +201,11 @@ "metadata": {}, "outputs": [], "source": [ - "ds_flat = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"flat\").isel(time=0, depth=0)\n", + "ds_flat = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"flat\")\n", "ds_flat[\"temperature\"] = ds_flat[\"U\"] + 20 # add temperature field of 20 deg\n", "ds_flat[\"U\"].data[:] = 1.0 # set U to 1 m/s\n", "ds_flat[\"V\"].data[:] = 1.0 # set V to 1 m/s\n", - "grid = parcels.XGrid.from_dataset(ds_flat, mesh=\"flat\")\n", - "U = parcels.Field(\"U\", ds_flat[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "V = parcels.Field(\"V\", ds_flat[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", - "UV = parcels.VectorField(\"UV\", U, V)\n", - "temperature = parcels.Field(\n", - " \"temperature\",\n", - " ds_flat[\"temperature\"],\n", - " grid,\n", - " interp_method=parcels.interpolators.XLinear,\n", - ")\n", - "fieldset_flat = parcels.FieldSet([U, V, UV, temperature])\n", + "fieldset_flat = parcels.FieldSet.from_sgrid_conventions(ds_flat, mesh=\"flat\")\n", "\n", "plt.pcolormesh(\n", " fieldset_flat.U.grid.lon,\n", diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 1477479d8..f9bcb795f 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -351,12 +351,11 @@ 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"]) if "W" in ds.data_vars: fields["W"] = Field("W", ds["W"], grid, XLinear) fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) - else: - fields["UV"] = VectorField("UV", fields["U"], fields["V"]) 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/_datasets/structured/generated.py b/src/parcels/_datasets/structured/generated.py index 8b7f8af08..20440717c 100644 --- a/src/parcels/_datasets/structured/generated.py +++ b/src/parcels/_datasets/structured/generated.py @@ -3,7 +3,13 @@ import numpy as np import xarray as xr +from parcels._core.utils.sgrid import ( + DimDimPadding, + Grid2DMetadata, + Padding, +) from parcels._core.utils.time import timedelta_to_float +from parcels._datasets.utils import _attach_sgrid_metadata def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh="spherical"): @@ -21,6 +27,18 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh="spherical"): "lat": (["YG"], np.linspace(-90, 90, 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( + _attach_sgrid_metadata, + Grid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + face_dimensions=( + DimDimPadding("XC", "XG", Padding.LOW), + DimDimPadding("YC", "YG", Padding.LOW), + ), + vertical_dimensions=(DimDimPadding("ZC", "depth", Padding.BOTH),), + ), ) diff --git a/src/parcels/_datasets/structured/generic.py b/src/parcels/_datasets/structured/generic.py index 1290d99d5..43302635d 100644 --- a/src/parcels/_datasets/structured/generic.py +++ b/src/parcels/_datasets/structured/generic.py @@ -4,12 +4,12 @@ from parcels._core.utils.sgrid import ( DimDimPadding, Grid2DMetadata, - Grid3DMetadata, Padding, ) from parcels._core.utils.sgrid import ( rename_dims as sgrid_rename_dims, ) +from parcels._datasets.utils import _attach_sgrid_metadata from . import T, X, Y, Z @@ -18,18 +18,6 @@ TIME = xr.date_range("2000", "2001", T) -def _attach_sgrid_metadata(ds, grid: Grid2DMetadata | Grid3DMetadata): - """Copies the dataset and attaches the SGRID metadata in 'grid' variable. Modifies 'conventions' attribute.""" - ds = ds.copy() - ds["grid"] = ( - [], - 0, - grid.to_attrs(), - ) - ds.attrs["Conventions"] = "SGRID" - return ds - - def _rotated_curvilinear_grid(): XG = np.arange(X) YG = np.arange(Y) diff --git a/src/parcels/_datasets/utils.py b/src/parcels/_datasets/utils.py index 15afeca46..0c475fac5 100644 --- a/src/parcels/_datasets/utils.py +++ b/src/parcels/_datasets/utils.py @@ -4,9 +4,26 @@ import numpy as np import xarray as xr +from parcels._core.utils.sgrid import ( + Grid2DMetadata, + Grid3DMetadata, +) + _SUPPORTED_ATTR_TYPES = int | float | str | np.ndarray +def _attach_sgrid_metadata(ds, grid: Grid2DMetadata | Grid3DMetadata): + """Copies the dataset and attaches the SGRID metadata in 'grid' variable. Modifies 'conventions' attribute.""" + ds = ds.copy() + ds["grid"] = ( + [], + 0, + grid.to_attrs(), + ) + ds.attrs["Conventions"] = "SGRID" + return ds + + def _print_mismatched_keys(d1: dict[Any, Any], d2: dict[Any, Any]) -> None: k1 = set(d1.keys()) k2 = set(d2.keys()) diff --git a/tests/test_advection.py b/tests/test_advection.py index 95b007e76..0a0b7484a 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -33,11 +33,7 @@ def test_advection_zonal(mesh, npart=10): """Particles at high latitude move geographically faster due to the pole correction in `GeographicPolar`.""" ds = simple_UV_dataset(mesh=mesh) ds["U"].data[:] = 1.0 - grid = XGrid.from_dataset(ds, mesh=mesh) - U = Field("U", ds["U"], grid, interp_method=XLinear) - V = Field("V", ds["V"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) - fieldset = FieldSet([U, V, UV]) + 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")) @@ -53,11 +49,7 @@ def test_advection_zonal_with_particlefile(tmp_store): npart = 10 ds = simple_UV_dataset(mesh="flat") ds["U"].data[:] = 1.0 - 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) - fieldset = FieldSet([U, V, UV]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") pset = ParticleSet(fieldset, lon=np.zeros(npart) + 20.0, lat=np.linspace(0, 80, npart)) pfile = ParticleFile(tmp_store, outputdt=np.timedelta64(30, "m")) @@ -85,11 +77,7 @@ def test_advection_zonal_periodic(): halo.XG.values = ds.XG.values[1] + 2 ds = xr.concat([ds, halo], dim="XG") - 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) - fieldset = FieldSet([U, V, UV]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") PeriodicParticle = Particle.add_variable(Variable("total_dlon", initial=0)) startlon = np.array([0.5, 0.4]) @@ -104,12 +92,8 @@ 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") ds["U"].data[:] = 1.0 - grid = XGrid.from_dataset(ds, mesh="flat") - U = Field("U", ds["U"], grid, interp_method=XLinear) - U.data[:, 0, :, :] = 0.0 # Set U to 0 at the surface - V = Field("V", ds["V"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) - fieldset = FieldSet([U, V, UV]) + ds["U"].data[:, 0, :, :] = 0.0 # Set U to 0 at the surface + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") 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")) @@ -122,15 +106,10 @@ def test_horizontal_advection_in_3D_flow(npart=10): @pytest.mark.parametrize("wErrorThroughSurface", [True, False]) def test_advection_3D_outofbounds(direction, wErrorThroughSurface): ds = simple_UV_dataset(mesh="flat") - grid = XGrid.from_dataset(ds, mesh="flat") - U = Field("U", ds["U"], grid, interp_method=XLinear) - U.data[:] = 0.01 # Set U to small value (to avoid horizontal out of bounds) - V = Field("V", ds["V"], grid, interp_method=XLinear) - W = Field("W", ds["V"], grid, interp_method=XLinear) # Use V as W for testing - W.data[:] = -1.0 if direction == "up" else 1.0 - UVW = VectorField("UVW", U, V, W) - UV = VectorField("UV", U, V) - fieldset = FieldSet([U, V, W, UVW, UV]) + ds["W"] = ds["V"].copy() # Just to have W field present + ds["U"].data[:] = 0.01 # Set U to small value (to avoid horizontal out of bounds) + ds["W"].data[:] = -1.0 if direction == "up" else 1.0 + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") def DeleteParticle(particles, fieldset): # pragma: no cover particles.state = np.where(particles.state == StatusCode.ErrorOutOfBounds, StatusCode.Delete, particles.state) diff --git a/tests/test_diffusion.py b/tests/test_diffusion.py index f06ceba27..d392c780f 100644 --- a/tests/test_diffusion.py +++ b/tests/test_diffusion.py @@ -5,7 +5,6 @@ from scipy import stats from parcels import ( - Field, FieldSet, GeographicPolarSquare, GeographicSquare, @@ -13,12 +12,9 @@ ParticleSet, Unity, Variable, - VectorField, - XGrid, ) from parcels._core.utils.time import timedelta_to_float from parcels._datasets.structured.generated import simple_UV_dataset -from parcels.interpolators import XLinear from parcels.kernels import AdvectionDiffusionEM, AdvectionDiffusionM1, DiffusionUniformKh from tests.utils import create_fieldset_zeros_conversion @@ -32,11 +28,7 @@ def test_fieldKh_Brownian(mesh): ds = simple_UV_dataset(dims=(2, 1, 2, 2), mesh=mesh) ds["lon"].data = np.array([-1e6, 1e6]) ds["lat"].data = np.array([-1e6, 1e6]) - grid = XGrid.from_dataset(ds, mesh=mesh) - U = Field("U", ds["U"], grid, interp_method=XLinear) - V = Field("V", ds["V"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) - fieldset = FieldSet([U, V, UV]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) fieldset.add_constant_field("Kh_zonal", kh_zonal, mesh=mesh) fieldset.add_constant_field("Kh_meridional", kh_meridional, mesh=mesh) @@ -74,9 +66,6 @@ def test_fieldKh_SpatiallyVaryingDiffusion(mesh, kernel): ds = simple_UV_dataset(dims=(2, 1, ydim, xdim), mesh=mesh) ds["lon"].data = np.linspace(-1e6, 1e6, xdim) ds["lat"].data = np.linspace(-1e6, 1e6, ydim) - grid = XGrid.from_dataset(ds, mesh=mesh) - U = Field("U", ds["U"], grid, interp_method=XLinear) - V = Field("V", ds["V"], grid, interp_method=XLinear) Kh = np.zeros((ydim, xdim), dtype=np.float32) for x in range(xdim): @@ -84,10 +73,7 @@ def test_fieldKh_SpatiallyVaryingDiffusion(mesh, kernel): ds["Kh_zonal"] = (["time", "depth", "YG", "XG"], np.full((2, 1, ydim, xdim), Kh)) ds["Kh_meridional"] = (["time", "depth", "YG", "XG"], np.full((2, 1, ydim, xdim), Kh)) - Kh_zonal = Field("Kh_zonal", ds["Kh_zonal"], grid=grid, interp_method=XLinear) - Kh_meridional = Field("Kh_meridional", ds["Kh_meridional"], grid=grid, interp_method=XLinear) - UV = VectorField("UV", U, V) - fieldset = FieldSet([U, V, UV, Kh_zonal, Kh_meridional]) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) fieldset.add_constant("dres", float(ds["lon"][1] - ds["lon"][0])) npart = 10000 diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index efb3d0092..131826647 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -317,7 +317,7 @@ def test_fieldset_from_copernicusmarine_with_W(caplog): assert "U" in fieldset.fields assert "V" in fieldset.fields assert "W" in fieldset.fields - assert "UV" not in fieldset.fields + assert "UV" in fieldset.fields assert "UVW" in fieldset.fields assert "renamed it to 'W'" in caplog.text diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index c282fe7f0..39b588f2b 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -166,23 +166,20 @@ def test_invdistland_interpolation(field, func, t, z, y, x, expected): def test_interpolation_mesh_type(mesh, npart=10): ds = simple_UV_dataset(mesh=mesh) ds["U"].data[:] = 1.0 - grid = XGrid.from_dataset(ds, mesh=mesh) - U = Field("U", ds["U"], grid, interp_method=XLinear) - V = Field("V", ds["V"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) lat = 30.0 time = 0.0 u_expected = 1.0 if mesh == "flat" else 1.0 / (1852 * 60 * np.cos(np.radians(lat))) - assert np.isclose(U[time, 0, lat, 0], u_expected, atol=1e-7) - assert V[time, 0, lat, 0] == 0.0 + assert np.isclose(fieldset.U[time, 0, lat, 0], u_expected, atol=1e-7) + assert fieldset.V[time, 0, lat, 0] == 0.0 - u, v = UV[time, 0, lat, 0] + u, v = fieldset.UV[time, 0, lat, 0] assert np.isclose(u, u_expected, atol=1e-7) assert v == 0.0 - assert U.eval(time, 0, lat, 0, applyConversion=False) == 1 + assert fieldset.U.eval(time, 0, lat, 0, applyConversion=False) == 1.0 interp_methods = { diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index c998677bb..b989356dc 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -55,11 +55,7 @@ def fieldset_no_time_interval() -> FieldSet: def zonal_flow_fieldset() -> FieldSet: ds = simple_UV_dataset(mesh="flat") ds["U"].data[:] = 1.0 - 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) - return FieldSet([U, V, UV]) + return FieldSet.from_sgrid_conventions(ds, mesh="flat") def test_pset_execute_invalid_arguments(fieldset, fieldset_no_time_interval): diff --git a/tests/utils.py b/tests/utils.py index 4ae8a083b..62420a638 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -10,10 +10,9 @@ import xarray as xr import parcels -from parcels import Field, FieldSet, Particle, Variable, VectorField +from parcels import FieldSet, Particle, Variable from parcels._core.xgrid import _FIELD_DATA_ORDERING, XGrid, get_axis_from_dim_name from parcels._datasets.structured.generated import simple_UV_dataset -from parcels.interpolators import XLinear PROJECT_ROOT = Path(__file__).resolve().parents[1] TEST_ROOT = PROJECT_ROOT / "tests" @@ -78,12 +77,7 @@ def create_fieldset_zeros_conversion(mesh="spherical", xdim=200, ydim=100) -> Fi ds = simple_UV_dataset(dims=(2, 1, ydim, xdim), mesh=mesh) ds["lon"].data = np.linspace(-1e6 * mesh_conversion, 1e6 * mesh_conversion, xdim) ds["lat"].data = np.linspace(-1e6 * mesh_conversion, 1e6 * mesh_conversion, ydim) - grid = XGrid.from_dataset(ds, mesh=mesh) - U = Field("U", ds["U"], grid, interp_method=XLinear) - V = Field("V", ds["V"], grid, interp_method=XLinear) - - UV = VectorField("UV", U, V) - return FieldSet([U, V, UV]) + return FieldSet.from_sgrid_conventions(ds, mesh=mesh) def create_simple_pset(n=1):