Merge pull request #2459 from Parcels-code/removing_unitconverters

Removing UnitConverters from v4

Author: erikvansebille

docs/user_guide/examples/explanation_kernelloop.md

@@ -43,6 +43,7 @@ Besides having commutable Kernels, the main advantage of this implementation is
 Below is a simple example of some particles at the surface of the ocean. We create an idealised zonal wind flow that will "push" a particle that is already affected by the surface currents. The Kernel loop ensures that these two forces act at the same time and location.
 
 ```{code-cell}
+:tags: [hide-output]
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray as xr
@@ -69,21 +70,18 @@ ds_fields["VWind"] = xr.DataArray(
 
 fieldset = parcels.FieldSet.from_copernicusmarine(ds_fields)
 
-# Set unit converters for custom wind fields
-fieldset.UWind.units = parcels.GeographicPolar()
-fieldset.VWind.units = parcels.Geographic()
+# Create a vecorfield for the wind
+windvector = parcels.VectorField("Wind", fieldset.UWind, fieldset.VWind)
+fieldset.add_field(windvector)
 ```
 
 Now we define a wind kernel that uses a forward Euler method to apply the wind forcing. Note that we update the `particles.dlon` and `particles.dlat` variables, rather than `particles.lon` and `particles.lat` directly.
 
 ```{code-cell}
 def wind_kernel(particles, fieldset):
-    particles.dlon += (
-        fieldset.UWind[particles] * particles.dt
-    )
-    particles.dlat += (
-        fieldset.VWind[particles] * particles.dt
-    )
+    uwind, vwind = fieldset.Wind[particles]
+    particles.dlon += uwind * particles.dt
+    particles.dlat += vwind * particles.dt
 ```
 
 First run a simulation where we apply kernels as `[AdvectionRK4, wind_kernel]`

docs/user_guide/examples/tutorial_unitconverters.ipynb

