functions.py

# A bunch of helpful functions for the parcels simulations

import parcels
from parcels import FieldSet, StatusCode
from parcels.tools.converters import Geographic, GeographicPolar
from datetime import timedelta
import os
import pandas as pd
import numpy as np
from glob import glob

import xarray as xr
import xgcm as xgcm
import xnemogcm as xn


def check_conda_environment(loaded_env, project_env):
    """ Function to ensure the defined conda environment is in use. """
    if loaded_env != project_env:
        raise ValueError(f"The incorrect conda environment is loaded: {loaded_env}. Please use the correct conda environment: {project_env}")


def create_dataset(derivatives, time_counter):
    """ Function to create the derivatives dataset. """
    ds = xr.Dataset(derivatives)
    ds["time_counter"] = time_counter
    ds = ds.assign_coords({"time_counter": time_counter})
    ds.attrs["Creation date"] = pd.Timestamp.now().strftime(
        "%Y-%m-%d %H:%M:%S")

    return ds


def export_dataset(ds, timestamp, outputdir, vars=['dudx', 'dudy', 'dvdx', 'dvdy']):
    """ Function to export the derivatives dataset to a netCDF file. """
    ds.to_netcdf(f"{outputdir}psy4v3r1_Lofoten_velocity_derivatives_{str(timestamp)}.nc",
                 encoding={var: {"zlib": True, "complevel": 1} for var in vars})


# Fieldset creation helper functions
def select_files(dirread, string_, date_current, dt_sim, i_date_s=-13, i_date_e=-3, dt_margin=8):
    # set dt_margin to e.g. 32 when dealing with monthly data, or 8 when dealing with weekly data
    yr0 = date_current.year

    time_start = date_current - timedelta(days=dt_margin)
    time_end = date_current + timedelta(days=dt_sim) + timedelta(days=dt_margin)
    yrEnd = time_end.year

    ragged_files = [glob(dirread + string_ % (yr0+i)) for i in range(-1, yrEnd-yr0+1)]
    possible_files = []

    for i in range(len(ragged_files)):
        for file_ in ragged_files[i]:
            possible_files.append(file_)

    possible_files = sorted(possible_files)

    indices_use = []
    for i1, file_ in enumerate(possible_files):
        str_date = os.path.basename(file_)[i_date_s:i_date_e]
        date_ = pd.Timestamp(str_date)

        if date_.date() > time_start and date_.date() < time_end:
            indices_use.append(i1)

    files_use = list(possible_files[i] for i in indices_use)
    return files_use



