🖥️ Spherical grids and velocity conversion#

In most applications, Parcels works with SphericalMesh() meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use meshes where longitude and latitude are given in meters (note that the dimensions are then still called longitude and latitude for consistency reasons).

In all cases, velocities are given in m s^-1. But for advection (such as particles.dy = v * particles.dt), the velocity needs to be in degrees s^-1.

Parcels seamlessly converts between meters s^-1 and degrees s^-1, under the hood. For transparency, this guide explains how this works.

Note

Parcels automatically detects whether the mesh is spherical based on the units attribute of the coordinates (degrees for spherical meshes, meters otherwise) and applies the appropriate unit conversion for velocities.

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^-1 and the temperature 20°C.

import matplotlib.pyplot as plt
import numpy as np

import parcels

Plotting the U field indeed shows a uniform 1 m s^-1 eastward flow.

plt.pcolormesh(
    fieldset.U.grid.lon,
    fieldset.U.grid.lat,
    fieldset.U.data[0, 0, :, :],
    vmin=0,
    vmax=1,
    shading="gouraud",
)
plt.ylabel("Latitude")
plt.xlabel("Longitude")
plt.colorbar()
plt.show()

../../_images/33e9c1ffe0f8c123034658c8bbd47ac143aa049c1d7b04a1bf0deb76825a2b2c.png

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

Sampling a velocity in Parcels is done by calling the fieldset.UV VectorField; see also the section “Sampling U and V separately” below.

time = np.array([0])
z = np.array([0])
lat = np.array([40])
lon = np.array([-5])
u, v = fieldset.UV[time, z, lat, lon]
temp = fieldset.temperature[time, z, lat, lon]

print(f"(u, v) = ({u}, {v})")
print(f"temperature = {temp}")

(u, v) = ([1.17477258e-05], [8.99928006e-06])
temperature = [20.]

While the temperature field indeed is 20°C, as we defined, these printed velocities are much smaller.

This is because Parcels converts under the hood from m s^-1 to degrees s^-1.

Indeed, if we multiply the value of u with 1852 * 60 * cos(lat) (the number of meters in 1 degree of longitude) and the value of v with 1852 * 60 (the number of meters in 1 degree of latitude), we get the expected 1 m s^-1 for u and v.

u, v = fieldset.UV[time, z, lat, lon]

u = u * 1852 * 60 * np.cos(np.deg2rad(lat))
v = v * 1852 * 60
print(f"(u, v) = ({u}, {v})")

assert np.isclose(u, 1.0)
assert np.isclose(v, 1.0)

(u, v) = ([1.], [1.])

Don’t sample U and V separately#

Sampling U and V separately will not convert to degrees s^-1, so these velocities should not 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.

for fld in [fieldset.U, fieldset.V]:
    print(f"{fld.name}: {fld.eval(time, z, lat, lon)}")

U: [1.]
V: [1.]

Unit conversion for other fields such as diffusivity#

For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a SphericalMesh(). For example, if we define a diffusivity field in m² s^-1, Parcels will not convert this to degrees² s^-1 under the hood.

If you want to work with diffusivity in degrees² s^-1 (for example to move particles using a random walk), you will have to convert this yourself in your kernel.

Note that for the built-in parcels.kernels.DiffusionUniformKh, parcels.kernels.AdvectionDiffusionM1 and parcels.kernels.AdvectionDiffusionEM, the conversion is done automatically.

Parcels on other planets#

Parcels can also be used to simulate particle movement on other planets. In such cases, the radius of the planet may differ from Earth’s radius, and this can affect the conversion between meters and degrees for spherical meshes. You can specify the planetary radius when creating a SphericalMesh() to ensure correct unit conversions for velocities.

radius_jupiter = 69911e3  # radius of Jupiter (in meters)
mesh_jupiter = parcels.SphericalMesh(radius=radius_jupiter)
fieldset_jupiter = parcels.FieldSet.from_sgrid_conventions(ds, mesh=mesh_jupiter)

u_earth, v_earth = fieldset.UV[time, z, lat, lon]

u_jupiter, v_jupiter = fieldset_jupiter.UV[time, z, lat, lon]
print(f"zonal velocity of 1 m/s on Jupiter: {u_jupiter[0]:.2e} degrees/s")
print(f"Zonal velocity of 1 m/s on Earth: {u_earth[0]:.2e} degrees/s")

assert u_jupiter[0] / u_earth[0] < 0.1

zonal velocity of 1 m/s on Jupiter: 1.07e-06 degrees/s
Zonal velocity of 1 m/s on Earth: 1.17e-05 degrees/s

As you see above, the same speed is 10 times smaller in degrees per second on Jupiter compared to Earth.