Skip to content

Commit a4cb13d

Browse files
erikvansebilleerikvansebille
authored
Merge pull request #2199 from OceanParcels/v4-argo-float-tutorial
Updating Argo tutorial to v4
2 parents 97cd0c3 + 8819aa2 commit a4cb13d

File tree

5 files changed

+243
-94
lines changed

5 files changed

+243
-94
lines changed

‎docs/community/v4-migration.md‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,10 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
1919

2020
- `interp_method` has to be an Interpolation function, instead of a string.
2121

22+
## Particle
23+
24+
- `Particle.add_variables()` has been replaced by `Particle.add_variable()`, which now also takes a list of `Variables`.
25+
2226
## ParticleSet
2327

2428
- `repeatdt` and `lonlatdepth_dtype` have been removed from the ParticleSet.
@@ -28,3 +32,4 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
2832
## ParticleFile
2933

3034
- Particlefiles should be created by `ParticleFile(...)` instead of `pset.ParticleFile(...)`
35+
- The `name` argument in `ParticleFile` has been replaced by `store` and can now be a string, a Path or a zarr store.

‎docs/examples/tutorial_Argofloats.ipynb‎

Lines changed: 171 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,9 @@
1313
"cell_type": "markdown",
1414
"metadata": {},
1515
"source": [
16-
"This tutorial shows how simple it is to construct a Kernel in Parcels that mimics the [vertical movement of Argo floats](https://www.aoml.noaa.gov/phod/argo/images/argo_float_mission.jpg).\n"
16+
"This tutorial shows how simple it is to construct a Kernel in Parcels that mimics the [vertical movement of Argo floats](https://www.aoml.noaa.gov/phod/argo/images/argo_float_mission.jpg).\n",
17+
"\n",
18+
"We first define the kernels for each phase of the Argo cycle."
1719
]
1820
},
1921
{
@@ -22,52 +24,90 @@
2224
"metadata": {},
2325
"outputs": [],
2426
"source": [
25-
"# Define the new Kernel that mimics Argo vertical movement\n",
26-
"def ArgoVerticalMovement(particle, fieldset, time):\n",
27-
" driftdepth = 1000 # maximum depth in m\n",
28-
" maxdepth = 2000 # maximum depth in m\n",
29-
" vertical_speed = 0.10 # sink and rise speed in m/s\n",
30-
" cycletime = 10 * 86400 # total time of cycle in seconds\n",
31-
" drifttime = 9 * 86400 # time of deep drift in seconds\n",
32-
"\n",
33-
" if particle.cycle_phase == 0:\n",
34-
" # Phase 0: Sinking with vertical_speed until depth is driftdepth\n",
35-
" particle_ddepth += vertical_speed * particle.dt\n",
36-
" if particle.depth + particle_ddepth >= driftdepth:\n",
37-
" particle_ddepth = driftdepth - particle.depth\n",
38-
" particle.cycle_phase = 1\n",
39-
"\n",
40-
" elif particle.cycle_phase == 1:\n",
41-
" # Phase 1: Drifting at depth for drifttime seconds\n",
42-
" particle.drift_age += particle.dt\n",
43-
" if particle.drift_age >= drifttime:\n",
44-
" particle.drift_age = 0 # reset drift_age for next cycle\n",
45-
" particle.cycle_phase = 2\n",
46-
"\n",
47-
" elif particle.cycle_phase == 2:\n",
48-
" # Phase 2: Sinking further to maxdepth\n",
49-
" particle_ddepth += vertical_speed * particle.dt\n",
50-
" if particle.depth + particle_ddepth >= maxdepth:\n",
51-
" particle_ddepth = maxdepth - particle.depth\n",
52-
" particle.cycle_phase = 3\n",
53-
"\n",
54-
" elif particle.cycle_phase == 3:\n",
55-
" # Phase 3: Rising with vertical_speed until at surface\n",
56-
" particle_ddepth -= vertical_speed * particle.dt\n",
57-
" # particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon] # if fieldset has temperature\n",
58-
" if particle.depth + particle_ddepth <= fieldset.mindepth:\n",
59-
" particle_ddepth = fieldset.mindepth - particle.depth\n",
60-
" # particle.temp = 0./0. # reset temperature to NaN at end of sampling cycle\n",
61-
" particle.cycle_phase = 4\n",
62-
"\n",
63-
" elif particle.cycle_phase == 4:\n",
64-
" # Phase 4: Transmitting at surface until cycletime is reached\n",
65-
" if particle.cycle_age > cycletime:\n",
66-
" particle.cycle_phase = 0\n",
67-
" particle.cycle_age = 0\n",
68-
"\n",
69-
" if particle.state == StatusCode.Evaluate:\n",
70-
" particle.cycle_age += particle.dt # update cycle_age"
27+
"import numpy as np\n",
28+
"\n",
29+
"# Define the new Kernels that mimic Argo vertical movement\n",
30+
"driftdepth = 1000 # maximum depth in m\n",
31+
"maxdepth = 2000 # maximum depth in m\n",
32+
"vertical_speed = 0.10 # sink and rise speed in m/s\n",
33+
"cycletime = (\n",
34+
" 10 * 86400\n",
35+
") # total time of cycle in seconds # TODO update to \"timedelta64[s]\"\n",
36+
"drifttime = 9 * 86400 # time of deep drift in seconds\n",
37+
"\n",
38+
"\n",
39+
"def ArgoPhase1(particles, fieldset):\n",
40+
" dt = particles.dt / np.timedelta64(1, \"s\") # convert dt to seconds\n",
41+
"\n",
42+
" def SinkingPhase(p):\n",
43+
" \"\"\"Phase 0: Sinking with vertical_speed until depth is driftdepth\"\"\"\n",
44+
" p.ddepth += vertical_speed * dt\n",
45+
" p.cycle_phase = np.where(p.depth + p.ddepth >= driftdepth, 1, p.cycle_phase)\n",
46+
" p.ddepth = np.where(\n",
47+
" p.depth + p.ddepth >= driftdepth, driftdepth - p.depth, p.ddepth\n",
48+
" )\n",
49+
"\n",
50+
" SinkingPhase(particles[particles.cycle_phase == 0])\n",
51+
"\n",
52+
"\n",
53+
"def ArgoPhase2(particles, fieldset):\n",
54+
" dt = particles.dt / np.timedelta64(1, \"s\") # convert dt to seconds\n",
55+
"\n",
56+
" def DriftingPhase(p):\n",
57+
" \"\"\"Phase 1: Drifting at depth for drifttime seconds\"\"\"\n",
58+
" p.drift_age += dt\n",
59+
" p.cycle_phase = np.where(p.drift_age >= drifttime, 2, p.cycle_phase)\n",
60+
" p.drift_age = np.where(p.drift_age >= drifttime, 0, p.drift_age)\n",
61+
"\n",
62+
" DriftingPhase(particles[particles.cycle_phase == 1])\n",
63+
"\n",
64+
"\n",
65+
"def ArgoPhase3(particles, fieldset):\n",
66+
" dt = particles.dt / np.timedelta64(1, \"s\") # convert dt to seconds\n",
67+
"\n",
68+
" def SecondSinkingPhase(p):\n",
69+
" \"\"\"Phase 2: Sinking further to maxdepth\"\"\"\n",
70+
" p.ddepth += vertical_speed * dt\n",
71+
" p.cycle_phase = np.where(p.depth + p.ddepth >= maxdepth, 3, p.cycle_phase)\n",
72+
" p.ddepth = np.where(\n",
73+
" p.depth + p.ddepth >= maxdepth, maxdepth - p.depth, p.ddepth\n",
74+
" )\n",
75+
"\n",
76+
" SecondSinkingPhase(particles[particles.cycle_phase == 2])\n",
77+
"\n",
78+
"\n",
79+
"def ArgoPhase4(particles, fieldset):\n",
80+
" dt = particles.dt / np.timedelta64(1, \"s\") # convert dt to seconds\n",
81+
"\n",
82+
" def RisingPhase(p):\n",
83+
" \"\"\"Phase 3: Rising with vertical_speed until at surface\"\"\"\n",
84+
" p.ddepth -= vertical_speed * dt\n",
85+
" p.temp = fieldset.thetao[p.time, p.depth, p.lat, p.lon]\n",
86+
" p.cycle_phase = np.where(\n",
87+
" p.depth + p.ddepth <= fieldset.mindepth, 4, p.cycle_phase\n",
88+
" )\n",
89+
" p.ddepth = np.where(\n",
90+
" p.depth + p.ddepth <= fieldset.mindepth,\n",
91+
" fieldset.mindepth - p.depth,\n",
92+
" p.ddepth,\n",
93+
" )\n",
94+
"\n",
95+
" RisingPhase(particles[particles.cycle_phase == 3])\n",
96+
"\n",
97+
"\n",
98+
"def ArgoPhase5(particles, fieldset):\n",
99+
" def TransmittingPhase(p):\n",
100+
" \"\"\"Phase 4: Transmitting at surface until cycletime is reached\"\"\"\n",
101+
" p.cycle_phase = np.where(p.cycle_age >= cycletime, 0, p.cycle_phase)\n",
102+
" p.cycle_age = np.where(p.cycle_age >= cycletime, 0, p.cycle_age)\n",
103+
" p.temp = np.nan # no temperature measurement when at surface\n",
104+
"\n",
105+
" TransmittingPhase(particles[particles.cycle_phase == 4])\n",
106+
"\n",
107+
"\n",
108+
"def ArgoPhase6(particles, fieldset):\n",
109+
" dt = particles.dt / np.timedelta64(1, \"s\") # convert dt to seconds\n",
110+
" particles.cycle_age += dt # update cycle_age"
71111
]
72112
},
73113
{
@@ -77,11 +117,7 @@
77117
"source": [
78118
"And then we can run Parcels with this 'custom kernel'.\n",
79119
"\n",
80-
"Note that below we use the two-dimensional velocity fields of GlobCurrent, as these are provided as example_data with Parcels.\n",
81-
"\n",
82-
"We therefore assume that the horizontal velocities are the same throughout the entire water column. However, the `ArgoVerticalMovement` kernel will work on any `FieldSet`, including from full three-dimensional hydrodynamic data.\n",
83-
"\n",
84-
"If the hydrodynamic data also has a Temperature Field, then uncommenting the lines about temperature will also simulate the sampling of temperature.\n"
120+
"Below we use the horizontal velocity fields of CopernicusMarine, which are provided as example_data with Parcels.\n"
85121
]
86122
},
87123
{
@@ -92,54 +128,58 @@
92128
"source": [
93129
"from datetime import timedelta\n",
94130
"\n",
95-
"import numpy as np\n",
131+
"import xarray as xr\n",
96132
"\n",
97133
"import parcels\n",
98134
"\n",
99-
"# Load the GlobCurrent data in the Agulhas region from the example_data\n",
100-
"example_dataset_folder = parcels.download_example_dataset(\"GlobCurrent_example_data\")\n",
101-
"filenames = {\n",
102-
" \"U\": f\"{example_dataset_folder}/20*.nc\",\n",
103-
" \"V\": f\"{example_dataset_folder}/20*.nc\",\n",
104-
"}\n",
105-
"variables = {\n",
106-
" \"U\": \"eastward_eulerian_current_velocity\",\n",
107-
" \"V\": \"northward_eulerian_current_velocity\",\n",
108-
"}\n",
109-
"dimensions = {\"lat\": \"lat\", \"lon\": \"lon\", \"time\": \"time\"}\n",
110-
"fieldset = parcels.FieldSet.from_netcdf(filenames, variables, dimensions)\n",
111-
"# uppermost layer in the hydrodynamic data\n",
112-
"fieldset.mindepth = fieldset.U.depth[0]\n",
135+
"# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n",
136+
"example_dataset_folder = parcels.download_example_dataset(\n",
137+
" \"CopernicusMarine_data_for_Argo_tutorial\"\n",
138+
")\n",
139+
"\n",
140+
"ds = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n",
141+
"\n",
142+
"# TODO check how we can get good performance without loading full dataset in memory\n",
143+
"ds.load() # load the dataset into memory\n",
113144
"\n",
145+
"fieldset = parcels.FieldSet.from_copernicusmarine(ds)\n",
146+
"fieldset.add_constant(\"mindepth\", 1.0)\n",
114147
"\n",
115148
"# Define a new Particle type including extra Variables\n",
116-
"ArgoParticle = parcels.Particle.add_variables(\n",
149+
"ArgoParticle = parcels.Particle.add_variable(\n",
117150
" [\n",
118-
" # Phase of cycle:\n",
119-
" # init_descend=0,\n",
120-
" # drift=1,\n",
121-
" # profile_descend=2,\n",
122-
" # profile_ascend=3,\n",
123-
" # transmit=4\n",
124151
" parcels.Variable(\"cycle_phase\", dtype=np.int32, initial=0.0),\n",
125-
" parcels.Variable(\"cycle_age\", dtype=np.float32, initial=0.0),\n",
152+
" parcels.Variable(\n",
153+
" \"cycle_age\", dtype=np.float32, initial=0.0\n",
154+
" ), # TODO update to \"timedelta64[s]\"\n",
126155
" parcels.Variable(\"drift_age\", dtype=np.float32, initial=0.0),\n",
127-
" # if fieldset has temperature\n",
128-
" # Variable('temp', dtype=np.float32, initial=np.nan),\n",
156+
" parcels.Variable(\"temp\", dtype=np.float32, initial=np.nan),\n",
129157
" ]\n",
130158
")\n",
131159
"\n",
132160
"# Initiate one Argo float in the Agulhas Current\n",
133161
"pset = parcels.ParticleSet(\n",
134-
" fieldset=fieldset, pclass=ArgoParticle, lon=[32], lat=[-31], depth=[0]\n",
162+
" fieldset=fieldset,\n",
163+
" pclass=ArgoParticle,\n",
164+
" lon=[32],\n",
165+
" lat=[-31],\n",
166+
" depth=[fieldset.mindepth],\n",
135167
")\n",
136168
"\n",
137169
"# combine Argo vertical movement kernel with built-in Advection kernel\n",
138-
"kernels = [ArgoVerticalMovement, parcels.AdvectionRK4]\n",
170+
"kernels = [\n",
171+
" ArgoPhase1,\n",
172+
" ArgoPhase2,\n",
173+
" ArgoPhase3,\n",
174+
" ArgoPhase4,\n",
175+
" ArgoPhase5,\n",
176+
" ArgoPhase6,\n",
177+
" parcels.AdvectionRK4,\n",
178+
"]\n",
139179
"\n",
140180
"# Create a ParticleFile object to store the output\n",
141-
"output_file = pset.ParticleFile(\n",
142-
" name=\"argo_float\",\n",
181+
"output_file = parcels.ParticleFile(\n",
182+
" store=\"argo_float.zarr\",\n",
143183
" outputdt=timedelta(minutes=15),\n",
144184
" chunks=(1, 500), # setting to write in chunks of 500 observations\n",
145185
")\n",
@@ -158,7 +198,25 @@
158198
"cell_type": "markdown",
159199
"metadata": {},
160200
"source": [
161-
"Now we can plot the trajectory of the Argo float with some simple calls to netCDF4 and matplotlib\n"
201+
"Now we can plot the trajectory of the Argo float with some simple calls to netCDF4 and matplotlib.\n",
202+
"\n",
203+
"First plot the depth as a function of time, with the temperature as color (only on the upcast)."
204+
]
205+
},
206+
{
207+
"cell_type": "code",
208+
"execution_count": null,
209+
"metadata": {},
210+
"outputs": [],
211+
"source": [
212+
"ds_out = xr.open_zarr(\n",
213+
" output_file.store, decode_times=False\n",
214+
") # TODO fix without using decode_times=False\n",
215+
"x = ds_out[\"lon\"][:].squeeze()\n",
216+
"y = ds_out[\"lat\"][:].squeeze()\n",
217+
"z = ds_out[\"z\"][:].squeeze()\n",
218+
"time = ds_out[\"time\"][:].squeeze() / 86400 # convert time to days\n",
219+
"temp = ds_out[\"temp\"][:].squeeze()"
162220
]
163221
},
164222
{
@@ -168,29 +226,49 @@
168226
"outputs": [],
169227
"source": [
170228
"import matplotlib.pyplot as plt\n",
171-
"import xarray as xr\n",
172-
"from mpl_toolkits.mplot3d import Axes3D\n",
173229
"\n",
174-
"ds = xr.open_zarr(\"argo_float.zarr\")\n",
175-
"x = ds[\"lon\"][:].squeeze()\n",
176-
"y = ds[\"lat\"][:].squeeze()\n",
177-
"z = ds[\"z\"][:].squeeze()\n",
178-
"ds.close()\n",
230+
"fig = plt.figure(figsize=(13, 6))\n",
231+
"ax = plt.axes()\n",
232+
"ax.plot(time, z, color=\"gray\")\n",
233+
"cb = ax.scatter(time, z, c=temp, s=20, marker=\"o\", zorder=2)\n",
234+
"ax.set_xlabel(\"Time [days]\")\n",
235+
"ax.set_ylabel(\"Depth (m)\")\n",
236+
"ax.invert_yaxis()\n",
237+
"fig.colorbar(cb, label=\"Temperature (°C)\")\n",
238+
"plt.show()"
239+
]
240+
},
241+
{
242+
"cell_type": "markdown",
243+
"metadata": {},
244+
"source": [
245+
"We can also make a 3D plot of the trajectory colored by temperature."
246+
]
247+
},
248+
{
249+
"cell_type": "code",
250+
"execution_count": null,
251+
"metadata": {},
252+
"outputs": [],
253+
"source": [
254+
"from mpl_toolkits.mplot3d import Axes3D\n",
179255
"\n",
180-
"fig = plt.figure(figsize=(13, 10))\n",
256+
"fig = plt.figure(figsize=(13, 8))\n",
181257
"ax = plt.axes(projection=\"3d\")\n",
182-
"cb = ax.scatter(x, y, z, c=z, s=20, marker=\"o\")\n",
258+
"ax.plot3D(x, y, z, color=\"gray\")\n",
259+
"cb = ax.scatter(x, y, z, c=temp, s=20, marker=\"o\", zorder=2)\n",
183260
"ax.set_xlabel(\"Longitude\")\n",
184261
"ax.set_ylabel(\"Latitude\")\n",
185262
"ax.set_zlabel(\"Depth (m)\")\n",
186263
"ax.set_zlim(np.max(z), 0)\n",
264+
"fig.colorbar(cb, label=\"Temperature (°C)\")\n",
187265
"plt.show()"
188266
]
189267
}
190268
],
191269
"metadata": {
192270
"kernelspec": {
193-
"display_name": "parcels",
271+
"display_name": "parcels-v4",
194272
"language": "python",
195273
"name": "python3"
196274
},
@@ -204,7 +282,7 @@
204282
"name": "python",
205283
"nbconvert_exporter": "python",
206284
"pygments_lexer": "ipython3",
207-
"version": "3.12.3"
285+
"version": "3.12.11"
208286
}
209287
},
210288
"nbformat": 4,

