Skip to content

Commit 85860a5

Browse files
erikvansebilleerikvansebille
committed
FInalising unstuck tutorial to v4
1 parent dee92eb commit 85860a5

1 file changed

Lines changed: 92 additions & 133 deletions

File tree

‎docs/user_guide/examples/documentation_unstuck_Agrid.ipynb‎

Lines changed: 92 additions & 133 deletions
Original file line numberDiff line numberDiff line change
@@ -42,15 +42,12 @@
4242
"metadata": {},
4343
"outputs": [],
4444
"source": [
45-
"from copy import copy\n",
4645
"from datetime import timedelta\n",
4746
"\n",
4847
"import matplotlib.gridspec as gridspec\n",
4948
"import matplotlib.pyplot as plt\n",
5049
"import numpy as np\n",
5150
"import xarray as xr\n",
52-
"from matplotlib.lines import Line2D\n",
53-
"from netCDF4 import Dataset\n",
5451
"from scipy import interpolate\n",
5552
"\n",
5653
"import parcels\n",
@@ -98,7 +95,12 @@
9895
"\n",
9996
"fields = {\"U\": ds[\"uo\"], \"V\": ds[\"vo\"]}\n",
10097
"ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)\n",
101-
"fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)"
98+
"ds_fset = ds_fset.fillna(\n",
99+
" 0\n",
100+
") # TODO remove when fillna done in convert.copernicusmarine_to_sgrid\n",
101+
"\n",
102+
"fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)\n",
103+
"fieldset = fieldset.to_windowed_arrays()"
102104
]
103105
},
104106
{
@@ -655,10 +657,6 @@
655657
"metadata": {},
656658
"outputs": [],
657659
"source": [
658-
"fieldset.models[0].data = fieldset.models[0].data.fillna(\n",
659-
" 0\n",
660-
") # TODO remove when fillna done in convert.copernicusmarine_to_sgrid\n",
661-
"\n",
662660
"pset = parcels.ParticleSet(fieldset=fieldset, x=X, y=Y)\n",
663661
"\n",
664662
"kernels = parcels.kernels.AdvectionRK2\n",
@@ -781,10 +779,28 @@
781779
")\n",
782780
"for j, traj in enumerate(df_stuck.partition_by(\"particle_id\", maintain_order=True)):\n",
783781
" label = \"stuck\" if j == 0 else None\n",
784-
" ax.plot(traj[\"x\"], traj[\"y\"], \".-\", linewidth=2, zorder=1, color=\"y\", label=label)\n",
782+
" ax.plot(\n",
783+
" traj[\"x\"],\n",
784+
" traj[\"y\"],\n",
785+
" \"-\",\n",
786+
" linewidth=2,\n",
787+
" zorder=3,\n",
788+
" color=\"y\",\n",
789+
" alpha=0.7,\n",
790+
" label=label,\n",
791+
" )\n",
785792
"for j, traj in enumerate(df_displace.partition_by(\"particle_id\", maintain_order=True)):\n",
786793
" label = \"displaced\" if j == 0 else None\n",
787-
" ax.plot(traj[\"x\"], traj[\"y\"], \".-\", linewidth=2, zorder=1, color=\"c\", label=label)\n",
794+
" ax.plot(\n",
795+
" traj[\"x\"],\n",
796+
" traj[\"y\"],\n",
797+
" \"-\",\n",
798+
" linewidth=2,\n",
799+
" zorder=2,\n",
800+
" color=\"c\",\n",
801+
" alpha=0.7,\n",
802+
" label=label,\n",
803+
" )\n",
788804
"ax.plot(\n",
789805
" X.flatten(), Y.flatten(), \"o\", color=\"w\", markersize=5, label=\"initial positions\"\n",
790806
")\n",
@@ -964,8 +980,6 @@
964980
"source": [
965981
"fieldset.UV.interp_method = parcels.interpolators.XPartialslip()\n",
966982
"\n",
967-
"fieldset.to_windowed_arrays() # TODO make partialslip also work in dask\n",
968-
"\n",
969983
"pset = parcels.ParticleSet(fieldset=fieldset, x=X, y=Y)\n",
970984
"\n",
971985
"output_file = parcels.ParticleFile(\n",
@@ -996,33 +1010,30 @@
9961010
"metadata": {},
9971011
"outputs": [],
9981012
"source": [
999-
"fieldset = parcels.FieldSet.from_netcdf(\n",
1000-
" filenames,\n",
1001-
" variables,\n",
1002-
" dimensions,\n",
1003-
" indices=indices, # TODO v4: Remove `indices` argument from this cell\n",
1004-
" interp_method={\n",
1005-
" \"U\": \"freeslip\",\n",
1006-
" \"V\": \"freeslip\",\n",
1007-
" }, # Setting the interpolation for U and V\n",
1008-
")\n",
1009-
"pset = parcels.ParticleSet(\n",
1010-
" fieldset=fieldset, pclass=parcels.Particle, lon=lons, lat=lats, time=time\n",
1011-
")\n",
1013+
"fieldset.UV.interp_method = parcels.interpolators.XFreeslip()\n",
10121014
"\n",
1013-
"kernels = pset.Kernel(parcels.kernels.AdvectionRK4)\n",
1015+
"pset = parcels.ParticleSet(fieldset=fieldset, x=X, y=Y)\n",
10141016
"\n",
1015-
"output_file = pset.ParticleFile(name=\"SMOC_freeslip.zarr\", outputdt=timedelta(hours=1))\n",
1017+
"output_file = parcels.ParticleFile(\n",
1018+
" \"freeslip_particles.parquet\",\n",
1019+
" outputdt=timedelta(hours=1),\n",
1020+
" mode=\"w\",\n",
1021+
")\n",
10161022
"\n",
1017-
"pset.execute(kernels, runtime=runtime, dt=dt, output_file=output_file)"
1023+
"pset.execute(\n",
1024+
" parcels.kernels.AdvectionRK2,\n",
1025+
" runtime=runtime,\n",
1026+
" dt=dt,\n",
1027+
" output_file=output_file,\n",
1028+
")"
10181029
]
10191030
},
10201031
{
10211032
"attachments": {},
10221033
"cell_type": "markdown",
10231034
"metadata": {},
10241035
"source": [
1025-
"Now we can load and plot the three different `interpolation_methods`\n"
1036+
"Now we can load and plot the two `interpolation_methods`\n"
10261037
]
10271038
},
10281039
{
@@ -1031,124 +1042,72 @@
10311042
"metadata": {},
10321043
"outputs": [],
10331044
"source": [
1034-
"ds_SMOC = xr.open_zarr(\"SMOC.zarr\")\n",
1035-
"ds_SMOC_part = xr.open_zarr(\"SMOC_partialslip.zarr\")\n",
1036-
"ds_SMOC_free = xr.open_zarr(\"SMOC_freeslip.zarr\")"
1037-
]
1038-
},
1039-
{
1040-
"cell_type": "code",
1041-
"execution_count": null,
1042-
"metadata": {},
1043-
"outputs": [],
1044-
"source": [
1045-
"fig = plt.figure(figsize=(18, 5), constrained_layout=True)\n",
1046-
"fig.suptitle(\"Figure 8. Solution comparison\", fontsize=18, y=1.06)\n",
1047-
"gs = gridspec.GridSpec(ncols=3, nrows=1, figure=fig)\n",
1045+
"df_freeslip = parcels.read_particlefile(\"freeslip_particles.parquet\")\n",
1046+
"df_partialslip = parcels.read_particlefile(\"partialslip_particles.parquet\")\n",
10481047
"\n",
1049-
"n_p = [[0, 1, 3, 4, 6, 7, 8], 0, 6]\n",
1050-
"\n",
1051-
"for i in range(3):\n",
1052-
" ax = fig.add_subplot(gs[0, i])\n",
1053-
" ax.set_title(chr(i + 65) + \") Trajectory \" + str(n_p[i]), fontsize=18)\n",
1054-
" land = ax.pcolormesh(\n",
1055-
" lon_vals[latmin : latmax + 1, lonmin : lonmax + 1],\n",
1056-
" lat_vals[latmin : latmax + 1, lonmin : lonmax + 1],\n",
1057-
" lmask.mask[latmin:latmax, lonmin:lonmax],\n",
1058-
" cmap=\"Reds_r\",\n",
1059-
" )\n",
1060-
" ax.scatter(\n",
1061-
" lons_plot,\n",
1062-
" lats_plot,\n",
1063-
" c=landmask[latmin:latmax, lonmin:lonmax],\n",
1064-
" s=50,\n",
1065-
" cmap=\"Reds_r\",\n",
1066-
" vmin=-0.05,\n",
1067-
" vmax=0.05,\n",
1068-
" edgecolors=\"k\",\n",
1069-
" )\n",
1048+
"fig, ax = plt.subplots(1, 1, figsize=(6, 5), constrained_layout=True)\n",
10701049
"\n",
1071-
" ax.scatter(\n",
1072-
" ds_SMOC[\"lon\"][n_p[i]],\n",
1073-
" ds_SMOC[\"lat\"][n_p[i]],\n",
1074-
" s=30,\n",
1075-
" color=\"limegreen\",\n",
1076-
" zorder=2,\n",
1077-
" )\n",
1050+
"fig.suptitle(\"Figure 7. Solution comparison\", fontsize=16)\n",
1051+
"ax.set_ylabel(\"Latitude [degrees]\")\n",
1052+
"ax.set_xlabel(\"Longitude [degrees]\")\n",
10781053
"\n",
1079-
" ax.scatter(\n",
1080-
" ds_SMOC_disp[\"lon\"][n_p[i]],\n",
1081-
" ds_SMOC_disp[\"lat\"][n_p[i]],\n",
1082-
" s=25,\n",
1083-
" color=\"tab:blue\",\n",
1084-
" zorder=2,\n",
1054+
"ax.pcolormesh(\n",
1055+
" x_centers,\n",
1056+
" y_centers,\n",
1057+
" lmask.mask,\n",
1058+
" cmap=\"Blues_r\",\n",
1059+
")\n",
1060+
"for j, traj in enumerate(df_stuck.partition_by(\"particle_id\", maintain_order=True)):\n",
1061+
" label = \"stuck\" if j == 0 else None\n",
1062+
" ax.plot(\n",
1063+
" traj[\"x\"],\n",
1064+
" traj[\"y\"],\n",
1065+
" \"-\",\n",
1066+
" linewidth=2,\n",
1067+
" zorder=3,\n",
1068+
" color=\"y\",\n",
1069+
" alpha=0.7,\n",
1070+
" label=label,\n",
10851071
" )\n",
1086-
"\n",
1087-
" ax.scatter(\n",
1088-
" ds_SMOC_part[\"lon\"][n_p[i]],\n",
1089-
" ds_SMOC_part[\"lat\"][n_p[i]],\n",
1090-
" s=20,\n",
1091-
" color=\"magenta\",\n",
1072+
"for j, traj in enumerate(\n",
1073+
" df_partialslip.partition_by(\"particle_id\", maintain_order=True)\n",
1074+
"):\n",
1075+
" label = \"partialslip\" if j == 0 else None\n",
1076+
" ax.plot(\n",
1077+
" traj[\"x\"],\n",
1078+
" traj[\"y\"],\n",
1079+
" \"-\",\n",
1080+
" linewidth=2,\n",
10921081
" zorder=2,\n",
1082+
" color=\"r\",\n",
1083+
" alpha=0.7,\n",
1084+
" label=label,\n",
10931085
" )\n",
1094-
"\n",
1095-
" ax.scatter(\n",
1096-
" ds_SMOC_free[\"lon\"][n_p[i]],\n",
1097-
" ds_SMOC_free[\"lat\"][n_p[i]],\n",
1098-
" s=15,\n",
1099-
" color=\"gold\",\n",
1100-
" zorder=2,\n",
1086+
"for j, traj in enumerate(df_freeslip.partition_by(\"particle_id\", maintain_order=True)):\n",
1087+
" label = \"freelslip\" if j == 0 else None\n",
1088+
" ax.plot(\n",
1089+
" traj[\"x\"],\n",
1090+
" traj[\"y\"],\n",
1091+
" \"-\",\n",
1092+
" linewidth=2,\n",
1093+
" zorder=1,\n",
1094+
" color=\"m\",\n",
1095+
" alpha=0.7,\n",
1096+
" label=label,\n",
11011097
" )\n",
1102-
"\n",
1103-
" ax.set_xlim(6.9, 7.6)\n",
1104-
" ax.set_ylim(53.4, 53.9)\n",
1105-
" ax.set_ylabel(\"Latitude [degrees]\")\n",
1106-
" ax.set_xlabel(\"Longitude [degrees]\")\n",
1107-
"\n",
1108-
" color_land = copy(plt.get_cmap(\"Reds\"))(0)\n",
1109-
" color_ocean = copy(plt.get_cmap(\"Reds\"))(128)\n",
1110-
"\n",
1111-
" custom_lines = [\n",
1112-
" Line2D([0], [0], c=\"limegreen\", marker=\"o\", markersize=10, lw=0),\n",
1113-
" Line2D([0], [0], c=\"tab:blue\", marker=\"o\", markersize=10, lw=0),\n",
1114-
" Line2D([0], [0], c=\"magenta\", marker=\"o\", markersize=10, lw=0),\n",
1115-
" Line2D([0], [0], c=\"gold\", marker=\"o\", markersize=10, lw=0),\n",
1116-
" Line2D(\n",
1117-
" [0],\n",
1118-
" [0],\n",
1119-
" c=color_ocean,\n",
1120-
" marker=\"o\",\n",
1121-
" markersize=10,\n",
1122-
" markeredgecolor=\"k\",\n",
1123-
" lw=0,\n",
1124-
" ),\n",
1125-
" Line2D(\n",
1126-
" [0], [0], c=color_land, marker=\"o\", markersize=10, markeredgecolor=\"k\", lw=0\n",
1127-
" ),\n",
1128-
" ]\n",
1129-
" ax.legend(\n",
1130-
" custom_lines,\n",
1131-
" [\n",
1132-
" \"basic RK4\",\n",
1133-
" \"displacement\",\n",
1134-
" \"partial slip\",\n",
1135-
" \"free slip\",\n",
1136-
" \"ocean point\",\n",
1137-
" \"land point\",\n",
1138-
" ],\n",
1139-
" bbox_to_anchor=(0.01, 0.8),\n",
1140-
" loc=\"center left\",\n",
1141-
" borderaxespad=0.0,\n",
1142-
" framealpha=1,\n",
1143-
" )"
1098+
"ax.plot(\n",
1099+
" X.flatten(), Y.flatten(), \"o\", color=\"w\", markersize=5, label=\"initial positions\"\n",
1100+
")\n",
1101+
"plt.legend()\n",
1102+
"plt.show()"
11441103
]
11451104
},
11461105
{
11471106
"attachments": {},
11481107
"cell_type": "markdown",
11491108
"metadata": {},
11501109
"source": [
1151-
"**Figure 8** shows the influence of the different solutions on the particle trajectories near the shore. **Subfigure 8B** shows how the different solutions make trajectory 0 move along the shore and around the corner of the model geometry. **Subfigure 8C** shows how trajectories are unaffected by the different interpolation scheme as long as they do not cross a coastal gridcell.\n"
1110+
"**Figure 7** shows the influence of the different solutions on the particle trajectories near the shore.\n"
11521111
]
11531112
}
11541113
],

0 commit comments

Comments
 (0)