def create_full_fieldset(settings):
    """ A constructor method to create a Parcels.Fieldset from hydrodynamic model data

    Parameters
    ----------
    settings :
        A dictionary of settings used to create the fieldset

    Returns
    -------
    fieldset
        A parcels.FieldSet object
    """

    # Location of hydrodynamic data
    dirread_model = os.path.join(settings['ocean']['directory'], settings['ocean']['filename_style'])

    # Start date and runtime of the simulation
    startdate = settings['simulation']['startdate']
    runtime = int(np.ceil(settings['simulation']['runtime'].total_seconds()/86400.))  # convert to days

    # Mesh masks
    ocean_mesh = os.path.join(settings['ocean']['directory'], settings['ocean']['ocean_mesh'])  # mesh_mask

    # Setup input for fieldset creation
    ufiles = select_files(dirread_model, 'U_%4i*.nc', startdate, runtime, dt_margin=3)
    vfiles = select_files(dirread_model, 'V_%4i*.nc', startdate, runtime, dt_margin=3)
    wfiles = select_files(dirread_model, 'W_%4i*.nc', startdate, runtime, dt_margin=3)
    tfiles = select_files(dirread_model, 'T_%4i*.nc', startdate, runtime, dt_margin=3)
    sfiles = select_files(dirread_model, 'S_%4i*.nc', startdate, runtime, dt_margin=3)
    kzfiles = select_files(dirread_model, 'KZ_%4i*.nc', startdate, runtime, dt_margin=3)

    #derivs_files = select_files("/storage/shared/oceanparcels/output_data/data_Michael/EddyPlasticEntrainment/data/psy4v3r1/derivatives/psy4v3r1_Lofoten_velocity_",
    #                            'derivatives_%4i*.nc', startdate, runtime, dt_margin=3)

    #derivs_files = sorted(glob("/storage/shared/oceanparcels/output_data/data_Michael/EddyPlasticEntrainment/data/psy4v3r1/derivatives/psy4v3r1_Lofoten_velocity_derivatives_2019-*.nc"))
    #derivs_mesh = "/storage/shared/oceanparcels/output_data/data_Michael/EddyPlasticEntrainment/data/psy4v3r1/derivatives/PSY4V3R1_Lofoten_mesh.nc"
    filenames = {'U': {'lon': ocean_mesh, 'lat': ocean_mesh, 'depth': wfiles[0], 'data': ufiles},
                 'V': {'lon': ocean_mesh, 'lat': ocean_mesh, 'depth': wfiles[0], 'data': vfiles},
                 'W': {'lon': ocean_mesh, 'lat': ocean_mesh, 'depth': wfiles[0], 'data': wfiles},
                 'conservative_temperature': {'lon': ocean_mesh, 'lat': ocean_mesh, 'depth': wfiles[0], 'data': tfiles},
                 'absolute_salinity': {'lon': ocean_mesh, 'lat': ocean_mesh, 'depth': wfiles[0], 'data': sfiles},
                 'mixing_kz' : {'lon': ocean_mesh, 'lat': ocean_mesh, 'depth': wfiles[0], 'data': kzfiles}}
    variables = settings['ocean']['variables']
    dimensions = settings['ocean']['dimensions']
    indices = settings['ocean']['indices']

    # Load the fieldset
    fieldset = FieldSet.from_nemo(filenames, variables, dimensions,
                                  indices=indices, allow_time_extrapolation=settings['allow_time_extrapolation'])

    # Add bathymetry
    fieldset.add_constant('z_start', 0.)
    bathymetry_variables = settings['ocean']['bathymetry_variables']
    bathymetry_dimensions = settings['ocean']['bathymetry_dimensions']
    bathymetry_mesh = os.path.join(settings['ocean']['directory'], settings['ocean']['bathymetry_mesh'])
    bathymetry_field = parcels.Field.from_netcdf(bathymetry_mesh, bathymetry_variables, bathymetry_dimensions, indices=indices)
    fieldset.add_field(bathymetry_field)

    return fieldset



def create_full_fieldset_copernicus(settings):
    """ A constructor method to create a Parcels.Fieldset from hydrodynamic model data
    downloaded from Copernicus Marine Store

    Parameters
    ----------
    settings :
        A dictionary of settings used to create the fieldset

    Returns
    -------
    fieldset
        A parcels.FieldSet object
    """
     # Location of hydrodynamic data
    dirread_model = settings['ocean']['directory']

    # Start date and runtime of the simulation
    startdate = settings['simulation']['startdate']
    runtime = int(np.ceil(settings['simulation']['runtime'].total_seconds()/86400.))  # convert to days

    filenames = select_files(dirread_model, settings['ocean']['filename_style'], startdate.to_pydatetime().date(), runtime, dt_margin=3)
    variables = settings['ocean']['variables']
    dimensions = settings['ocean']['dimensions']
    
    # Load the fieldset
    fieldset = parcels.FieldSet.from_netcdf(filenames, variables, dimensions,
                                            allow_time_extrapolation=settings['allow_time_extrapolation'],
                                            interp_method={
                                                "U": "freeslip",
                                                "V": "freeslip",
                                                'conservative_temperature': "linear_invdist_land_tracer",
                                                "absolute_salinity": "linear_invdist_land_tracer"}
                                            )


    # Add bathymetry
    bathy_filename = dirread_model + settings['ocean']['bathymetry_mesh']
    bathy_variables = settings['ocean']['bathymetry_variables']
    bathy_dimensions = settings['ocean']['bathymetry_dimensions']
    bathy_fieldset = parcels.FieldSet.from_netcdf(bathy_filename, bathy_variables, bathy_dimensions)
    fieldset.add_field(bathy_fieldset.bathymetry)


    # Add constants to the fieldset
    fieldset.add_constant('use_mixing', settings['use_mixing'])
    fieldset.add_constant('use_biofouling', settings['use_biofouling'])
    fieldset.add_constant('use_stokes', settings['use_stokes'])
    fieldset.add_constant('use_wind', settings['use_wind'])
    fieldset.add_constant('G', 9.81)  # Gravitational constant [m s-1]
    fieldset.add_constant('use_3D', settings['use_3D'])
    
    fieldset.add_periodic_halo(zonal=True)

    return fieldset