‎parcels/particleset.py‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import numpy as np
88
import xarray as xr
99
from tqdm import tqdm
10+
from zarr.storage import DirectoryStore
1011

1112
from parcels._core.utils.time import TimeInterval, maybe_convert_python_timedelta_to_numpy
1213
from parcels._reprs import particleset_repr
@@ -522,7 +523,7 @@ def execute(
522523

523524
# Set up pbar
524525
if output_file:
525-
logger.info(f"Output files are stored in {output_file.store}.")
526+
logger.info(f"Output files are stored in {_format_output_location(output_file.store)}")
526527

527528
if verbose_progress:
528529
pbar = tqdm(total=(end_time - start_time) / np.timedelta64(1, "s"), file=sys.stdout)
@@ -636,3 +637,9 @@ def _get_start_time(first_release_time, time_interval, sign_dt, runtime):
636637

637638
start_time = first_release_time if not np.isnat(first_release_time) else fieldset_start
638639
return start_time
640+
641+
642+
def _format_output_location(zarr_obj):
643+
if isinstance(zarr_obj, DirectoryStore):
644+
return zarr_obj.path
645+
return repr(zarr_obj)

‎parcels/tools/exampledata_utils.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,10 @@
4545
f"{date.strftime('%Y%m%d')}000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc"
4646
for date in ([datetime(2002, 1, 1) + timedelta(days=x) for x in range(0, 365)] + [datetime(2003, 1, 1)])
4747
],
48+
"CopernicusMarine_data_for_Argo_tutorial": [
49+
"cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m_uo-vo_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc",
50+
"cmems_mod_glo_phy-thetao_anfc_0.083deg_P1D-m_thetao_31.00E-33.00E_33.00S-30.00S_0.49-2225.08m_2024-01-01-2024-02-01.nc",
51+
],
4852
"DecayingMovingEddy_data": [
4953
"decaying_moving_eddyU.nc",
5054
"decaying_moving_eddyV.nc",

0 commit comments

Comments
 (0)