|
5 | 5 | "cell_type": "markdown", |
6 | 6 | "metadata": {}, |
7 | 7 | "source": [ |
8 | | - "# The homepage animation" |
| 8 | + "# 🎓 The homepage animation" |
| 9 | + ] |
| 10 | + }, |
| 11 | + { |
| 12 | + "cell_type": "markdown", |
| 13 | + "metadata": {}, |
| 14 | + "source": [ |
| 15 | + "This notebook shows how to create the animation on the homepage of the Parcels documentation. It uses a dataset of virtual particles at the surface of the global ocean simulation that can be retrieved from the [Copernicus Marine Service](https://marine.copernicus.eu/)." |
9 | 16 | ] |
10 | 17 | }, |
11 | 18 | { |
12 | 19 | "cell_type": "code", |
13 | 20 | "execution_count": null, |
14 | | - "metadata": { |
15 | | - "tags": [ |
16 | | - "raises-exception" |
17 | | - ] |
18 | | - }, |
| 21 | + "metadata": {}, |
19 | 22 | "outputs": [], |
20 | 23 | "source": [ |
21 | | - "%matplotlib qt\n", |
22 | | - "import copy\n", |
23 | | - "\n", |
24 | 24 | "import cartopy\n", |
25 | 25 | "import cartopy.crs as ccrs\n", |
26 | | - "import matplotlib.animation as animation\n", |
| 26 | + "import copernicusmarine\n", |
27 | 27 | "import matplotlib.gridspec as gridspec\n", |
28 | | - "import matplotlib.image as mpimg\n", |
29 | 28 | "import matplotlib.pyplot as plt\n", |
30 | 29 | "import numpy as np\n", |
| 30 | + "import polars as pl\n", |
31 | 31 | "import xarray as xr\n", |
32 | | - "from matplotlib import colors\n", |
33 | | - "from matplotlib.animation import FuncAnimation, PillowWriter, writers" |
| 32 | + "from matplotlib.animation import FuncAnimation, PillowWriter\n", |
| 33 | + "\n", |
| 34 | + "import parcels" |
34 | 35 | ] |
35 | 36 | }, |
36 | 37 | { |
37 | | - "attachments": {}, |
38 | 38 | "cell_type": "markdown", |
39 | 39 | "metadata": {}, |
40 | 40 | "source": [ |
41 | | - "Load in particle data, define the time range at which to plot (**`plottimes`**) and select the indices of the first timestep in the variable `'b'`.\n" |
| 41 | + "The cell below provides the code needed to run this simulation - but because it is a time-consuming run (~20 minutes) and requires a login on the Copernicus Marine Service, we also provide the output file for download directly at TODO" |
42 | 42 | ] |
43 | 43 | }, |
44 | 44 | { |
45 | 45 | "cell_type": "code", |
46 | 46 | "execution_count": null, |
47 | | - "metadata": { |
48 | | - "tags": [ |
49 | | - "raises-exception" |
50 | | - ] |
51 | | - }, |
| 47 | + "metadata": {}, |
52 | 48 | "outputs": [], |
53 | 49 | "source": [ |
54 | | - "filename = \"medusarun.nc\"\n", |
55 | | - "pfile = xr.open_dataset(str(filename), decode_cf=True)\n", |
56 | | - "lon = np.ma.filled(pfile.variables[\"lon\"], np.nan)\n", |
57 | | - "lat = np.ma.filled(pfile.variables[\"lat\"], np.nan)\n", |
58 | | - "time = np.ma.filled(pfile.variables[\"time\"], np.nan)\n", |
| 50 | + "particle_filename = \"copernicusmarine_globalsurface.parquet\"\n", |
| 51 | + "\n", |
| 52 | + "\n", |
| 53 | + "def run_global_copernicusmarine():\n", |
| 54 | + " copernicusmarine.login()\n", |
59 | 55 | "\n", |
60 | | - "pfile.close()\n", |
| 56 | + " ds = copernicusmarine.open_dataset(\n", |
| 57 | + " dataset_id=\"cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m\",\n", |
| 58 | + " variables=[\"uo\", \"vo\"],\n", |
| 59 | + " start_datetime=\"2024-01-01\",\n", |
| 60 | + " end_datetime=\"2024-12-31\",\n", |
| 61 | + " minimum_depth=0.5,\n", |
| 62 | + " maximum_depth=0.5,\n", |
| 63 | + " service=\"arco-geo-series\",\n", |
| 64 | + " chunk_size_limit=1,\n", |
| 65 | + " )\n", |
| 66 | + " ds = ds.fillna(0)\n", |
61 | 67 | "\n", |
62 | | - "plottimes = np.arange(time[0, 0], np.nanmax(time), np.timedelta64(10, \"D\"))\n", |
63 | | - "starttime = 20\n", |
64 | | - "b = time == plottimes[0 + starttime]" |
| 68 | + " ds_wrap = xr.concat(\n", |
| 69 | + " [\n", |
| 70 | + " ds,\n", |
| 71 | + " ds.isel(longitude=slice(0, 1)).assign_coords(\n", |
| 72 | + " longitude=ds.longitude.isel(longitude=slice(0, 1)) + 360\n", |
| 73 | + " ),\n", |
| 74 | + " ],\n", |
| 75 | + " dim=\"longitude\",\n", |
| 76 | + " )\n", |
| 77 | + " ds = parcels.convert.copernicusmarine_to_sgrid(\n", |
| 78 | + " fields={\"U\": ds_wrap[\"uo\"], \"V\": ds_wrap[\"vo\"]}\n", |
| 79 | + " )\n", |
| 80 | + " fieldset = parcels.FieldSet.from_sgrid_conventions(ds)\n", |
| 81 | + " # fieldset.UV.interp_method = parcels.interpolators.XFreeslip()\n", |
| 82 | + " fieldset.to_windowed_arrays()\n", |
| 83 | + "\n", |
| 84 | + " lon, lat = np.meshgrid(np.arange(-179, 180, 2), np.arange(-89, 90, 2))\n", |
| 85 | + " lon = lon.flatten()\n", |
| 86 | + " lat = lat.flatten()\n", |
| 87 | + "\n", |
| 88 | + " # Filter out particles that are not in the ocean (i.e. where speed is zero)\n", |
| 89 | + " u, v = fieldset.UV[np.zeros_like(lon), np.zeros_like(lat), lat, lon]\n", |
| 90 | + " speed = np.sqrt(u**2 + v**2)\n", |
| 91 | + " inocean = speed > 0\n", |
| 92 | + "\n", |
| 93 | + " pset = parcels.ParticleSet(fieldset, x=lon[inocean], y=lat[inocean])\n", |
| 94 | + "\n", |
| 95 | + " oufile = parcels.ParticleFile(\n", |
| 96 | + " particle_filename,\n", |
| 97 | + " outputdt=np.timedelta64(5, \"D\"),\n", |
| 98 | + " )\n", |
| 99 | + "\n", |
| 100 | + " def AdvectionRK2_periodic(particles, fieldset): # pragma: no cover\n", |
| 101 | + " \"\"\"Advection of particles using second-order Runge-Kutta integration,\n", |
| 102 | + " keeping particles within the periodic domain [-180, 180].\n", |
| 103 | + " \"\"\"\n", |
| 104 | + " (u1, v1) = fieldset.UV[particles]\n", |
| 105 | + " x1 = particles.x + u1 * 0.5 * particles.dt\n", |
| 106 | + " x1 = ((x1 + 180) % 360) - 180\n", |
| 107 | + " y1 = particles.y + v1 * 0.5 * particles.dt\n", |
| 108 | + " (u2, v2) = fieldset.UV[\n", |
| 109 | + " particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles\n", |
| 110 | + " ]\n", |
| 111 | + " particles.dx += u2 * particles.dt\n", |
| 112 | + " particles.dx = ((particles.dx + particles.x + 180) % 360) - (particles.x + 180)\n", |
| 113 | + " particles.dy += v2 * particles.dt\n", |
| 114 | + "\n", |
| 115 | + " pset.execute(\n", |
| 116 | + " [AdvectionRK2_periodic],\n", |
| 117 | + " dt=np.timedelta64(1, \"h\"),\n", |
| 118 | + " runtime=np.timedelta64(365, \"D\"),\n", |
| 119 | + " output_file=oufile,\n", |
| 120 | + " )\n", |
| 121 | + "\n", |
| 122 | + "\n", |
| 123 | + "# run_global_copernicusmarine()" |
65 | 124 | ] |
66 | 125 | }, |
67 | 126 | { |
68 | 127 | "attachments": {}, |
69 | 128 | "cell_type": "markdown", |
70 | 129 | "metadata": {}, |
71 | 130 | "source": [ |
72 | | - "The particle data in `medusarun.nc` is available at https://surfdrive.surf.nl/files/index.php/s/lwjW5w05jtHuYz9\n" |
| 131 | + "Now load in the particle data, and select which time indices to plot\n" |
| 132 | + ] |
| 133 | + }, |
| 134 | + { |
| 135 | + "cell_type": "code", |
| 136 | + "execution_count": null, |
| 137 | + "metadata": {}, |
| 138 | + "outputs": [], |
| 139 | + "source": [ |
| 140 | + "df = parcels.read_particlefile(particle_filename)\n", |
| 141 | + "\n", |
| 142 | + "startframe = 20\n", |
| 143 | + "endframe = 41\n", |
| 144 | + "plottimes = sorted(np.unique(df[\"t\"]))[startframe:endframe]" |
73 | 145 | ] |
74 | 146 | }, |
75 | 147 | { |
76 | 148 | "attachments": {}, |
77 | 149 | "cell_type": "markdown", |
78 | 150 | "metadata": {}, |
79 | 151 | "source": [ |
80 | | - "The animation consists of three figures: the northern hemisphere, the southern hemisphere and the oceanparcels logo. To organize their locations we use [matplotlib.gridspec](https://matplotlib.org/tutorials/intermediate/gridspec.html). The animation spans 12 frames and updates the particle positions based on the timestep in `plottimes`.\n" |
| 152 | + "The animation consists of two figures: the northern hemisphere and the southern hemisphere, using the [Cartopy package](https://cartopy.readthedocs.io/stable/) for map projections.\n" |
81 | 153 | ] |
82 | 154 | }, |
83 | 155 | { |
|
93 | 165 | "fig = plt.figure(figsize=(8, 4))\n", |
94 | 166 | "gs = gridspec.GridSpec(ncols=8, nrows=4, figure=fig)\n", |
95 | 167 | "\n", |
96 | | - "### Northern Hemisphere\n", |
97 | | - "ax1 = fig.add_subplot(\n", |
98 | | - " gs[:, :4],\n", |
99 | | - " projection=ccrs.NearsidePerspective(\n", |
100 | | - " central_latitude=90, central_longitude=-30, satellite_height=15000000\n", |
101 | | - " ),\n", |
102 | | - ")\n", |
103 | | - "ax1.add_feature(cartopy.feature.LAND, zorder=1)\n", |
104 | | - "ax1.add_feature(cartopy.feature.OCEAN, zorder=1)\n", |
105 | | - "ax1.coastlines()\n", |
106 | | - "scat1 = ax1.scatter(\n", |
107 | | - " lon[b],\n", |
108 | | - " lat[b],\n", |
109 | | - " marker=\".\",\n", |
110 | | - " s=25,\n", |
111 | | - " c=\"#AB2200\",\n", |
112 | | - " edgecolor=\"white\",\n", |
113 | | - " linewidth=0.15,\n", |
114 | | - " transform=ccrs.PlateCarree(),\n", |
115 | | - ")\n", |
116 | | - "\n", |
117 | | - "### Southern Hemisphere\n", |
118 | | - "ax2 = fig.add_subplot(\n", |
119 | | - " gs[:, 4:],\n", |
120 | | - " projection=ccrs.NearsidePerspective(\n", |
121 | | - " central_latitude=-90, central_longitude=-30, satellite_height=15000000\n", |
122 | | - " ),\n", |
123 | | - ")\n", |
124 | | - "ax2.add_feature(cartopy.feature.LAND, zorder=1)\n", |
125 | | - "ax2.add_feature(cartopy.feature.OCEAN, zorder=1)\n", |
126 | | - "ax2.coastlines()\n", |
127 | | - "scat2 = ax2.scatter(\n", |
128 | | - " lon[b],\n", |
129 | | - " lat[b],\n", |
130 | | - " marker=\".\",\n", |
131 | | - " s=25,\n", |
132 | | - " c=\"#AB2200\",\n", |
133 | | - " edgecolor=\"white\",\n", |
134 | | - " linewidth=0.15,\n", |
135 | | - " transform=ccrs.PlateCarree(),\n", |
136 | | - ")\n", |
137 | | - "\n", |
138 | | - "frames = np.arange(0, 20)\n", |
139 | | - "\n", |
140 | | - "\n", |
141 | | - "def animate(t):\n", |
142 | | - " b = time == plottimes[t + starttime]\n", |
143 | | - " scat1.set_offsets(np.vstack((lon[b], lat[b])).transpose())\n", |
144 | | - " scat2.set_offsets(np.vstack((lon[b], lat[b])).transpose())\n", |
145 | | - " return scat1, scat2\n", |
146 | | - "\n", |
147 | | - "\n", |
148 | | - "anim = animation.FuncAnimation(fig, animate, frames=frames, interval=150, blit=True)\n", |
149 | | - "anim\n", |
150 | | - "\n", |
151 | | - "# needed for tight_layout to work with cartopy\n", |
152 | | - "fig.canvas.draw()\n", |
153 | | - "plt.tight_layout()\n", |
154 | | - "# writergif = PillowWriter(fps=6)\n", |
155 | | - "# anim.save('homepageshort.gif', writer=writergif, savefig_kwargs={\"transparent\": True})" |
| 168 | + "scat = [None, None]\n", |
| 169 | + "particles = df.filter(pl.col(\"t\") == pl.lit(plottimes[0]))\n", |
| 170 | + "\n", |
| 171 | + "for i, central_latitude in enumerate([90, -90]):\n", |
| 172 | + " ax = fig.add_subplot(\n", |
| 173 | + " gs[:, :4] if central_latitude == 90 else gs[:, 4:],\n", |
| 174 | + " projection=ccrs.NearsidePerspective(\n", |
| 175 | + " central_latitude=central_latitude,\n", |
| 176 | + " central_longitude=-30,\n", |
| 177 | + " satellite_height=15000000,\n", |
| 178 | + " ),\n", |
| 179 | + " )\n", |
| 180 | + " ax.add_feature(cartopy.feature.LAND, zorder=1)\n", |
| 181 | + " ax.add_feature(cartopy.feature.OCEAN, zorder=1)\n", |
| 182 | + " ax.coastlines()\n", |
| 183 | + " scat[i] = ax.scatter(\n", |
| 184 | + " particles[\"x\"],\n", |
| 185 | + " particles[\"y\"],\n", |
| 186 | + " marker=\".\",\n", |
| 187 | + " s=25,\n", |
| 188 | + " c=\"#AB2200\",\n", |
| 189 | + " edgecolor=\"white\",\n", |
| 190 | + " linewidth=0.15,\n", |
| 191 | + " transform=ccrs.PlateCarree(),\n", |
| 192 | + " )\n", |
| 193 | + "\n", |
| 194 | + "\n", |
| 195 | + "def animate(i):\n", |
| 196 | + " particles = df.filter(pl.col(\"t\") == pl.lit(plottimes[i]))\n", |
| 197 | + " scat[0].set_offsets(np.c_[particles[\"x\"], particles[\"y\"]])\n", |
| 198 | + " scat[1].set_offsets(np.c_[particles[\"x\"], particles[\"y\"]])\n", |
| 199 | + " return scat[0], scat[1]\n", |
| 200 | + "\n", |
| 201 | + "\n", |
| 202 | + "plt.rcParams[\"animation.html\"] = \"jshtml\"\n", |
| 203 | + "anim = FuncAnimation(fig, animate, frames=len(plottimes), interval=150, blit=True)\n", |
| 204 | + "plt.close(fig)\n", |
| 205 | + "anim.save(\n", |
| 206 | + " particle_filename.replace(\".parquet\", \".gif\"),\n", |
| 207 | + " writer=PillowWriter(fps=6),\n", |
| 208 | + " savefig_kwargs={\"transparent\": True},\n", |
| 209 | + ")" |
156 | 210 | ] |
157 | 211 | }, |
158 | 212 | { |
|
162 | 216 | "source": [ |
163 | 217 | "The resulting animation is then\n", |
164 | 218 | "\n", |
165 | | - "" |
| 219 | + "" |
166 | 220 | ] |
167 | 221 | } |
168 | 222 | ], |
169 | 223 | "metadata": { |
170 | 224 | "celltoolbar": "Metagegevens bewerken", |
171 | 225 | "kernelspec": { |
172 | | - "display_name": "Python 3", |
| 226 | + "display_name": "Parcels:docs (3.14.6)", |
173 | 227 | "language": "python", |
174 | 228 | "name": "python3" |
175 | 229 | }, |
|
183 | 237 | "name": "python", |
184 | 238 | "nbconvert_exporter": "python", |
185 | 239 | "pygments_lexer": "ipython3", |
186 | | - "version": "3.6.11" |
| 240 | + "version": "3.14.6" |
187 | 241 | } |
188 | 242 | }, |
189 | 243 | "nbformat": 4, |
|
0 commit comments