Skip to content
Merged
Show file tree
Hide file tree
Changes from 10 commits
Commits
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
1 change: 0 additions & 1 deletion docs/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,6 @@
"**.ipynb_checkpoints",
"user_guide/examples_v3",
".jupyter_cache",
"user_guide/examples/explanation_kernelloop.md", # TODO v4: https://github.com/Parcels-code/Parcels/issues/2695
]

# The reST default role (used for this markup: `text`) to use for all
Expand Down
22 changes: 11 additions & 11 deletions docs/getting_started/explanation_concepts.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ Here, we will explain the most important classes and functions. This overview ca
A Parcels simulation is generally built up from four different components:

1. [**FieldSet**](#1-fieldset). The input dataset of gridded fields (e.g. ocean current velocity, temperature) in which virtual particles are defined.
2. [**ParticleSet**](#2-particleset). The dataset of virtual particles. These always contain time, z, lat, and lon, for which initial values must be defined. The ParticleSet may also contain other, custom variables.
2. [**ParticleSet**](#2-particleset). The dataset of virtual particles. These always contain time, z, y, and x, for which initial values must be defined. The ParticleSet may also contain other, custom variables.
3. [**Kernels**](#3-kernels). Kernels perform some specific operation on the particles every time step (e.g. advect the particles with the three-dimensional flow; or interpolate the temperature field to the particle location).
4. [**Execute**](#4-execute). Execute the simulation. The core method which integrates the operations defined in Kernels for a given runtime and timestep, and writes output to a ParticleFile.

Expand Down Expand Up @@ -82,16 +82,16 @@ Once the environment has a `parcels.FieldSet` object, you can start defining you

1. The `parcels.FieldSet` object in which the particles will be released.
2. The type of `parcels.Particle`: A default `Particle` or a custom `Particle`-type with additional `Variable`s (see the [custom kernel example](custom-kernel)).
3. Initial conditions for each `Variable` defined in the `Particle`, most notably the release coordinates of `time`, `z`, `lat` and `lon`.
3. Initial conditions for each `Variable` defined in the `Particle`, most notably the release coordinates of `time`, `z`, `y` and `x`.

```python
time = np.array([0])
z = np.array([0])
lat = np.array([0])
lon = np.array([0])
y = np.array([0])
x = np.array([0])

# Create a ParticleSet
pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, time=time, z=z, lat=lat, lon=lon)
pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, time=time, z=z, y=y, x=x)
```

```{admonition} 🖥️ Learn more about how to create ParticleSets
Expand All @@ -103,7 +103,7 @@ pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, time=time

A **`parcels.Kernel`** object is a little snippet of code, which is applied to the particles in the `ParticleSet`, for every time step during a simulation. Kernels define the computation or numerical integration done by Parcels, and can represent many processes such as advection, ageing, growth, or simply the sampling of a field.

Advection of a particle by the flow, the change in position $\mathbf{x}(t) = (lon(t), lat(t))$ at time $t$, can be described by the equation:
Advection of a particle by the flow, the change in position $\mathbf{x}(t) = (x(t), y(t))$ at time $t$, can be described by the equation:

$$
\begin{aligned}
Expand All @@ -113,20 +113,20 @@ $$

where $\mathbf{v}(\mathbf{x},t) = (u(\mathbf{x},t), v(\mathbf{x},t))$ describes the ocean velocity field at position $\mathbf{x}$ at time $t$.

In Parcels, we can write a kernel function which integrates this equation at each timestep `particles.dt`. To do so, we need the ocean velocity field `fieldset.UV` at the `particles` location, and compute the change in position, `particles.dlon` and `particles.dlat`.
In Parcels, we can write a kernel function which integrates this equation at each timestep `particles.dt`. To do so, we need the ocean velocity field `fieldset.UV` at the `particles` location, and compute the change in position, `particles.dx` and `particles.dy`.

```python
def AdvectionEE(particles, fieldset):
"""Advection of particles using Explicit Euler (aka Euler Forward) integration."""
(u1, v1) = fieldset.UV[particles]
particles.dlon += u1 * particles.dt
particles.dlat += v1 * particles.dt
particles.dx += u1 * particles.dt
particles.dy += v1 * particles.dt
```

Basic kernels are included in Parcels to compute advection and diffusion. The standard advection kernel is `parcels.kernels.AdvectionRK2`, a [second-order Runge-Kutta integrator](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method) of the advection function.

```{warning}
It is advised _not_ to update the particle coordinates (`particles.time`, `particles.z`, `particles.lat`, or `particles.lon`) directly within a Kernel, as that can negatively interfere with the way that particle movements by different kernels are vectorially added. Use a change in the coordinates: `particles.dlon`, `particles.dlat` and/or `particles.dz`. Read the [kernel loop tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html) to understand why.
It is advised _not_ to update the particle coordinates (`particles.time`, `particles.z`, `particles.y`, or `particles.x`) directly within a Kernel, as that can negatively interfere with the way that particle movements by different kernels are vectorially added. Use a change in the coordinates: `particles.dy`, `particles.dx` and/or `particles.dz`. Read the [kernel loop tutorial](../user_guide/examples/explanation_kernelloop.md) to understand why.
```

(custom-kernel)=
Expand Down Expand Up @@ -186,7 +186,7 @@ pset.execute(kernels=kernels, dt=dt, runtime=runtime)

### Output

To analyse the particle data generated in the simulation, we need to define a `parcels.ParticleFile` and add it as an argument to `parcels.ParticleSet.execute()`. The output will be written in a [parquet format](https://parquet.apache.org/), which can be opened as a `polars.DataFrame`. The dataset will contain the particle data with at least `time`, `z`, `lat` and `lon`, for each particle at timesteps defined by the `outputdt` argument.
To analyse the particle data generated in the simulation, we need to define a `parcels.ParticleFile` and add it as an argument to `parcels.ParticleSet.execute()`. The output will be written in a [parquet format](https://parquet.apache.org/), which can be opened as a `polars.DataFrame`. The dataset will contain the particle data with at least `time`, `z`, `y` and `x`, for each particle at timesteps defined by the `outputdt` argument.

There are many ways to analyze particle output, and although we provide [a short tutorial to get started](./tutorial_output.ipynb), we recommend writing your own analysis code and checking out [related Lagrangian analysis projects in our community page](../community/index.md#analysis-code).

Expand Down
32 changes: 16 additions & 16 deletions docs/getting_started/tutorial_output.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -77,11 +77,11 @@
"npart = 10 # number of particles to be released\n",
"lon = 32 * np.ones(npart)\n",
"lat = np.linspace(-32.5, -30.5, npart, dtype=np.float32)\n",
"time = ds_fields.time.values[0] + np.arange(0, npart) * np.timedelta64(2, \"h\")\n",
"z = np.repeat(ds_fields.depth.values[0], npart)\n",
"time = ds_fields.time.values[0] + np.arange(0, npart) * np.timedelta64(2, \"h\")\n",
"\n",
"pset = parcels.ParticleSet(\n",
" fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z\n",
" fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, z=z, time=time\n",
")\n",
"\n",
"output_file = parcels.ParticleFile(\"output.parquet\", outputdt=np.timedelta64(2, \"h\"))"
Expand Down Expand Up @@ -145,9 +145,9 @@
"```{note}\n",
"As of Parcels v4, the default output format is [`parquet`](https://parquet.apache.org/) (instead of `zarr`). The `parquet` output format is a tabular format, in which every row corresponds to an observation of a particle trajectory. The `zarr` output format is a multidimensional array format, in which the data is stored in a 2D array with dimensions `traj` and `obs`. The `parquet` format is more compact and faster to read.\n",
"\n",
"However, the `parquet` format does not support the [CF-convention for trajectories data](http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_multidimensional_array_representation_of_trajectories) implemented with the [NCEI trajectory template](https://www.ncei.noaa.gov/data/oceans/ncei/formats/netcdf/v2.0/trajectoryIncomplete.cdl). We are working on efficient tooling to convert the parcels `parquet` output into a CF-compliant format.\n",
"However, the `parquet` format does not support the [CF-convention for trajectories data](http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_multidimensional_array_representation_of_trajectories) implemented with the [NCEI trajectory template](https://www.ncei.noaa.gov/data/oceans/ncei/formats/netcdf/v2.0/trajectoryIncomplete.cdl). We have implemented a `parcels.read_particlefile()` function to facilitate reading `parquet` output files, see more information below.\n",
"\n",
"TODO: Add link to tracking issue on github for this tooling.\n",
"TODO: Add information on conversion functions once https://github.com/Parcels-code/Parcels/issues/2600 is resolved.\n",
"```"
]
},
Expand Down Expand Up @@ -263,7 +263,7 @@
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(5, 3))\n",
"ax.plot(df_particles[\"lon\"], df_particles[\"lat\"], \".-\")\n",
"ax.plot(df_particles[\"x\"], df_particles[\"y\"], \".-\")\n",
"plt.show()"
]
},
Expand All @@ -283,7 +283,7 @@
"fig, ax = plt.subplots(figsize=(5, 3))\n",
"for traj in df_particles.partition_by(\"particle_id\", maintain_order=True):\n",
" traj_id = traj[\"particle_id\"][0]\n",
" ax.plot(traj[\"lon\"], traj[\"lat\"], \".-\", label=f\"P{traj_id}\")\n",
" ax.plot(traj[\"x\"], traj[\"y\"], \".-\", label=f\"P{traj_id}\")\n",
"ax.legend(loc=\"center left\", bbox_to_anchor=(1.02, 0.5), borderaxespad=0.0)\n",
"plt.tight_layout()\n",
"plt.show()"
Expand All @@ -307,7 +307,7 @@
"particles = df_particles.filter(pl.col(\"time\") == pl.lit(time_to_plot))\n",
"\n",
"fig, ax = plt.subplots(figsize=(5, 3))\n",
"ax.plot(particles[\"lon\"], particles[\"lat\"], \"o\")\n",
"ax.plot(particles[\"x\"], particles[\"y\"], \"o\")\n",
"title_time = pd.to_datetime(time_to_plot).strftime(\"%Y-%m-%d %H:%M:%S\")\n",
"ax.set_title(f\"Particle locations at {title_time}\")\n",
"plt.show()"
Expand All @@ -332,7 +332,7 @@
")\n",
"\n",
"fig, ax = plt.subplots(figsize=(5, 3))\n",
"ax.plot(particles[\"lon\"], particles[\"lat\"], \"o\")\n",
"ax.plot(particles[\"x\"], particles[\"y\"], \"o\")\n",
"ax.set_title(f\"Particle locations {time_step} after their release\")\n",
"plt.show()"
]
Expand All @@ -354,7 +354,7 @@
"\n",
"for traj in df_particles.partition_by(\"particle_id\", maintain_order=True):\n",
" distance = np.sqrt(\n",
" (traj[\"lon\"] - traj[\"lon\"][0]) ** 2 + (traj[\"lat\"] - traj[\"lat\"][0]) ** 2\n",
" (traj[\"x\"] - traj[\"x\"][0]) ** 2 + (traj[\"y\"] - traj[\"y\"][0]) ** 2\n",
" )\n",
" ax[0].plot(traj[\"time\"], distance, \".-\", label=f\"P{traj['particle_id'][0]}\")\n",
"\n",
Expand Down Expand Up @@ -439,8 +439,8 @@
"# --> plot first timestep\n",
"particles = df_particles.filter(pl.col(\"time\") == pl.lit(timerange[0]))\n",
"scatter = ax.scatter(\n",
" particles[\"lon\"],\n",
" particles[\"lat\"],\n",
" particles[\"x\"],\n",
" particles[\"y\"],\n",
" s=10,\n",
" c=[trajectory_to_color[p] for p in particles[\"particle_id\"]],\n",
")\n",
Expand All @@ -464,7 +464,7 @@
" particles = df_particles.filter(pl.col(\"time\") == pl.lit(timerange[i]))\n",
"\n",
" if len(particles) > 0:\n",
" scatter.set_offsets(np.c_[particles[\"lon\"], particles[\"lat\"]])\n",
" scatter.set_offsets(np.c_[particles[\"x\"], particles[\"y\"]])\n",
" scatter.set_color([trajectory_to_color[p] for p in particles[\"particle_id\"]])\n",
"\n",
" # --> reset trails\n",
Expand All @@ -481,8 +481,8 @@
" )\n",
" if len(traj_trail) > 1:\n",
" (trail,) = ax.plot(\n",
" traj_trail[\"lon\"],\n",
" traj_trail[\"lat\"],\n",
" traj_trail[\"x\"],\n",
" traj_trail[\"y\"],\n",
" color=trajectory_to_color[traj],\n",
" linewidth=0.6,\n",
" alpha=0.6,\n",
Expand All @@ -502,7 +502,7 @@
"metadata": {
"celltoolbar": "Metagegevens bewerken",
"kernelspec": {
"display_name": "docs",
"display_name": "Parcels:docs (3.14.6)",
"language": "python",
"name": "python3"
},
Expand All @@ -516,7 +516,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.4"
"version": "3.14.6"
}
},
"nbformat": 4,
Expand Down
16 changes: 8 additions & 8 deletions docs/getting_started/tutorial_quickstart.md
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ the virtual particles for which we will calculate the trajectories.

