parcels/examples/tutorial_Agulhasparticles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tutorial showing how to create Parcels in Agulhas animated gif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This brief tutorial shows how to recreate the [animated gif](http://oceanparcels.org/animated-gifs/globcurrent_fullyseeded.gif) showing particles in the Agulhas region south of Africa."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with importing the relevant modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4, ErrorCode\n",
    "from datetime import timedelta\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now load the Globcurrent fields from the `GlobCurrent_example_data` directory (note that unlike in the main Parcels tutorial we don't use a dictionary for the filenames here; as they are the same for all variables, we don't need to)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Casting lon data to np.float32\n",
      "WARNING: Casting lat data to np.float32\n",
      "WARNING: Casting depth data to np.float32\n"
     ]
    }
   ],
   "source": [
    "filenames = \"GlobCurrent_example_data/20*.nc\"\n",
    "variables = {'U': 'eastward_eulerian_current_velocity',\n",
    "             'V': 'northward_eulerian_current_velocity'}\n",
    "dimensions = {'lat': 'lat',\n",
    "              'lon': 'lon',\n",
    "              'time': 'time'}\n",
    "fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create vectors of Longitude and Latitude starting locations on a regular mesh, and use these to initialise a `ParticleSet` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lons, lats = np.meshgrid(range(15, 35), range(-40, -30))\n",
    "pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=lons, lat=lats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to advect the particles. However, the Globcurrent data that we loaded in is only for a limited, regional domain and particles might be able to leave this domain. We therefore need to tell Parcels that particles that leave the domain need to be deleted. We do that using a `Recovery Kernel`, which will be invoked when a particle encounters an `ErrorOutOfBounds` error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DeleteParticle(particle, fieldset, time):\n",
    "    particle.delete()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can advect the particles. Note that we do this inside a `for`-loop, so we can save a plot every six hours (which is the value of `runtime`). See the [plotting tutorial](http://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/examples/tutorial_plotting.ipynb) for more information on the `pset.show()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/erik/miniconda2/envs/py2_parcels_cartopy/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py:1744: RuntimeWarning: invalid value encountered in less\n",
      "  u, v = self.projection.transform_vectors(t, x, y, u, v)\n",
      "/Users/erik/miniconda2/envs/py2_parcels_cartopy/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py:1744: RuntimeWarning: invalid value encountered in greater\n",
      "  u, v = self.projection.transform_vectors(t, x, y, u, v)\n",
      "INFO: Plot saved to particles00.png\n",
      "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gr/T/parcels-504/de5d5c9b14f1551d894b622aa4ca7898.so\n",
      "100% (21600.0 of 21600.0) |##############| Elapsed Time: 0:00:00 Time:  0:00:00\n",
      "INFO: Plot saved to particles01.png\n",
      "100% (21600.0 of 21600.0) |##############| Elapsed Time: 0:00:00 Time:  0:00:00\n",
      "INFO: Plot saved to particles02.png\n",
      "100% (21600.0 of 21600.0) |##############| Elapsed Time: 0:00:00 Time:  0:00:00\n"
     ]
    }
   ],
   "source": [
    "for cnt in range(3):\n",
    "    # First plot the particles\n",
    "    pset.show(savefile='particles'+str(cnt).zfill(2), field='vector', land=True, vmax=2.0)\n",
    "\n",
    "    # Then advect the particles for 6 hours\n",
    "    pset.execute(AdvectionRK4,\n",
    "                 runtime=timedelta(hours=6),  # runtime controls the interval of the plots\n",
    "                 dt=timedelta(minutes=5),\n",
    "                 recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})  # the recovery kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This now has created 3 plots. Note that the original animated gif contained 20 plots, but to keep running of this notebook fast we have reduced the number here. Of course, it is trivial to increase the number of plots by changing the value in the `range()` in the cell above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a final step, you can use [ImageMagick](http://www.imagemagick.org/script/index.php) or an online tool to stitch these individual plots together in an animated gif."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}