.ipynb .pdf

repository open issue suggest edit

Binder

Contents

Particles tracking

The OceanParcels project develops Parcels (Probably A Really Computationally Efficient Lagrangian Simulator), a set of Python classes and methods to create customisable particle tracking simulations using output from Ocean Circulation models such as the eReefs one.

Note

Parcels can be used to track passive and active particulates such as:

• water,

• plankton,

• plastic and

• fish.

In this notebook, we will first cover how to run a set of particles from exported eReefs data. Then we will show how to use particles to sample a field such as temperature and how to write a kernel that tracks the distance travelled by the particles.

See also

The OceanParcels team has designed a series of Jupyter notebooks to help people to get started with Parcels. You can find the tutorials on this link!

Load the required Python libraries

First of all, load the necessary libraries:

import os
import numpy as np
import xarray as xr

import cmocean

import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
cartopy.config['data_dir'] = os.getenv('CARTOPY_DIR', cartopy.config.get('data_dir'))

from parcels import FieldSet, Field, ParticleSet, Variable, JITParticle
from parcels import AdvectionRK4, plotTrajectoriesFile, ErrorCode

import math
from datetime import timedelta
from operator import attrgetter

from matplotlib import pyplot as plt
#%config InlineBackend.figure_format = 'retina'
plt.ion()  # To trigger the interactive inline mode

INFO: Compiled ParcelsRandom ==> /tmp/parcels-1001/libparcels_random_4cae6ef5-3eba-46fd-ab4f-692a1b6b577c.so

Build multi-file dataset

We will use the open_mfdataset function from xArray to open multiple netCDF files into a single xarray Dataset.

We will query load the GBR4km dataset from the AIMS server, so let’s first define the base URL:

# For the hydro dataset
base_url = "http://thredds.ereefs.aims.gov.au/thredds/dodsC/s3://aims-ereefs-public-prod/derived/ncaggregate/ereefs/gbr4_v2/daily-monthly/EREEFS_AIMS-CSIRO_gbr4_v2_hydro_daily-monthly-"

For the sake of the demonstration, we will only use 1 month:

month_st = 1   # Starting month 
month_ed = 1   # Ending month 
year = 2018    # Year

# Based on the server the file naming convention 
hydrofiles = [f"{base_url}{year}-{month:02}.nc" for month in range(month_st, month_ed+1)]

Loading dataset into xArray

Using xArray, we open these files into a Dataset:

ds_hydro = xr.open_mfdataset(hydrofiles)

Clip the Dataset

To reduce the Dataset size we will clip the spatial extent based on longitudinal and latitudinal values.

This is easely done using the sel function with the slice method.

min_lon = 146     # lower left longitude
min_lat = -21     # lower left latitude
max_lon = 150     # upper right longitude
max_lat = -16     # upper right latitude

# Defining the boundaries
lon_bnds = [min_lon, max_lon]
lat_bnds = [min_lat, max_lat]

# Performing the reduction and only taking the surface dataset (k=-1)
ds_hydro_clip = ds_hydro.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds), k=-1)
ds_hydro_clip

