Updating tutorials that combine fields (#2740)

* Updating tutorials that combine fields

The particle release was not very good; too close to the boundary. Also made a few other fixes, such as using RK2 for wind advection and improving wind field

* Update tutorial_manipulating_field_data.ipynb

* Remove time and depth dimension from field

Author: erikvansebille

docs/user_guide/examples/explanation_kernelloop.md

@@ -62,7 +62,7 @@ ds_fields.load()  # load the dataset into memory
 # Create an idealised wind field and add it to the dataset
 tdim, ydim, xdim = (len(ds_fields.time),len(ds_fields.latitude), len(ds_fields.longitude))
 ds_fields["UWind"] = xr.DataArray(
-    data=np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values)[None, :, None],
+    data=0.5 * np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[None, :, None],
     coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude])
 
 ds_fields["VWind"] = xr.DataArray(
@@ -88,10 +88,12 @@ fieldset = parcels.FieldSet.from_sgrid_conventions(
 Now we define a wind Kernel that uses a forward Euler method to apply the wind forcing. Note that we update the `particles.dx` and `particles.dy` variables, rather than `particles.x` and `particles.y` directly.
 
 ```{code-cell}
-def wind_kernel(particles, fieldset):
-    uwind, vwind = fieldset.Wind[particles]
-    particles.dx += uwind * particles.dt
-    particles.dy += vwind * particles.dt
+def windRK2(particles, fieldset):
+    (u1, v1) = fieldset.Wind[particles]
+    x1, y1 = (particles.x + u1 * 0.5 * particles.dt, particles.y + v1 * 0.5 * particles.dt)
+    (u2, v2) = fieldset.Wind[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles]
+    particles.dx += u2 * particles.dt
+    particles.dy += v2 * particles.dt
 ```
 
 First run a simulation where we apply Kernels as `[AdvectionRK2, wind_kernel]`
@@ -100,22 +102,22 @@ First run a simulation where we apply Kernels as `[AdvectionRK2, wind_kernel]`
 :tags: [hide-output]
 npart = 10
 z = np.repeat(ds_fields.depth[0].values, npart)
-lons = np.repeat(31, npart)
+lons = np.repeat(32.2, npart)
 lats = np.linspace(-32.5, -30.5, npart)
 
 pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons)
 output_file = parcels.ParticleFile(
     path="advection_then_wind.parquet", outputdt=np.timedelta64(6,'h')
 )
 pset.execute(
-    [parcels.kernels.AdvectionRK2, wind_kernel],
+    [parcels.kernels.AdvectionRK2, windRK2],
     runtime=np.timedelta64(5,'D'),
     dt=np.timedelta64(1,'h'),
     output_file=output_file,
 )
 ```
 
-Then also run a simulation where we apply the Kernels in the reverse order as `[wind_kernel, AdvectionRK2]`
+Then also run a simulation where we apply the Kernels in the reverse order as `[windRK2, AdvectionRK2]`
 
 ```{code-cell}
 :tags: [hide-output]
@@ -126,7 +128,7 @@ output_file_reverse = parcels.ParticleFile(
     path="wind_then_advection.parquet", outputdt=np.timedelta64(6,"h")
 )
 pset_reverse.execute(
-    [wind_kernel, parcels.kernels.AdvectionRK2],
+    [windRK2, parcels.kernels.AdvectionRK2],
     runtime=np.timedelta64(5,"D"),
     dt=np.timedelta64(1,"h"),
     output_file=output_file_reverse,

docs/user_guide/examples/tutorial_manipulating_field_data.ipynb

@@ -52,19 +52,17 @@
     "ds_fields.load()  # load the dataset into memory\n",
     "\n",
     "# Create an idealised wind field and add it to the dataset\n",
-    "tdim, ydim, xdim = (\n",
-    "    len(ds_fields.time),\n",
-    "    len(ds_fields.latitude),\n",
-    "    len(ds_fields.longitude),\n",
-    ")\n",
+    "ydim, xdim = len(ds_fields.latitude), len(ds_fields.longitude)\n",
     "ds_fields[\"UWind\"] = xr.DataArray(\n",
-    "    data=np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values)[None, :, None],\n",
-    "    coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n",
+    "    data=0.5\n",
+    "    * np.ones((ydim, xdim))\n",
+    "    * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[:, None],\n",
+    "    coords=[ds_fields.latitude, ds_fields.longitude],\n",
     ")\n",
     "\n",
     "ds_fields[\"VWind\"] = xr.DataArray(\n",
-    "    data=np.zeros((tdim, ydim, xdim)),\n",
-    "    coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n",
+    "    data=np.zeros((ydim, xdim)),\n",
+    "    coords=[ds_fields.latitude, ds_fields.longitude],\n",
     ")"
    ]
   },
@@ -126,12 +124,13 @@
     "\n",
     "npart = 10\n",
     "z = np.repeat(ds_fields.depth[0].values, npart)\n",
-    "lons = np.repeat(31, npart)\n",
+    "lons = np.repeat(32.2, npart)\n",
     "lats = np.linspace(-32.5, -30.5, npart)\n",
     "\n",
     "pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons)\n",
     "output_file = parcels.ParticleFile(\n",
-    "    path=\"summed_advection_wind.parquet\", outputdt=np.timedelta64(6, \"h\")\n",
+    "    path=\"summed_advection_wind.parquet\",\n",
+    "    outputdt=np.timedelta64(6, \"h\"),\n",
     ")\n",
     "pset.execute(\n",
     "    [parcels.kernels.AdvectionRK2],\n",
@@ -158,8 +157,9 @@
    "source": [
     "# Plot the resulting particle trajectories overlapped for both cases\n",
     "summed_advection_wind = parcels.read_particlefile(\"summed_advection_wind.parquet\")\n",
+    "fig, ax = plt.subplots(figsize=(5, 3))\n",
     "for traj in summed_advection_wind.partition_by(\"particle_id\", maintain_order=True):\n",
-    "    plt.plot(traj[\"x\"], traj[\"y\"], \"-\")\n",
+    "    ax.plot(traj[\"x\"], traj[\"y\"], \"-\")\n",
     "plt.show()"
    ]
   }