@@ -5,7 +5,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# 🖥️ Spherical grids and unit converters"
+    "# 🖥️ Spherical grids and velocity conversion"
    ]
   },
   {
@@ -23,7 +23,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "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"
+    "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<sup>-1</sup> and the temperature 20°C.\n"
    ]
   },
   {
@@ -67,7 +67,7 @@
     "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"
+    "Plotting the `U` field indeed shows a uniform 1 m s<sup>-1</sup> eastward flow.\n"
    ]
   },
   {
@@ -104,7 +104,11 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "However, printing the velocites directly shows something perhaps surprising. Here, we use the square-bracket field-interpolation notation to print the field value at (5W, 40N, 0m depth) at time 0. _Note that sampling a velocity in Parcels is done by calling the `fieldset.UV` VectorField; see the [Field Sampling tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_sampling.html#Sampling-velocity-fields) for more information._\n"
+    "However, printing the velocites directly shows something perhaps surprising. Here, we use the square-bracket field-interpolation notation to print the field value at (5W, 40N, 0m depth) at time 0. \n",
+    "\n",
+    "```{note}\n",
+    "Sampling a velocity in Parcels is done by calling the `fieldset.UV` VectorField; see also the section \"Sampling U and V separately\" below.\n",
+    "```\n"
    ]
   },
   {
@@ -117,38 +121,29 @@
     "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])"
+    "u, v = fieldset.UV[time, z, lat, lon]\n",
+    "temp = fieldset.temperature[time, z, lat, lon]\n",
+    "\n",
+    "print(f\"(u, v) = ({u}, {v})\")\n",
+    "print(f\"temperature = {temp}\")"
    ]
   },
   {
    "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "While the temperature field indeed is 20C, as we defined, these printed velocities are much smaller.\n",
+    "While the temperature field indeed is 20°C, 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 `parcels.converters` object, which is stored in the `.units` attribute of each Field. Below, we print these\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "for fld in [fieldset.U, fieldset.V, fieldset.temperature]:\n",
-    "    print(f\"{fld.name}: {fld.units}\")"
+    "This is because Parcels converts under the hood from m s<sup>-1</sup> to degrees s<sup>-1</sup>."
    ]
   },
   {
    "attachments": {},
    "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 `Unity` 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"
+    "Indeed, if we multiply the value of the U field with 1852 \\* 60 \\* cos(lat) (the number of meters in 1 degree of longitude), we get the expected 1 m s<sup>-1</sup>.\n"
    ]
   },
   {
@@ -158,15 +153,30 @@
    "outputs": [],
    "source": [
     "u, v = fieldset.UV[time, z, lat, lon]\n",
-    "print(v * 1852 * 60)"
+    "\n",
+    "u = u * 1852 * 60 * np.cos(np.deg2rad(lat))\n",
+    "v = v * 1852 * 60 * np.cos(np.deg2rad(lat))\n",
+    "print(f\"(u, v) = ({u}, {v})\")\n",
+    "\n",
+    "assert np.isclose(u, 1.0)\n",
+    "assert np.isclose(v, 1.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "It may be surprising that the conversion factor depends on latitude also for the V velocity component. This is because Parcels assumes a flux-based interpolation, so that the velocity components are defined on the faces of the grid cells, which vary in size with latitude.\n",
+    "```"
    ]
   },
   {
    "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Note that you can also interpolate the Field without a unit conversion, by using the `eval()` method and setting `applyConversion=False`, as below\n"
+    "You can also interpolate the Field to these values directly by using the `eval()` method and setting `applyConversion=False`, as below\n"
    ]
   },
   {
@@ -187,134 +197,17 @@
    ]
   },
   {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## UnitConverters for `mesh='flat'`\n"
-   ]
-  },
-  {
-   "attachments": {},
-   "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"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "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.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(\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}\")"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Indeed, in this case all Fields have the same default `Unity` object.\n"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## UnitConverters for Diffusion fields\n"
-   ]
-  },
-  {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The units for Brownian diffusion are in $m^2/s$. If (and only if!) the diffusion fields are called \"Kh_zonal\" and \"Kh_meridional\", Parcels will automatically assign the correct Unitconverter objects to these fields.\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "kh_zonal = 100  # in m^2/s\n",
-    "kh_meridional = 100  # in m^2/s\n",
-    "\n",
-    "ds[\"Kh_zonal\"] = xr.DataArray(\n",
-    "    data=kh_zonal * np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\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=fieldset.U.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",
-    "    val = fld[time, z, lat, lon]\n",
-    "    print(f\"{fld.name}: {val} {fld.units}\")"
+    "## Don't sample U and V separately"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Here, the unitconverters are `GeographicPolarSquare` and `GeographicSquare`, respectively.\n",
-    "\n",
-    "Indeed, multiplying with $(1852\\cdot60)^2$ returns the original value\n"
+    "Sampling `U` and `V` separately will _not_ convert to degrees s<sup>-1</sup>, so these velocities cannot be used directly for advection on spherical coordinates. This is one of the main reasons to always use the `UV` VectorField for velocity sampling in Parcels.\n"
    ]
   },
   {
@@ -323,76 +216,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "deg_to_m = 1852 * 60\n",
-    "print(fieldset.Kh_meridional[time, z, lat, lon] * deg_to_m**2)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Adding a UnitConverter object to a Field\n"
+    "for fld in [fieldset.U, fieldset.V]:\n",
+    "    print(f\"{fld.name}: {fld.eval(time, z, lat, lon)}\")"
    ]
   },
   {
-   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "So, to summarise, here is a table with all the conversions\n",
-    "\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   | `Unity`                 | 1                                                         | 1                            |\n",
+    "## Unit conversion for other fields such as diffusivity\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",
+    "For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a `spherical` mesh. For example, if we define a diffusivity field with value 10 m<sup>2</sup> s<sup>-1</sup>, Parcels will not convert this to degrees<sup>2</sup> s<sup>-1</sup> under the hood. \n",
     "\n",
-    "Fortunately, you can always add a UnitConverter later, as explained below:\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "ds[\"Ustokes\"] = xr.DataArray(\n",
-    "    data=np.ones((nlat, nlon), dtype=np.float32), dims=[\"YG\", \"XG\"]\n",
-    ")\n",
+    "If you want to work with diffusivity in degrees<sup>2</sup> s<sup>-1</sup> (for example to move particles using a random walk), you will have to convert this yourself in your kernel. \n",
     "\n",
-    "fieldset.add_field(\n",
-    "    parcels.Field(\n",
-    "        \"Ustokes\",\n",
-    "        ds[\"Ustokes\"],\n",
-    "        grid=fieldset.U.grid,\n",
-    "        interp_method=parcels.interpolators.XLinear,\n",
-    "    )\n",
-    ")\n",
-    "print(fieldset.Ustokes[time, z, lat, lon])"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This value for `Ustokes` of course is not as expected, since the mesh is spherical and hence this would mean 1 degree/s velocity. Assigning the correct `GeographicPolar` Unitconverter gives\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "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))"
+    "Note that for the built-in `DiffusionUniformKh`, `AdvectionDiffusionM1` and `AdvectionDiffusionEM`, the conversion is done automatically."
    ]
   }
  ],

docs/user_guide/examples_v3/example_moving_eddies.py

@@ -25,8 +25,7 @@ def moving_eddies_fieldset(xdim=200, ydim=350, mesh="flat"):
     xdim :
         Vertical dimension of the generated fieldset (Default value = 200)
     mesh : str
-        String indicating the type of mesh coordinates and
-        units used during velocity interpolation:
+        String indicating the type of mesh coordinates used during velocity interpolation:
 
         1. spherical: Lat and lon in degree, with a
            correction for zonal velocity U near the poles.

docs/user_guide/examples_v3/example_peninsula.py

@@ -24,8 +24,7 @@ def peninsula_fieldset(xdim, ydim, mesh="flat", grid_type="A"):
     xdim :
         Vertical dimension of the generated fieldset
     mesh : str
-        String indicating the type of mesh coordinates and
-        units used during velocity interpolation:
+        String indicating the type of mesh coordinates used during velocity interpolation:
 
         1. spherical: Lat and lon in degree, with a
            correction for zonal velocity U near the poles.