Skip to content

Removing UnitConverters from v4 #2459

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 10 commits into
base: unitconversion_for_cgridvelocity
Choose a base branch
from
+186 −448
Open

Removing UnitConverters from v4 #2459

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
c7b9799
Removing Geographic and GeographicPolar unit converters
erikvansebille Jan 9, 2026
ae080ac
Removing Geographis(Polar)Square converters
erikvansebille Jan 9, 2026
124c032
Removing UnitConverter library altogether
erikvansebille Jan 9, 2026
6254524
Also removing _convert_to_flat_array
erikvansebille Jan 9, 2026
cd190f7
Cleaning up docstrings
erikvansebille Jan 9, 2026
892ed42
Updating tutorials
erikvansebille Jan 9, 2026
c9f63cc
Update explanation_kernelloop.md
erikvansebille Jan 9, 2026
aa16654
Further clean-up of tutorial
erikvansebille Jan 9, 2026
b678ccd
Removing applyConversion option from field.eval
erikvansebille Jan 10, 2026
aeaf27e
Updating field.eval docstrings
erikvansebille Jan 10, 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
16 changes: 7 additions & 9 deletions docs/user_guide/examples/explanation_kernelloop.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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]`
Expand Down
244 changes: 41 additions & 203 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 @@ -5,7 +5,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# 🖥️ Spherical grids and unit converters"
"# 🖥️ Spherical grids and velocity conversion"
]
},
{
Expand All @@ -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"
]
},
{
Expand Down Expand Up @@ -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"
]
},
{
Expand Down Expand Up @@ -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"
]
},
{
Expand All @@ -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"
]
},
{
Expand All @@ -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"
]
},
{
Expand All @@ -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"
]
},
{
Expand All @@ -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."
]
}
],
Expand Down
3 changes: 1 addition & 2 deletions docs/user_guide/examples_v3/example_moving_eddies.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
3 changes: 1 addition & 2 deletions docs/user_guide/examples_v3/example_peninsula.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
Loading
Loading