Particles Start Out Of Bounds

[ingwag]

Question

Question

Hi, I am trying to set up some simple particle tracking with an analytically calculated velocity field. But I am getting an issue that some particles seem to be initialized out of bounds. I just don't understand why. The coordinates of the particles are very much within the coordinates specified in the grid.

Some of the executions with a small number of particles work. But often they don't. The problematic point that crashes the execution according to the error message is not further away from the center than other points.

Supporting code/error messages

Error:

FieldOutOfBoundError: Field sampled out-of-bound, at (depth=0.0, lat=19589.552734375, lon=5321.353515625)  

fieldset.gridset.grids:

[RectilinearZGrid(lon=array([-200000.00, -195000.00, -190000.00, ...,  190000.00,  195000.00,  200000.00], shape=(81,), dtype=float32), 
lat=array([-200000.00, -195000.00, -190000.00, ...,  190000.00,  195000.00,  200000.00], shape=(81,), dtype=float32), 
time=array([ 0.00,  300.00,  600.00, ...,  863400.00,  863700.00,  864000.00], shape=(2881,)), time_origin=0, mesh='spherical')]

minimal example

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from IPython.display import HTML
import parcels
from datetime import timedelta

def run_parcels():
    #grid definition
    dx = 5e3
    dy = 5e3
    xrange = 2e5
    yrange = 2e5
    x = np.arange(-xrange,xrange+1,dx)
    y = np.arange(-yrange,yrange+1,dx)
    X,Y = np.meshgrid(x,y)

    U = np.zeros(X.shape)
    V = np.zeros(Y.shape)

    
    # time array
    dt = 300
    nt = 86400*10//dt + 1#86400/300 
    time = np.arange(0,nt*dt,dt)

    U = np.repeat(U[np.newaxis,:,:],nt,axis=0) # repeat U nt times in along a third dimension
    V = np.repeat(V[np.newaxis,:,:],nt,axis=0) # repeat V nt times in along a third dimension

    # create xarray
    ds = xr.Dataset(data_vars={'U':(['time','y','x'],U/10),'V':(['time','y','x'],V/10)},
                    coords={'x':(['x'],x),'y':(['y'],y),'time':(['time'],time)})

    #set up the fieldset
    variables = {'U':'U','V':'V'}
    dimensions = {'U': {'lon':'x','lat':'y','time':'time'},
                    'V': {'lon':'x','lat':'y','time':'time'}}

    fieldset = parcels.FieldSet.from_xarray_dataset(ds, variables, dimensions)                

    print(fieldset.gridset.grids)

    #generate particles
    xlocs = np.random.uniform(-2e4,2e4,10)
    ylocs = np.random.uniform(-2e4,2e4,10)

    print(xlocs, ylocs)
    
    # set up the particle set
    pset = parcels.ParticleSet.from_list(
    fieldset=fieldset,
    pclass=parcels.JITParticle,
    lon = xlocs,  # particle longitudes
    lat = ylocs,  # particle latitudes
    )
    output_file = pset.ParticleFile(
    name="testParticles.zarr", outputdt=timedelta(hours=1)
    )   

    pset.execute(
        [parcels.AdvectionRK4],
        runtime=timedelta(days=10),
        dt=timedelta(seconds=300),
        output_file=output_file,
    )

if __name__=='__main__':
    run_parcels()

---

Thank you so much if you take your time to have a look at my problem.

best,
Ingo

[erikvansebille]

Thanks for reporting @ingwag. I see in the output of your print-statement that your Grid is assumed to be spherical (i.e. in degrees longitude and latitude) instead of flat (i.e. in km).

I'm not entirely sure if this causes the issue, but could you try adding the extra argument mesh='flat' when you create your fieldset?

Also, does the problem appear also when you start a particle at lon=0, lat=0? Or does it only happen close to the boundaries?

[ingwag]

Ah, yes thank you so much. I certainly missed that setting. It is supposed to be flat and setting it to flat does fix the issue!

I learned another thing and won't forget to specify the mesh type again :)