Particles initial locations out of bounds on run with HYCOM

[molkjam]

I am trying to run PARCELS with the output of the HYCOM hydrodynamic model.
I am able to initialise the fieldset with:

fname = 'uv3z.nc4'
filenames = {'U': fname, 'V': fname}
variables = {'U': 'water_u', 'V': 'water_v'}
dimensions = {'lat': 'lat', 'lon': 'lon', 'time': 'time','depth':'depth'} # In the GlobCurrent data the dimensions are also called 'lon', 'lat' and 'time'
             
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions,vmin=-2000)

and define the initial locations of my particles with:

#Define the particles type and initial conditions in a ParticleSet object
sitelocs = pandas.read_csv('sites.csv')

lat = sitelocs.lat.tolist()
lon =sitelocs.lon.tolist()
sites = len(lat)
n = 5 #number of replicate releases from each site

print('releasing ',n*sites,'particles from ',sites,'sites (',n,'particles per site)')
lat = list(itertools.chain.from_iterable(itertools.repeat(x, n) for x in lat))
lon = list(itertools.chain.from_iterable(itertools.repeat(x, n) for x in lon))

pset = ParticleSet(fieldset=fieldset,   
                   pclass=JITParticle,  
                   lon=lon,              # release longitude
                   lat=lat)             # release latitude


pset.show(field=fieldset.U)

However, when I try to execute the kernel I fall at the first hurdle:

#Execute an advection kernel on the given fieldset

output_file = pset.ParticleFile(name="reverse5.nc", outputdt=timedelta(hours=1)) 
pset.execute(AdvectionRK4,                 
             runtime=timedelta(days=1),    
             dt=-timedelta(minutes=5),      
             output_file=output_file)


INFO: Compiled JITParticleAdvectionRK4 ==> /tmp/parcels-15572/75f11f043bd3a0fa57b03bbe68452681_0.so

---------------------------------------------------------------------------
OutOfBoundsError                          Traceback (most recent call last)
<ipython-input-4-d68316a1e96a> in <module>()
      9              runtime=timedelta(days=1),    # the total length of the run
     10              dt=-timedelta(minutes=5),      # the timestep of the kernel (negative timestep for backwards run)
---> 11              output_file=output_file)
     12 

