Skip to content

Commit e51196c

Browse files
reint-fischerreint-fischer
authored andcommitted
complete writing basic quickstart tutorial
1 parent a135b45 commit e51196c

File tree

1 file changed

+73
-28
lines changed

1 file changed

+73
-28
lines changed

docs/getting_started/tutorial_quickstart.md

Lines changed: 73 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,6 @@ kernelspec:
66

77
# Quickstart tutorial
88

9-
```{note}
10-
TODO: Completely rewrite examples/parcels_tutorial.ipynb to be this quickstart tutorial. Decide which file format and notebook testing to do so this file is checked, which is not in the "examples" folder
11-
```
12-
139
Welcome to the **Parcels** quickstart tutorial, in which we will go through all the necessary steps to run a simulation. The code in this notebook can be used as a starting point to run Parcels in your own environment. Along the way we will familiarize ourselves with some specific classes and methods. If you are ever confused about one of these and want to read more, we have a [concepts overview](concepts_overview.md) discussing them in more detail. Let's dive in!
1410

1511
## Imports
@@ -22,7 +18,7 @@ import xarray as xr
2218
import parcels
2319
```
2420

25-
## Input 1: `FieldSet`
21+
## Input flow fields: `FieldSet`
2622

2723
A Parcels simulation of Lagrangian trajectories of virtual particles requires two inputs; the first is a set of hydrodynamics fields in which the particles are tracked. This set of vectorfields, with `U`, `V` (and `W`) flow velocities, can be read in to a `parcels.FieldSet` object from many types of models or observations. Here we provide an example using a subset of the [Global Ocean Physics Reanalysis](https://doi.org/10.48670/moi-00021) from the Copernicus Marine Service.
2824

@@ -33,81 +29,130 @@ example_dataset_folder = parcels.download_example_dataset(
3329
3430
ds_in = xr.open_mfdataset(f"{example_dataset_folder}/*.nc", combine="by_coords")
3531
ds_in.load() # load the dataset into memory
36-
37-
fieldset = parcels.FieldSet.from_copernicusmarine(ds_in)
32+
ds_in
3833
```
3934

40-
The reanalysis data files contain `U`, `V`, potential temperature (`thetao`) and salinity (`so`) fields:
35+
As we can see, the reanalysis dataset contains eastward velocity `uo`, northward velocity `vo`, potential temperature (`thetao`) and salinity (`so`) fields. To load the dataset into parcels we create a `FieldSet`, which recognizes the standard names of a velocity field:
4136

4237
```{code-cell}
43-
ds_in
38+
fieldset = parcels.FieldSet.from_copernicusmarine(ds_in)
4439
```
4540

4641
The subset contains a region of the Agulhas current along the southeastern coast of Africa:
4742

4843
```{code-cell}
49-
import matplotlib.pyplot as plt
50-
import cartopy.crs as ccrs
51-
52-
temperature = ds_in.isel(time=0, depth=0).thetao.plot(cmap="magma",subplot_kws=dict(projection=ccrs.PlateCarree()))
44+
temperature = ds_in.isel(time=0, depth=0).thetao.plot(cmap="magma")
5345
velocity = ds_in.isel(time=0, depth=0).plot.quiver(x="longitude", y="latitude", u="uo", v="vo")
54-
55-
temperature.axes.set_extent([30, 34, -34, -29])
56-
gl = temperature.axes.gridlines(draw_labels=True)
57-
gl.top_labels = False
58-
gl.right_labels = False
59-
temperature.axes.coastlines()
6046
```
6147

62-
## Input 2: `ParticleSet`
48+
## Input virtual particles: `ParticleSet`
6349

6450
Now that we have read in the hydrodynamic data, we need to provide our second input: the virtual particles for which we will calculate the trajectories. We need to define the initial time and position and read those into a `parcels.ParticleSet` object, which also needs to know about the `FieldSet` in which the particles "live". Finally, we need to specify the type of `parcels.Particle` we want to use. The default particles have `time`, `lon`, `lat`, and `z` to keep track of, but with Parcels you can easily build your own particles to mimic plastic or an [ARGO float](../user_guide/tutorial_Argofloats.ipynb), adding variables such as size, temperature, or age.
6551

6652
```{code-cell}
6753
# Particle locations and initial time
6854
npart = 10 # number of particles to be released
6955
lon = np.repeat(32, npart)
70-
lat = np.linspace(-32.5, -30.5, npart)
71-
time = np.repeat(ds.time.values[0], npart)
72-
z = np.repeat(ds.depth.values[0], npart)
56+
lat = np.linspace(-32.5, -30.5, npart) # release particles in a line along a meridian
57+
time = np.repeat(ds_in.time.values[0], npart) # at initial time of input data
58+
z = np.repeat(ds_in.depth.values[0], npart) # at the first depth (surface)
7359
7460
pset = parcels.ParticleSet(
7561
fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z
7662
)
7763
```
7864

65+
```{code-cell}
66+
temperature = ds_in.isel(time=0, depth=0).thetao.plot(cmap="magma")
67+
velocity = ds_in.isel(time=0, depth=0).plot.quiver(x="longitude", y="latitude", u="uo", v="vo")
68+
ax = temperature.axes
69+
ax.scatter(lon,lat,s=40,c='w',edgecolors='r')
70+
```
71+
7972
## Compute: `Kernel`
8073

