Summed Field Indexing Error

[JonathanC33]

Hi, I'm trying to run ocean parcels using a summed windage and ocean velocity field and the initial particle set locations are determined from real-world drifter starting location and times. Every time I start the pset.execute command though I get an immediate error "index 990 is out of bounds for axis 0 with size 990" where 990 is one less than the number of the rows in the windage velocity field. I've also checked the starting location of the drifter and it is within the bounds of both the windage and ocean velocity fields so I'm not sure where the error comes. I've included my code below.

`
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation, writers, PillowWriter

import matplotlib
matplotlib.interactive(False)
from parcels import Field, NestedField, FieldSet, ParticleSet, Variable, JITParticle, AdvectionRK4, plotTrajectoriesFile, ErrorCode
import numpy as np
import math
from datetime import timedelta
from operator import attrgetter
import xarray as xr

import cartopy.crs as ccrs
import cartopy

#Import the Bay of Fundy large domain U and V velocity fields
filepath = '/home/joc000/Documents/Data/ocean_drifters/'
filepath_const = '/home/stt001/gpfs4/CONSTANTS/Fundy500/'
coord = filepath_const + 'Bathymetry_Fundy500.nc'
mesh_mask = filepath_const + 'mesh-mask_Fundy500.nc'
filepath_depth = filepath + 'FUNDY500_HC_1h_grid_U_20191225-20191225.nc'

filepath_ocean_u = filepath+'FUNDY500_HC_1h_grid_U_20*.nc'
filepath_ocean_v = filepath+'FUNDY500_HC_1h_grid_V_20*.nc'

#Names for the U and V files, opens multiple
filenames = {
'U': {'lon': coord, 'lat': coord, 'depth':filepath_depth, 'data': filepath_ocean_u},
'V': {'lon': coord, 'lat': coord, 'depth':filepath_depth, 'data': filepath_ocean_v}
}

#Define the dictionary of the variables and the dimensions, no depth in this case
variables = {
'U': 'uo',
'V': 'vo'
}

dimensions = {
'lat': 'nav_lat',
'lon': 'nav_lon',
'depth': 'depthu',
'time': 'time_instant'
}

#Take the depth at 15m (approximately, it's at 14.68m)
indices = {
'U': {'depth': [0]},
'V': {'depth': [0]}
}

#Fieldset using FieldSet.from_netcdf using the filenames, variables, and dimensions
fieldset_ocean = FieldSet.from_netcdf(
filenames,
variables,
dimensions,
indices=indices,
gridindexingtype='nemo',
chunksize=False,
allow_time_extrapolation=True)

#Import the windage
filepath_windg = '/home/joc000/Documents/Data/windage_drifters/'
filepath_wind = filepath_windg+'HRDPS_OPPeast_ps2.5km_y20*.nc'
filepath_coord = filepath_windg+'HRDPS_OPPeast_ps2.5km_y2019m12d20.nc'

#Account for the difference in longitude coordinates
#pfile = xr.open_dataset(filepath_coord)
#pfile['nav_lon'] = pfile['nav_lon'] - 360
#pfile.to_netcdf(filepath_windg+'windage_coords.nc')
filepath_coord = filepath_windg+'windage_coords.nc'

pfile = xr.open_dataset(filepath_windg+'HRDPS_OPPeast_ps2.5km_y2020m01d20.nc')

filenames = {
'U': {'lon': filepath_coord, 'lat': filepath_coord, 'data': filepath_wind},
'V': {'lon': filepath_coord, 'lat': filepath_coord, 'data': filepath_wind}
}

#Define the dictionary of the variables and the dimensions, no depth in this case
variables = {
'U': 'u_wind',
'V': 'v_wind'
}

dimensions = {
'lat': 'nav_lat',
'lon': 'nav_lon',
'time': 'time_counter'
}

#Fieldset using FieldSet.from_netcdf using the filenames, variables, and dimensions
fieldset_wind = FieldSet.from_netcdf(
filenames,
variables,
dimensions,
chunksize=False,
gridindexingtype='nemo',
allow_time_extrapolation=True)

#Scale the wind some percentage
percentage_wind = 0.5
fieldset_wind.U.set_scaling_factor(percentage_wind/100)
fieldset_wind.V.set_scaling_factor(percentage_wind/100)

#Add the two together, wind and ocean
fieldset_sum = FieldSet(
U=fieldset_wind.U+fieldset_ocean.U,
V=fieldset_wind.V+fieldset_ocean.V)
#,fields={'umask': fieldset_ocean.umask})

#Import the Davis drifter data
drifter_number = ['00','01','02','03','04','05','06','07','08','09','10','11','12','13']
drifter_dates = ['20191220-20200224','20191220-20200111','20200106-20200112',
'20200106-20200224','20200120-20200414','20200120-20200221',
'20200203-20200228','20200203-20200414','20200218-20200323',
'20200218-20200414','20200302-20200405','20200302-20200406',
'20200319-20200414','20200319-20200414']

#Run through each of the Davis drifters
#Run for the entire length of the drifter
for i in np.arange(np.shape(drifter_number)[0]):

#Import the relevant data
filepath_drifter = '/gpfs/fs4/dfo/dpnm/dfo_dpnm/sdfo600/OBSERVATIONS/OBS-EAST/FERRY/'
filepath_drifter = filepath_drifter+'SJHDrifters-Davis-Winter1920_'+\
drifter_number[i]+'_'+\
drifter_dates[i]+'.nc'

pfile = xr.open_dataset(filepath_drifter)
drifter_time = pfile.time.values
drifter_lon = pfile.longitude.values
drifter_lat = pfile.latitude.values
pfile.close()

#Set the location of the OP particles using the initial drifter location,0.0002 spacing in between
drifter_start = 2
lonmin, lonmax = drifter_lon[drifter_start]-0.001,drifter_lon[drifter_start]+0.001
latmin, latmax = drifter_lat[drifter_start]-0.001,drifter_lat[drifter_start]+0.001
lons,lats = np.meshgrid(np.arange(lonmin,lonmax,0.0002),np.arange(latmin,latmax,0.0002))

#Define the starting time, lons, lats in the pset
pset = ParticleSet.from_list(fieldset=fieldset_sum, pclass=JITParticle,
lon=lons.flatten(),
lat=lats.flatten(),
time=drifter_time[drifter_start])

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

#Define the starting time of the execute
output_file = pset.ParticleFile(name='/home/joc000/Documents/OceanParcelsOutput/Drifter_OP_execute/DavisDrifter_'+\
                                drifter_number[i]+'_0pt5percwindage.nc',
                                outputdt=timedelta(minutes=60))
    
kernels = pset.Kernel(AdvectionRK4) #+ LandWind
pset.execute(kernels,
runtime=timedelta(hours=24*10),
dt=timedelta(minutes=1),
recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle},
output_file=output_file)
output_file.export()`