<xarray.Dataset>
Dimensions:      (latitude: 167, longitude: 134, time: 31)
Coordinates:
    zc           float64 dask.array<chunksize=(), meta=np.ndarray>
  * time         (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:...
  * latitude     (latitude) float64 -20.99 -20.96 -20.93 ... -16.04 -16.01
  * longitude    (longitude) float64 146.0 146.0 146.1 ... 149.9 150.0 150.0
Data variables:
    mean_cur     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    salt         (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    temp         (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    u            (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    v            (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    mean_wspeed  (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    eta          (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_u     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    wspeed_v     (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time

We will now drop all unnecessary variables.

Basically we will only need the current velocities (u and v) and the variable we want to track with the parcels (here temp):

surf_data = ds_hydro_clip.drop(['zc','mean_wspeed','salt',
                    'eta','wspeed_u','wspeed_v', 
                    'mean_cur'])

surf_data

<xarray.Dataset>
Dimensions:    (latitude: 167, longitude: 134, time: 31)
Coordinates:
  * time       (time) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30T14:00:00
  * latitude   (latitude) float64 -20.99 -20.96 -20.93 ... -16.07 -16.04 -16.01
  * longitude  (longitude) float64 146.0 146.0 146.1 146.1 ... 149.9 150.0 150.0
Data variables:
    temp       (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    u          (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
    v          (time, latitude, longitude) float32 dask.array<chunksize=(31, 167, 134), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                     CF-1.0
    Run_ID:                          2
    _CoordSysBuilder:                ucar.nc2.dataset.conv.CF1Convention
    aims_ncaggregate_buildDate:      2020-08-21T12:50:07+10:00
    aims_ncaggregate_datasetId:      products__ncaggregate__ereefs__gbr4_v2__...
    aims_ncaggregate_firstDate:      2018-01-01T00:00:00+10:00
    ...                              ...
    paramhead:                       GBR 4km resolution grid
    shoc_version:                    v1.1 rev(5620)
    technical_guide_link:            https://eatlas.org.au/pydio/public/aims-...
    technical_guide_publish_date:    2020-08-18
    title:                           eReefs AIMS-CSIRO GBR4 Hydrodynamic v2 d...
    DODS_EXTRA.Unlimited_Dimension:  time

OceanParcels will need to read the file locally and we save it using the Xarray funciton to_netcdf. Otherwise you could use the from_xarray_dataset function to directly read the dataset in Parcels.

Note

Here the cell has been commented to make the notebook run faster, but you will need to uncomment it for running your how experiment…

data_name ='ereefdata.nc'
# try:
#     os.remove(data_name)
# except OSError:
#     pass

# surf_data.to_netcdf(path=data_name)

Reading velocities into Parcels

As we used the NetCDF format, it is fairly easy to load the velocity fields into the FieldSet.from_netcdf() function available in Parcels:

filenames = {'U': data_name,
             'V': data_name}

Then, define a dictionary of the variables (U and V) and dimensions (lon, lat and time; note that in this case there is no depth because we only took the surface variables from the eReefs datastet:

variables = {'U': 'u',
             'V': 'v'}

dimensions = {'lat': 'latitude',
              'lon': 'longitude',
              'time': 'time'}

Finally, we read in the fieldset using the FieldSet.from_netcdf function with the above-defined filenames, variables and dimensions

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 define a ParticleSet, in this case with 5 particles starting on a line between (147E, 18.5S) and (148E, 17.5S) using the ParticleSet.from_line constructor method:

pset = ParticleSet.from_line(fieldset=fieldset, pclass=JITParticle,
                             size=5,               # releasing 5 particles
                             start=(147, -18.5),   # releasing on a line: the start longitude and latitude
                             finish=(148, -17.5))  # releasing on a line: the end longitude and latitude

Note

Another approach if you want to list a series of initial positions consists in using the from_list function which takes:

• lon – List of initial longitude values for particles

• lat – List of initial latitude values for particles

Now we want to advect the particles. However, the eReefs 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:

def DeleteParticle(particle, fieldset, time):
    particle.delete()

Now we can advect the particles by executing the ParticleSet for 30 days using 4th order Runge-Kutta (AdvectionRK4):

output_nc = 'CurrentParticles.nc'
try:
    os.remove(output_nc)
except OSError:
    pass

output_file = pset.ParticleFile(name=output_nc, 
                                outputdt=timedelta(hours=6))

pset.execute(AdvectionRK4,
             runtime=timedelta(days=30),
             dt=timedelta(minutes=5),
             output_file=output_file,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-1001/0afc2149946c974dbfdf33957dcb4b68_0.so

We now export the particles recorded every 6h in our output file:

output_file.export()

Plotting outputs

We open the particles file with Xarray:

parcels = xr.open_dataset(output_nc)
parcels

<xarray.Dataset>
Dimensions:     (obs: 121, traj: 5)
Dimensions without coordinates: obs, traj
Data variables:
    trajectory  (traj, obs) float64 0.0 0.0 0.0 0.0 0.0 ... 4.0 4.0 4.0 4.0 4.0
    time        (traj, obs) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30...
    lat         (traj, obs) float32 -18.5 -18.54 -18.57 ... -18.79 -18.8 -18.82
    lon         (traj, obs) float32 147.0 147.0 147.1 ... 146.8 146.8 146.9
    z           (traj, obs) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
    feature_type:           trajectory
    Conventions:            CF-1.6/CF-1.7
    ncei_template_version:  NCEI_NetCDF_Trajectory_Template_v2.0
    parcels_version:        2.2.2
    parcels_mesh:           spherical

We will now use the Xarray plot function as we have done in the past…

See also

Check out this link to see how to use Xarray’s convenient matplotlib-backed plotting interface to visualize your datasets!

# Figure size
size = (9, 10)

# Color from cmocean
color = cmocean.cm.speed

# Defining the figure
fig = plt.figure(figsize=size, facecolor='w', edgecolor='k')

# Axes with Cartopy projection
ax = plt.axes(projection=ccrs.PlateCarree())
# and extent
ax.set_extent([min_lon, max_lon, min_lat, max_lat], ccrs.PlateCarree())

# Plotting using Matplotlib 
# We plot the PH at the surface at the final recorded time interval
cf = ds_hydro_clip.mean_cur.isel(time=10).plot( 
    transform=ccrs.PlateCarree(), cmap=color,
    vmin = 0.1, vmax = 1.0, alpha=0.2, 
    add_colorbar=False
)

# Color bar
cbar = fig.colorbar(cf, ax=ax, fraction=0.027, pad=0.045, 
                    orientation="horizontal")
cbar.set_label(ds_hydro_clip.mean_cur.long_name+' '+ds_hydro_clip.mean_cur.units, rotation=0, 
               labelpad=5, fontsize=10)
cbar.ax.tick_params(labelsize=8)

# Title
plt.title('Parcels evolution',
          fontsize=13
         )

# Plot lat/lon grid 
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=0.1, color='k', alpha=1, 
                  linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 8}
gl.ylabel_style = {'size': 8} 

# Add map features with Cartopy 
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '10m', 
                                            edgecolor='face', 
                                            facecolor='lightgray'))
ax.coastlines(linewidth=1)

for k in range(parcels.lon.shape[0]):
    ax.scatter(parcels.lon.isel(traj=k), parcels.lat.isel(traj=k), s=40, edgecolors='w', 
               linewidth=0.2, transform=ccrs.PlateCarree()).set_zorder(11)

plt.tight_layout()
plt.show()
fig.clear()
plt.close(fig)
plt.clf()

<Figure size 432x288 with 0 Axes>

Sampling a Field with Particles

One typical use case of particle simulations is to sample a Field (such as temperature, salinity, surface hight) along a particle trajectory. In Parcels, this is very easy to do, with a custom Kernel.

Let’s define the FieldSet as above:

filenames = {'U': data_name,
             'V': data_name,
            }

variables = {'U': 'u',
             'V': 'v',
            }

dimensions = {'lat': 'latitude',
              'lon': 'longitude',
              'time': 'time'}

Finally, we read in the fieldset using the FieldSet.from_netcdf function with the above-defined filenames, variables and dimensions

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)

We now add the field temperature to the FieldSet

filenames = {'lat': data_name,
             'lon': data_name,
             'data': data_name}

variable = ('T', 'temp')

dimensions = {'lat': 'latitude',
              'lon': 'longitude',
              'time': 'time'}

field = Field.from_netcdf(filenames, variable, dimensions)
fieldset.add_field(field)

Now define a new Particle class that has an extra Variable: the temperature. We initialise this by sampling the fieldset2.T field.

class SampleParticle(JITParticle):          # Define a new particle class
    t = Variable('t', initial=fieldset.T)  # Variable 't' initialised by sampling the temperature

Now define a ParticleSet using the from_line method also used above. Plot the pset and print their temperature values t:

pset = ParticleSet.from_line(fieldset=fieldset, pclass=SampleParticle, 
                             size=5,               # releasing 5 particles
                             start=(147, -18.5),   # releasing on a line: the start longitude and latitude
                             finish=(148, -17.5),  # releasing on a line: the end longitude and latitude
                             time=0)
print('t values before execution:', [p.t for p in pset])

WARNING: Particle initialisation from field can be very slow as it is computed in scipy mode.

t values before execution: [28.90722, 29.133469, 29.18559, 29.126036, 29.002943]

Now create a custom function that samples the fieldset.T field at the particle location. Cast this function to a Kernel.

def SampleT(particle, fieldset, time):  # Custom function that samples fieldset2.T at particle location
    particle.t = fieldset.T[time, particle.depth, particle.lat, particle.lon]

k_sample = pset.Kernel(SampleT)    # Casting the SampleT function to a kernel.

Finally, execute the pset with a combination of the AdvectionRK4 and SampleT kernels, plot the pset and print their new temperature values t:

output_nc_temp = 'CurrentParticlesTemp.nc'
try:
    os.remove(output_nc_temp)
except OSError:
    pass

output_file_temp = pset.ParticleFile(name=output_nc_temp, 
                                outputdt=timedelta(hours=6))

pset.execute(AdvectionRK4 + k_sample,    # Add kernels using the + operator.
             runtime=timedelta(days=30),
             dt=timedelta(minutes=5),
             output_file=output_file_temp,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

output_file_temp.export()

INFO: Compiled SampleParticleAdvectionRK4SampleT ==> /tmp/parcels-1001/10e9ccb780cd6a66e3739935a7c2464b_0.so

Plotting the result

We can extract the netcdf file values…

parcels_temp = xr.open_dataset(output_nc_temp)
parcels_temp

<xarray.Dataset>
Dimensions:     (obs: 121, traj: 5)
Dimensions without coordinates: obs, traj
Data variables:
    trajectory  (traj, obs) float64 5.0 5.0 5.0 5.0 5.0 ... 9.0 9.0 9.0 9.0 9.0
    time        (traj, obs) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30...
    lat         (traj, obs) float32 -18.5 -18.54 -18.57 ... -18.79 -18.8 -18.82
    lon         (traj, obs) float32 147.0 147.0 147.1 ... 146.8 146.8 146.9
    z           (traj, obs) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    t           (traj, obs) float32 28.91 28.92 28.91 ... 28.85 28.95 29.04
Attributes:
    feature_type:           trajectory
    Conventions:            CF-1.6/CF-1.7
    ncei_template_version:  NCEI_NetCDF_Trajectory_Template_v2.0
    parcels_version:        2.2.2
    parcels_mesh:           spherical

And plot the resulting parcels output on a map:

# Figure size
size = (9, 10)

# Color from cmocean
color = cmocean.cm.speed

# Defining the figure
fig = plt.figure(figsize=size, facecolor='w', edgecolor='k')

# Axes with Cartopy projection
ax = plt.axes(projection=ccrs.PlateCarree())
# and extent
ax.set_extent([min_lon, max_lon, min_lat, max_lat], ccrs.PlateCarree())

# Plotting using Matplotlib 
# We plot the PH at the surface at the final recorded time interval
cf = ds_hydro_clip.mean_cur.isel(time=10).plot( 
    transform=ccrs.PlateCarree(), cmap=color,
    vmin = 0.1, vmax = 1.0, alpha=0.2, 
    add_colorbar=False
)

# Color bar
cbar = fig.colorbar(cf, ax=ax, fraction=0.027, pad=0.045, 
                    orientation="horizontal")
cbar.set_label(ds_hydro_clip.mean_cur.long_name+' '+ds_hydro_clip.mean_cur.units, rotation=0, 
               labelpad=5, fontsize=10)
cbar.ax.tick_params(labelsize=8)

# Title
plt.title('Parcels evolution coloured by surface temperature',
          fontsize=13
         )

# Plot lat/lon grid 
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=0.1, color='k', alpha=1, 
                  linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 8}
gl.ylabel_style = {'size': 8} 

# Add map features with Cartopy 
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '10m', 
                                            edgecolor='face', 
                                            facecolor='lightgray'))
ax.coastlines(linewidth=1)

distmin = parcels_temp.t.min().item()
distmax = parcels_temp.t.max().item()

for k in range(parcels_temp.lon.shape[0]):
    sc = plt.scatter(parcels_temp.lon.isel(traj=k), parcels_temp.lat.isel(traj=k), s=40, 
               c=parcels_temp.t.isel(traj=k), edgecolors='k', 
               cmap=cmocean.cm.balance, vmin=distmin, vmax=distmax, 
               linewidth=0.2, transform=ccrs.PlateCarree()).set_zorder(11)

# Color bar
cbar2 = plt.colorbar(sc, ax=ax, fraction=0.027, pad=0.045)
cbar2.set_label('Surface temperature in degree C', rotation=90, 
               labelpad=5, fontsize=10)
cbar2.ax.tick_params(labelsize=8)

plt.tight_layout()
plt.show()
fig.clear()
plt.close(fig)
plt.clf()

<Figure size 432x288 with 0 Axes>

Calculating distance travelled

As a second example of what custom kernels can do, we will now show how to create a kernel that logs the total distance that particles have travelled.

First, we need to create a new Particle class that includes three extra variables. The distance variable will be written to output, but the auxiliary variables prev_lon and prev_lat won’t be written to output (can be controlled using the to_write keyword)

class DistParticle(JITParticle):  # Define a new particle class that contains three extra variables
    
    distance = Variable('distance', initial=0., dtype=np.float32)  # the distance travelled
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
                        initial=attrgetter('lon'))  # the previous longitude
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
                        initial=attrgetter('lat'))  # the previous latitude.

Now define a new function TotalDistance that calculates the sum of Euclidean distances between the old and new locations in each RK4 step:

def TotalDistance(particle, fieldset, time):
    
    # Calculate the distance in latitudinal direction (using 1.11e2 kilometer per degree latitude)
    lat_dist = (particle.lat - particle.prev_lat) * 1.11e2
    # Calculate the distance in longitudinal direction, using cosine(latitude) - spherical earth
    lon_dist = (particle.lon - particle.prev_lon) * 1.11e2 * math.cos(particle.lat * math.pi / 180)
    # Calculate the total Euclidean distance travelled by the particle
    particle.distance += math.sqrt(math.pow(lon_dist, 2) + math.pow(lat_dist, 2))

    particle.prev_lon = particle.lon  # Set the stored values for next iteration.
    particle.prev_lat = particle.lat

Note

Here it is assumed that the latitude and longitude are measured in degrees North and East, respectively. However, some datasets give them measured in (kilo)meters, in which case we must not include the factor 1.11e2.

We will run the TotalDistance function on a ParticleSet containing the five particles within the eReefs fieldset similar to the one we did above.

Note that pclass=DistParticle in this case.

filenames = {'U': data_name,
             'V': data_name,
            }

variables = {'U': 'u',
             'V': 'v'}

dimensions = {'lat': 'latitude',
              'lon': 'longitude',
              'time': 'time'}

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)

pset = ParticleSet.from_line(fieldset=fieldset, 
                             pclass=DistParticle,
                             size=5,               # releasing 5 particles
                             start=(147, -18.5),   # releasing on a line: the start longitude and latitude
                             finish=(148, -17.5),  # releasing on a line: the end longitude and latitude
                             )

Again we define a new kernel to include the function written above and execute the ParticleSet.

output_nc_dist = 'CurrentParticlesDist.nc'
try:
    os.remove(output_nc_dist)
except OSError:
    pass

file_dist = pset.ParticleFile(name=output_nc_dist, 
                                outputdt=timedelta(hours=1))

k_dist = pset.Kernel(TotalDistance)  # Casting the TotalDistance function to a kernel.

pset.execute(AdvectionRK4 + k_dist,  # Add kernels using the + operator.
             runtime=timedelta(days=30),
             dt=timedelta(minutes=5),
             output_file=file_dist,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

file_dist.export()

INFO: Compiled DistParticleAdvectionRK4TotalDistance ==> /tmp/parcels-1001/db33ec2ecce7d6d05221b0c0912c7346_0.so

We can now print the distance in km that each particle has travelled:

print([p.distance for p in pset]) # the distances in km travelled by the particles

[278.97974, 459.3084, 507.74344, 508.9902]

Plotting the result

We can extract the netcdf file values…

parcels_dist = xr.open_dataset(output_nc_dist)
parcels_dist

<xarray.Dataset>
Dimensions:     (obs: 721, traj: 5)
Dimensions without coordinates: obs, traj
Data variables:
    trajectory  (traj, obs) float64 10.0 10.0 10.0 10.0 ... 14.0 14.0 14.0 14.0
    time        (traj, obs) datetime64[ns] 2017-12-31T14:00:00 ... 2018-01-30...
    lat         (traj, obs) float32 -18.5 -18.51 -18.51 ... -18.81 -18.81 -18.82
    lon         (traj, obs) float32 147.0 147.0 147.0 ... 146.8 146.8 146.9
    z           (traj, obs) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    distance    (traj, obs) float32 0.0 0.9391 1.884 2.838 ... 508.3 508.6 509.0
Attributes:
    feature_type:           trajectory
    Conventions:            CF-1.6/CF-1.7
    ncei_template_version:  NCEI_NetCDF_Trajectory_Template_v2.0
    parcels_version:        2.2.2
    parcels_mesh:           spherical

And plot the resulting parcels output on a map:

# Figure size
size = (9, 10)

# Color from cmocean
color = cmocean.cm.speed

# Defining the figure
fig = plt.figure(figsize=size, facecolor='w', edgecolor='k')

# Axes with Cartopy projection
ax = plt.axes(projection=ccrs.PlateCarree())
# and extent
ax.set_extent([min_lon, max_lon, min_lat, max_lat], ccrs.PlateCarree())

# Plotting using Matplotlib 
# We plot the PH at the surface at the final recorded time interval
cf = ds_hydro_clip.mean_cur.isel(time=10).plot( 
    transform=ccrs.PlateCarree(), cmap=color,
    vmin = 0.1, vmax = 1.0, alpha=0.2, 
    add_colorbar=False
)

# Color bar
cbar = fig.colorbar(cf, ax=ax, fraction=0.027, pad=0.045, 
                    orientation="horizontal")
cbar.set_label(ds_hydro_clip.mean_cur.long_name+' '+ds_hydro_clip.mean_cur.units, rotation=0, 
               labelpad=5, fontsize=10)
cbar.ax.tick_params(labelsize=8)

# Title
plt.title('Parcels evolution coloured by travelled distance',
          fontsize=13
         )

# Plot lat/lon grid 
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=0.1, color='k', alpha=1, 
                  linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 8}
gl.ylabel_style = {'size': 8} 

# Add map features with Cartopy 
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '10m', 
                                            edgecolor='face', 
                                            facecolor='lightgray'))
ax.coastlines(linewidth=1)

distmin = parcels_dist.distance.min().item()
distmax = parcels_dist.distance.max().item()

for k in range(parcels_dist.lon.shape[0]):
    sc = plt.scatter(parcels_dist.lon.isel(traj=k), parcels_dist.lat.isel(traj=k), s=40, 
               c=parcels_dist.distance.isel(traj=k), edgecolors='w', 
               cmap='jet', vmin=distmin, vmax=distmax, 
               linewidth=0.2, transform=ccrs.PlateCarree()).set_zorder(11)

# Color bar
cbar2 = plt.colorbar(sc, ax=ax, fraction=0.027, pad=0.045)
cbar2.set_label('Travelled distance in km', rotation=90, 
               labelpad=5, fontsize=10)
cbar2.ax.tick_params(labelsize=8)

plt.tight_layout()
plt.show()
fig.clear()
plt.close(fig)
plt.clf()

<Figure size 432x288 with 0 Axes>

See also

Tutorial on how to analyse Parcels output

PyViz Application Other libraries of interest