81-
After setting up the input data, we need to specify how to calculate the advection of the particles. These calculations, or numerical integrations, will be performed by what we call a `Kernel`, operating on each particle in the `ParticleSet`. The most common calculation is the advection of particles through the velocity field. Parcels comes with a number of standard kernels, from which we will use the Euler-forward advection kernel:
74+
After setting up the input data, we need to specify how to calculate the advection of the particles. These calculations, or numerical integrations, will be performed by what we call a `Kernel`, operating on each particle in the `ParticleSet`. The most common calculation is the advection of particles through the velocity field. Parcels comes with a number of standard kernels, from which we will use the Euler-forward advection kernel `AdvectionEE`:
8275

8376
```{note}
84-
TODO: link to list of included kernels
77+
TODO: link to a list of included kernels
8578
```
8679

8780
```{code-cell}
8881
kernels = parcels.kernels.AdvectionEE
8982
```
9083

9184
## Prepare output: `ParticleFile`
85+
Before starting the simulation, we must define where and how frequent we want to write the output of our simulation. We can define this in a `ParticleFile` object:
9286

9387
```{code-cell}
9488
output_file = parcels.ParticleFile("Output.zarr", outputdt=np.timedelta64(1, "h"))
9589
```
9690

91+
The output files are in `.zarr` [format]([format](https://zarr.readthedocs.io/en/stable/).), which can be read by `xarray`. See the [Parcels output tutorial](../user_guide/examples/tutorial_output.ipynb) for more information on the zarr format. We want to choose the `outputdt` argument such they capture the smallest timescales of our interest.
92+
9793
## Run Simulation: `ParticleSet.execute()`
94+
Finally, we can run the simulation by *executing* the `ParticleSet` using the specified `kernels`. Additionally, we need to specify:
95+
- the `runtime`: for how long we want to simulate particles.
96+
- the `dt`: the timestep with which to perform the numerical integration in the `kernels`. Depending on the numerical integration scheme, the accuracy of our simulation will depend on `dt`. Read [this notebook](https://github.com/Parcels-code/10year-anniversary-session2/blob/8931ef69577dbf00273a5ab4b7cf522667e146c5/advection_and_windage.ipynb) to learn more about numerical accuracy.
97+
98+
```{note}
99+
TODO: add Michaels 10-years Parcels notebook to the user guide
100+
```
98101

99102
```{code-cell}
103+
:tags: [hide-output]
100104
pset.execute(
101105
kernels,
102-
runtime=np.timedelta64(5, "h"),
106+
runtime=np.timedelta64(1, "D"),
103107
dt=np.timedelta64(5, "m"),
104108
output_file=output_file,
105109
)
106110
```
107111

108112
## Read output
113+
To start analyzing the trajectories computed by **Parcels**, we can open the `ParticleFile` using `xarray`:
109114

110115
```{code-cell}
111-
data_xarray = xr.open_zarr("Output.zarr")
112-
data_xarray
116+
ds_out = xr.open_zarr("Output.zarr")
117+
ds_out
113118
```
119+
120+
The 10 particle trajectories are stored along the `trajectory` dimension, and each trajectory contains 25 observations (initial values + 24 hourly timesteps) along the `obs` dimension. The [Working with Parcels output tutorial](../user_guide/examples/tutorial_output.ipynb) provides more detail about the dataset and how to analyse it. Let's verify that Parcels has computed the advection of the virtual particles!
121+
122+
```{code-cell}
123+
import matplotlib.pyplot as plt
124+
125+
# plot positions and color particles by number of observation
126+
plt.scatter(ds_out.lon.T, ds_out.lat.T, c=np.repeat(ds_out.obs.values,npart))
127+
plt.xlabel("Longitude [deg E]")
128+
plt.ylabel("Latitude [deg N]")
129+
plt.show()
130+
```
131+
That looks good! The virtual particles released in a line along the 32nd meridian (dark blue) have been advected by the flow field.
132+
133+
## Running a simulation backwards in time
134+
Now that we know how to run a simulation, we can easily run another and change one of the settings. We can trace back the particles from their current to their original position by running the simulation backwards in time. To do so, we can simply make `dt` < 0.
135+
136+
```{code-cell}
137+
:tags: [hide-output]
138+
# set up output file
139+
output_file = parcels.ParticleFile("Output-backwards.zarr", outputdt=np.timedelta64(1, "h"))
140+
141+
# execute simulation in backwards time
142+
pset.execute(
143+
kernels,
144+
runtime=np.timedelta64(1, "D"),
145+
dt=-np.timedelta64(5, "m"),
146+
output_file=output_file,
147+
)
148+
```
149+
When we check the output, we can see that the particles have returned to their original position!
150+
151+
```{code-cell}
152+
ds_out = xr.open_zarr("Output-backwards.zarr")
153+
154+
plt.scatter(ds_out.lon.T, ds_out.lat.T, c=np.repeat(ds_out.obs.values,npart))
155+
plt.xlabel("Longitude [deg E]")
156+
plt.ylabel("Latitude [deg N]")
157+
plt.show()
158+
```

0 commit comments

Comments
 (0)