[erikvansebille]

Thank you for asking this question, @JonathanC33. It's hard to find out what's happening without seeing the full error message. However, could it be that your drifter longitudes are defined between 0 and 360 degrees, while your fieldset longitude is defined between -180 and 180 degrees (or the other way around)? That kind of mismatch in coordinate systems often gives errors on the western hemisphere

[JonathanC33]

Thank you for the response. I've double checked both velocity fields and the drifter data and all are defined using longitudes from -180 to 180. I also checked the summed velocity field and the longitudes are also -180 to 180. I'll inlcude the full error message below:

`INFO: Compiled ArrayJITParticleAdvectionRK4 ==> /tmp/parcels-68372/lib754ae70b005a15a472c88ba65e93b0c9_0.so
Traceback (most recent call last):

Input In [42] in
pset.execute(kernels,

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/particleset/baseparticleset.py:433 in execute
next_input = self.fieldset.computeTimeChunk(time, np.sign(dt)) if self.fieldset is not None else np.inf

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldset.py:1072 in computeTimeChunk
data = f.computeTimeChunk(data, tind)

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/field.py:1375 in computeTimeChunk
buffer_data = filebuffer.data

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldfilebuffer.py:159 in data
return self.data_access()

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldfilebuffer.py:164 in data_access
data = self._apply_indices(data, ti)

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/parcels/fieldfilebuffer.py:149 in _apply_indices
data = data[ti, self.indices['lat'], self.indices['lon']]

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/dataarray.py:745 in getitem
return self.isel(indexers=self._item_key_to_dict(key))

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/dataarray.py:1200 in isel
ds = self._to_temp_dataset()._isel_fancy(

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/dataset.py:2419 in _isel_fancy
new_var = var.isel(indexers=var_indexers)

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/variable.py:1193 in isel
return self[key]

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/variable.py:786 in getitem
data = as_indexable(self._data)[indexer]

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/indexing.py:555 in getitem
return type(self)(_wrap_numpy_scalars(self.array[key]))

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/indexing.py:524 in getitem
return type(self)(_wrap_numpy_scalars(self.array[key]))

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/indexing.py:431 in getitem
return type(self)(self.array, self._updated_key(indexer))

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/indexing.py:403 in _updated_key
full_key.append(_index_indexer_1d(k, next(iter_new_key), size))

File ~/miniconda3/envs/py3_parcels/lib/python3.9/site-packages/xarray/core/indexing.py:176 in _index_indexer_1d
indexer = _expand_slice(old_indexer, size)[applied_indexer]

IndexError: index 990 is out of bounds for axis 0 with size 990`

[erikvansebille]

Hmm, I don't immediately see what's wrong, based on this error message and your code. The error is somewhere in the loading of the FieldSet, but not sure where and why.
I suggest that you carefully 'debug' your code in the sense that you try to strip it down to its simplest version and then one by one add extra complexity. For example, does the code work if you don't use SummedFields (i.e. only use the ocean flow or only use the wind)? If it still breaks with either of these non-summed FieldSets, that would be a hint...

[JonathanC33]

I did what you suggested and found that the windage data I'm using changes size at a specific date and that was causing the error. Thank you for your help!