🖥️ Using the built-in Interpolators#
Parcels comes with a number of different interpolation methods for fields on structured (X) and unstructured (Ux) grids. Here, we will look at a few common parcels.interpolators for tracer fields, and how to configure them in an idealised example. For more guidance on the sampling of such fields, check out the sampling tutorial.
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
import parcels
Interpolators on structured grids#
We will first look at interpolation schemes which work on tracer fields defined on parcels.XGrid objects.
We create a small 2D structured grid, where P is a tracer that we want to interpolate. In each grid cell, P has a random value between 0.1 and 1.1. We then set P[1,1] to 0, which for Parcels specifies that this is a land cell
from parcels._datasets.structured.generated import simple_UV_dataset
ds = simple_UV_dataset(n_times=1, nz=1, ny=5, nx=4, mesh_str="flat")
ds["lat"][:] = np.linspace(0.0, 1.0, len(ds.YG))
ds["lon"][:] = np.linspace(0.0, 1.0, len(ds.XG))
dx, dy = 1.0 / len(ds.XG), 1.0 / len(ds.YG)
ds["P"] = ds["U"] + np.random.rand(1, 1, 5, 4) + 0.1
ds["P"][:, :, 1, 1] = 0
ds
<xarray.Dataset> Size: 720B
Dimensions: (time: 1, depth: 1, YG: 5, XG: 4, YC: 5, XC: 4)
Coordinates:
* time (time) datetime64[ns] 8B 2000-01-01
* depth (depth) float64 8B 0.0
* YG (YG) int64 40B 0 1 2 3 4
lat (YG) float64 40B 0.0 0.25 0.5 0.75 1.0
* XG (XG) int64 32B 0 1 2 3
lon (XG) float64 32B 0.0 0.3333 0.6667 1.0
* YC (YC) float64 40B 0.5 1.5 2.5 3.5 4.5
* XC (XC) float64 32B 0.5 1.5 2.5 3.5
Data variables:
U (time, depth, YG, XG) float64 160B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
V (time, depth, YG, XG) float64 160B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
grid int64 8B 0
P (time, depth, YG, XG) float64 160B 0.437 0.4078 ... 0.6681 0.727
Attributes:
Conventions: SGRIDFrom this dataset we create a parcels.FieldSet using the parcels.FieldSet.from_sgrid_conventions() constructor, which automatically sets up the grid and fields according to Parcels’ s-grid conventions.
The default interpolator for fields on structured A-grid grids is (tri)linear interpolation, implemented in parcels.interpolators.XLinear.
fieldset = parcels.FieldSet.from_sgrid_conventions(ds)
assert isinstance(fieldset.P.interp_method, parcels.interpolators.XLinear)
We create a Particle class that can sample this field P.
SampleParticle = parcels.Particle.add_variable(
parcels.Variable("p", dtype=np.float32, initial=np.nan)
)
def SampleP(particles, fieldset):
particles.p = fieldset.P[particles]
Now, we perform four different interpolations on P, which we can control by setting fieldset.P.interp_method. Note that this can always be done after the FieldSet creation. We store the results of each interpolation method in an entry in the dictionary pset.
pset = {}
from parcels.interpolators import (
CGrid_Tracer,
XLinear,
XLinearInvdistLandTracer,
XNearest,
)
interp_methods = [XLinear(), XLinearInvdistLandTracer(), XNearest(), CGrid_Tracer()]
for p_interp in interp_methods:
fieldset.P.interp_method = (
p_interp # setting the interpolation method for fieldset.P
)
fieldset.describe()
xv, yv = np.meshgrid(np.linspace(0, 1, 8), np.linspace(0, 1, 8))
pset[p_interp.__class__.__name__] = parcels.ParticleSet(
fieldset, pclass=SampleParticle, x=xv.flatten(), y=yv.flatten()
)
pset[p_interp.__class__.__name__].execute(
SampleP,
runtime=np.timedelta64(1, "D"),
dt=np.timedelta64(1, "D"),
verbose_progress=False,
)
And then we can show each of the four interpolation methods, by plotting the interpolated values on the Particles locations (circles) on top of the Field values (background colours)
fig, ax = plt.subplots(1, 4, figsize=(18, 6))
for i, p in enumerate(pset.keys()):
data = np.copy(fieldset.P.data[0, 0, :, :])
data[1, 1] = np.nan
x = np.linspace(-dx / 2, 1 + dx / 2, len(ds.XG) + 1)
y = np.linspace(-dy / 2, 1 + dy / 2, len(ds.YG) + 1)
if p == "CGrid_Tracer":
for lat in fieldset.P.grid.lat:
ax[i].axhline(lat, color="k", linestyle="--")
for lon in fieldset.P.grid.lon:
ax[i].axvline(lon, color="k", linestyle="--")
pc = ax[i].pcolormesh(x, y, data, vmin=0.1, vmax=1.1)
ax[i].scatter(
pset[p].x, pset[p].y, c=pset[p].p, edgecolors="k", s=50, vmin=0.1, vmax=1.1
)
xp, yp = np.meshgrid(fieldset.P.grid.lon, fieldset.P.grid.lat)
ax[i].plot(xp, yp, "kx")
ax[i].set_title(f"Using interp_method='{p}'")
plt.colorbar(pc, ax=ax, orientation="horizontal", shrink=0.5)
plt.show()
The white box is here the ‘land’ point where the tracer is set to zero and the crosses are the locations of the grid points. As you see, the interpolated value is always equal to the field value if the particle is exactly on the grid point (circles on crosses).
For interp_method='XNearest', the particle values are the same for all particles in a grid cell. They are also the same for interp_method='CGrid_Tracer', but the grid cells have then shifted. That is because in a C-grid, the tracer is defined on the corners of the velocity grid (black dashed lines in right-most panel).
For interp_method='XLinearInvdistLandTracer', we see that values are the same as interp_method='XLinear' for grid cells that don’t border the land point. For grid cells that do border the land cell, the XLinearInvdistLandTracer interpolation method gives higher values, as also shown in the difference plot below.
plt.scatter(
pset["XLinear"].x,
pset["XLinear"].y,
c=pset["XLinearInvdistLandTracer"].p - pset["XLinear"].p,
edgecolors="k",
s=50,
cmap=cm.bwr,
vmin=-0.05,
vmax=0.05,
)
xp, yp = np.meshgrid(fieldset.P.grid.lon, fieldset.P.grid.lat)
plt.plot(xp, yp, "kx")
plt.colorbar()
plt.title(
"Difference between 'interp_method=XLinear' "
"and 'interp_method=XLinearInvdistLandTracer'"
)
plt.show()
So in summary, Parcels has four different interpolation schemes for tracers on structured grids:
interp_method=parcels.interpolators.XLinear: compute linear interpolationinterp_method=parcels.interpolators.XLinearInvdistLandTracer: compute linear interpolation except near land (where field value is zero). In that case, inverse distance weighting interpolation is computed, weighting by squares of the distance.interp_method=parcels.interpolators.XNearest: return nearest field valueinterp_method=parcels.interpolators.CGridTracer: return nearest field value supposing C cells
In the special case where a parcels.Field is constant in time and space, such as implementing constant diffusion on a spherical mesh, we can use:
interp_method=parcels.interpolators.XConstantField: return single value of a Constant Field
[..] API reference
All available interpolation methods are listed here: parcels.interpolators
Interpolation at boundaries#
In some cases, we need to implement specific boundary conditions, for example to prevent particles from
getting “stuck” near land. This guide describes
how to implement this in Parcels using parcels.interpolators.XFreeslip and parcels.interpolators.XPartialslip.
Interpolators on unstructured grids#
Parcels v4 supports the use of general circulation model output that is defined on unstructured grids. We include basic interpolators to help you get started, including
UxConstantFaceConstantZC- this interpolator implements piecewise constant interpolation in both the lateral and vertical directions. It is appropriate for data that is registered to the face centres of the unstructured grid and centred on vertical layers.UxLinearNodeConstantZC- this interpolator implements barycentric interpolation in the lateral direction and piecewise constant in the vertical direction. It is appropriate for data that is registered to the corner vertices of the unstructured grid faces and centred on vertical layers.UxLinearNodeLinearZF- this interpolator implements barycentric interpolation in the lateral direction and piecewise linear interpolation in the vertical direction. It is appropriate for data that is registered to the corner vertices of the unstructured grid faces and on vertical layer interfaces.UxConstantFaceLinearZF- this interpolator implements piecewise constant interpolation in the lateral direction and piecewise linear interpolation in the vertical direction. It is appropriate for data that is registered to the face centres of the unstructured grid and centred on vertical layers
To get started, we use a very simple generated UxArray.UxDataset that is included with Parcels.
from parcels._datasets.unstructured.generated import simple_small_delaunay
ds = simple_small_delaunay(nx=4, ny=4)
ds
<xarray.UxDataset> Size: 916B
Dimensions: (time: 1, zc: 1, n_face: 16, zf: 2, n_node: 16)
Coordinates:
* time (time) datetime64[ns] 8B 2000-01-01
* zc (zc) float32 4B 500.0
* zf (zf) float32 8B 0.0 1e+03
Dimensions without coordinates: n_face, n_node
Data variables:
U (time, zc, n_face) float64 128B -0.1699 -0.03851 ... 0.4288 0.09719
V (time, zc, n_face) float64 128B -0.6425 -0.9843 ... -0.2431 -0.3724
W (time, zf, n_node) float64 256B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
p (time, zc, n_face) float64 128B 0.01428 0.02188 ... 0.03833 0.05872
T_face (time, zc, n_face) float64 128B -0.4924 -0.7544 ... 0.604 0.9254
T_node (time, zc, n_node) float64 128B 0.0 0.0 ... -3.812e-08 -1.421e-07Next, we create theFieldset object that will be used later in advancing particles. When using the FieldSet.from_ugrid_conventions() method, Parcels takes care of automatically assigning the appropriate interpolator based on the coordinates of each data variable.
import numpy as np
import parcels
fieldset = parcels.FieldSet.from_ugrid_conventions(ds)
We can now create particles to sample data using either the node registered or face registered data.
SampleParticle = parcels.Particle.add_variable(
parcels.Variable("tracer", dtype=np.float32, initial=np.nan)
)
def SampleTracer_Node(particles, fieldset):
particles.tracer = fieldset.T_node[particles]
def SampleTracer_Face(particles, fieldset):
particles.tracer = fieldset.T_face[particles]
Now, we can create a ParticleSet for each interpolation experiment. In one of the ParticleSet objects we use the SampleTracer_Node kernel to showcase interpolating with the UxLinearNodeLinearZF interpolator for node-registered data. In the other, we use the SampleTracer_Face to showcase interpolating with the UxConstantFaceConstantZC interpolator for face-registered data.
pset = {}
xv, yv = np.meshgrid(np.linspace(0.2, 0.8, 8), np.linspace(0.2, 0.9, 8))
pset["node"] = parcels.ParticleSet(
fieldset, pclass=SampleParticle, x=xv.flatten(), y=yv.flatten()
)
pset["node"].execute(
SampleTracer_Node,
runtime=np.timedelta64(1, "D"),
dt=np.timedelta64(1, "D"),
verbose_progress=False,
)
pset["face"] = parcels.ParticleSet(
fieldset, pclass=SampleParticle, x=xv.flatten(), y=yv.flatten()
)
pset["face"].execute(
SampleTracer_Face,
runtime=np.timedelta64(1, "D"),
dt=np.timedelta64(1, "D"),
verbose_progress=False,
)
Now, we can plot the node-registered (T_node) and face-registered (T_face) fields with the values interpolated onto the particles for each field.
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
fig, ax = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
x = ds.uxgrid.node_lon.values
y = ds.uxgrid.node_lat.values
tri = ds.uxgrid.face_node_connectivity.values
triang = mtri.Triangulation(x, y, triangles=tri)
# Shade the triangle faces using smooth shading between the node registered data
T_node = np.squeeze(ds["T_node"][0, 0, :].values)
tpc = ax[0].tripcolor(
triang,
T_node,
shading="gouraud", # smooth shading
cmap="viridis",
vmin=-1.0,
vmax=1.0,
)
# Plot the element edges
ax[0].triplot(triang, color="black", linewidth=0.5)
# Plot the node registered data
sc1 = ax[0].scatter(
x, y, c=T_node, cmap="viridis", s=150, edgecolors="black", vmin=-1.0, vmax=1.0
)
ax[0].scatter(
pset["node"].x,
pset["node"].y,
c=pset["node"].tracer,
cmap="viridis",
edgecolors="k",
s=50,
vmin=-1.0,
vmax=1.0,
)
cbar = fig.colorbar(sc1, ax=ax[0])
cbar.set_label("T")
ax[0].set_aspect("equal", "box")
ax[0].set_xlabel("x")
ax[0].set_ylabel("y")
ax[0].set_title("Node Registered Data")
# # Plot the face registered data
T_face = np.squeeze(ds["T_face"][0, 0, :].values)
assert triang.triangles.shape[0] == T_face.shape[0]
tpc2 = ax[1].tripcolor(
triang,
facecolors=T_face,
shading="flat",
cmap="viridis",
edgecolors="k",
linewidth=0.5,
vmin=-1.0,
vmax=1.0,
)
# Plot the element edges
ax[1].triplot(triang, color="black", linewidth=0.5)
# Plot the face registered data
xf = ds.uxgrid.face_lon.values
yf = ds.uxgrid.face_lat.values
ax[1].scatter(
pset["face"].x,
pset["face"].y,
c=pset["face"].tracer,
cmap="viridis",
edgecolors="k",
s=50,
vmin=-1.0,
vmax=1.0,
)
cbar = fig.colorbar(tpc, ax=ax[1])
cbar.set_label("T")
ax[1].set_aspect("equal", "box")
ax[1].set_xlabel("x")
ax[1].set_ylabel("y")
ax[1].set_title("Face Registered Data")
plt.show()
In both plots the black lines show the edges that define the boundaries between the faces in the unstructured grid. For the node registered field, the background colouring is done using smooth shading based on the value of T_node at the corner nodes. For the face registered fields, each face is coloured according to the value of T_face at the face centre. The colour of the particles is the value of the function that the particles take on via the corresponding interpolation method.
Write your own interpolator#
The interpolators included in Parcels are designed for common interpolation schemes in Parcels simulations. If we want to add a custom interpolation method, we can write a custom function (e.g. spline) using the interpolator API.