We need to create a {py:obj}`parcels.ParticleSet` object with the particles' initial time and position. The `parcels.ParticleSet`
object also needs to know about the `FieldSet` in which the particles "live". Finally, we need to specify the type of
{py:obj}`parcels.ParticleClass` we want to use. The default particles have `time`, `z`, `lat`, and `lon`, but you can easily add
{py:obj}`parcels.ParticleClass` we want to use. The default particles have `time`, `z`, `y`, and `x`, but you can easily add
other {py:obj}`parcels.Variable`s such as size, temperature, or age to create your own particles to mimic plastic or an [ARGO float](../user_guide/examples/tutorial_Argofloats.ipynb).

```{code-cell}
Expand All @@ -86,7 +86,7 @@ time = np.repeat(ds_fields.time.values[0], npart) # at initial time of input dat
z = np.repeat(ds_fields.depth.values[0], npart) # at the first depth (surface)

pset = parcels.ParticleSet(
fieldset=fieldset, pclass=parcels.Particle, time=time, z=z, lat=lat, lon=lon
fieldset=fieldset, pclass=parcels.Particle, time=time, z=z, y=lat, x=lon
)
```

Expand Down Expand Up @@ -168,8 +168,8 @@ Let's verify that Parcels has computed the advection of the virtual particles!
import matplotlib.pyplot as plt

