diff --git a/docs/documentation/index.md b/docs/documentation/index.md index b9f99adf9..d57007514 100644 --- a/docs/documentation/index.md +++ b/docs/documentation/index.md @@ -34,7 +34,7 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4. :caption: Creating ParticleSets :name: tutorial-particlesets - +../examples/tutorial_delaystart.ipynb ``` ```{nbgallery} @@ -43,6 +43,8 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4. ../examples/tutorial_sampling.ipynb +../examples/tutorial_Argofloats.ipynb +../examples/tutorial_gsw_density.ipynb @@ -66,7 +68,7 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4. :caption: Worked examples :name: tutorial-examples -../examples/tutorial_Argofloats.ipynb + ``` diff --git a/docs/examples/tutorial_Argofloats.ipynb b/docs/examples/tutorial_Argofloats.ipynb index 11927c44b..883340fcb 100644 --- a/docs/examples/tutorial_Argofloats.ipynb +++ b/docs/examples/tutorial_Argofloats.ipynb @@ -274,7 +274,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.11" + "version": "3.11.0" } }, "nbformat": 4, diff --git a/docs/examples/tutorial_delaystart.ipynb b/docs/examples/tutorial_delaystart.ipynb new file mode 100644 index 000000000..23e446eb6 --- /dev/null +++ b/docs/examples/tutorial_delaystart.ipynb @@ -0,0 +1,452 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Delayed starts\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In many applications, it is needed to 'delay' the start of particle advection. For example because particles need to be released at different times throughout an experiment. Or because particles need to be released at a constant rate from the same set of locations.\n", + "\n", + "This tutorial will show how this can be done. We start with importing the relevant modules.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import timedelta\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import xarray as xr\n", + "from IPython.display import HTML\n", + "from matplotlib.animation import FuncAnimation\n", + "\n", + "import parcels" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First import a `FieldSet` (from the Argo example, in this case)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n", + "example_dataset_folder = parcels.download_example_dataset(\n", + " \"CopernicusMarine_data_for_Argo_tutorial\"\n", + ")\n", + "\n", + "ds = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n", + "ds.load() # load the dataset into memory\n", + "\n", + "fieldset = parcels.FieldSet.from_copernicusmarine(ds)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining initial `time` as particle release\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Defining the initial times of particles is done when the `ParticleSet` is defined. Although `time` and `z` are optional arguments (with FieldSet t=0 and z=0 as defaults), it is good practice to define them explicitly to ensure expected behavior. The simplest way to delay the start of a particle is to use the `time` argument for each particle.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "npart = 10 # number of particles to be released\n", + "lon = 32 * np.ones(npart)\n", + "lat = np.linspace(-31.5, -30.5, npart, dtype=np.float32)\n", + "# release every particle one hour later from the initial fieldset time\n", + "time = ds.time.values[0] + np.arange(0, npart) * np.timedelta64(1, \"h\")\n", + "z = np.repeat(ds.depth.values[0], npart)\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we can execute the `pset` as usual\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "output_file = parcels.ParticleFile(\n", + " \"DelayParticle_time.zarr\", outputdt=np.timedelta64(1, \"h\")\n", + ")\n", + "\n", + "pset.execute(\n", + " parcels.kernels.AdvectionRK4,\n", + " runtime=np.timedelta64(24, \"h\"),\n", + " dt=np.timedelta64(5, \"m\"),\n", + " output_file=output_file,\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And then finally, we can show a movie of the particles. Note that the southern-most particles start to move first.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "\n", + "ds_out = xr.open_zarr(\"DelayParticle_time.zarr\")\n", + "\n", + "fig = plt.figure(figsize=(7, 5), constrained_layout=True)\n", + "ax = fig.add_subplot()\n", + "\n", + "ax.set_ylabel(\"Latitude [deg N]\")\n", + "ax.set_xlabel(\"Longitude [deg E]\")\n", + "ax.set_xlim(31, 33)\n", + "ax.set_ylim(-32, -30)\n", + "\n", + "timerange = np.unique(ds_out[\"time\"].values[np.isfinite(ds_out[\"time\"])])\n", + "\n", + "# Indices of the data where time = 0\n", + "time_id = np.where(ds_out[\"time\"] == timerange[0])\n", + "\n", + "sc = ax.scatter(ds_out[\"lon\"].values[time_id], ds_out[\"lat\"].values[time_id])\n", + "\n", + "t = timerange[0].astype(\"datetime64[h]\")\n", + "title = ax.set_title(f\"Particles at t = {t}\")\n", + "\n", + "\n", + "def animate(i):\n", + " t = timerange[i].astype(\"datetime64[h]\")\n", + " title.set_text(f\"Particles at t = {t}\")\n", + "\n", + " time_id = np.where(ds_out[\"time\"] == timerange[i])\n", + " sc.set_offsets(np.c_[ds_out[\"lon\"].values[time_id], ds_out[\"lat\"].values[time_id]])\n", + "\n", + "\n", + "anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "HTML(anim.to_jshtml())" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Release particles repeatedly at a set frequency using `np.broadcast_to`\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When you want to repeatedly release particles from the same set of locations, you can expand the initial array of particle release locations. Here we show how to release particles at a set frequency `repeatdt`, for a fixed number of releases `nrepeat`. The total number of particles released is then **nrepeat** x **npart**." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Note that the `repeatdt` argument in `parcels.ParticleSet` used in previous versions of parcels is no longer supported as of v4\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "npart = 10 # number of particles to be released\n", + "\n", + "lon_i = 32 * np.ones(npart)\n", + "lat_i = np.linspace(-31.5, -30.5, npart, dtype=np.float32)\n", + "time_i = np.repeat(ds.time.values[0], npart)\n", + "z_i = np.repeat(ds.depth.values[0], npart)\n", + "\n", + "# repeat release at frequency `repeatdt` for `nrepeat` different releases\n", + "repeatdt = np.timedelta64(6, \"h\")\n", + "nrepeat = 3\n", + "\n", + "lon = np.broadcast_to(lon_i, (nrepeat, npart))\n", + "lat = np.broadcast_to(lat_i, (nrepeat, npart))\n", + "time = (\n", + " np.broadcast_to(time_i, (nrepeat, npart))\n", + " + np.arange(0, nrepeat)[:, np.newaxis] * repeatdt\n", + ")\n", + "print(f\"pset.time shape = {time.shape}\")\n", + "z = np.broadcast_to(z_i, (nrepeat, npart))\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we again define an output file and execute the `pset` as usual.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "output_file = parcels.ParticleFile(\n", + " \"DelayParticle_releasedt.zarr\", outputdt=timedelta(hours=1)\n", + ")\n", + "\n", + "pset.execute(\n", + " parcels.kernels.AdvectionRK4,\n", + " runtime=timedelta(hours=24),\n", + " dt=timedelta(minutes=5),\n", + " output_file=output_file,\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And we get an animation where a new particle is released every 6 hours from each start location\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "\n", + "ds_out = xr.open_zarr(\"DelayParticle_releasedt.zarr\")\n", + "\n", + "fig = plt.figure(figsize=(7, 5), constrained_layout=True)\n", + "ax = fig.add_subplot()\n", + "\n", + "ax.set_ylabel(\"Latitude [deg N]\")\n", + "ax.set_xlabel(\"Longitude [deg E]\")\n", + "ax.set_xlim(31, 33)\n", + "ax.set_ylim(-32, -30)\n", + "\n", + "timerange = np.unique(ds_out[\"time\"].values[np.isfinite(ds_out[\"time\"])])\n", + "\n", + "# Indices of the data where time = 0\n", + "time_id = np.where(ds_out[\"time\"] == timerange[0])\n", + "\n", + "sc = ax.scatter(ds_out[\"lon\"].values[time_id], ds_out[\"lat\"].values[time_id])\n", + "\n", + "t = timerange[0].astype(\"datetime64[h]\")\n", + "title = ax.set_title(f\"Particles at t = {t}\")\n", + "\n", + "\n", + "def animate(i):\n", + " t = timerange[i].astype(\"datetime64[h]\")\n", + " title.set_text(f\"Particles at t = {t}\")\n", + "\n", + " time_id = np.where(ds_out[\"time\"] == timerange[i])\n", + " sc.set_offsets(np.c_[ds_out[\"lon\"].values[time_id], ds_out[\"lat\"].values[time_id]])\n", + "\n", + "\n", + "anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "HTML(anim.to_jshtml())" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Synced `time` in the output file\n", + "\n", + "Note that, because the `outputdt` variable controls the Kernel-loop, all particles are written _at the same time_, even when they start at a non-multiple of `outputdt`.\n", + "\n", + "For example, if your particles start at `time=[0, 1, 2]` and `outputdt=2`, then the times written (for `dt=1` and `endtime=4`) will be\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "outtime_expected = np.array(\n", + " [[0, 2, 4], [2, 4, np.datetime64(\"NaT\")], [2, 4, np.datetime64(\"NaT\")]],\n", + " dtype=\"timedelta64[h]\",\n", + ")\n", + "print(outtime_expected)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "outfilepath = \"DelayParticle_nonmatchingtime.zarr\"\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " pclass=parcels.Particle,\n", + " lat=[-31] * 3,\n", + " lon=[32] * 3,\n", + " time=ds.time.values[0] + np.arange(3) * np.timedelta64(1, \"h\"),\n", + " z=[0.5] * 3,\n", + ")\n", + "\n", + "output_file = parcels.ParticleFile(outfilepath, outputdt=np.timedelta64(2, \"h\"))\n", + "pset.execute(\n", + " parcels.kernels.AdvectionRK4,\n", + " endtime=ds.time.values[0] + np.timedelta64(4, \"h\"),\n", + " dt=np.timedelta64(1, \"h\"),\n", + " output_file=output_file,\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And indeed, the `time` values in the NetCDF output file are as expected\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "outtime_infile = (\n", + " xr.open_zarr(outfilepath).time.values[:] - ds.time.values[0]\n", + ") # subtract initial time to convert from datetime64 to timedelta64\n", + "print(outtime_infile.astype(\"timedelta64[h]\"))\n", + "\n", + "assert (\n", + " outtime_expected[np.isfinite(outtime_expected)]\n", + " == outtime_infile[np.isfinite(outtime_infile)]\n", + ").all()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, for some applications, this behavior may be undesirable; for example when particles need to be analyzed at a same age (instead of at a same time). In that case, we recommend either changing `outputdt` so that it is a common divisor of all start times; or doing multiple Parcels runs with subsets of the original `ParticleSet` (e.g., in the example above, one run with the Particles that start at `time=[0, 2]` and one with the Particle at `time=[1]`). In that case, you will get two files:\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for times in np.array([[0, 2], [1, 3]]).astype(\"timedelta64[h]\"):\n", + " pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " pclass=parcels.Particle,\n", + " lat=[-31] * len(times),\n", + " lon=[32] * len(times),\n", + " time=ds.time.values[0] + times,\n", + " z=[0.5] * len(times),\n", + " )\n", + " output_file = parcels.ParticleFile(outfilepath, outputdt=np.timedelta64(2, \"h\"))\n", + " pset.execute(\n", + " parcels.kernels.AdvectionRK4,\n", + " runtime=np.timedelta64(4, \"h\"),\n", + " dt=np.timedelta64(1, \"h\"),\n", + " output_file=output_file,\n", + " )\n", + " outtime_infile = xr.open_zarr(outfilepath).time.values[:] - ds.time.values[0]\n", + " print(outtime_infile.astype(\"timedelta64[h]\"))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "test-notebooks", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/docs/examples_v3/tutorial_output.ipynb b/docs/examples/tutorial_output.ipynb similarity index 62% rename from docs/examples_v3/tutorial_output.ipynb rename to docs/examples/tutorial_output.ipynb index f562f4e90..13778feca 100644 --- a/docs/examples_v3/tutorial_output.ipynb +++ b/docs/examples/tutorial_output.ipynb @@ -21,7 +21,7 @@ "- [**Plotting**](#Plotting)\n", "- [**Animations**](#Animations)\n", "\n", - "First we need to create some parcels output to analyze. We simulate a set of particles using the setup described in the [Delay start tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_delaystart.html). We will also add some user defined metadata to the output file.\n" + "For more advanced reading and tutorials on the analysis of Lagrangian trajectories, we recommend checking out the [Lagrangian Diagnostics](https://lagrangian-diags.readthedocs.io/en/latest/index.html) project. The [TrajAn package]() can be used to read and plot datasets of Lagrangian trajectories." ] }, { @@ -33,47 +33,60 @@ "from datetime import datetime, timedelta\n", "\n", "import numpy as np\n", + "import xarray as xr\n", "\n", "import parcels" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First we need to create some parcels output to analyze. We simulate a set of particles using the setup described in the [Delay start tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_delaystart.html). We will also add some user defined metadata to the output file." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "example_dataset_folder = parcels.download_example_dataset(\"Peninsula_data\")\n", - "example_dataset_folder = parcels.download_example_dataset(\"Peninsula_data\")\n", - "filenames = {\n", - " \"U\": str(example_dataset_folder / \"peninsulaU.nc\"),\n", - " \"V\": str(example_dataset_folder / \"peninsulaV.nc\"),\n", - "}\n", - "variables = {\"U\": \"vozocrtx\", \"V\": \"vomecrty\"}\n", - "dimensions = {\"lon\": \"nav_lon\", \"lat\": \"nav_lat\", \"time\": \"time_counter\"}\n", - "fieldset = parcels.FieldSet.from_netcdf(\n", - " filenames, variables, dimensions, allow_time_extrapolation=True\n", + "# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n", + "example_dataset_folder = parcels.download_example_dataset(\n", + " \"CopernicusMarine_data_for_Argo_tutorial\"\n", ")\n", "\n", - "npart = 10 # number of particles to be released\n", - "lon = 3e3 * np.ones(npart)\n", - "lat = np.linspace(3e3, 45e3, npart, dtype=np.float32)\n", + "ds = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n", + "ds.load() # load the dataset into memory\n", "\n", - "# release every particle two hours later\n", - "time = np.arange(npart) * timedelta(hours=2).total_seconds()\n", + "fieldset = parcels.FieldSet.from_copernicusmarine(ds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Particle locations and initial time\n", + "npart = 10 # number of particles to be released\n", + "lon = 32 * np.ones(npart)\n", + "lat = np.linspace(-32.5, -30.5, npart, dtype=np.float32)\n", + "time = ds.time.values[0] + np.arange(0, npart) * np.timedelta64(2, \"h\")\n", + "z = np.repeat(ds.depth.values[0], npart)\n", "\n", "pset = parcels.ParticleSet(\n", - " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time\n", + " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z\n", ")\n", "\n", - "output_file = pset.ParticleFile(name=\"Output.zarr\", outputdt=timedelta(hours=2))" + "output_file = parcels.ParticleFile(\"Output.zarr\", outputdt=np.timedelta64(2, \"h\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Parcels saves some metadata in the output file with every simulation (Parcels version, CF convention information, etc.). This metadata is just a dictionary which is propogated to `xr.Dataset(attrs=...)` and is stored in the `.metadata` attribute. The user is free to manipulate this dictionary to add any custom, xarray-compatible metadata relevant to their simulation. Here we add a custom metadata field `date_created` to the output file." + "Parcels saves some metadata in the output file with every simulation (Parcels version, CF convention information, etc.). This metadata is just a dictionary which is propogated to `xr.Dataset(attrs=...)` and is stored in the `.metadata` attribute. We are free to manipulate this dictionary to add any custom, xarray-compatible metadata relevant to their simulation. Here we add a custom metadata field `date_created` to the output file." ] }, { @@ -82,7 +95,15 @@ "metadata": {}, "outputs": [], "source": [ - "output_file.metadata[\"date_created\"] = datetime.now().isoformat()" + "output_file.metadata[\"date_created\"] = datetime.now().isoformat()\n", + "output_file.metadata" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To write the metadata to the output_file, we need to add it before running `pset.execute()` which writes the particleset including the metadata to the output_file." ] }, { @@ -93,8 +114,8 @@ "source": [ "pset.execute(\n", " parcels.kernels.AdvectionRK4,\n", - " runtime=timedelta(hours=24),\n", - " dt=timedelta(minutes=5),\n", + " runtime=np.timedelta64(48, \"h\"),\n", + " dt=np.timedelta64(5, \"m\"),\n", " output_file=output_file,\n", ")" ] @@ -106,7 +127,7 @@ "source": [ "## Reading the output file\n", "\n", - "Parcels exports output trajectories in `zarr` [format](https://zarr.readthedocs.io/en/stable/). Files in `zarr` are typically _much_ smaller in size than netcdf, although may be slightly more challenging to handle (but `xarray` has a fairly seamless `open_zarr()` method). Note when when we display the dataset we cam see our custom metadata field `date_created`.\n" + "Parcels exports output trajectories in `zarr` [format](https://zarr.readthedocs.io/en/stable/). Files in `zarr` are typically _much_ smaller in size than netcdf, although may be slightly more challenging to handle (but `xarray` has a fairly seamless `open_zarr()` method). Note when we display the dataset we can see our custom metadata field `date_created`.\n" ] }, { @@ -115,21 +136,11 @@ "metadata": {}, "outputs": [], "source": [ - "import xarray as xr\n", - "\n", "data_xarray = xr.open_zarr(\"Output.zarr\")\n", + "\n", "print(data_xarray)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(data_xarray[\"trajectory\"])" - ] - }, { "attachments": {}, "cell_type": "markdown", @@ -160,10 +171,9 @@ "source": [ "np.set_printoptions(linewidth=160)\n", "one_hour = np.timedelta64(1, \"h\") # Define timedelta object to help with conversion\n", + "time_from_start = data_xarray[\"time\"].values - fieldset.time_interval.left\n", "\n", - "print(\n", - " data_xarray[\"time\"].data.compute() / one_hour\n", - ") # timedelta / timedelta -> float number of hours" + "print(time_from_start / one_hour) # timedelta / timedelta -> float number of hours" ] }, { @@ -198,9 +208,9 @@ " np.sqrt(np.square(np.diff(x)) + np.square(np.diff(y))), axis=1\n", ") # d = (dx^2 + dy^2)^(1/2)\n", "\n", - "real_time = data_xarray[\"time\"] / one_hour # convert time to hours\n", + "real_time = time_from_start / one_hour # convert time to hours\n", "time_since_release = (\n", - " real_time.values.transpose() - real_time.values[:, 0]\n", + " real_time.transpose() - real_time[:, 0]\n", ") # substract the initial time from each timeseries" ] }, @@ -361,7 +371,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Trajectory plots like the ones above can become very cluttered for large sets of particles. To better see patterns, it's a good idea to create an animation in time and space. To do this, matplotlib offers an [animation package](https://matplotlib.org/stable/api/animation_api.html). Here we show how to use the [**FuncAnimation**](https://matplotlib.org/3.3.2/api/_as_gen/matplotlib.animation.FuncAnimation.html#matplotlib.animation.FuncAnimation) class to animate parcels trajectory data.\n", + "Trajectory plots like the ones above can become very cluttered for large sets of particles. To better see patterns, it's a good idea to create an animation in time and space. To do this, matplotlib offers an [animation package](https://matplotlib.org/stable/api/animation_api.html). Here we show how to use the [**FuncAnimation**](https://matplotlib.org/3.3.2/api/_as_gen/matplotlib.animation.FuncAnimation.html#matplotlib.animation.FuncAnimation) class to animate parcels trajectory data, based on [this visualisation tutorial](https://github.com/Parcels-code/10year-anniversary-session5/blob/eaf7ac35f43c222280fa5577858be81dc346c06b/animations_tutorial.ipynb) from 10-years Parcels. \n", "\n", "To correctly reveal the patterns in time we must remember that the `obs` dimension does not necessarily correspond to the `time` variable ([see the section of Trajectory data structure above](#Trajectory-data-structure)). In the animation of the particles, we usually want to draw the points at each consecutive moment in time, not necessarily at each moment since the start of the trajectory. To do this we must [select the correct data](#Conditional-selection) in each rendering.\n" ] @@ -372,17 +382,10 @@ "metadata": {}, "outputs": [], "source": [ - "from IPython.display import HTML\n", - "from matplotlib.animation import FuncAnimation\n", - "\n", - "outputdt = timedelta(hours=2)\n", - "\n", - "# timerange in nanoseconds\n", - "timerange = np.arange(\n", - " np.nanmin(data_xarray[\"time\"].values),\n", - " np.nanmax(data_xarray[\"time\"].values) + np.timedelta64(outputdt),\n", - " outputdt,\n", - ")" + "import cartopy.crs as ccrs\n", + "import cartopy.feature as cfeature\n", + "import matplotlib\n", + "from matplotlib.animation import FuncAnimation" ] }, { @@ -391,38 +394,149 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", - "fig = plt.figure(figsize=(5, 5), constrained_layout=True)\n", - "ax = fig.add_subplot()\n", - "\n", - "ax.set_ylabel(\"Meridional distance [m]\")\n", - "ax.set_xlabel(\"Zonal distance [m]\")\n", - "ax.set_xlim(0, 90000)\n", - "ax.set_ylim(0, 50000)\n", - "plt.xticks(rotation=45)\n", - "\n", - "# Indices of the data where time = 0\n", - "time_id = np.where(data_xarray[\"time\"] == timerange[0])\n", - "\n", - "scatter = ax.scatter(\n", - " data_xarray[\"lon\"].values[time_id], data_xarray[\"lat\"].values[time_id]\n", - ")\n", - "\n", - "t = str(timerange[0].astype(\"timedelta64[h]\"))\n", - "title = ax.set_title(\"Particles at t = \" + t)\n", - "\n", - "\n", - "def animate(i):\n", - " t = str(timerange[i].astype(\"timedelta64[h]\"))\n", - " title.set_text(\"Particles at t = \" + t)\n", - "\n", - " time_id = np.where(data_xarray[\"time\"] == timerange[i])\n", - " scatter.set_offsets(\n", - " np.c_[data_xarray[\"lon\"].values[time_id], data_xarray[\"lat\"].values[time_id]]\n", + "# for interactive display of animation\n", + "plt.rcParams[\"animation.html\"] = \"jshtml\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Number of timesteps to animate\n", + "nframes = 25 # use less frames for testing purposes\n", + "nreducedtrails = 1 # every 10th particle will have a trail (if 1, all particles have trails. Adjust for faster performance)\n", + "\n", + "\n", + "# Set up the colors and associated trajectories:\n", + "# get release times for each particle (first valide obs for each trajectory)\n", + "release_times = data_xarray[\"time\"].min(dim=\"obs\", skipna=True).values\n", + "\n", + "# get unique release times and assign colors\n", + "unique_release_times = np.unique(release_times[~np.isnat(release_times)])\n", + "n_release_times = len(unique_release_times)\n", + "print(f\"Number of unique release times: {n_release_times}\")\n", + "\n", + "# choose a continuous colormap\n", + "colormap = matplotlib.colormaps[\"tab20b\"]\n", + "\n", + "# set up a unique color for each release time\n", + "release_time_to_color = {}\n", + "for i, release_time in enumerate(unique_release_times):\n", + " release_time_to_color[release_time] = colormap(i / max(n_release_times - 1, 1))\n", + "\n", + "\n", + "# --> Store data for all timeframes (this is needed for faster performance)\n", + "print(\"Pre-computing all particle positions...\")\n", + "all_particles_data = []\n", + "for i, target_time in enumerate(timerange):\n", + " time_id = np.where(data_xarray[\"time\"] == target_time)\n", + " lons = data_xarray[\"lon\"].values[time_id]\n", + " lats = data_xarray[\"lat\"].values[time_id]\n", + " particle_indices = time_id[0]\n", + " valid = ~np.isnan(lons) & ~np.isnan(lats)\n", + "\n", + " all_particles_data.append(\n", + " {\n", + " \"lons\": lons[valid],\n", + " \"lats\": lats[valid],\n", + " \"particle_indices\": particle_indices[valid],\n", + " \"valid_count\": np.sum(valid),\n", + " }\n", " )\n", "\n", "\n", - "anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)" + "# figure setup\n", + "fig, ax = plt.subplots(figsize=(6, 5), subplot_kw={\"projection\": ccrs.PlateCarree()})\n", + "ax.set_xlim(30, 33)\n", + "ax.set_xticks(np.arange(30, 33.5, 0.5))\n", + "ax.set_xlabel(\"Longitude (deg E)\")\n", + "ax.set_ylim(-33, -30)\n", + "ax.set_yticks(ticks=np.arange(-33, -29.5, 0.5))\n", + "ax.set_yticklabels(np.arange(33, 29.5, -0.5).astype(str))\n", + "ax.set_ylabel(\"Latitude (deg S)\")\n", + "ax.coastlines(color=\"saddlebrown\")\n", + "ax.add_feature(cfeature.LAND, alpha=0.5, facecolor=\"saddlebrown\")\n", + "\n", + "# --> Use pre-computed data for initial setup\n", + "initial_data = all_particles_data[0]\n", + "initial_colors = []\n", + "for particle_idx in initial_data[\"particle_indices\"]:\n", + " rt = release_times[particle_idx]\n", + " if rt in release_time_to_color:\n", + " initial_colors.append(release_time_to_color[rt])\n", + " else:\n", + " initial_colors.append(\"blue\")\n", + "\n", + "# --> plot first timestep\n", + "scatter = ax.scatter(initial_data[\"lons\"], initial_data[\"lats\"], s=10, c=initial_colors)\n", + "\n", + "# --> initialize trails\n", + "trail_plot = []\n", + "\n", + "# Set initial title\n", + "t_str = str(timerange[0])[:19] # Format datetime nicely\n", + "title = ax.set_title(f\"Particles at t = {t_str}\")\n", + "\n", + "\n", + "# loop over for animation\n", + "def animate(i):\n", + " print(f\"Animating frame {i + 1}/{len(timerange)} at time {timerange[i]}\")\n", + " t_str = str(timerange[i])[:19]\n", + " title.set_text(f\"Particles at t = {t_str}\")\n", + "\n", + " # Find particles at current time\n", + " current_data = all_particles_data[i]\n", + "\n", + " if current_data[\"valid_count\"] > 0:\n", + " current_colors = []\n", + " for particle_idx in current_data[\"particle_indices\"]:\n", + " rt = release_times[particle_idx]\n", + " current_colors.append(release_time_to_color[rt])\n", + "\n", + " scatter.set_offsets(np.c_[current_data[\"lons\"], current_data[\"lats\"]])\n", + " scatter.set_color(current_colors)\n", + "\n", + " # --> add trails\n", + "\n", + " for trail in trail_plot:\n", + " trail.remove()\n", + " trail_plot.clear()\n", + "\n", + " trail_length = min(10, i) # trails will have max length of 10 time steps\n", + "\n", + " if trail_length > 0:\n", + " sampled_particles = current_data[\"particle_indices\"][\n", + " ::nreducedtrails\n", + " ] # use all or sample if you want faster computation\n", + "\n", + " for particle_idx in sampled_particles:\n", + " trail_lons = []\n", + " trail_lats = []\n", + " for j in range(i - trail_length, i + 1):\n", + " past_data = all_particles_data[j]\n", + " if particle_idx in past_data[\"particle_indices\"]:\n", + " idx = np.where(past_data[\"particle_indices\"] == particle_idx)[\n", + " 0\n", + " ][0]\n", + " trail_lons.append(past_data[\"lons\"][idx])\n", + " trail_lats.append(past_data[\"lats\"][idx])\n", + " if len(trail_lons) > 1:\n", + " rt = release_times[particle_idx]\n", + " color = release_time_to_color[rt]\n", + " (trail,) = ax.plot(\n", + " trail_lons, trail_lats, color=color, linewidth=0.6, alpha=0.6\n", + " )\n", + " trail_plot.append(trail)\n", + "\n", + " else:\n", + " scatter.set_offsets(np.empty((0, 2)))\n", + "\n", + "\n", + "# Create animation\n", + "anim = FuncAnimation(fig, animate, frames=nframes, interval=100)\n", + "anim" ] }, { @@ -430,15 +544,13 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "HTML(anim.to_jshtml())" - ] + "source": [] } ], "metadata": { "celltoolbar": "Metagegevens bewerken", "kernelspec": { - "display_name": "parcels", + "display_name": "test-notebooks", "language": "python", "name": "python3" }, @@ -452,7 +564,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.3" + "version": "3.11.0" } }, "nbformat": 4, diff --git a/docs/examples_v3/tutorial_delaystart.ipynb b/docs/examples_v3/tutorial_delaystart.ipynb deleted file mode 100644 index 6584dfc15..000000000 --- a/docs/examples_v3/tutorial_delaystart.ipynb +++ /dev/null @@ -1,619 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Delayed starts\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In many applications, it is needed to 'delay' the start of particle advection. For example because particles need to be released at different times throughout an experiment. Or because particles need to be released at a constant rate from the same set of locations.\n", - "\n", - "This tutorial will show how this can be done. We start with importing the relevant modules.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datetime import timedelta\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import xarray as xr\n", - "from IPython.display import HTML\n", - "from matplotlib.animation import FuncAnimation\n", - "\n", - "import parcels" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "First import a `FieldSet` (from the Peninsula example, in this case)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "example_dataset_folder = parcels.download_example_dataset(\"Peninsula_data\")\n", - "filenames = {\n", - " \"U\": str(example_dataset_folder / \"peninsulaU.nc\"),\n", - " \"V\": str(example_dataset_folder / \"peninsulaV.nc\"),\n", - "}\n", - "variables = {\"U\": \"vozocrtx\", \"V\": \"vomecrty\"}\n", - "dimensions = {\"lon\": \"nav_lon\", \"lat\": \"nav_lat\", \"time\": \"time_counter\"}\n", - "fieldset = parcels.FieldSet.from_netcdf(\n", - " filenames, variables, dimensions, allow_time_extrapolation=True\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, there are two ways to delay the start of particles. Either by defining the whole `ParticleSet` at initialisation and giving each particle its own `time`. Or by using the `repeatdt` argument. We will show both options here\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Assigning each particle its own `time`\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The simplest way to delay the start of a particle is to use the `time` argument for each particle\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "npart = 10 # number of particles to be released\n", - "lon = 3e3 * np.ones(npart)\n", - "lat = np.linspace(3e3, 45e3, npart, dtype=np.float32)\n", - "\n", - "# release every particle one hour later\n", - "time = np.arange(0, npart) * timedelta(hours=1).total_seconds()\n", - "\n", - "pset = parcels.ParticleSet(\n", - " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Then we can execute the `pset` as usual\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "output_file = pset.ParticleFile(\n", - " name=\"DelayParticle_time.zarr\", outputdt=timedelta(hours=1)\n", - ")\n", - "\n", - "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", - " runtime=timedelta(hours=24),\n", - " dt=timedelta(minutes=5),\n", - " output_file=output_file,\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And then finally, we can show a movie of the particles. Note that the southern-most particles start to move first.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%capture\n", - "\n", - "ds = xr.open_zarr(\"DelayParticle_time.zarr\")\n", - "\n", - "fig = plt.figure(figsize=(7, 5), constrained_layout=True)\n", - "ax = fig.add_subplot()\n", - "\n", - "ax.set_ylabel(\"Meridional distance [m]\")\n", - "ax.set_xlabel(\"Zonal distance [m]\")\n", - "ax.set_xlim(0, 9e4)\n", - "ax.set_ylim(0, 5e4)\n", - "\n", - "timerange = np.unique(ds[\"time\"].values[np.isfinite(ds[\"time\"])])\n", - "\n", - "# Indices of the data where time = 0\n", - "time_id = np.where(ds[\"time\"] == timerange[0])\n", - "\n", - "sc = ax.scatter(ds[\"lon\"].values[time_id], ds[\"lat\"].values[time_id])\n", - "\n", - "t = str(timerange[0].astype(\"timedelta64[h]\"))\n", - "title = ax.set_title(f\"Particles at t = {t}\")\n", - "\n", - "\n", - "def animate(i):\n", - " t = str(timerange[i].astype(\"timedelta64[h]\"))\n", - " title.set_text(f\"Particles at t = {t}\")\n", - "\n", - " time_id = np.where(ds[\"time\"] == timerange[i])\n", - " sc.set_offsets(np.c_[ds[\"lon\"].values[time_id], ds[\"lat\"].values[time_id]])\n", - "\n", - "\n", - "anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "HTML(anim.to_jshtml())" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using the `repeatdt` argument\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The second method to delay the start of particle releases is to use the `repeatdt` argument when constructing a `ParticleSet`. This is especially useful if you want to repeatedly release particles from the same set of locations.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "npart = 10 # number of particles to be released\n", - "lon = 3e3 * np.ones(npart)\n", - "lat = np.linspace(3e3, 45e3, npart, dtype=np.float32)\n", - "repeatdt = timedelta(hours=3) # release from the same set of locations every 3h\n", - "\n", - "pset = parcels.ParticleSet(\n", - " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, repeatdt=repeatdt\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we again define an output file and execute the `pset` as usual.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "output_file = pset.ParticleFile(\n", - " name=\"DelayParticle_releasedt\", outputdt=timedelta(hours=1)\n", - ")\n", - "\n", - "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", - " runtime=timedelta(hours=24),\n", - " dt=timedelta(minutes=5),\n", - " output_file=output_file,\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And we get an animation where a new particle is released every 3 hours from each start location\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%capture\n", - "\n", - "ds = xr.open_zarr(\"DelayParticle_releasedt.zarr\")\n", - "\n", - "fig = plt.figure(figsize=(7, 5), constrained_layout=True)\n", - "ax = fig.add_subplot()\n", - "\n", - "ax.set_ylabel(\"Meridional distance [m]\")\n", - "ax.set_xlabel(\"Zonal distance [m]\")\n", - "ax.set_xlim(0, 9e4)\n", - "ax.set_ylim(0, 5e4)\n", - "\n", - "timerange = np.unique(ds[\"time\"].values[np.isfinite(ds[\"time\"])])\n", - "\n", - "# Indices of the data where time = 0\n", - "time_id = np.where(ds[\"time\"] == timerange[0])\n", - "\n", - "sc = ax.scatter(ds[\"lon\"].values[time_id], ds[\"lat\"].values[time_id])\n", - "\n", - "t = str(timerange[0].astype(\"timedelta64[h]\"))\n", - "title = ax.set_title(f\"Particles at t = {t}\")\n", - "\n", - "\n", - "def animate(i):\n", - " t = str(timerange[i].astype(\"timedelta64[h]\"))\n", - " title.set_text(f\"Particles at t = {t}\")\n", - "\n", - " time_id = np.where(ds[\"time\"] == timerange[i])\n", - " sc.set_offsets(np.c_[ds[\"lon\"].values[time_id], ds[\"lat\"].values[time_id]])\n", - "\n", - "\n", - "anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "HTML(anim.to_jshtml())" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that, if you want to if you want to at some point stop the repeatdt, the easiest implementation is to use two calls to `pset.execute()`. For example, if in the above example you only want four releases of the pset, you could do the following\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "pset = parcels.ParticleSet(\n", - " fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, repeatdt=repeatdt\n", - ")\n", - "output_file = pset.ParticleFile(\n", - " name=\"DelayParticle_releasedt_9hrs\", outputdt=timedelta(hours=1)\n", - ")\n", - "\n", - "# first run for 3 * 3 hrs\n", - "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", - " runtime=timedelta(hours=9),\n", - " dt=timedelta(minutes=5),\n", - " output_file=output_file,\n", - ")\n", - "\n", - "# now stop the repeated release\n", - "pset.repeatdt = None\n", - "\n", - "# now continue running for the remaining 15 hours\n", - "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", - " runtime=timedelta(hours=15),\n", - " dt=timedelta(minutes=5),\n", - " output_file=output_file,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%capture\n", - "\n", - "ds = xr.open_zarr(\"DelayParticle_releasedt_9hrs.zarr\")\n", - "\n", - "fig = plt.figure(figsize=(7, 5), constrained_layout=True)\n", - "ax = fig.add_subplot()\n", - "\n", - "ax.set_ylabel(\"Meridional distance [m]\")\n", - "ax.set_xlabel(\"Zonal distance [m]\")\n", - "ax.set_xlim(0, 9e4)\n", - "ax.set_ylim(0, 5e4)\n", - "\n", - "timerange = np.unique(ds[\"time\"].values[np.isfinite(ds[\"time\"])])\n", - "\n", - "# Indices of the data where time = 0\n", - "time_id = np.where(ds[\"time\"] == timerange[0])\n", - "\n", - "sc = ax.scatter(ds[\"lon\"].values[time_id], ds[\"lat\"].values[time_id])\n", - "\n", - "t = str(timerange[0].astype(\"timedelta64[h]\"))\n", - "title = ax.set_title(f\"Particles at t = {t}\")\n", - "\n", - "\n", - "def animate(i):\n", - " t = str(timerange[i].astype(\"timedelta64[h]\"))\n", - " title.set_text(f\"Particles at t = {t}\")\n", - "\n", - " time_id = np.where(ds[\"time\"] == timerange[i])\n", - " sc.set_offsets(np.c_[ds[\"lon\"].values[time_id], ds[\"lat\"].values[time_id]])\n", - "\n", - "\n", - "anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "HTML(anim.to_jshtml())" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Synced `time` in the output file\n", - "\n", - "Note that, because the `outputdt` variable controls the Kernel-loop, all particles are written _at the same time_, even when they start at a non-multiple of `outputdt`.\n", - "\n", - "For example, if your particles start at `time=[0, 1, 2]` and `outputdt=2`, then the times written (for `dt=1` and `endtime=4`) will be\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "outtime_expected = np.array(\n", - " [[0, 2, 4], [2, 4, np.datetime64(\"NaT\")], [2, 4, np.datetime64(\"NaT\")]],\n", - " dtype=\"timedelta64[s]\",\n", - ")\n", - "print(outtime_expected)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "outfilepath = \"DelayParticle_nonmatchingtime.zarr\"\n", - "\n", - "pset = parcels.ParticleSet(\n", - " fieldset=fieldset,\n", - " pclass=parcels.Particle,\n", - " lat=[3e3] * 3,\n", - " lon=[3e3] * 3,\n", - " time=[0, 1, 2],\n", - ")\n", - "\n", - "output_file = pset.ParticleFile(name=outfilepath, outputdt=2)\n", - "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", - " endtime=4,\n", - " dt=1,\n", - " output_file=output_file,\n", - ")\n", - "\n", - "# Note that we also need to write the final time to the file\n", - "output_file.write_latest_locations(pset, 4)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And indeed, the `time` values in the NetCDF output file are as expected\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "outtime_infile = xr.open_zarr(outfilepath).time.values[:]\n", - "print(outtime_infile.astype(\"timedelta64[s]\"))\n", - "\n", - "assert (\n", - " outtime_expected[np.isfinite(outtime_expected)]\n", - " == outtime_infile[np.isfinite(outtime_infile)]\n", - ").all()" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, for some applications, this behavior may be undesirable; for example when particles need to be analyzed at a same age (instead of at a same time). In that case, we recommend either changing `outputdt` so that it is a common divisor of all start times; or doing multiple Parcels runs with subsets of the original `ParticleSet` (e.g., in the example above, one run with the Particles that start at `time=[0, 2]` and one with the Particle at `time=[1]`). In that case, you will get two files:\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for times in [[0, 2], [1]]:\n", - " pset = parcels.ParticleSet(\n", - " fieldset=fieldset,\n", - " pclass=parcels.Particle,\n", - " lat=[3e3] * len(times),\n", - " lon=[3e3] * len(times),\n", - " time=times,\n", - " )\n", - " output_file = pset.ParticleFile(name=outfilepath, outputdt=2)\n", - " pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", - " endtime=4,\n", - " dt=1,\n", - " output_file=output_file,\n", - " )\n", - " # Note that we also need to write the final time to the file\n", - " output_file.write_latest_locations(pset, 4)\n", - " print(xr.open_zarr(outfilepath).time.values[:].astype(\"timedelta64[s]\"))" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Adding new particles to a ParticleSet during runtime\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the examples above, all particles were defined at the start of the simulation. There are use-cases, though, where it is important to be able to add particles 'on-the-fly', during the runtime of a Parcels simulation.\n", - "\n", - "Unfortuantely, Parcels does not (yet) support adding new particles _within_ a `Kernel`. A workaround is to temporarily leave the `execution()` call to add particles via the `ParticleSet.add()` method, before continuing with `execution()`.\n", - "\n", - "See the example below, where we add 'mass' to a particle each timestep, based on a probablistic condition, and then split the particle once its 'mass' is larger than 5\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "GrowingParticle = parcels.Particle.add_variables(\n", - " [\n", - " parcels.Variable(\"mass\", initial=0),\n", - " parcels.Variable(\"splittime\", initial=-1),\n", - " parcels.Variable(\"splitmass\", initial=0),\n", - " ]\n", - ")\n", - "\n", - "\n", - "def GrowParticles(particle, fieldset, time):\n", - " import random\n", - "\n", - " # 25% chance per timestep for particle to grow\n", - " if random.random() < 0.25:\n", - " particle.mass += 1.0\n", - " if (particle.mass >= 5.0) and (particle.splittime < 0):\n", - " particle.splittime = time\n", - " particle.splitmass = particle.mass / 2.0\n", - " particle.mass = particle.mass / 2.0\n", - "\n", - "\n", - "pset = parcels.ParticleSet(fieldset=fieldset, pclass=GrowingParticle, lon=0, lat=0)\n", - "outfile = parcels.ParticleFile(\"growingparticles.zarr\", pset, outputdt=1)\n", - "\n", - "for t in range(40):\n", - " pset.execute(\n", - " GrowParticles, runtime=1, dt=1, output_file=outfile, verbose_progress=False\n", - " )\n", - " for p in pset:\n", - " if p.splittime > 0:\n", - " pset.add(\n", - " parcels.ParticleSet(\n", - " fieldset=fieldset,\n", - " pclass=GrowingParticle,\n", - " lon=0,\n", - " lat=0,\n", - " time=p.splittime,\n", - " mass=p.splitmass,\n", - " )\n", - " )\n", - " p.splittime = -1 # reset splittime" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The 'trick' is that we place the `pset.execute()` call in a for-loop, so that we leave the Kernel-loop and can add Particles to the ParticleSet.\n", - "\n", - "Indeed, if we plot the mass of particles as a function of time, we see that new particles are created every time a particle reaches a mass of 5.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ds = xr.open_zarr(\"growingparticles.zarr\")\n", - "plt.plot(ds.time.values[:].astype(\"timedelta64[s]\").T, ds.mass.T)\n", - "plt.grid()\n", - "plt.xlabel(\"Time\")\n", - "plt.ylabel(\"Particle 'mass'\")\n", - "plt.show()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "parcels", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.3" - } - }, - "nbformat": 4, - "nbformat_minor": 1 -} diff --git a/docs/examples_v3/tutorial_splitparticles.ipynb b/docs/examples_v3/tutorial_splitparticles.ipynb new file mode 100644 index 000000000..02a0d26c9 --- /dev/null +++ b/docs/examples_v3/tutorial_splitparticles.ipynb @@ -0,0 +1,158 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding new particles to a ParticleSet during runtime\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are use-cases, where it is important to be able to add particles 'on-the-fly', during the runtime of a Parcels simulation.\n", + "\n", + "Unfortuantely, Parcels does not (yet) support adding new particles _within_ a `Kernel`. A workaround is to temporarily leave the `execution()` call to add particles via the `ParticleSet.add()` method, before continuing with `execution()`.\n", + "\n", + "See the example below, where we add 'mass' to a particle each timestep, based on a probablistic condition, and then split the particle once its 'mass' is larger than 5\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from datetime import timedelta\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import xarray as xr\n", + "from IPython.display import HTML\n", + "from matplotlib.animation import FuncAnimation\n", + "\n", + "import parcels" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n", + "example_dataset_folder = parcels.download_example_dataset(\n", + " \"CopernicusMarine_data_for_Argo_tutorial\"\n", + ")\n", + "\n", + "ds = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n", + "ds.load() # load the dataset into memory\n", + "\n", + "fieldset = parcels.FieldSet.from_copernicusmarine(ds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GrowingParticle = parcels.Particle.add_variable([parcels.Variable(\"mass\", initial=0)])\n", + "\n", + "\n", + "def GrowParticles(particles, fieldset):\n", + " # 25% chance per timestep for particle to grow\n", + " growing_particles = np.random.random_sample(len(particles)) < 0.25\n", + " particles[growing_particles].mass += 1.0\n", + "\n", + "\n", + "def SplitParticles(particles, fieldset):\n", + " splitting_particles = particles.mass >= 5.0\n", + " particles[splitting_particles].mass = particles[splitting_particles].mass / 2.0\n", + "\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " pclass=GrowingParticle,\n", + " lon=0,\n", + " lat=0,\n", + " time=fieldset.time_interval.left,\n", + ")\n", + "outfile = parcels.ParticleFile(\"growingparticles.zarr\", outputdt=np.timedelta64(1, \"h\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pset.execute(\n", + " GrowParticles,\n", + " runtime=np.timedelta64(40, \"h\"),\n", + " dt=np.timedelta64(1, \"h\"),\n", + " output_file=outfile,\n", + " verbose_progress=False,\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The 'trick' is that we place the `pset.execute()` call in a for-loop, so that we leave the Kernel-loop and can add Particles to the ParticleSet.\n", + "\n", + "Indeed, if we plot the mass of particles as a function of time, we see that new particles are created every time a particle reaches a mass of 5.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds_out = xr.open_zarr(\"growingparticles.zarr\")\n", + "plt.plot(\n", + " (ds_out.time.values[:].T - ds.time.values[0]).astype(\"timedelta64[h]\"),\n", + " ds_out.mass.T,\n", + ")\n", + "plt.grid()\n", + "plt.xlabel(\"Time [h]\")\n", + "plt.ylabel(\"Particle 'mass'\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "test-notebooks", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}