Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
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
Original file line number Diff line number Diff line change
Expand Up @@ -40,20 +40,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 `parcels.FieldSet`. Parcels requires an interpolation method to be set for each `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 `parcels.interpolators.XLinear`."
"From this dataset we create a `parcels.FieldSet`. The default interpolator for fields on structured A-grid grids is (tri)linear interpolation, implemented in `parcels.interpolators.XLinear`."
]
},
{
Expand All @@ -62,12 +62,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
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
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
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
18 changes: 18 additions & 0 deletions src/parcels/_datasets/structured/generated.py
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 @@ -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.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
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
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
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
13 changes: 5 additions & 8 deletions tests/test_interpolation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


interp_methods = {
Expand Down
Loading
Loading