# plot positions and color particles by time
scatter = plt.scatter(df['lon'], df['lat'], c=df['time'])
plt.scatter(df['lon'][:npart], df['lat'][:npart], facecolors="none", edgecolors='r') # starting positions
scatter = plt.scatter(df['x'], df['y'], c=df['time'])
plt.scatter(df['x'][:npart], df['y'][:npart], facecolors="none", edgecolors='r') # starting positions
plt.scatter(lon, lat, facecolors="none", edgecolors='r') # starting positions
plt.xlim(31,33)
plt.ylabel("Latitude [deg N]")
Expand Down Expand Up @@ -211,9 +211,9 @@ When we check the output, we can see that the particles have returned to their o
```{code-cell}
df_back = parcels.read_particlefile("output-backwards.parquet")

scatter = plt.scatter(df_back['lon'], df_back['lat'], c=df_back['time'])
scatter = plt.scatter(df_back['x'], df_back['y'], c=df_back['time'])
particles_at_max_time = df_back.filter(pl.col("time") == df_back["time"].max())
plt.scatter(particles_at_max_time['lon'], particles_at_max_time['lat'], facecolors="none", edgecolors='r') # starting positions
plt.scatter(particles_at_max_time['x'], particles_at_max_time['y'], facecolors="none", edgecolors='r') # starting positions
plt.xlabel("Longitude [deg E]")
plt.xlim(31,33)
plt.ylabel("Latitude [deg N]")
Expand All @@ -227,6 +227,6 @@ Using Euler forward advection, the final positions are equal to the original pos
```{code-cell}
# testing that final location == original location
particles_at_min_time = df_back.filter(pl.col("time") == df_back["time"].min())
np.testing.assert_almost_equal(particles_at_min_time["lat"], lat, 2)
np.testing.assert_almost_equal(particles_at_min_time['lon'], lon, 2)
np.testing.assert_almost_equal(particles_at_min_time["y"], lat, 2)
np.testing.assert_almost_equal(particles_at_min_time['x'], lon, 2)
```
28 changes: 18 additions & 10 deletions docs/user_guide/examples/explanation_interpolation.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,6 @@
Interpolation is an important functionality of Parcels. On this page we will discuss the way it is
implemented in **Parcels** and how to write a custom interpolator function.

```{note}
TODO: expand explanation (similar to Kernel loop explanation)
```

When we want to know the state of particles in an environmental field, such as temperature or velocity,
we _evaluate_ the `parcels.Field` at the particles real position in time and space (`t`, `z`, `lat`, `lon`).
In Parcels we can do this using square brackets:
Expand All @@ -26,7 +22,7 @@ If you want to sample at a different location, or time, that is not necessarily
```python
particles.temperature = fieldset.temperature[time, depth, lat, lon]
```
but this could be slower for curvilinear and unstructured because the entire grid needs to be searched.
but this could be slower for curvilinear and unstructured Grids because the entire grid needs to be searched.
````

The values of the `temperature` field at the particles' positions are determined using an interpolation
Expand All @@ -36,15 +32,23 @@ relate to the value at any point within a grid cell.
Each `parcels.Field` is defined on a (structured) `parcels.XGrid` or (unstructured) `parcels.UXGrid`.
The interpolation function takes information about the particles position relative to this grid (`grid_positions`),
as well as the values of the grid points of the `parcels.Field` in time and space, to calculate
the requested value at the particles location. Note that all grid values are available so that higher-order interpolation is possible.
the requested value at the particles location.

## Interpolator API

The interpolators included in Parcels are designed for common interpolation schemes in Parcels simulations.
If we want to add a custom interpolation method, we need to look at the interpolator API:
The interpolators included in Parcels are designed for common interpolation schemes in Parcels simulations; see the [Using the built-in interpolators tutorial](./tutorial_interpolation.ipynb).

We can write an interpolator function that takes a `parcels.Field` (or `parcels.VectorField`), a dictionary with the `particle_positions`
in real space and time, and a dictionary with the `grid_positions`.
If we want to create a custom interpolation method, we need to look at the interpolator API. Each interpolator is a class that inherits from either the `ScalarInterpolator` or `VectorInterpolator` class. The `ScalarInterpolator` class is used for scalar fields, such as temperature or salinity, while the `VectorInterpolator` class is used for vector fields, such as velocity. An interpolator class than has to have a `.interp()` method with the following signature:

```python
def interp(
self,
particle_positions: dict[str, float | np.ndarray],
grid_positions: dict[ptyping.XgridAxis, dict[str, int | float | np.ndarray]],
field: Field,
):
...
```

The `particle_positions` dictionary contains:

Expand Down Expand Up @@ -74,3 +78,7 @@ grid_positions = {
"FACE": {"index": fi, "bcoord": bcoord}
}
```

The `.interp()` method should return a float (in the case o a `ScalarInterpolator` or a tuple of three floats `(u, v, w)` in the case of a `VectorInterpolator`).

Writing custom interpolators is not trivial, so we recommend that you have a look at the built-in interpolators in `parcels.interpolators._xinterpolators` or `parcels.interpolators._uxinterpolators` to see how they are implemented.
Loading
Loading