Add UnitConversion guide v4

docs/conf.py

@@ -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"

docs/user_guide/examples_v3/tutorial_unitconverters.ipynb → docs/user_guide/examples/tutorial_unitconverters.ipynb

@@ -5,7 +5,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Spherical grids and UnitConversion\n"
+    "# Spherical grids and unit converters\n"
    ]
   },
   {
@@ -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 guide 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"
    ]
   },
   {
@@ -34,6 +34,7 @@
    "source": [
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
+    "import xarray as xr\n",
     "\n",
     "import parcels"
    ]
@@ -44,26 +45,29 @@
    "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"
+    "Plotting the `U` field indeed shows a uniform 1 m/s eastward flow.\n"
    ]
   },
   {
@@ -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=\"gouraud\",\n",
+    ")\n",
+    "plt.ylabel(\"Latitude\")\n",
+    "plt.xlabel(\"Longitude\")\n",
     "plt.colorbar()\n",
     "plt.show()"
    ]
@@ -92,8 +113,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(fieldset.UV[0, 0, 40, -5])\n",
-    "print(fieldset.temp[0, 0, 40, -5])"
+    "time = np.array([0])\n",
+    "z = np.array([0])\n",
+    "lat = np.array([40])\n",
+    "lon = np.array([-5])\n",
+    "print(fieldset.UV[time, z, lat, lon])\n",
+    "print(fieldset.temperature[time, z, lat, lon])"
    ]
   },
   {
@@ -103,7 +128,7 @@
    "source": [
     "While the temperature field indeed is 20C, as we defined, these printed velocities are much smaller.\n",
     "\n",
-    "This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a `UnitConverter` object, which is stored in the `.units` attribute of each Field. Below, we print these\n"
+    "This is because Parcels converts under the hood from m/s to degrees/s. This conversion is done with a `parcels.converters` object, which is stored in the `.units` attribute of each Field. Below, we print these\n"
    ]
   },
   {
@@ -112,7 +137,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}\")"
    ]
   },
@@ -121,7 +146,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` Unitconverter and the `temp` field has a `UnitConverter` object.\n",
+    "So the U field has a `GeographicPolar` UnitConverter object, the V field has a `Geographic` UnitConverter and the `temp` field has a `UnitConverter` object.\n",
     "\n",
     "Indeed, if we multiply the value of the V field with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m/s.\n"
    ]
@@ -132,7 +157,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "u, v = fieldset.UV[0, 0, 40, -5]\n",
+    "u, v = fieldset.UV[time, z, lat, lon]\n",
     "print(v * 1852 * 60)"
    ]
   },
@@ -150,7 +175,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(fieldset.UV.eval(0, 0, 40, -5, applyConversion=False))"
+    "print(\n",
+    "    fieldset.UV.eval(\n",
+    "        time,\n",
+    "        z,\n",
+    "        lat,\n",
+    "        lon,\n",
+    "        applyConversion=False,\n",
+    "    )\n",
+    ")"
    ]
   },
   {
@@ -166,7 +199,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"
    ]
   },
   {
@@ -175,20 +208,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=\"gouraud\",\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[time, z, lat, lon],\n",
+    ")\n",
+    "for fld in [fieldset_flat.U, fieldset_flat.V, fieldset_flat.temperature]:\n",
     "    print(f\"{fld.name}: {fld.units}\")"
    ]
   },
@@ -225,23 +276,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[time, z, lat, lon]\n",
+    "    print(f\"{fld.name}: {val} {fld.units}\")"
    ]
   },
   {
@@ -261,7 +324,7 @@
    "outputs": [],
    "source": [
     "deg_to_m = 1852 * 60\n",
-    "print(fieldset.Kh_meridional[0, 0, 40, -5] * deg_to_m**2)"
+    "print(fieldset.Kh_meridional[time, z, lat, lon] * deg_to_m**2)"
    ]
   },
   {
@@ -281,11 +344,11 @@
     "\n",
     "| Field name       | Converter object        | Conversion for `mesh='spherical'`                         | Conversion for `mesh='flat'` |\n",
     "| ---------------- | ----------------------- | --------------------------------------------------------- | ---------------------------- |\n",
-    "| 'U'              | `GeographicPolar`       | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$     | 1                            |\n",
-    "| 'V'              | `Geographic`            | $1852 \\cdot 60$                                           | 1                            |\n",
-    "| 'Kh_zonal'       | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1                            |\n",
-    "| 'Kh_meridional'  | `GeographicSquare`      | $(1852 \\cdot 60)^2$                                       | 1                            |\n",
-    "| All other fields | `UnitConverter`         | 1                                                         | 1                            |\n",
+    "| `\"U\"`              | `GeographicPolar`       | $1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180})$     | 1                            |\n",
+    "| `\"V\"`              | `Geographic`            | $1852 \\cdot 60$                                           | 1                            |\n",
+    "| `\"Kh_zonal\"`       | `GeographicPolarSquare` | $(1852 \\cdot 60 \\cdot \\cos(lat \\cdot \\frac{\\pi}{180}))^2$ | 1                            |\n",
+    "| `\"Kh_meridional\"`  | `GeographicSquare`      | $(1852 \\cdot 60)^2$                                       | 1                            |\n",
+    "| All other fields   | `UnitConverter`         | 1                                                         | 1                            |\n",
     "\n",
     "Only four Field names are recognised and assigned an automatic UnitConverter object. This means that things might go very wrong when e.g. a velocity field is not called `U` or `V`.\n",
     "\n",
@@ -298,12 +361,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[time, z, lat, lon])"
    ]
   },
   {
@@ -320,17 +390,15 @@
    "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[time, z, lat, lon])\n",
+    "print(fieldset.Ustokes[time, z, lat, lon] * 1852 * 60 * np.cos(40 * np.pi / 180))"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "test-notebooks",
    "language": "python",
    "name": "python3"
   },
@@ -344,7 +412,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.0"
   }
  },
  "nbformat": 4,

docs/user_guide/index.md

@@ -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 -->
@@ -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