Skip to content

Removing UnitConverters from v4 (attempt 2) #2465

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 11 commits into
base: v4-dev
Choose a base branch
from
Open

Removing UnitConverters from v4 (attempt 2) #2465

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
be42132
Removing Geographic and GeographicPolar unit converters
erikvansebille Jan 9, 2026
80fb1b3
Removing Geographis(Polar)Square converters
erikvansebille Jan 9, 2026
c23c5f8
Removing UnitConverter library altogether
erikvansebille Jan 9, 2026
aa15d10
Also removing _convert_to_flat_array
erikvansebille Jan 9, 2026
3a760ac
Cleaning up docstrings
erikvansebille Jan 9, 2026
b0f6ec4
Updating tutorials
erikvansebille Jan 9, 2026
c7555d7
Update explanation_kernelloop.md
erikvansebille Jan 9, 2026
b1da81e
Further clean-up of tutorial
erikvansebille Jan 9, 2026
8cd02cb
Removing applyConversion option from field.eval
erikvansebille Jan 10, 2026
1bab706
Updating field.eval docstrings
erikvansebille Jan 10, 2026
7fccebf
Fixes from merge conflicts
erikvansebille Jan 12, 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
234 changes: 41 additions & 193 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 @@ -63,7 +63,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"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 @@ -91,7 +91,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 @@ -104,38 +108,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 @@ -145,15 +140,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 `apply_conversion=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 @@ -174,19 +184,17 @@
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## UnitConverters for `mesh='flat'`\n"
"## Don't sample U and V separately"
]
},
{
"attachments": {},
"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"
"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 @@ -195,181 +203,21 @@
"metadata": {},
"outputs": [],
"source": [
"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",
"fieldset_flat = parcels.FieldSet.from_sgrid_conventions(ds_flat, mesh=\"flat\")\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"
"for fld in [fieldset.U, fieldset.V]:\n",
" print(f\"{fld.name}: {fld.eval(time, z, lat, lon)}\")"
]
},
{
"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",
"## Unit conversion for other fields such as diffusivity\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",
"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",
"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",
"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(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}\")"
]
},
{
"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"
]
},
{
"cell_type": "code",
"execution_count": null,
"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"
]
},
{
"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",
"\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",
"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",
"\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
3 changes: 2 additions & 1 deletion docs/user_guide/v4-migration.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,12 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
- `Field.eval()` returns an array of floats instead of a single float (related to the vectorization)
- `Field.eval()` does not throw OutOfBounds or other errors
- The `NestedField` class has been removed. See the Nested Grids how-to guide for how to set up Nested Grids in v4.
- `applyConversion` has been renamed to `apply_conversion` and only works for VectorFields. Conversion of units should be handled in Kernels.

## GridSet

- `GridSet` is now a list, so change `fieldset.gridset.grids[0]` to `fieldset.gridset[0]`.

## UnitConverters

- The default `UnitConverter` is now called `Unity()`
- UnitConverters have been removed. Instead, Interpolation functions should handle unit conversion internally, based on the value of `grid._mesh` ("spherical" or "flat").
Loading
Loading