Skip to content

Using FieldSet.from_sgrid_conventions() in tutorials #2437

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

Open
erikvansebille wants to merge 9 commits into
base: v4-dev
Choose a base branch
from
+73 −139
Open

Using FieldSet.from_sgrid_conventions() in tutorials #2437

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
cc68a74
Using parcels.FieldSet.from_sgrid_conventions() in diffusion tutorial
erikvansebille Dec 16, 2025
88d32c6
Updating simple_UV_dataset to return SGRID compliant dataset
erikvansebille Dec 16, 2025
23dcbe0
Updating unit tests to use simple_UV_dataset with SGRID convention
erikvansebille Dec 16, 2025
5631b23
Using FieldSet.from_sgrid_conventions() for simple_UV_dataset in all …
erikvansebille Dec 16, 2025
692d2dc
Merge branch 'v4-dev' into using_fieldset_from_sgrid_in_tutorials
erikvansebille Dec 16, 2025
ddb2f55
Fixing review comments
erikvansebille Dec 17, 2025
44f7a46
Updating simple_UV_dataset to 2D topology
erikvansebille Dec 17, 2025
d672fba
Fixing type in simple_UV_dataset.lat
erikvansebille Dec 17, 2025
1fcb075
Merge branch 'v4-dev' into using_fieldset_from_sgrid_in_tutorials
erikvansebille Jan 6, 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
17 changes: 2 additions & 15 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,9 +192,8 @@
"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",
"ds"
]
Expand All @@ -205,20 +204,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
17 changes: 7 additions & 10 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\"] = ds[\"U\"] + np.random.rand(1, 1, 5, 4) + 0.1\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`. The default interpolator for fields on structured A-grid grids is (tri)linear interpolation, implemented in `parcels.interpolators.XLinear`."
]
},
{
Expand All @@ -68,12 +68,9 @@
"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\")\n",
"\n",
"assert fieldset.P.interp_method == parcels.interpolators.XLinear"
]
},
{
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
35 changes: 6 additions & 29 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,24 +49,20 @@
"\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",
"ds"
"display(ds)\n",
"\n",
"fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh=\"spherical\")"
]
},
{
"attachments": {},
"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",
"\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",
"```\n",
"\n",
"Plotting the `U` field indeed shows a uniform 1 m/s eastward flow.\n"
]
},
Expand All @@ -76,15 +72,6 @@
"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",
"\n",
"plt.pcolormesh(\n",
" fieldset.U.grid.lon,\n",
" fieldset.U.grid.lat,\n",
Expand Down Expand Up @@ -208,21 +195,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
12 changes: 12 additions & 0 deletions src/parcels/_core/utils/sgrid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -432,6 +432,18 @@ def _metadata_rename_dims(grid: Grid2DMetadata, dims_dict: dict[str, str]) -> Gr
def _metadata_rename_dims(grid: Grid3DMetadata, dims_dict: dict[str, str]) -> Grid3DMetadata: ...


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 _metadata_rename_dims(grid, dims_dict):
"""
Renames dimensions in SGrid metadata.
Expand Down
20 changes: 19 additions & 1 deletion src/parcels/_datasets/structured/generated.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,12 @@
import numpy as np
import xarray as xr

from parcels._core.utils.sgrid import (
DimDimPadding,
Grid2DMetadata,
Padding,
_attach_sgrid_metadata,
)
from parcels._core.utils.time import timedelta_to_float


Expand All @@ -18,9 +24,21 @@ 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(-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.HIGH),
DimDimPadding("YC", "YG", Padding.HIGH),
),
vertical_dimensions=(DimDimPadding("depth", "depth", Padding.NONE),),
),
)


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,8 +4,8 @@
from parcels._core.utils.sgrid import (
DimDimPadding,
Grid2DMetadata,
Grid3DMetadata,
Padding,
_attach_sgrid_metadata,
)
from parcels._core.utils.sgrid import (
rename_dims as sgrid_rename_dims,
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
39 changes: 10 additions & 29 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])
fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat")
fieldset.U.data[:, 0, :, :] = 0.0 # Set U to 0 at the surface

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,12 @@ 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
fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat")
fieldset.U.data[:] = 0.01 # Set U to small value (to avoid horizontal out of bounds)
W = Field("W", ds["V"], fieldset.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])
UVW = VectorField("UVW", fieldset.U, fieldset.V, W)
fieldset = FieldSet([fieldset.U, fieldset.V, W, UVW, fieldset.UV])

def DeleteParticle(particles, fieldset): # pragma: no cover
particles.state = np.where(particles.state == StatusCode.ErrorOutOfBounds, StatusCode.Delete, particles.state)
Expand Down
20 changes: 6 additions & 14 deletions tests/test_diffusion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,6 @@
ParticleSet,
Unity,
Variable,
VectorField,
XGrid,
)
from parcels._core.utils.time import timedelta_to_float
from parcels._datasets.structured.generated import simple_UV_dataset
Expand All @@ -32,11 +30,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)

Expand Down Expand Up @@ -74,20 +68,18 @@ 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)
fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh)

Kh = np.zeros((ydim, xdim), dtype=np.float32)
for x in range(xdim):
Kh[:, x] = np.tanh(ds["lon"][x] / ds["lon"][-1] * 10.0) * xdim / 2.0 + xdim / 2.0 + 100.0

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])
Kh_zonal = Field("Kh_zonal", ds["Kh_zonal"], grid=fieldset.U.grid, interp_method=XLinear)
Kh_meridional = Field("Kh_meridional", ds["Kh_meridional"], grid=fieldset.V.grid, interp_method=XLinear)
fieldset.add_field(Kh_zonal)
fieldset.add_field(Kh_meridional)
fieldset.add_constant("dres", float(ds["lon"][1] - ds["lon"][0]))

npart = 10000
Expand Down
Loading
Loading