def add_additional_fields(fieldset, settings):
    ### Based on create_fieldset from plasticparcels
    """A constructor method to create a `Parcels.Fieldset` with all fields
    necessary for a plasticparcels simulation (e.g., a hydrodynamic model
    velocity field, biogeochemical model variable fields, a wind field, etc.).

    Parameters
    ----------
    settings :
        A dictionary of settings that contains an ocean model information,
        biogeochemical model information, wind model information,
        and other optional settings.

    Returns
    -------
    fieldset
        A `parcels.FieldSet` object.
    """
    # First create the hydrodynamic fieldset
    #fieldset = create_hydrodynamic_fieldset(settings)

    # Now add the other fields
    # Start date and runtime of the simulation
    startdate = settings['simulation']['startdate']
    runtime = int(np.ceil(settings['simulation']['runtime'].total_seconds()/86400.))  # convert to days

    if fieldset.use_biofouling:
        # MOi glossary: https://www.mercator-ocean.eu/wp-content/uploads/2021/11/Glossary.pdf
        # and https://catalogue.marine.copernicus.eu/documents/PUM/CMEMS-GLO-PUM-001-028.pdf

        # Add BGC constants to current fieldset
        for key in settings['bgc']['constants']:
            fieldset.add_constant(key, settings['bgc']['constants'][key])

        # Create a fieldset with BGC data
        dirread_bgc = os.path.join(settings['bgc']['directory'], settings['bgc']['filename_style'])
        bgc_mesh = os.path.join(settings['bgc']['directory'], settings['bgc']['bgc_mesh'])  # mesh_mask_4th

        dirread_model = os.path.join(settings['ocean']['directory'], settings['ocean']['filename_style'])
        wfiles = select_files(dirread_model, 'W_%4i*.nc', startdate, runtime, dt_margin=3)

        ppfiles = select_files(dirread_bgc, 'nppv_%4i*.nc', startdate, runtime, dt_margin=8)
        phy1files = select_files(dirread_bgc, 'phy_%4i*.nc', startdate, runtime, dt_margin=8)
        phy2files = select_files(dirread_bgc, 'phy2_%4i*.nc', startdate, runtime, dt_margin=8)

        filenames_bio = {'pp_phyto': {'lon': bgc_mesh, 'lat': bgc_mesh, 'depth': wfiles[0], 'data': ppfiles},
                         'bio_nanophy': {'lon': bgc_mesh, 'lat': bgc_mesh, 'depth': wfiles[0], 'data': phy1files},
                         'bio_diatom': {'lon': bgc_mesh, 'lat': bgc_mesh, 'depth': wfiles[0], 'data': phy2files}}

        variables_bio = settings['bgc']['variables']
        dimensions_bio = settings['bgc']['dimensions']
        indices_bio = settings['bgc']['indices']

        # Create the BGC fieldset
        bio_fieldset = FieldSet.from_nemo(filenames_bio, variables_bio, dimensions_bio, indices=indices_bio)

        # Add the fields to the main fieldset
        fieldset.add_field(bio_fieldset.pp_phyto)  # phytoplankton primary productivity
        fieldset.add_field(bio_fieldset.bio_nanophy)  # nanopyhtoplankton concentration [mmol C m-3]
        fieldset.add_field(bio_fieldset.bio_diatom)  # diatom concentration [mmol C m-3]

    if fieldset.use_stokes:
        dirread_Stokes = os.path.join(settings['stokes']['directory'], settings['stokes']['filename_style'])
        wavesfiles = select_files(dirread_Stokes, '%4i*.nc', startdate, runtime, dt_margin=32)

        filenames_Stokes = {'Stokes_U': wavesfiles,
                            'Stokes_V': wavesfiles,
                            'wave_Tp': wavesfiles}

        variables_Stokes = settings['stokes']['variables']
        dimensions_Stokes = settings['stokes']['dimensions']
        indices_Stokes = settings['stokes']['indices']

        fieldset_Stokes = FieldSet.from_netcdf(filenames_Stokes, variables_Stokes, dimensions_Stokes, indices=indices_Stokes, mesh='spherical')
        fieldset_Stokes.Stokes_U.units = GeographicPolar()
        fieldset_Stokes.Stokes_V.units = Geographic()
        fieldset_Stokes.add_periodic_halo(zonal=True)

        fieldset.add_field(fieldset_Stokes.Stokes_U)
        fieldset.add_field(fieldset_Stokes.Stokes_V)
        fieldset.add_field(fieldset_Stokes.wave_Tp)

    if fieldset.use_wind:
        dirread_wind = os.path.join(settings['wind']['directory'], settings['wind']['filename_style'])
        windfiles = select_files(dirread_wind, '%4i*.nc', startdate, runtime, dt_margin=32)

        filenames_wind = {'Wind_U': windfiles,
                          'Wind_V': windfiles}

        variables_wind = settings['wind']['variables']
        dimensions_wind = settings['wind']['dimensions']
        indices_wind = settings['wind']['indices']

        fieldset_wind = FieldSet.from_netcdf(filenames_wind, variables_wind, dimensions_wind, indices=indices_wind, mesh='spherical')
        fieldset_wind.Wind_U.units = GeographicPolar()
        fieldset_wind.Wind_V.units = Geographic()
        fieldset_wind.add_periodic_halo(zonal=True)

        fieldset.add_field(fieldset_wind.Wind_U)
        fieldset.add_field(fieldset_wind.Wind_V)

    # Apply unbeaching currents when Stokes/Wind can push particles into land cells
    # if fieldset.use_stokes: # or fieldset.use_wind > 0:
    #     unbeachfiles = settings['unbeaching']['filename']
    #     filenames_unbeach = {'unbeach_U': unbeachfiles,
    #                          'unbeach_V': unbeachfiles}

    #     variables_unbeach = settings['unbeaching']['variables']

    #     dimensions_unbeach = settings['unbeaching']['dimensions']

    #     indices_unbeach = settings['unbeaching']['indices']

    #     fieldset_unbeach = FieldSet.from_netcdf(filenames_unbeach, variables_unbeach, dimensions_unbeach, indices=indices_unbeach, mesh='spherical')
    #     fieldset_unbeach.unbeach_U.units = GeographicPolar()
    #     fieldset_unbeach.unbeach_V.units = Geographic()

    #     fieldset.add_field(fieldset_unbeach.unbeach_U)
    #     fieldset.add_field(fieldset_unbeach.unbeach_V)

    fieldset.add_constant('verbose_delete', settings['verbose_delete'])

    return fieldset


