Skip to content

Using from_sgrid_conventions for simple_UV_dataset throughout #2453

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
erikvansebille merged 6 commits into from
Jan 9, 2026
Merged

Using from_sgrid_conventions for simple_UV_dataset throughout #2453

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
295cfef
Cleaning AdvectEE_NestedGrids
erikvansebille Jan 7, 2026
99c550b
Always attach a UV Field
erikvansebille Jan 2, 2026
45b54df
Fixing test
erikvansebille Jan 2, 2026
a43573c
Using from_sgrid_conventions for fieldset creation
erikvansebille Jan 7, 2026
f9053b3
Update tutorial_unitconverters.ipynb
erikvansebille Jan 7, 2026
3627866
Merge branch 'v4-dev' into use_sgrid_for_simple_UV_dataset
erikvansebille Jan 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 2 additions & 14 deletions docs/user_guide/examples/tutorial_diffusion.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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)"
]
},
Expand Down
13 changes: 4 additions & 9 deletions docs/user_guide/examples/tutorial_interpolation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,20 +46,20 @@
"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"
]
},
{
"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."
]
},
{
Expand All @@ -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\")"
]
},
{
Expand Down
24 changes: 3 additions & 21 deletions docs/user_guide/examples/tutorial_nestedgrids.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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",
Expand Down
9 changes: 2 additions & 7 deletions docs/user_guide/examples/tutorial_statuscodes.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`:
Expand Down
29 changes: 6 additions & 23 deletions docs/user_guide/examples/tutorial_unitconverters.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -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"
]
},
{
Expand All @@ -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",
Expand Down
3 changes: 1 addition & 2 deletions src/parcels/_core/fieldset.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
18 changes: 18 additions & 0 deletions src/parcels/_datasets/structured/generated.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"):
Expand All @@ -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),),
),
)


Expand Down
14 changes: 1 addition & 13 deletions src/parcels/_datasets/structured/generic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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)
Expand Down
17 changes: 17 additions & 0 deletions src/parcels/_datasets/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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())
Expand Down
39 changes: 9 additions & 30 deletions tests/test_advection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand All @@ -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"))
Expand Down Expand Up @@ -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])
Expand All @@ -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"))
Expand All @@ -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)
Expand Down
Loading
Loading