|
15 | 15 | "metadata": {}, |
16 | 16 | "outputs": [], |
17 | 17 | "source": [ |
18 | | - "from parcels import FieldSet, ParticleSet, JITParticle, Variable, AdvectionRK4, ScipyParticle, StatusCode\n", |
| 18 | + "from parcels import (\n", |
| 19 | + " FieldSet,\n", |
| 20 | + " ParticleSet,\n", |
| 21 | + " JITParticle,\n", |
| 22 | + " Variable,\n", |
| 23 | + " AdvectionRK4,\n", |
| 24 | + " ScipyParticle,\n", |
| 25 | + " StatusCode,\n", |
| 26 | + ")\n", |
19 | 27 | "import numpy as np\n", |
20 | 28 | "import math\n", |
21 | 29 | "import xarray as xr\n", |
|
62 | 70 | "outputs": [], |
63 | 71 | "source": [ |
64 | 72 | "# Execute the code below to create what's in Parcels called a `FieldSet` object\n", |
65 | | - "filename = 'global-analysis-forecast-phy-001-024_1641284158025.nc'\n", |
66 | | - "variables = {'U': 'uo', # dictionary with names of the variables in the netcdf file\n", |
67 | | - " 'V': 'vo',\n", |
68 | | - " 'S': 'so',\n", |
69 | | - " 'T': 'thetao'}\n", |
70 | | - "dimensions = {'lon': 'longitude', # dictionary with names of the dimensions in the netcdf file\n", |
71 | | - " 'lat': 'latitude',\n", |
72 | | - " 'depth': 'depth',\n", |
73 | | - " 'time': 'time'}\n", |
| 73 | + "filename = \"global-analysis-forecast-phy-001-024_1641284158025.nc\"\n", |
| 74 | + "variables = {\n", |
| 75 | + " \"U\": \"uo\", # dictionary with names of the variables in the netcdf file\n", |
| 76 | + " \"V\": \"vo\",\n", |
| 77 | + " \"S\": \"so\",\n", |
| 78 | + " \"T\": \"thetao\",\n", |
| 79 | + "}\n", |
| 80 | + "dimensions = {\n", |
| 81 | + " \"lon\": \"longitude\", # dictionary with names of the dimensions in the netcdf file\n", |
| 82 | + " \"lat\": \"latitude\",\n", |
| 83 | + " \"depth\": \"depth\",\n", |
| 84 | + " \"time\": \"time\",\n", |
| 85 | + "}\n", |
74 | 86 | "\n", |
75 | 87 | "fieldset = FieldSet.from_netcdf(filename, variables, dimensions)\n", |
76 | 88 | "fieldset.mindepth = fieldset.U.depth[0] # uppermost layer in the hydrodynamic data" |
|
111 | 123 | } |
112 | 124 | ], |
113 | 125 | "source": [ |
114 | | - "fieldset.computeTimeChunk() # first load the data in memory\n", |
| 126 | + "fieldset.computeTimeChunk() # first load the data in memory\n", |
115 | 127 | "\n", |
116 | 128 | "fig, ax = plt.subplots(1, 2, sharex=True, sharey=True)\n", |
117 | | - "p0 = ax[0].pcolormesh(fieldset.U.lon, fieldset.U.lat, fieldset.U.data[0,0,:,:], vmin=-1.5, vmax=1.5, cmap=cmocean.cm.balance)\n", |
118 | | - "p1 = ax[1].pcolormesh(fieldset.V.lon, fieldset.V.lat, fieldset.V.data[0,0,:,:], vmin=-1.5, vmax=1.5, cmap=cmocean.cm.balance)\n", |
119 | | - "ax[0].set_title('zonal velocity')\n", |
120 | | - "ax[1].set_title('meridional velocity')\n", |
121 | | - "ax[0].set_ylabel('latitude [°N]')\n", |
122 | | - "ax[0].set_xlabel('longitude [°E]')\n", |
123 | | - "ax[1].set_xlabel('longitude [°E]')\n", |
| 129 | + "p0 = ax[0].pcolormesh(\n", |
| 130 | + " fieldset.U.lon,\n", |
| 131 | + " fieldset.U.lat,\n", |
| 132 | + " fieldset.U.data[0, 0, :, :],\n", |
| 133 | + " vmin=-1.5,\n", |
| 134 | + " vmax=1.5,\n", |
| 135 | + " cmap=cmocean.cm.balance,\n", |
| 136 | + ")\n", |
| 137 | + "p1 = ax[1].pcolormesh(\n", |
| 138 | + " fieldset.V.lon,\n", |
| 139 | + " fieldset.V.lat,\n", |
| 140 | + " fieldset.V.data[0, 0, :, :],\n", |
| 141 | + " vmin=-1.5,\n", |
| 142 | + " vmax=1.5,\n", |
| 143 | + " cmap=cmocean.cm.balance,\n", |
| 144 | + ")\n", |
| 145 | + "ax[0].set_title(\"zonal velocity\")\n", |
| 146 | + "ax[1].set_title(\"meridional velocity\")\n", |
| 147 | + "ax[0].set_ylabel(\"latitude [°N]\")\n", |
| 148 | + "ax[0].set_xlabel(\"longitude [°E]\")\n", |
| 149 | + "ax[1].set_xlabel(\"longitude [°E]\")\n", |
124 | 150 | "ax[1].yaxis.set_tick_params(left=False, right=True, labelleft=False, labelright=True)\n", |
125 | 151 | "ax[0].grid()\n", |
126 | 152 | "ax[1].grid()\n", |
127 | | - "cb = fig.colorbar(p1, ax=ax, orientation='horizontal')\n", |
128 | | - "cb.set_label('[m/s]')\n", |
| 153 | + "cb = fig.colorbar(p1, ax=ax, orientation=\"horizontal\")\n", |
| 154 | + "cb.set_label(\"[m/s]\")\n", |
129 | 155 | "gd = fieldset.U.grid\n", |
130 | | - "fig.suptitle(f\"velocity at depth {gd.depth[0]:.2f} m, {gd.timeslices[0][0].astype('datetime64[m]')}\")\n", |
| 156 | + "fig.suptitle(\n", |
| 157 | + " f\"velocity at depth {gd.depth[0]:.2f} m, {gd.timeslices[0][0].astype('datetime64[m]')}\"\n", |
| 158 | + ")\n", |
131 | 159 | "plt.show()\n", |
132 | 160 | "\n", |
133 | | - "print(f\"The dimensions of the data are:\\nlon: {gd.xdim}\\nlat: {gd.ydim}\\ndepth: {gd.zdim}\\ntime: {len(gd.time_full)} ({gd.tdim} loaded)\")" |
| 161 | + "print(\n", |
| 162 | + " f\"The dimensions of the data are:\\nlon: {gd.xdim}\\nlat: {gd.ydim}\\ndepth: {gd.zdim}\\ntime: {len(gd.time_full)} ({gd.tdim} loaded)\"\n", |
| 163 | + ")" |
134 | 164 | ] |
135 | 165 | }, |
136 | 166 | { |
|
183 | 213 | " elif particle.cycle_phase == 3:\n", |
184 | 214 | " # Phase 3: Rising with vertical_speed until at surface\n", |
185 | 215 | " particle_ddepth -= vertical_speed * particle.dt\n", |
186 | | - " particle.cycle_age += particle.dt # solve issue of not updating cycle_age during ascent\n", |
| 216 | + " particle.cycle_age += (\n", |
| 217 | + " particle.dt\n", |
| 218 | + " ) # solve issue of not updating cycle_age during ascent\n", |
187 | 219 | " if particle.depth <= fieldset.mindepth:\n", |
188 | 220 | " particle.depth = fieldset.mindepth\n", |
189 | | - " particle.temp = math.nan # reset temperature to NaN at end of sampling cycle\n", |
| 221 | + " particle.temp = (\n", |
| 222 | + " math.nan\n", |
| 223 | + " ) # reset temperature to NaN at end of sampling cycle\n", |
190 | 224 | " particle.salt = math.nan # idem\n", |
191 | 225 | " particle.cycle_phase = 4\n", |
192 | 226 | " else:\n", |
|
226 | 260 | "source": [ |
227 | 261 | "# Define the ArgoParticle class\n", |
228 | 262 | "variables = [\n", |
229 | | - " Variable('cycle_phase', dtype=np.int32, initial=0.0),\n", |
230 | | - " Variable('cycle_age', dtype=np.float32, initial=0.0),\n", |
231 | | - " Variable('drift_age', dtype=np.float32, initial=0.0),\n", |
232 | | - " Variable('temp', dtype=np.float32, initial=np.nan),\n", |
233 | | - " Variable('salt', dtype=np.float32, initial=np.nan),\n", |
| 263 | + " Variable(\"cycle_phase\", dtype=np.int32, initial=0.0),\n", |
| 264 | + " Variable(\"cycle_age\", dtype=np.float32, initial=0.0),\n", |
| 265 | + " Variable(\"drift_age\", dtype=np.float32, initial=0.0),\n", |
| 266 | + " Variable(\"temp\", dtype=np.float32, initial=np.nan),\n", |
| 267 | + " Variable(\"salt\", dtype=np.float32, initial=np.nan),\n", |
234 | 268 | "]\n", |
235 | 269 | "ArgoParticle = JITParticle.add_variables(variables)" |
236 | 270 | ] |
|
248 | 282 | "metadata": {}, |
249 | 283 | "outputs": [], |
250 | 284 | "source": [ |
251 | | - "# Define the ParticleSet \n", |
| 285 | + "# Define the ParticleSet\n", |
252 | 286 | "argoset = ParticleSet(\n", |
253 | | - " fieldset=fieldset, pclass=ArgoParticle, lon=[-37.09499], lat=[6.73472],\n", |
254 | | - " depth=fieldset.mindepth, time=np.datetime64('2019-07-25T11:14:00')\n", |
| 287 | + " fieldset=fieldset,\n", |
| 288 | + " pclass=ArgoParticle,\n", |
| 289 | + " lon=[-37.09499],\n", |
| 290 | + " lat=[6.73472],\n", |
| 291 | + " depth=fieldset.mindepth,\n", |
| 292 | + " time=np.datetime64(\"2019-07-25T11:14:00\"),\n", |
255 | 293 | ")" |
256 | 294 | ] |
257 | 295 | }, |
|
271 | 309 | "outputs": [], |
272 | 310 | "source": [ |
273 | 311 | "argoset.execute(\n", |
274 | | - " [ArgoVerticalMovement, AdvectionRK4, KeepAtSurface], # list of kernels to be executed\n", |
| 312 | + " [\n", |
| 313 | + " ArgoVerticalMovement,\n", |
| 314 | + " AdvectionRK4,\n", |
| 315 | + " KeepAtSurface,\n", |
| 316 | + " ], # list of kernels to be executed\n", |
275 | 317 | " runtime=timedelta(days=31),\n", |
276 | 318 | " dt=timedelta(minutes=5),\n", |
277 | | - " output_file=argoset.ParticleFile(name=\"argo_float\", outputdt=timedelta(minutes=5), chunks=(1,500)),\n", |
| 319 | + " output_file=argoset.ParticleFile(\n", |
| 320 | + " name=\"argo_float\", outputdt=timedelta(minutes=5), chunks=(1, 500)\n", |
| 321 | + " ),\n", |
278 | 322 | ")" |
279 | 323 | ] |
280 | 324 | }, |
|
303 | 347 | ], |
304 | 348 | "source": [ |
305 | 349 | "# Plot the trajectory of the Argo float you've just simulated\n", |
306 | | - "ds_argo = xr.open_zarr('argo_float.zarr')\n", |
| 350 | + "ds_argo = xr.open_zarr(\"argo_float.zarr\")\n", |
307 | 351 | "\n", |
308 | 352 | "x = ds_argo[\"lon\"][:].squeeze()\n", |
309 | 353 | "y = ds_argo[\"lat\"][:].squeeze()\n", |
310 | 354 | "z = ds_argo[\"z\"][:].squeeze()\n", |
311 | 355 | "\n", |
312 | | - "fig = plt.figure(figsize=(13,10))\n", |
313 | | - "ax = plt.axes(projection='3d')\n", |
| 356 | + "fig = plt.figure(figsize=(13, 10))\n", |
| 357 | + "ax = plt.axes(projection=\"3d\")\n", |
314 | 358 | "ax.invert_zaxis()\n", |
315 | | - "p = ax.scatter(x,y,z, c=z, s=20, marker=\"o\", cmap='viridis_r', depthshade=False)\n", |
| 359 | + "p = ax.scatter(x, y, z, c=z, s=20, marker=\"o\", cmap=\"viridis_r\", depthshade=False)\n", |
316 | 360 | "cb = plt.colorbar(p, shrink=0.5)\n", |
317 | 361 | "cb.ax.invert_yaxis()\n", |
318 | | - "cb.set_label('depth [m]')" |
| 362 | + "cb.set_label(\"depth [m]\")" |
319 | 363 | ] |
320 | 364 | }, |
321 | 365 | { |
|
352 | 396 | "source": [ |
353 | 397 | "# Answer\n", |
354 | 398 | "\n", |
355 | | - "def ArgoVerticalMovement2(particle, fieldset, time): # -- change kernel name\n", |
| 399 | + "\n", |
| 400 | + "def ArgoVerticalMovement2(particle, fieldset, time): # -- change kernel name\n", |
356 | 401 | " driftdepth = 1000 # maximum depth in m\n", |
357 | 402 | " maxdepth = 2000 # maximum depth in m\n", |
358 | 403 | " vertical_speed = 0.10 # sink and rise speed in m/s\n", |
|
380 | 425 | "\n", |
381 | 426 | " elif particle.cycle_phase == 3:\n", |
382 | 427 | " # Phase 3: Rising with vertical_speed until at surface\n", |
383 | | - " particle_ddepth -= vertical_speed * particle.dt/ 2.0 # ------------- make twice as slow\n", |
384 | | - " particle.cycle_age += particle.dt \n", |
| 428 | + " particle_ddepth -= (\n", |
| 429 | + " vertical_speed * particle.dt / 2.0\n", |
| 430 | + " ) # ------------- make twice as slow\n", |
| 431 | + " particle.cycle_age += particle.dt\n", |
385 | 432 | " particle.temp = fieldset.T[time, particle.depth, particle.lat, particle.lon]\n", |
386 | 433 | " particle.salt = fieldset.S[time, particle.depth, particle.lat, particle.lon]\n", |
387 | 434 | " if particle.depth <= fieldset.mindepth:\n", |
388 | 435 | " particle.depth = fieldset.mindepth\n", |
389 | | - " particle.temp = math.nan # reset temperature to NaN at end of sampling cycle\n", |
| 436 | + " particle.temp = (\n", |
| 437 | + " math.nan\n", |
| 438 | + " ) # reset temperature to NaN at end of sampling cycle\n", |
390 | 439 | " particle.salt = math.nan # idem\n", |
391 | 440 | " particle.cycle_phase = 4\n", |
392 | | - " \n", |
| 441 | + "\n", |
393 | 442 | " elif particle.cycle_phase == 4:\n", |
394 | 443 | " # Phase 4: Transmitting at surface until cycletime is reached\n", |
395 | 444 | " if particle.cycle_age > cycletime:\n", |
|
407 | 456 | "\n", |
408 | 457 | "\n", |
409 | 458 | "argoset = ParticleSet( # ------------ redefine argoset, otherwise simulation continues from previous run\n", |
410 | | - " fieldset=fieldset, pclass=ArgoParticle, lon=[-37.09499], lat=[6.73472],\n", |
411 | | - " depth=fieldset.mindepth, time=np.datetime64('2019-07-25T11:14:00')\n", |
| 459 | + " fieldset=fieldset,\n", |
| 460 | + " pclass=ArgoParticle,\n", |
| 461 | + " lon=[-37.09499],\n", |
| 462 | + " lat=[6.73472],\n", |
| 463 | + " depth=fieldset.mindepth,\n", |
| 464 | + " time=np.datetime64(\"2019-07-25T11:14:00\"),\n", |
412 | 465 | ")\n", |
413 | 466 | "argoset.execute(\n", |
414 | 467 | " [ArgoVerticalMovement2, AdvectionRK4, KeepAtSurface], # ----- changed kernel\n", |
415 | 468 | " runtime=timedelta(days=31),\n", |
416 | 469 | " dt=timedelta(minutes=5),\n", |
417 | | - " output_file=argoset.ParticleFile(name=\"argo_float2\", outputdt=timedelta(minutes=30), chunks=(1,500)), # -- changed filename\n", |
| 470 | + " output_file=argoset.ParticleFile(\n", |
| 471 | + " name=\"argo_float2\", outputdt=timedelta(minutes=30), chunks=(1, 500)\n", |
| 472 | + " ), # -- changed filename\n", |
418 | 473 | ")\n", |
419 | 474 | "\n", |
420 | 475 | "# Plot depth as a function of time for both the regular and slowly-ascending floats to check if your code works\n", |
421 | | - "ds_argo2 = xr.open_zarr('argo_float2.zarr')\n", |
| 476 | + "ds_argo2 = xr.open_zarr(\"argo_float2.zarr\")\n", |
422 | 477 | "fig = plt.figure()\n", |
423 | | - "plt.plot(ds_argo.time[0,:], -ds_argo.z[0,:], 'k', label='original')\n", |
424 | | - "plt.plot(ds_argo2.time[0,:], -ds_argo2.z[0,:], 'b:',label='edited')\n", |
| 478 | + "plt.plot(ds_argo.time[0, :], -ds_argo.z[0, :], \"k\", label=\"original\")\n", |
| 479 | + "plt.plot(ds_argo2.time[0, :], -ds_argo2.z[0, :], \"b:\", label=\"edited\")\n", |
425 | 480 | "plt.legend()\n", |
426 | 481 | "plt.grid()\n", |
427 | | - "plt.ylabel('depth [m]')\n", |
| 482 | + "plt.ylabel(\"depth [m]\")\n", |
428 | 483 | "plt.show()" |
429 | 484 | ] |
430 | 485 | }, |
|
446 | 501 | ], |
447 | 502 | "source": [ |
448 | 503 | "# Plots of the vertical temperature and salinity profile, which were recorded during the ascents to the surface.\n", |
449 | | - "fig = plt.figure(figsize=(12,6))\n", |
450 | | - "fig.suptitle(f'Temperature, salinity and time during phase 3')\n", |
| 504 | + "fig = plt.figure(figsize=(12, 6))\n", |
| 505 | + "fig.suptitle(f\"Temperature, salinity and time during phase 3\")\n", |
451 | 506 | "\n", |
452 | | - "days = (ds_argo.time[0,:]-ds_argo.time[0,0]).data.astype('timedelta64[m]').astype('int')/1440\n", |
453 | | - "tids3 = ~ds_argo.temp[0,:].isnull().compute() # time indices during cycle phase 3 (ascent)\n", |
| 507 | + "days = (ds_argo.time[0, :] - ds_argo.time[0, 0]).data.astype(\"timedelta64[m]\").astype(\n", |
| 508 | + " \"int\"\n", |
| 509 | + ") / 1440\n", |
| 510 | + "tids3 = (\n", |
| 511 | + " ~ds_argo.temp[0, :].isnull().compute()\n", |
| 512 | + ") # time indices during cycle phase 3 (ascent)\n", |
454 | 513 | "\n", |
455 | 514 | "ax1 = fig.add_subplot(131)\n", |
456 | | - "ax1.plot(ds_argo.temp[0,:],-ds_argo.z[0,:], 'k')\n", |
| 515 | + "ax1.plot(ds_argo.temp[0, :], -ds_argo.z[0, :], \"k\")\n", |
457 | 516 | "ax1.grid()\n", |
458 | | - "ax1.set_xlabel('temperature [°C]')\n", |
459 | | - "ax1.set_ylabel('depth [m]')\n", |
| 517 | + "ax1.set_xlabel(\"temperature [°C]\")\n", |
| 518 | + "ax1.set_ylabel(\"depth [m]\")\n", |
460 | 519 | "\n", |
461 | 520 | "ax2 = fig.add_subplot(132, sharey=ax1)\n", |
462 | | - "ax2.plot(ds_argo.salt[0,:], -ds_argo.z[0,:], 'k')\n", |
| 521 | + "ax2.plot(ds_argo.salt[0, :], -ds_argo.z[0, :], \"k\")\n", |
463 | 522 | "ax2.grid()\n", |
464 | 523 | "ax2.yaxis.set_tick_params(left=False, labelleft=False)\n", |
465 | | - "ax2.set_xlabel('salinity [g/kg]')\n", |
| 524 | + "ax2.set_xlabel(\"salinity [g/kg]\")\n", |
466 | 525 | "\n", |
467 | 526 | "ax3 = fig.add_subplot(133, sharey=ax1)\n", |
468 | | - "ax3.plot(days[tids3], -ds_argo.z[0,tids3], 'k-o', ms=4)\n", |
| 527 | + "ax3.plot(days[tids3], -ds_argo.z[0, tids3], \"k-o\", ms=4)\n", |
469 | 528 | "ax3.yaxis.tick_right()\n", |
470 | | - "ax3.yaxis.label_position = 'right'\n", |
| 529 | + "ax3.yaxis.label_position = \"right\"\n", |
471 | 530 | "ax3.grid()\n", |
472 | | - "ax3.set_xlabel('days from launch')\n", |
473 | | - "ax3.set_ylabel('depth [m]', rotation=270)\n", |
| 531 | + "ax3.set_xlabel(\"days from launch\")\n", |
| 532 | + "ax3.set_ylabel(\"depth [m]\", rotation=270)\n", |
474 | 533 | "plt.show()" |
475 | 534 | ] |
476 | 535 | } |
|
0 commit comments