def add_additional_fields_copernicus(fieldset, settings):
    ### Based on create_fieldset from plasticparcels
    """A constructor method to create a `Parcels.Fieldset` with all fields
    necessary for a plasticparcels simulation (e.g., a hydrodynamic model
    velocity field, biogeochemical model variable fields, a wind field, etc.),
    using output downloaded from the Copernicus Marine Store.

    Parameters
    ----------
    settings :
        A dictionary of settings that contains an ocean model information,
        biogeochemical model information, wind model information,
        and other optional settings.

    Returns
    -------
    fieldset
        A `parcels.FieldSet` object.
    """
    
    # Start date and runtime of the simulation
    startdate = settings['simulation']['startdate']
    runtime = int(np.ceil(settings['simulation']['runtime'].total_seconds()/86400.))  # convert to days

    if fieldset.use_biofouling:
        # Add BGC constants to current fieldset
        for key in settings['bgc']['constants']:
            fieldset.add_constant(key, settings['bgc']['constants'][key])

        # Create a fieldset with BGC data
        dirread_bgc = os.path.join(settings['bgc']['directory'], settings['bgc']['filename_style'])

        ppfiles = select_files(dirread_bgc, 'nppv_%4i*.nc', startdate.to_pydatetime().date(), runtime, dt_margin=3) # daily nppv
        phycfiles = select_files(dirread_bgc, 'phyc_%4i*.nc', startdate.to_pydatetime().date(), runtime, dt_margin=32) # monthly phyc

        variables_bio = settings['bgc']['variables']
        dimensions_bio = settings['bgc']['dimensions']

        pp_phyto_fieldset = parcels.FieldSet.from_netcdf(ppfiles, variables_bio, dimensions_bio)
        pp_phyto_fieldset.add_periodic_halo(zonal=True)
        phyc_fieldset = parcels.FieldSet.from_netcdf(phycfiles, variables_bio, dimensions_bio)
        phyc_fieldset.add_periodic_halo(zonal=True)

        fieldset.add_field(pp_phyto_fieldset.pp_phyto)
        fieldset.add_field(phyc_fieldset.bio_nanophy)

    if fieldset.use_stokes:
        dirread_Stokes = settings['stokes']['directory']
        wavesfiles = select_files(dirread_Stokes, settings['stokes']['filename_style'], startdate.to_pydatetime().date(), runtime, dt_margin=3) # daily files

        filenames_Stokes = {'Stokes_U': wavesfiles,
                            'Stokes_V': wavesfiles,
                            'wave_Tp': wavesfiles}

        variables_Stokes = settings['stokes']['variables']
        dimensions_Stokes = settings['stokes']['dimensions']
        indices_Stokes = settings['stokes']['indices']

        fieldset_Stokes = parcels.FieldSet.from_netcdf(filenames_Stokes, variables_Stokes, dimensions_Stokes, indices=indices_Stokes, mesh='spherical')
        fieldset_Stokes.Stokes_U.units = GeographicPolar()
        fieldset_Stokes.Stokes_V.units = Geographic()
        fieldset_Stokes.add_periodic_halo(zonal=True)

        fieldset.add_field(fieldset_Stokes.Stokes_U)
        fieldset.add_field(fieldset_Stokes.Stokes_V)
        fieldset.add_field(fieldset_Stokes.wave_Tp)
    
    fieldset.add_constant('verbose_delete', settings['verbose_delete'])
    
    return fieldset

