FieldOutOfBoundError at beginning of integration

[RyMole]

Question

Question

Hi all 👋
Similar to a number of other questions, I am also running in to a FildOutOfBoundError at the beginning of the simulation, despite particles being within bounds of the field. The solution given in #2071 doesn't work as my grid is in lat/lon and the mesh is spherical.

I am working with SOSE 5-day output (it 155, year 2022) available here which covers longitude 0 to 360 (stored as 0 to 360) and latitude -78 to -30. To begin with I was initialising a large number of particles but in the minimum working example below I still have the problem when initialising a single particle.

Through attempting to solve the problem I have found that I can get simulation to run if I only load part of the dataset using the indices keyword but that this is highly dependent on the combination of initial position and chosen indices. An example of this is given after the MWE below.

I'm a bit stuck at how to proceed with this and would appreciate any help to get things going in the right direction again! Thanks for your time.

Ryan

---

Supporting code/error messages

parcels.__version__:

'3.1.4'

Error:

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

fieldset.gridset.grids:

[CurvilinearZGrid(lon=array([[ 157.50,  157.67,  157.83, ...,  0.00,  0.00,  0.00],
        [ 157.50,  157.67,  157.83, ...,  0.00,  0.00,  0.00],
        [ 157.50,  157.67,  157.83, ...,  0.00,  0.00,  0.00],
        ...,
        [ 157.50,  157.67,  157.83, ...,  359.50,  359.67,  359.83],
        [ 157.50,  157.67,  157.83, ...,  359.50,  359.67,  359.83],
        [ 157.50,  157.67,  157.83, ...,  359.50,  359.67,  359.83]], dtype=float32), lat=array([[-78.00, -78.00, -78.00, ...,  0.00,  0.00,  0.00],
        [-77.97, -77.97, -77.97, ...,  0.00,  0.00,  0.00],
        [-77.93, -77.93, -77.93, ...,  0.00,  0.00,  0.00],
        ...,
        [-30.09, -30.09, -30.09, ..., -30.09, -30.09, -30.09],
        [-29.95, -29.95, -29.95, ..., -29.95, -29.95, -29.95],
        [-29.80, -29.80, -29.80, ..., -29.80, -29.80, -29.80]], dtype=float32), time=array([ 432000.00,  864000.00]), time_origin=2022-01-03T00:00:00.000000000, mesh='spherical')]

Minimum Working Example

import matplotlib.pyplot as plt
import xarray as xr
import parcels
from datetime import timedelta


# set up fieldset

grid = 'grid.nc'
filenames = {
    'U': {'data': 'Uvel_bsoseI155_2022_5day.nc',
          'lon': grid,
          'lat': grid,
          'depth': grid
          },
    'V': {'data': 'Vvel_bsoseI155_2022_5day.nc',
          'lon': grid,
          'lat': grid,
          'depth': grid
          },
}
variables = {'U': 'UVEL', 'V': 'VVEL'}
dimensions = {'U': {'lon': 'XG', 'lat': 'YG', 'time': 'time', 'depth': 'RC'},
              'V': {'lon': 'XG', 'lat': 'YG', 'time': 'time', 'depth': 'RC'}}


fieldset = parcels.FieldSet.from_mitgcm(filenames,
                                        variables,
                                        dimensions,
                                        indices={'depth': [0],
                                                 # 'lon': range(945, 2160)
                                                 }
                                        )

# set up particles

lon = [226]
lat = [-60]
time = [0]

pset = parcels.ParticleSet.from_list(fieldset=fieldset,
                                     pclass=parcels.JITParticle,
                                     lon=lon,
                                     lat=lat,
                                     time=time)

name = 'mwe_parcels.zarr'
outputdt = timedelta(hours=1)
output_file = pset.ParticleFile(name=name, outputdt=outputdt)


# Plot the field and initial position

fieldset.computeTimeChunk()
plt.pcolormesh(fieldset.U.grid.lon,
               fieldset.U.grid.lat,
               fieldset.U.data[0, :, :])
for p in pset:
    plt.scatter(p.lon, p.lat, c='r')
plt.show()


# set kernels and integration settings

kernels = [parcels.AdvectionRK4]
runtime = timedelta(days=6)
dt = timedelta(minutes=5)

# execute the integration

pset.execute(kernels, runtime=runtime, dt=dt, output_file=output_file)


# plot the trajectory
ds = xr.open_zarr(output_file.fname)
plt.plot(ds.lon.T, ds.lat.T, ".-")
plt.show()

Further Information

Output of the First Plot

Controlling Field with Indices

Different combinations of loaded longitudes and particle initialisation longitudes can be tested by uncommenting the line in the fieldset creation. longitude indices range(945, 2160) and initial position lon = [226] will run successfully.

100%|██████████| 518400.0/518400.0 [00:02<00:00, 191789.94it/s]

Running with longitude indices range(944, 2160) and initial position lon = [226] will fail.

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

If I maintain the same longitude indices range(944, 2160) but change the initial position to lon = [225] it will again run successfuly.

100%|██████████| 518400.0/518400.0 [00:02<00:00, 184119.86it/s]

[michaeldenes]

Hi @RyMole, If I recall correctly from last year, @sruehs was using the same (or similar) dataset, and the problem was not with parcels, but the grid provided with the data. See this comment from last year #1685 (comment)

You can see the issue when you plot your data using pcolormesh, there are grid cells defined near [0,0] that shouldn't exist. To solve the problem you will need to find or correct the grid you use when creating your fieldset.

EDIT: I realise the link to the other issue seems unrelated, but I do remember the same issue popped up in subsequent discussions, hence the link.

[RyMole]

Hi @michaeldenes. Was just coming to comment that some of the lat/lon data for land cells are masked to 0 in the dataset. Given that the particle is not initialised near any of these cells and that the simulation can still run with some of these cells loaded, I have been assuming that this was not part of the problem.

Thanks for providing the link. I'll look through the discussion and try fixing the grid as a first step.

[RyMole]

I should add that the lat/lon data of land cells are masked to 0 and this results in the artefacts in the figures.

There also does not appear to be anything particularly wrong with the data along longitude index 944:

[michaeldenes]

I think the issue may be caused in the grid search part of the code, where perhaps when iterating through the grid to find which cell the particle is in, it encounters one of these [0,0] cells, and ultimately can't find the right cell and throws an out of bounds error. This is me speculating though, if you can try a corrected grid (or manufacture a grid?) to see if that works?

[RyMole]

Yes I think you might be right. As a quick test I have tried the simulation controlling which latitude bands are loaded via the indices argument. If I exclude the southerly bands with a lot of zero-masked cells then the simulation runs successfully (with all longitudes loaded). I'll generate a corrected grid next and see if I can run the full domain.

[RyMole]

Thanks for the thoughts and speculation @michaeldenes as this seems to have been the root of the problem. The program was able to run after I corrected the grid data to remove the masking. It's not clear to me why I was still able to run it after only loading a partial grid which still included masked locations, but I guess it comes down to how the grid search has been implemented.

I wasn't able to find a method to overwrite the grid values after the fieldset was already created, so ended up having to generate an additional netcdf file to use as input.