~/envs/parcels/lib/python3.6/site-packages/parcels/particleset.py in execute(self, pyfunc, endtime, runtime, dt, moviedt, recovery, output_file, movie_background_field, verbose_progress, postIterationCallbacks, callbackdt)
    767                 time = max(next_prelease, next_input, next_output, next_movie, next_callback, endtime)
    768             self.kernel.execute(self, endtime=time, dt=dt, recovery=recovery, output_file=output_file,
--> 769                                 execute_once=execute_once)
    770             if abs(time-next_prelease) < tol:
    771                 pset_new = ParticleSet(fieldset=self.fieldset, time=time, lon=self.repeatlon,

~/envs/parcels/lib/python3.6/site-packages/parcels/kernel.py in execute(self, pset, endtime, dt, recovery, output_file, execute_once)
    424                     recovery_kernel = recovery_map[particles.state]
    425                     particles.set_state(StateCode.Success)
--> 426                     recovery_kernel(particles, self.fieldset, particles.time)
    427                     if particles.state == StateCode.Success:
    428                         particles.set_state(StateCode.Evaluate)

~/envs/parcels/lib/python3.6/site-packages/parcels/tools/statuscodes.py in recovery_kernel_out_of_bounds(particle, fieldset, time)
    184         # TODO: JIT does not yet provide the context that created
    185         # the exception. We need to pass that info back from C.
--> 186         raise OutOfBoundsError(particle, fieldset)
    187     else:
    188         error = particle.exception

OutOfBoundsError: 0
Particle P[0](lon=-61.644100, lat=17.032516, depth=0.000000, time=10443600.000000)
Time: 2019-12-31T09:00:00.000000000,	timestep dt: -300.000000
Out-of-bounds sampling by particle at (-61.644100, 17.032516, 0.000000)

Its not too clear to me why this is happening as my plot seems to suggest that my particles are well within the bounds of my model. Any insight would be greatly appreciated.

[erikvansebille]

Hi @molkjam. Might it be that your model domain longitude is defined between 0 and 360 degrees? Because from the error I see that your initial longitude is set at -61 degrees. That would explain the out-of-bounds error. The solution would also be easy: just add 360 to all your initial longitudes

The reason why you the plotting works is that cartopy takes this 360 periodic boundary condition into account. We might thing of displaying a warning in pset.show() when the particle is offset by 360 degrees

Let us know if this possible solution works

[molkjam]

Hi @erikvansebille thanks for the quick reply.

This solution has indeed solved my issue. Many thanks

[AshleighWomack]

I have the same issue but with curvilinear data from a netcdf file (OSI-SAF ice drift data -polar stereographic).
I followed the same methodology as above.
The dataset lons go from -179 to 180 (but with y, x dimensions).
Everything seems to work until i run pset.execute.
Also when I try: pset.show(field='vector',land=True, domain={'N':-40, 'S':70, 'E':60 ,'W':-170})
I get the error: FieldOutOfBoundError: U sampled out-of-bound, at (-170, 70, 0).
Is there something wrong with my gridding?

[erikvansebille]

Thanks for reporting, @AshleighWomack. I've seen before that Parcels has difficulty with PolarStereographic projections, so it might be that you've just hit a limitation of Parcels. It there a way to convert the grid to curvilinear?

I'm afraid I can't help more without access to a minimal reproducible example of your code, including a snippet of your data.

[AshleighWomack]

Changing the grid worked. Thanks so much.

[Xyx105]

Hi @erikvansebille, I am struggling with this error as well:

`---------------------------------------------------------------------------
OutOfBoundsError Traceback (most recent call last)
Cell In[86], line 82
78 kernels = pset.Kernel(AdvectionRK4)
80 pset.show(fieldset.U)
---> 82 pset.execute(kernels, runtime=delta(days=1), dt=delta(hours=6), output_file=pfile)

File ~\Anaconda3\envs\parcels\Lib\site-packages\parcels\particleset\baseparticleset.py:548, in BaseParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, moviedt, recovery, output_file, movie_background_field, verbose_progress, postIterationCallbacks, callbackdt)
546 # If we don't perform interaction, only execute the normal kernel efficiently.
547 if self.interaction_kernel is None:
--> 548 self.kernel.execute(self, endtime=next_time, dt=dt, recovery=recovery, output_file=output_file,
549 execute_once=execute_once)
550 # Interaction: interleave the interaction and non-interaction kernel for each time step.
551 # E.g. Inter -> Normal -> Inter -> Normal if endtime-time == 2*dt
552 else:
553 cur_time = time

File ~\Anaconda3\envs\parcels\Lib\site-packages\parcels\kernel\kernelsoa.py:213, in KernelSOA.execute(self, pset, endtime, dt, recovery, output_file, execute_once)
211 recovery_kernel = recovery_map[p.state]
212 p.set_state(StateCode.Success)
--> 213 recovery_kernel(p, self.fieldset, p.time)
214 if p.isComputed():
215 p.reset_state()

File ~\Anaconda3\envs\parcels\Lib\site-packages\parcels\tools\statuscodes.py:182, in recovery_kernel_out_of_bounds(particle, fieldset, time)
178 """Default sampling error kernel that throws OutOfBoundsError."""
179 if particle.exception is None:
180 # TODO: JIT does not yet provide the context that created
181 # the exception. We need to pass that info back from C.
--> 182 raise OutOfBoundsError(particle, fieldset)
183 else:
184 error = particle.exception

OutOfBoundsError: 0
Particle P[732](lon=4.700000, lat=53.700000, depth=0.000000, time=0.000000)
Time: 2022-01-01T00:00:00.000000000
timestep dt: 7200.0
Out-of-bounds sampling by particle at (4.7, 53.7, 0.0)`

Here is my code:
`filenames = {
"U": {
"lon": f"output.nc",
"lat": f"output.nc",
"data": f"NC Files/NEATL*.nc",
},
"V": {
"lon": f"output.nc",
"lat": f"output.nc",
"data": f"NC Files/NEATL*.nc",
},
}

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

dimensions = {"lon": "lon_u", "lat": "lat_u", "time": "ocean_time"}

fieldset = FieldSet.from_nemo(
filenames, variables, dimensions, allow_time_extrapolation=True,
)

fieldset.computeTimeChunk(fieldset.U.grid.time[0], 6)
fieldset.computeTimeChunk(fieldset.V.grid.time[0], 6)

npart = 1
lonp = 4.7
latp = 53.7

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

pset = ParticleSet.from_list(field_set, JITParticle, lon=lonp, lat=latp, )
print(pset)

pfile = ParticleFile("nemo_particles", pset, outputdt=delta(days=1))
kernels = pset.Kernel(AdvectionRK4)

pset.show(fieldset.U)

pset.execute(kernels, runtime=delta(days=1), dt=delta(hours=6), output_file=pfile)`

I have tried adding 360 to my longitude but I still get the same error. Any help would be greatly appreciated!