def deleteParticle(particle, fieldset, time):
    """ Kernel for deleting particles if they throw an error other than through the surface.
    """
    if particle.state >= 50 and particle.state != StatusCode.ErrorThroughSurface:
        particle.inside_eddy = 0.
        particle.delete()

def checkErrorThroughSurface(particle, fieldset, time):
    """ Kernel to set the particle depth to the particle surface if it goes through the surface.
    """
    if particle.state == StatusCode.ErrorThroughSurface:
        particle_ddepth = - particle.depth # reset the depth to the surface
        particle.state = StatusCode.Success

def distance(origin, destination):
    """ Haversine function to compute the great-circle distance between two (lat/lon) points in km
    """
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371 # km

    dlat = np.radians(lat2-lat1)
    dlon = np.radians(lon2-lon1)
    a = np.sin(dlat/2) * np.sin(dlat/2) + np.cos(np.radians(lat1)) \
        * np.cos(np.radians(lat2)) * np.sin(dlon/2) * np.sin(dlon/2)
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
    d = radius * c

    return d

def displacement(inlon1, inlon2, inlat1, inlat2):
    """ Similar to the distance function, but read in lists of lons and lats seperately.

    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)

    ## Returns the distance in km
    """

    #Convert decimal degrees to Radians:
    lon1 = np.radians(inlon1)
    lat1 = np.radians(inlat1)
    lon2 = np.radians(inlon2)
    lat2 = np.radians(inlat2)
    radius = 6371
    # if np.isnan([lon1, lon2, lat1, lat2]).any():
    #     return mask_value
    #Implementing Haversine Formula: 
    dlon = np.subtract(lon2, lon1)
    dlat = np.subtract(lat2, lat1)

    a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2),  
                          np.multiply(np.cos(lat1), 
                                      np.multiply(np.cos(lat2), 
                                                  np.power(np.sin(np.divide(dlon, 2)), 2))))
    c = np.multiply(2, np.arcsin(np.sqrt(a)))
    d = radius * c
    return d

