Tutorial showing how to create Parcels in Agulhas animated gif

This brief tutorial shows how to recreate the animated gif showing particles in the Agulhas region south of Africa.

We start with importing the relevant modules

In [1]:
from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4, ErrorCode
from datetime import timedelta
import numpy as np

Now load the Globcurrent fields from the GlobCurrent_example_data directory

In [2]:
filenames = {'U': "GlobCurrent_example_data/20*.nc",
             'V': "GlobCurrent_example_data/20*.nc"}
variables = {'U': 'eastward_eulerian_current_velocity',
             'V': 'northward_eulerian_current_velocity'}
dimensions = {'lat': 'lat',
              'lon': 'lon',
              'time': 'time'}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
WARNING: Casting lon data to np.float32
WARNING: Casting lat data to np.float32
WARNING: Casting depth data to np.float32

Now create vectors of Longitude and Latitude starting locations on a regular mesh, and use these to initialise a ParticleSet object.

In [3]:
lons, lats = np.meshgrid(range(15, 35), range(-40, -30))
pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=lons, lat=lats)

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:

In [4]:
def DeleteParticle(particle, fieldset, time, dt):

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 for more information on the pset.show() method.

In [5]:
for cnt in range(3):
    # First plot the particles
    pset.show(savefile='particles'+str(cnt).zfill(2), field='vector', land=True, vmax=2.0)

    # Then advect the particles for 6 hours
                 starttime=pset[0].time,  # note that we have to set the starttime to the particle time
                 runtime=timedelta(hours=6),  # runtime controls the interval of the plots
                 recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})  # the recovery kernel
INFO: Plot saved to particles00.png
INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gr/T/parcels-504/27805ff3aa34ba12ddb373f3f2cb1d1b.so
INFO: Plot saved to particles01.png
INFO: Plot saved to particles02.png

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.

As a final step, you can use ImageMagick or an online tool to stitch these individual plots together in an animated gif.