Skip to content

Commit 770d424

Browse files
reint-fischerreint-fischer
authored
Merge pull request #2317 from reint-fischer/sampling-tutorial
Add Field sampling tutorial v4
2 parents a20e7e2 + 4493522 commit 770d424

40 files changed

+93
-71
lines changed

‎docs/documentation/index.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4.
4242
:name: tutorial-kernels
4343
4444
<!-- ../examples/tutorial_diffusion.ipynb -->
45-
<!-- ../examples/tutorial_sampling.ipynb -->
45+
../examples/tutorial_sampling.ipynb
4646
<!-- ../examples/tutorial_particle_field_interaction.ipynb -->
4747
<!-- ../examples/tutorial_interaction.ipynb -->
4848
<!-- ../examples/tutorial_analyticaladvection.ipynb -->

‎docs/examples/tutorial_Argofloats.ipynb‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -268,7 +268,7 @@
268268
],
269269
"metadata": {
270270
"kernelspec": {
271-
"display_name": "parcels-v4",
271+
"display_name": "test-notebooks",
272272
"language": "python",
273273
"name": "python3"
274274
},

‎docs/examples/tutorial_sampling.ipynb‎

Lines changed: 90 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@
5656
"cell_type": "markdown",
5757
"metadata": {},
5858
"source": [
59-
"Suppose we want to study the environmental temperature for plankton drifting around a peninsula. We have a dataset with surface ocean velocities and the corresponding sea surface temperature stored in netcdf files in the folder `\"Peninsula_data\"`. Besides the velocity fields, we load the temperature field using `extra_fields={'T': 'T'}`. The particles are released on the left hand side of the domain.\n"
59+
"Suppose we want to study the environmental temperature for plankton drifting in the Agulhas current. We have a CopernicusMarine dataset with surface ocean velocities and the corresponding potential temperature stored in netcdf files in the [parcels example dataset repository](https://github.com/OceanParcels/parcels-data). Loading in the FieldSet, parcels detects U and V because they have CF standard names and tells us that they are assumed as the velocity fields to be used in the simulation.\n"
6060
]
6161
},
6262
{
@@ -65,33 +65,49 @@
6565
"metadata": {},
6666
"outputs": [],
6767
"source": [
68-
"# Velocity and temperature fields\n",
69-
"example_dataset_folder = parcels.download_example_dataset(\"Peninsula_data\")\n",
70-
"filenames = {\n",
71-
" \"U\": str(example_dataset_folder / \"peninsulaU.nc\"),\n",
72-
" \"V\": str(example_dataset_folder / \"peninsulaV.nc\"),\n",
73-
" \"T\": str(example_dataset_folder / \"peninsulaT.nc\"),\n",
74-
"}\n",
75-
"variables = {\"U\": \"vozocrtx\", \"V\": \"vomecrty\", \"T\": \"T\"}\n",
76-
"dimensions = {\"lon\": \"nav_lon\", \"lat\": \"nav_lat\", \"time\": \"time_counter\"}\n",
77-
"fieldset = parcels.FieldSet.from_netcdf(\n",
78-
" filenames, variables, dimensions, allow_time_extrapolation=True\n",
68+
"# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n",
69+
"example_dataset_folder = parcels.download_example_dataset(\n",
70+
" \"CopernicusMarine_data_for_Argo_tutorial\"\n",
7971
")\n",
8072
"\n",
73+
"ds = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n",
74+
"ds.load() # load the dataset into memory\n",
75+
"\n",
76+
"fieldset = parcels.FieldSet.from_copernicusmarine(ds)"
77+
]
78+
},
79+
{
80+
"cell_type": "markdown",
81+
"metadata": {},
82+
"source": [
83+
" Ten particles are initialized at the surface in the center of our domain, at the initial time step."
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {},
90+
"outputs": [],
91+
"source": [
8192
"# Particle locations and initial time\n",
8293
"npart = 10 # number of particles to be released\n",
83-
"lon = 3e3 * np.ones(npart)\n",
84-
"lat = np.linspace(3e3, 45e3, npart, dtype=np.float32)\n",
85-
"time = (\n",
86-
" np.arange(0, npart) * timedelta(hours=2).total_seconds()\n",
87-
") # release each particle two hours later\n",
94+
"lon = 32 * np.ones(npart)\n",
95+
"lat = np.linspace(-32.5, -30.5, npart, dtype=np.float32)\n",
96+
"time = np.repeat(\n",
97+
" ds.time.values[0], npart\n",
98+
") # release all particles at the start time of the fieldset\n",
99+
"depth = np.repeat(ds.depth.values[0], npart)\n",
88100
"\n",
89101
"# Plot temperature field and initial particle locations\n",
90-
"T_data = xr.open_dataset(f\"{example_dataset_folder}/peninsulaT.nc\")\n",
91102
"plt.figure()\n",
92103
"ax = plt.axes()\n",
93104
"T_contour = ax.contourf(\n",
94-
" T_data.x.values, T_data.y.values, T_data.T.values[0, 0], cmap=plt.cm.inferno\n",
105+
" ds.longitude.values,\n",
106+
" ds.latitude.values,\n",
107+
" ds.thetao.values[0, 0],\n",
108+
" cmap=plt.cm.inferno,\n",
109+
" vmin=20,\n",
110+
" vmax=27,\n",
95111
")\n",
96112
"ax.scatter(lon, lat, c=\"w\")\n",
97113
"plt.colorbar(T_contour, label=r\"T [$^{\\circ} C$]\")\n",
@@ -112,15 +128,15 @@
112128
"metadata": {},
113129
"outputs": [],
114130
"source": [
115-
"SampleParticle = parcels.Particle.add_variable(\"temperature\")\n",
116-
"\n",
117-
"pset = parcels.ParticleSet(\n",
118-
" fieldset=fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=time\n",
131+
"SampleParticle = parcels.Particle.add_variable(\n",
132+
" parcels.Variable(\"temperature\", dtype=np.float32, initial=np.nan)\n",
119133
")\n",
120134
"\n",
121135
"\n",
122-
"def SampleT(particle, fieldset, time):\n",
123-
" particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]"
136+
"def SampleT(particles, fieldset):\n",
137+
" particles.temperature = fieldset.thetao[\n",
138+
" particles.time, particles.depth, particles.lat, particles.lon\n",
139+
" ]"
124140
]
125141
},
126142
{
@@ -138,10 +154,10 @@
138154
"outputs": [],
139155
"source": [
140156
"pset = parcels.ParticleSet(\n",
141-
" fieldset=fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=time\n",
157+
" fieldset=fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=time, depth=depth\n",
142158
")\n",
143159
"\n",
144-
"output_file = pset.ParticleFile(name=\"SampleTemp.zarr\", outputdt=timedelta(hours=1))\n",
160+
"output_file = parcels.ParticleFile(\"SampleTemp.zarr\", outputdt=timedelta(hours=1))\n",
145161
"\n",
146162
"pset.execute(\n",
147163
" [parcels.kernels.AdvectionRK4, SampleT],\n",
@@ -165,21 +181,19 @@
165181
"metadata": {},
166182
"outputs": [],
167183
"source": [
168-
"Particle_data = xr.open_zarr(\"SampleTemp.zarr\")\n",
184+
"Particle_data = xr.open_zarr(\"SampleTemp.zarr\", decode_times=False)\n",
169185
"\n",
170186
"plt.figure()\n",
171187
"ax = plt.axes()\n",
172-
"ax.set_ylabel(\"Y\")\n",
173-
"ax.set_xlabel(\"X\")\n",
174-
"ax.set_ylim(1000, 49000)\n",
175-
"ax.set_xlim(1000, 99000)\n",
188+
"ax.set_ylabel(\"Latitude\")\n",
189+
"ax.set_xlabel(\"Longitude\")\n",
176190
"ax.plot(Particle_data.lon.transpose(), Particle_data.lat.transpose(), c=\"k\", zorder=1)\n",
177191
"T_scatter = ax.scatter(\n",
178192
" Particle_data.lon,\n",
179193
" Particle_data.lat,\n",
180194
" c=Particle_data.temperature,\n",
181195
" cmap=plt.cm.inferno,\n",
182-
" norm=mpl.colors.Normalize(vmin=0.0, vmax=20.0),\n",
196+
" norm=mpl.colors.Normalize(vmin=22.0, vmax=24.0),\n",
183197
" edgecolor=\"k\",\n",
184198
" zorder=2,\n",
185199
")\n",
@@ -209,15 +223,19 @@
209223
"metadata": {},
210224
"outputs": [],
211225
"source": [
212-
"def SampleVel_wrong(particle, fieldset, time):\n",
213-
" u = fieldset.U[particle]\n",
226+
"def SampleVel_wrong(particles, fieldset):\n",
227+
" u = fieldset.U[particles]\n",
214228
"\n",
215229
"\n",
216230
"pset = parcels.ParticleSet(\n",
217-
" fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time\n",
231+
" fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, depth=depth\n",
218232
")\n",
219233
"\n",
220-
"pset.execute(SampleVel_wrong)"
234+
"pset.execute(\n",
235+
" SampleVel_wrong,\n",
236+
" runtime=timedelta(hours=30),\n",
237+
" dt=timedelta(minutes=5),\n",
238+
")"
221239
]
222240
},
223241
{
@@ -234,15 +252,15 @@
234252
"metadata": {},
235253
"outputs": [],
236254
"source": [
237-
"def SampleVel_correct(particle, fieldset, time):\n",
238-
" u, v = fieldset.UV[particle]\n",
255+
"def SampleVel_correct(particles, fieldset):\n",
256+
" u, v = fieldset.UV[particles]\n",
239257
"\n",
240258
"\n",
241259
"pset = parcels.ParticleSet(\n",
242-
" fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time\n",
260+
" fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, depth=depth\n",
243261
")\n",
244262
"\n",
245-
"pset.execute(SampleVel_correct)"
263+
"pset.execute(SampleVel_correct, runtime=timedelta(hours=30), dt=timedelta(minutes=5))"
246264
]
247265
},
248266
{
@@ -259,19 +277,17 @@
259277
"metadata": {},
260278
"outputs": [],
261279
"source": [
262-
"SampleParticle = parcels.Particle.add_variables(\n",
280+
"SampleParticle = parcels.Particle.add_variable(\n",
263281
" [\n",
264282
" parcels.Variable(\"U\", dtype=np.float32, initial=np.nan),\n",
265283
" parcels.Variable(\"V\", dtype=np.float32, initial=np.nan),\n",
266284
" ]\n",
267285
")\n",
268286
"\n",
269287
"\n",
270-
"def SampleVel_correct(particle, fieldset, time):\n",
288+
"def SampleVel_correct(particles, fieldset):\n",
271289
" # attention: samples particle velocity in units of the mesh (deg/s or m/s)\n",
272-
" particle.U, particle.V = fieldset.UV[\n",
273-
" time, particle.depth, particle.lat, particle.lon, particle\n",
274-
" ]"
290+
" particles.U, particles.V = fieldset.UV[particles]"
275291
]
276292
},
277293
{
@@ -306,11 +322,16 @@
306322
"outputs": [],
307323
"source": [
308324
"SampleParticleOnce = parcels.Particle.add_variable(\n",
309-
" \"temperature\", initial=0, to_write=\"once\"\n",
325+
" parcels.Variable(\"temperature\", initial=np.nan, to_write=\"once\")\n",
310326
")\n",
311327
"\n",
312328
"pset = parcels.ParticleSet(\n",
313-
" fieldset=fieldset, pclass=SampleParticleOnce, lon=lon, lat=lat, time=time\n",
329+
" fieldset=fieldset,\n",
330+
" pclass=SampleParticleOnce,\n",
331+
" lon=lon,\n",
332+
" lat=lat,\n",
333+
" time=time,\n",
334+
" depth=depth,\n",
314335
")"
315336
]
316337
},
@@ -320,7 +341,7 @@
320341
"metadata": {},
321342
"outputs": [],
322343
"source": [
323-
"output_file = pset.ParticleFile(name=\"WriteOnce.zarr\", outputdt=timedelta(hours=1))\n",
344+
"output_file = parcels.ParticleFile(\"WriteOnce.zarr\", outputdt=timedelta(hours=1))\n",
324345
"\n",
325346
"pset.execute(\n",
326347
" [parcels.kernels.AdvectionRK4, SampleT],\n",
@@ -331,11 +352,10 @@
331352
]
332353
},
333354
{
334-
"attachments": {},
335355
"cell_type": "markdown",
336356
"metadata": {},
337357
"source": [
338-
"Since all the particles are released at the same x-position and the temperature field is invariant in the y-direction, all particles have an initial temperature of 0.4$^\\circ$C\n"
358+
"We can compare the output where only the initial value is written to output, with the original simulation, where the temperature at each outputdt is written:"
339359
]
340360
},
341361
{
@@ -344,33 +364,35 @@
344364
"metadata": {},
345365
"outputs": [],
346366
"source": [
347-
"Particle_data = xr.open_zarr(\"WriteOnce.zarr\")\n",
367+
"Particle_data_once = xr.open_zarr(\"WriteOnce.zarr\", decode_times=False)\n",
348368
"\n",
349369
"plt.figure()\n",
350370
"ax = plt.axes()\n",
351-
"ax.set_ylabel(\"Y\")\n",
352-
"ax.set_xlabel(\"X\")\n",
353-
"ax.set_ylim(1000, 49000)\n",
354-
"ax.set_xlim(1000, 99000)\n",
355-
"ax.plot(Particle_data.lon.transpose(), Particle_data.lat.transpose(), c=\"k\", zorder=1)\n",
356-
"T_scatter = ax.scatter(\n",
357-
" Particle_data.lon,\n",
358-
" Particle_data.lat,\n",
359-
" c=np.tile(Particle_data.temperature, (Particle_data.lon.shape[1], 1)).T,\n",
360-
" cmap=plt.cm.inferno,\n",
361-
" norm=mpl.colors.Normalize(vmin=0.0, vmax=1.0),\n",
362-
" edgecolor=\"k\",\n",
363-
" zorder=2,\n",
371+
"ax.set_ylabel(\"Temperature [$^{\\\\circ}$C]\")\n",
372+
"ax.set_xlabel(\"Observation Number (-)\")\n",
373+
"ax.set_ylim(22.2, 24.5)\n",
374+
"l1 = ax.plot(Particle_data.obs, Particle_data.temperature.T, color=\"red\")\n",
375+
"l2 = ax.plot(\n",
376+
" Particle_data_once.obs,\n",
377+
" np.tile(Particle_data_once.temperature, (Particle_data_once.lon.shape[1], 1)),\n",
378+
" color=\"tab:blue\",\n",
364379
")\n",
365-
"plt.colorbar(T_scatter, label=r\"Initial T [$^{\\circ} C$]\")\n",
380+
"ax.legend([l1[0], l2[0]], [\"Write every outputdt\", \"Write once\"])\n",
366381
"plt.show()"
367382
]
383+
},
384+
{
385+
"cell_type": "code",
386+
"execution_count": null,
387+
"metadata": {},
388+
"outputs": [],
389+
"source": []
368390
}
369391
],
370392
"metadata": {
371393
"celltoolbar": "Raw-celnotatie",
372394
"kernelspec": {
373-
"display_name": "parcels",
395+
"display_name": "test-notebooks",
374396
"language": "python",
375397
"name": "python3"
376398
},
@@ -384,7 +406,7 @@
384406
"name": "python",
385407
"nbconvert_exporter": "python",
386408
"pygments_lexer": "ipython3",
387-
"version": "3.12.3"
409+
"version": "3.12.11"
388410
},
389411
"pycharm": {
390412
"stem_cell": {

docs/examples/documentation_LargeRunsOutput.ipynb renamed to docs/examples_v3/documentation_LargeRunsOutput.ipynb

File renamed without changes.

docs/examples/documentation_MPI.ipynb renamed to docs/examples_v3/documentation_MPI.ipynb

File renamed without changes.

docs/examples/documentation_advanced_zarr.ipynb renamed to docs/examples_v3/documentation_advanced_zarr.ipynb

File renamed without changes.

docs/examples/documentation_geospatial.ipynb renamed to docs/examples_v3/documentation_geospatial.ipynb

File renamed without changes.

docs/examples/documentation_homepage_animation.ipynb renamed to docs/examples_v3/documentation_homepage_animation.ipynb

File renamed without changes.

docs/examples/documentation_indexing.ipynb renamed to docs/examples_v3/documentation_indexing.ipynb

File renamed without changes.

docs/examples/documentation_stuck_particles.ipynb renamed to docs/examples_v3/documentation_stuck_particles.ipynb

File renamed without changes.

0 commit comments

Comments
 (0)