class PlasticParticle(parcels.JITParticle):
    """ A class to define the properties of plastic particles.
    """
    seawater_density = parcels.Variable('seawater_density', dtype=np.float32, to_write=False)
    wind_coefficient = parcels.Variable('wind_coefficient', dtype=np.float32, initial=0., to_write=False)

    settling_velocity = parcels.Variable('settling_velocity', dtype=np.float64, initial=0., to_write=False)
    absolute_salinity = parcels.Variable('absolute_salinity', dtype=np.float64, initial=np.nan, to_write=False)
    algae_amount = parcels.Variable('algae_amount', dtype=np.float64, initial=0., to_write=False)

    plastic_diameter = parcels.Variable('plastic_diameter', dtype=np.float32, to_write='once')
    plastic_density = parcels.Variable('plastic_density', dtype=np.float32, initial=1010., to_write='once')
    plastic_amount = parcels.Variable('plastic_amount', dtype=np.float64, initial=0., to_write='once')
    release_class = parcels.Variable('release_class', dtype=np.int32, to_write='once') # 0 = river, 1 = coast, 2 = fisheries
    release_id = parcels.Variable('release_id', dtype=np.int32, to_write='once') # unique ID for each release location from release map files
    # NOTE max value of an int32 is 2,147,483,647 - we are well within this range for the release_id


