Skip to content
Merged
Show file tree
Hide file tree
Changes from 4 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
2 changes: 1 addition & 1 deletion docs/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -521,7 +521,7 @@ def linkcode_resolve(domain, info):
# -- Options for MyST parser ----------------------------------------------
myst_heading_anchors = 3

myst_enable_extensions = ["substitution"]
myst_enable_extensions = ["substitution", "amsmath", "dollarmath"]

# -- Options for MyST-nb --------------------------------------------------
nb_execution_mode = "cache"
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,15 @@
"source": [
"In most applications, Parcels works with `spherical` meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use `flat` meshes, where longitude and latitude are given in meters (note that the dimensions are then still called `longitude` and `latitude` for consistency reasons).\n",
"\n",
"In all cases, velocities are given in m/s. So Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this tutorial explain how this works.\n"
"In all cases, velocities are given in m/s. Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this tutorial explains how this works.\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's first import the relevant modules, and create dictionaries for the `U`, `V` and `temp` data arrays, with the velocities 1 m/s and the temperature 20C.\n"
"Let's first import the relevant modules, and generate a simple dataset on a 2D spherical mesh, with `U`, `V` and `temperature` data arrays, with the velocities 1 m/s and the temperature 20C.\n"
]
},
{
Expand All @@ -34,6 +34,7 @@
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import xarray as xr\n",
"\n",
"import parcels"
]
Expand All @@ -44,24 +45,27 @@
"metadata": {},
"outputs": [],
"source": [
"xdim, ydim = (10, 20)\n",
"data = {\n",
" \"U\": np.ones((ydim, xdim), dtype=np.float32),\n",
" \"V\": np.ones((ydim, xdim), dtype=np.float32),\n",
" \"temp\": 20 * np.ones((ydim, xdim), dtype=np.float32),\n",
"}\n",
"dims = {\n",
" \"lon\": np.linspace(-15, 5, xdim, dtype=np.float32),\n",
" \"lat\": np.linspace(35, 60, ydim, dtype=np.float32),\n",
"}"
"from parcels._datasets.structured.generated import simple_UV_dataset\n",
"\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[\"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"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We can convert these data and dims to a FieldSet object using `FieldSet.from_data`. We add the argument `mesh='spherical'` (this is the default option) to signal that all longitudes and latitudes are in degrees.\n",
"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 1m/s eastward flow.\n"
]
Expand All @@ -72,8 +76,25 @@
"metadata": {},
"outputs": [],
"source": [
"fieldset = parcels.FieldSet.from_data(data, dims, mesh=\"spherical\")\n",
"plt.pcolormesh(fieldset.U.lon, fieldset.U.lat, fieldset.U.data[0, :, :], vmin=0, vmax=1)\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",
"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",
" fieldset.U.data[0, 0, :, :],\n",
" vmin=0,\n",
" vmax=1,\n",
" shading=\"nearest\",\n",
")\n",
"plt.ylabel(\"Latitude\")\n",
"plt.xlabel(\"Longitude\")\n",
"plt.colorbar()\n",
"plt.show()"
]
Expand All @@ -92,8 +113,10 @@
"metadata": {},
"outputs": [],
"source": [
"print(fieldset.UV[0, 0, 40, -5])\n",
"print(fieldset.temp[0, 0, 40, -5])"
"print(fieldset.UV[np.array([0]), np.array([0]), np.array([40]), np.array([-5])])\n",
"print(\n",
" fieldset.temperature[np.array([0]), np.array([0]), np.array([40]), np.array([-5])]\n",
")"
]
},
{
Expand All @@ -112,7 +135,7 @@
"metadata": {},
"outputs": [],
"source": [
"for fld in [fieldset.U, fieldset.V, fieldset.temp]:\n",
"for fld in [fieldset.U, fieldset.V, fieldset.temperature]:\n",
" print(f\"{fld.name}: {fld.units}\")"
]
},
Expand All @@ -132,7 +155,7 @@
"metadata": {},
"outputs": [],
"source": [
"u, v = fieldset.UV[0, 0, 40, -5]\n",
"u, v = fieldset.UV[np.array([0]), np.array([0]), np.array([40]), np.array([-5])]\n",
"print(v * 1852 * 60)"
]
},
Expand All @@ -150,7 +173,15 @@
"metadata": {},
"outputs": [],
"source": [
"print(fieldset.UV.eval(0, 0, 40, -5, applyConversion=False))"
"print(\n",
" fieldset.UV.eval(\n",
" np.array([0]),\n",
" np.array([0]),\n",
" np.array([40]),\n",
" np.array([-5]),\n",
" applyConversion=False,\n",
" )\n",
")"
]
},
{
Expand All @@ -166,7 +197,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the FieldSet object.\n"
"If longitudes and latitudes are given in meters, rather than degrees, simply add `mesh='flat'` when creating the XGrid object.\n"
]
},
{
Expand All @@ -175,20 +206,38 @@
"metadata": {},
"outputs": [],
"source": [
"fieldset_flat = parcels.FieldSet.from_data(data, dims, mesh=\"flat\")\n",
"ds_flat = simple_UV_dataset(dims=(1, 1, nlat, nlon), mesh=\"flat\").isel(time=0, depth=0)\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",
"\n",
"plt.pcolormesh(\n",
" fieldset_flat.U.lon,\n",
" fieldset_flat.U.lat,\n",
" fieldset_flat.U.data[0, :, :],\n",
" fieldset_flat.U.grid.lon,\n",
" fieldset_flat.U.grid.lat,\n",
" fieldset_flat.U.data[0, 0, :, :],\n",
" vmin=0,\n",
" vmax=1,\n",
" shading=\"nearest\",\n",
")\n",
"plt.colorbar()\n",
"plt.show()\n",
"\n",
"print(\"Velocities:\", fieldset_flat.UV[0, 0, 40, -5])\n",
"for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temp]:\n",
"print(\n",
" \"Velocities:\",\n",
" fieldset_flat.UV[np.array([0]), np.array([0]), np.array([40]), np.array([-5])],\n",
")\n",
"for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temperature]:\n",
" print(f\"{fld.name}: {fld.units}\")"
]
},
Expand Down Expand Up @@ -225,23 +274,35 @@
"kh_zonal = 100 # in m^2/s\n",
"kh_meridional = 100 # in m^2/s\n",
"\n",
"fieldset.add_field(\n",
" parcels.Field(\n",
" \"Kh_zonal\",\n",
" kh_zonal * np.ones((ydim, xdim), dtype=np.float32),\n",
" grid=fieldset.U.grid,\n",
" )\n",
"ds[\"Kh_zonal\"] = xr.DataArray(\n",
" data=kh_zonal * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
")\n",
"fieldset.add_field(\n",
" parcels.Field(\n",
" \"Kh_meridional\",\n",
" kh_meridional * np.ones((ydim, xdim), dtype=np.float32),\n",
" grid=fieldset.U.grid,\n",
" )\n",
"\n",
"kh_zonal_field = parcels.Field(\n",
" \"Kh_zonal\",\n",
" ds[\"Kh_zonal\"],\n",
" grid=fieldset.U.grid,\n",
" interp_method=parcels.interpolators.XLinear,\n",
")\n",
"\n",
"fieldset.add_field(kh_zonal_field)\n",
"\n",
"ds[\"Kh_meridional\"] = xr.DataArray(\n",
" data=kh_meridional * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
")\n",
"\n",
"kh_meridional_field = parcels.Field(\n",
" \"Kh_meridional\",\n",
" ds[\"Kh_meridional\"],\n",
" grid=grid,\n",
" interp_method=parcels.interpolators.XLinear,\n",
")\n",
"\n",
"fieldset.add_field(kh_meridional_field)\n",
"\n",
"for fld in [fieldset.Kh_zonal, fieldset.Kh_meridional]:\n",
" print(f\"{fld.name}: {fld[0, 0, 40, -5]:e} {fld.units}\")"
" val = fld[np.array([0]), np.array([0]), np.array([40]), np.array([-5])]\n",
" print(f\"{fld.name}: {val} {fld.units}\")"
]
},
{
Expand All @@ -261,7 +322,10 @@
"outputs": [],
"source": [
"deg_to_m = 1852 * 60\n",
"print(fieldset.Kh_meridional[0, 0, 40, -5] * deg_to_m**2)"
"print(\n",
" fieldset.Kh_meridional[np.array([0]), np.array([0]), np.array([40]), np.array([-5])]\n",
" * deg_to_m**2\n",
")"
]
},
{
Expand Down Expand Up @@ -298,12 +362,19 @@
"metadata": {},
"outputs": [],
"source": [
"ds[\"Ustokes\"] = xr.DataArray(\n",
" data=np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
")\n",
"\n",
"fieldset.add_field(\n",
" parcels.Field(\n",
" \"Ustokes\", np.ones((ydim, xdim), dtype=np.float32), grid=fieldset.U.grid\n",
" \"Ustokes\",\n",
" ds[\"Ustokes\"],\n",
" grid=fieldset.U.grid,\n",
" interp_method=parcels.interpolators.XLinear,\n",
" )\n",
")\n",
"print(fieldset.Ustokes[0, 0, 40, -5])"
"print(fieldset.Ustokes[np.array([0]), np.array([0]), np.array([40]), np.array([-5])])"
]
},
{
Expand All @@ -320,17 +391,20 @@
"metadata": {},
"outputs": [],
"source": [
"from parcels.tools.converters import GeographicPolar\n",
"\n",
"fieldset.Ustokes.units = GeographicPolar()\n",
"print(fieldset.Ustokes[0, 0, 40, -5])\n",
"print(fieldset.Ustokes[0, 0, 40, -5] * 1852 * 60 * np.cos(40 * np.pi / 180))"
"fieldset.Ustokes.units = parcels.GeographicPolar()\n",
"print(fieldset.Ustokes[np.array([0]), np.array([0]), np.array([40]), np.array([-5])])\n",
"print(\n",
" fieldset.Ustokes[np.array([0]), np.array([0]), np.array([40]), np.array([-5])]\n",
" * 1852\n",
" * 60\n",
" * np.cos(40 * np.pi / 180)\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "test-notebooks",
"language": "python",
"name": "python3"
},
Expand All @@ -344,7 +418,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
"version": "3.11.0"
}
},
"nbformat": 4,
Expand Down
2 changes: 1 addition & 1 deletion docs/user_guide/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4.
:titlesonly:

examples/tutorial_nemo_curvilinear.ipynb
examples/tutorial_unitconverters.ipynb
```

<!-- examples/documentation_indexing.ipynb -->
Expand All @@ -32,7 +33,6 @@ examples/tutorial_nemo_curvilinear.ipynb
<!-- examples/tutorial_timevaryingdepthdimensions.ipynb -->
<!-- examples/tutorial_periodic_boundaries.ipynb -->
<!-- examples/tutorial_interpolation.ipynb -->
<!-- examples/tutorial_unitconverters.ipynb -->

```{toctree}
:caption: Create ParticleSets
Expand Down