OutOfBoundError: I have this code from parcel but it looks like if the particle trajects at shallow depths, it produces an OutOfBounds Error message. I would like to either delete the particles that falls out of the domain or have it bounced back to its original position. I included the 'deleteparticle' function and also produces a different error message. Any help will be much appreciated. Thanks about parcels HOT 1 CLOSED

#This script shows the trajectory of a single particle through the ocean
#************************************************************************

import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc4
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as colors
from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4_3D, ErrorCode, ScipyParticle, plotTrajectoriesFile
from glob import glob
from datetime import timedelta as delta
#from parcels import logger, XarrayDecodedFilter
#logger.addFilter(XarrayDecodedFilter())  # Add a filter for the xarray decoding warning


location  = '/home/bejamino/Desktop/DataFiles/'

evel_Data = nc4.Dataset(location + 'EVELMASS/Monthly_Avg_Files/EVEL_01_avg.nc')
wvel_Data = nc4.Dataset(location + 'WVELMASS/Monthly_Avg_Files/WVEL_01_avg.nc')
nvel_Data = nc4.Dataset(location + 'NVELMASS/Monthly_Avg_Files/NVEL_01_avg.nc')

#loading the coordinates of the data
lonGrid       = evel_Data.variables['longitude']
latGrid       = evel_Data.variables['latitude'] 

#Making the streamplot at the surface in January
EVEL      = evel_Data.variables['EVELMASS'][0,0,:,:]
NVEL      = nvel_Data.variables['NVELMASS'][0,0,:,:]


ufiles = sorted(glob(location + 'EVELMASS/Monthly_Avg_Files/EVEL_*.nc'))
vfiles = sorted(glob(location + 'NVELMASS/Monthly_Avg_Files/NVEL_*.nc'))
wfiles = sorted(glob(location + 'WVELMASS/Monthly_Avg_Files/WVEL_*.nc'))


filenames = {'U': ufiles, 'V': vfiles,'W': wfiles}

variables = {'U': 'EVELMASS', 'V': 'NVELMASS','W': 'WVELMASS'}

dimensions = {'U': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'Z',  'time': 'time'},
              'V': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'Z',  'time': 'time'},
              'W': {'lon': 'longitude', 'lat': 'latitude', 'depth': 'Zl', 'time': 'time'}}

fset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True, time_periodic=False)


pset = ParticleSet.from_line(fieldset=fset, pclass=JITParticle,
                             size   = 1, #the number of particles that are released
                             start  = (-30, -40), #defining the starting location (lon,lat)
                             finish = (-30, 30), #defining the finishing location (lon,lat)
                             depth  =100) #the depth 

#here i defined the particle delete function
#def DeleteParticle(particle, fieldset, time):
 #   particle.delete()
    
    
#defining the particle locations for 8 particle sets (somewhere in the Atlantic)
#p_lons=np.array([-30,-30,-30,-30,-30,-30,-30,-30])
#p_lats=np.array([-40,-30,-20,-10,0,10,20,30])
#p_dep=np.array([100,100,100,100,100,100,100,100])

#pset = ParticleSet(fieldset=fset, pclass=JITParticle, lon=p_lons, lat=p_lats, depth=p_dep)

#kernels = pset.Kernel(AdvectionRK4_3D)
#pset.execute(kernels, runtime=delta(days=30), dt=delta(hours=6))

depth_level = 8
print("Level[%d] depth is: [%g %g]" % (depth_level, fset.W.grid.depth[depth_level],
                                       fset.W.grid.depth[depth_level+1]))


def colorGradient(val, minval, maxval):
   # min_depth = 100
    #max_depth = 1000
        
    if val < minval:
        return (0,0,0) #returns black particle
    
    if val > maxval:
        return (1,0,0) #returns red particle
    
#interpolating continuously the depth
    v = (val - minval) / (maxval - minval)
    return (v,0,0)

#Generating the trajectories for the particles
nParticles = len(pset)
nDays = 40
nSteps = 20 #this show how often the particle is propagated for nDays (if it goes in total beyong 360 is crushes)
kernels = pset.Kernel(AdvectionRK4_3D)
plocations = [[] for _ in range(nParticles)]

#recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle}

for i in range(nSteps):
    pset.execute(kernels, runtime=delta(days=nDays), dt=delta(hours=6))
    for k, p in enumerate(pset):
        if k == 0:
            print(f"Step {i}, particle {k}:")
            print(p)
            print(f"lon: {p.lon}, lat: {p.lat}, depth: {p.depth}")
        plocations[k].append((p.lon, p.lat, p.depth))
        
        
#making the plot
plt.ion()
fig = plt.figure(figsize=(12,8))

# miller projection with Basemaps
map = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90,llcrnrlon=-180,urcrnrlon=180)


# plot coastlines, draw label meridians and parallels.
map.drawmeridians(np.arange(0,360,60),labels=[False,False,False,True],linewidth=(0.0))
map.drawparallels(np.arange(-80,80,20),labels=[True,False,False,False],linewidth=(0.0))

# Fill the continents 
map.fillcontinents(color='0.7',lake_color='0.7')
map.drawmapboundary(fill_color='0.7')
map.fillcontinents(color='#ffe2ab',lake_color='0.7')
map.drawcoastlines()

magnitude  = (EVEL ** 2 + NVEL ** 2) ** 0.5       

#Plotting the stream pattern
map.streamplot(x=lonGrid, y=latGrid, u=EVEL, v=NVEL, linewidth=2, density=2, color=magnitude)
        

#plt.title('Put the title here', fontsize=8)

plt.xlabel('Longitude [$^\circ$]', labelpad=20, fontsize=12)
plt.ylabel('Latitude [$^\circ$]', labelpad=40, fontsize=12)
    
    
# Plot trajectories    
min_depth, max_depth = 100, 1000

#colorBy = "depth"
colorBy = "step"
for j, pos in enumerate(plocations):
    for i in range(nSteps):
        lon, lat, depth = pos[i]
        if colorBy == "depth":
            map.scatter(lon, lat, s=4, color=colorGradient(depth, min_depth, max_depth))
        elif colorBy == "step":
            map.scatter(lon, lat, s=4, color=colorGradient(i, 0, nSteps))