def create_particleset_from_release_files(startrelease, endrelease, release_folder, settings, fieldset):
    """ Function to create a ParticleSet from release files."""
    # Plastic sizes, 1000 mm = 1 m -> 1 mm = 0.001 m
    # 6 classes
    plastic_diameters = np.array([0.0001,    # 0.1 mm   -> < 0.3 mm
                            0.00035,   # 3.5 mm   -> 0.3-0.4 mm
                            0.0005,    # 0.5 mm   -> 0.4-0.6 mm
                            0.00105,   # 1.05 mm  -> 0.6-1.5 mm
                            0.002,     # 2 mm     -> 1.5-2.5 mm
                            0.00375])  # 3.75 mm  -> 2.5-5 mm

    plastic_density = 1010 # Density of plastic [kg m-3]

    release_files = np.array(os.listdir(release_folder))
    startrelease = np.datetime64(startrelease)
    endrelease = np.datetime64(endrelease)

    release_dates = np.array([np.datetime64(f.split('_')[2][:-4]) for f in release_files])
    valid_release_file_ids = np.intersect1d(np.where(release_dates >= startrelease), np.where(release_dates <= endrelease))
    valid_release_files = release_files[valid_release_file_ids]


    coastal_ds = pd.read_csv(settings['release_maps']['coastal'])
    rivers_ds = pd.read_csv(settings['release_maps']['rivers'])
    fisheries_ds = pd.read_csv(settings['release_maps']['fisheries'])

    # Arrays that will contain all the release information
    release_lons = []
    release_lats = []
    release_amounts = []
    release_plastic_diameters = []
    release_plastic_density = []
    release_class = []
    release_times = []
    release_ids = []

    # Loop over the release files and constuct a release array - noting we are only doing 1 density class
    for i, release_file in enumerate(valid_release_files):
        release_date = release_dates[valid_release_file_ids][i]
        # Load the release_file information
        release_info = np.load(release_folder+release_file, allow_pickle=True)
        rivers_ids = release_info['rivers_ids']
        coastal_ids = release_info['coastal_ids']
        fisheries_ids = release_info['fisheries_ids']

        # Extract just the relevant release info from the datasets
        rivers_info = rivers_ds.query('release_id in @rivers_ids')
        coastal_info = coastal_ds.query('release_id in @coastal_ids')
        fisheries_info = fisheries_ds.query('release_id in @fisheries_ids')

        # Construct the arrays for the rivers release location
        rivers_lons = np.array(rivers_info['Longitude'].values)
        rivers_lats = np.array(rivers_info['Latitude'].values)
        rivers_amounts = np.array(rivers_info['Emissions'].values)

        rivers_release_lons = np.repeat(rivers_lons, len(plastic_diameters))
        rivers_release_lats = np.repeat(rivers_lats, len(plastic_diameters))
        rivers_release_amounts = np.repeat(rivers_amounts, len(plastic_diameters))

        rivers_plastic_diameters = np.tile(plastic_diameters, len(rivers_lons))
        rivers_plastic_density = np.repeat(plastic_density, len(rivers_release_lons))
        rivers_release_class = np.zeros_like(rivers_release_lons)

        # Construct the arrays for the coastal release location
        coastal_lons = np.array(coastal_info['Longitude'].values)
        coastal_lats = np.array(coastal_info['Latitude'].values)
        coastal_amounts = np.array(coastal_info['MPW_Cell'].values)

        coastal_release_lons = np.repeat(coastal_lons, len(plastic_diameters))
        coastal_release_lats = np.repeat(coastal_lats, len(plastic_diameters))
        coastal_release_amounts = np.repeat(coastal_amounts, len(plastic_diameters))

        coastal_plastic_diameters = np.tile(plastic_diameters, len(coastal_lons))
        coastal_plastic_density = np.repeat(plastic_density, len(coastal_release_lons))
        coastal_release_class = np.ones_like(coastal_release_lons)

        # Construct the arrays for the fisheries release location
        fisheries_lons = np.array(fisheries_info['Longitude'].values)
        fisheries_lats = np.array(fisheries_info['Latitude'].values)
        fisheries_amounts = np.array(fisheries_info['fishing_hours'].values)

        fisheries_release_lons = np.repeat(fisheries_lons, len(plastic_diameters))
        fisheries_release_lats = np.repeat(fisheries_lats, len(plastic_diameters))
        fisheries_release_amounts = np.repeat(fisheries_amounts, len(plastic_diameters))

        fisheries_plastic_diameters = np.tile(plastic_diameters, len(fisheries_lons))
        fisheries_plastic_density = np.repeat(plastic_density, len(fisheries_release_lons))
        fisheries_release_class = np.full_like(fisheries_release_lons, 2)

        # Combine all the release locations
        combined_release_lons = np.concatenate((rivers_release_lons, coastal_release_lons, fisheries_release_lons))
        combined_release_lats = np.concatenate((rivers_release_lats, coastal_release_lats, fisheries_release_lats))
        combined_release_amounts = np.concatenate((rivers_release_amounts, coastal_release_amounts, fisheries_release_amounts))
        combined_plastic_diameters = np.concatenate((rivers_plastic_diameters, coastal_plastic_diameters, fisheries_plastic_diameters))
        combined_plastic_density = np.concatenate((rivers_plastic_density, coastal_plastic_density, fisheries_plastic_density))
        combined_release_class = np.concatenate((rivers_release_class, coastal_release_class, fisheries_release_class))
        combined_release_ids = np.concatenate((np.repeat(rivers_ids, len(plastic_diameters)), np.repeat(coastal_ids, len(plastic_diameters)), np.repeat(fisheries_ids, len(plastic_diameters))))

        # Construct an array of the release times
        combined_release_times = np.repeat(release_date, len(combined_release_lons))

        # Extend to the complete arrays
        release_lons.extend(combined_release_lons)
        release_lats.extend(combined_release_lats)
        release_amounts.extend(combined_release_amounts)
        release_plastic_diameters.extend(combined_plastic_diameters)
        release_plastic_density.extend(combined_plastic_density)
        release_class.extend(combined_release_class)
        release_times.extend(combined_release_times)
        release_ids.extend(combined_release_ids)


    release_lons = np.array(release_lons)
    release_lats = np.array(release_lats)
    release_amounts = np.array(release_amounts)
    release_plastic_diameters = np.array(release_plastic_diameters)
    release_plastic_density = np.array(release_plastic_density)
    release_class = np.array(release_class)
    release_times = np.array(release_times)
    release_ids = np.array(release_ids)


    pset = parcels.ParticleSet(fieldset=fieldset,
                            pclass=PlasticParticle,
                            lon=release_lons,
                            lat=release_lats,
                            time=release_times,
                            plastic_diameter=release_plastic_diameters,
                            plastic_density=release_plastic_density,
                            plastic_amount=release_amounts,
                            release_class=release_class,
                            release_id=release_ids,
                            )

    return pset