TypeError in ChunkCachedArray: '<' not supported between instances of 'slice' and 'int'

[erikvansebille]

Parcels version

v4.0.0

Description

When trying simulation on a Delft3D fieldset with to_chunk_cached_arrays(), I ran into the following error. I'm afraid it won't be easy to find a minimal breaking example - but perhaps the full error log can point to what is going on?

Code sample

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[9], line 15
     11 def DeleteOnAnyError(particles, fieldset):
     12     any_error = particles.state >= 50  # This captures all Errors
     13     particles[any_error].state = parcels.StatusCode.Delete
     14 
---> 15 pset.execute(
     16     [parcels.kernels.AdvectionRK2, DeleteOnAnyError],
     17     endtime=fieldset.time_interval.right,
     18     dt=np.timedelta64(60, 's'),

File ~/parcels/src/parcels/_core/particleset.py:454, in ParticleSet.execute(self, kernels, dt, endtime, runtime, output_file, verbose_progress)
    451 else:
    452     next_time = end_time
--> 454 self._kernel.execute(self, endtime=next_time, dt=dt)
    456 if next_output is not None:
    457     if np.abs(next_time - next_output) < 0.001:

File ~/parcels/src/parcels/_core/kernel.py:209, in Kernel.execute(self, pset, endtime, dt)
    206 with warnings.catch_warnings():
    207     warnings.simplefilter("ignore", FieldEvalWarning)
--> 209     f(pset[evaluate_particles], self._fieldset)
    211     # check for particles that have to be repeated
    212     repeat_particles = pset.state == StatusCode.Repeat

File ~/parcels/src/parcels/kernels/_advection.py:22, in AdvectionRK2(particles, fieldset)
     20 def AdvectionRK2(particles, fieldset):  # pragma: no cover
     21     """Advection of particles using second-order Runge-Kutta integration."""
---> 22     (u1, v1) = fieldset.UV[particles]
     23     x1 = particles.x + u1 * 0.5 * particles.dt
     24     y1 = particles.y + v1 * 0.5 * particles.dt

File ~/parcels/src/parcels/_core/field.py:300, in VectorField.__getitem__(self, key)
    298 try:
    299     if isinstance(key, (ParticleSetView, ParticleSet)):
--> 300         return self.eval(key.t, key.z, key.y, key.x, key)
    301     else:
    302         return self.eval(*key)

File ~/parcels/src/parcels/_core/field.py:286, in VectorField.eval(self, t, z, y, x, particles)
    282 x = np.atleast_1d(x)
    284 particle_positions, grid_positions = _get_positions(self.U, t, z, y, x, particles, _ei)
--> 286 (u, v, w) = self._interp_method.interp(particle_positions, grid_positions, self)
    288 for vel in (u, v, w):
    289     _update_particle_states_interp_value(particles, vel)

File ~/parcels/src/parcels/interpolators/_xinterpolators.py:489, in XFreeslip.interp(self, particle_positions, grid_positions, vectorfield)
    482 def interp(
    483     self,
    484     particle_positions: dict[str, float | np.ndarray],
    485     grid_positions: dict[ptyping.XgridAxis, dict[str, int | float | np.ndarray]],
    486     vectorfield: VectorField,
    487 ):
    488     """Free-slip boundary condition interpolation for velocity fields."""
--> 489     return _Spatialslip(particle_positions, grid_positions, vectorfield, a=1.0, b=0.0)

File ~/parcels/src/parcels/interpolators/_xinterpolators.py:405, in _Spatialslip(particle_positions, grid_positions, vectorfield, a, b)
    402 npart = len(xsi)
    404 _xlinear = XLinear()
--> 405 u = _xlinear.interp(particle_positions, grid_positions, vectorfield.U)
    406 v = _xlinear.interp(particle_positions, grid_positions, vectorfield.V)
    407 if vectorfield.W:

File ~/parcels/src/parcels/interpolators/_xinterpolators.py:133, in XLinear.interp(self, particle_positions, grid_positions, field)
    130 lenT = 2 if np.any(tau > 0) else 1
    131 lenZ = 2 if np.any(zeta > 0) else 1
--> 133 corner_data = _get_corner_data_Agrid(data, ti, zi, yi, xi, lenT, lenZ, len(xsi), axis_dim)
    135 if lenT == 2:
    136     tau = tau[np.newaxis, :]

File ~/parcels/src/parcels/interpolators/_xinterpolators.py:96, in _get_corner_data_Agrid(data, ti, zi, yi, xi, lenT, lenZ, npart, axis_dim)
     89 """Helper function to get the corner data for a given A-grid field and position."""
     90 levels: dict[ptyping.XgcmAxisDirection, tuple[np.ndarray, ...]] = {
     91     "T": (ti,) if lenT == 1 else (ti, np.clip(ti + 1, 0, data.shape[0] - 1)),
     92     "Z": (zi,) if lenZ == 1 else (zi, np.clip(zi + 1, 0, data.shape[1] - 1)),
     93     "Y": (yi, np.clip(yi + 1, 0, data.shape[2] - 1)),
     94     "X": (xi, np.clip(xi + 1, 0, data.shape[3] - 1)),
     95 }
---> 96 return _gather_corners(data, axis_dim, levels, npart)

File ~/parcels/src/parcels/interpolators/_xinterpolators.py:75, in _gather_corners(data, axis_dim, levels, npart)
     72         in_slot_i = np.expand_dims(stacked, other_slots)
     73         selection_dict[dims[axis]] = xr.DataArray(np.broadcast_to(in_slot_i, shape).reshape(-1), dims="points")
---> 75 return data.isel(selection_dict).data.reshape(shape)

File ~/PortOfRotterdam/.pixi/envs/default/lib/python3.14/site-packages/xarray/core/dataarray.py:1586, in DataArray.isel(self, indexers, drop, missing_dims, **indexers_kwargs)
   1583 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
   1585 if any(is_fancy_indexer(idx) for idx in indexers.values()):
-> 1586     ds = self._to_temp_dataset()._isel_fancy(
   1587         indexers, drop=drop, missing_dims=missing_dims
   1588     )
   1589     return self._from_temp_dataset(ds)
   1591 # Much faster algorithm for when all indexers are ints, slices, one-dimensional
   1592 # lists, or zero or one-dimensional np.ndarray's

File ~/PortOfRotterdam/.pixi/envs/default/lib/python3.14/site-packages/xarray/core/dataset.py:2976, in Dataset._isel_fancy(self, indexers, drop, missing_dims)
   2972                 var_indexers = {
   2973                     k: v for k, v in valid_indexers.items() if k in var.dims
   2974                 }
   2975                 if var_indexers:
-> 2976                     new_var = var.isel(indexers=var_indexers)
   2977                     # drop scalar coordinates
   2978                     # https://github.com/pydata/xarray/issues/6554
   2979                     if name in self.coords and drop and new_var.ndim == 0:

File ~/PortOfRotterdam/.pixi/envs/default/lib/python3.14/site-packages/xarray/core/variable.py:1144, in Variable.isel(self, indexers, missing_dims, **indexers_kwargs)
   1141 indexers = drop_dims_from_indexers(indexers, self.dims, missing_dims)
   1143 key = tuple(indexers.get(dim, slice(None)) for dim in self.dims)
-> 1144 return self[key]

File ~/PortOfRotterdam/.pixi/envs/default/lib/python3.14/site-packages/xarray/core/variable.py:831, in Variable.__getitem__(self, key)
    828 dims, indexer, new_order = self._broadcast_indexes(key)
    829 indexable = as_indexable(self._data)
--> 831 data = indexing.apply_indexer(indexable, indexer)
    833 if new_order:
    834     data = duck_array_ops.moveaxis(data, range(len(new_order)), new_order)

File ~/PortOfRotterdam/.pixi/envs/default/lib/python3.14/site-packages/xarray/core/indexing.py:1188, in apply_indexer(indexable, indexer)
   1186 """Apply an indexer to an indexable object."""
   1187 if isinstance(indexer, VectorizedIndexer):
-> 1188     return indexable.vindex[indexer]
   1189 elif isinstance(indexer, OuterIndexer):
   1190     return indexable.oindex[indexer]

File ~/PortOfRotterdam/.pixi/envs/default/lib/python3.14/site-packages/xarray/core/indexing.py:472, in IndexCallable.__getitem__(self, key)
    471 def __getitem__(self, key: Any) -> Any:
--> 472     return self.getter(key)

File ~/parcels/src/parcels/_chunk_cached_array/core.py:143, in ChunkCachedArray._vindex_get(self, indexer)
    141 def _vindex_get(self, indexer: VectorizedIndexer):
    142     key = indexer.tuple
--> 143     return self._raw_vindex(*key)

File ~/parcels/src/parcels/_chunk_cached_array/core.py:99, in ChunkCachedArray._raw_vindex(self, *indices)
     94 n_points = len(indices[0])
     96 # Step 1: Map global indices to chunk coords and local indices.
     97 # Normalize negative indices (e.g. -1 → last element) to positive,
     98 # matching standard numpy fancy-indexing semantics.
---> 99 indices = tuple(np.where(idx < 0, idx + self.array.shape[d], idx) for d, idx in enumerate(indices))
    100 chunk_ids = np.empty((ndim, n_points), dtype=np.intp)
    101 local_indices = np.empty((ndim, n_points), dtype=np.intp)

File ~/parcels/src/parcels/_chunk_cached_array/core.py:99, in <genexpr>(.0)
     94 n_points = len(indices[0])
     96 # Step 1: Map global indices to chunk coords and local indices.
     97 # Normalize negative indices (e.g. -1 → last element) to positive,
     98 # matching standard numpy fancy-indexing semantics.
---> 99 indices = tuple(np.where(idx < 0, idx + self.array.shape[d], idx) for d, idx in enumerate(indices))
    100 chunk_ids = np.empty((ndim, n_points), dtype=np.intp)
    101 local_indices = np.empty((ndim, n_points), dtype=np.intp)

TypeError: '<' not supported between instances of 'slice' and 'int'

The dataset is created with

ds = xr.open_mfdataset(
    "simona/fews_rmm_osr_*.nc",
    combine="by_coords",
    preprocess=lambda dataset: dataset.squeeze("analysis_time", drop=True),
    chunks={"time": 1, "row": 256, "col": 256},
)
ds = ds.rename({'lon': 'XZETA', 'lat': 'YZETA', 'row': 'N', 'col': 'M'})
ds = ds.isel(N=slice(800, 1500), M=slice(100, 200))
ds

which gives the following display(ds)

[VeckoTheGecko]

If you put a breakpoint at src/parcels/interpolators/_xinterpolators.py:96 , what are the values of data, axis_dim, levels, npart

And for src/parcels/interpolators/_xinterpolators.py:75 what is the value of selection_dict?

[erikvansebille]

CoPilot/Claude came up with the two-line solution in https://github.com/Parcels-code/Parcels/pull/2898/changes. It fixes the error in my case. Do you think this is indeed a robust solution?

[VeckoTheGecko]

Do you think this is indeed a robust solution?

No, I don't think so.

I don't think falling back to Dask is a viable solution

[erikvansebille]

OK thanks; I've closed #2898. I'll try to get you the values of these arrays next

[erikvansebille]

If you put a breakpoint at src/parcels/interpolators/_xinterpolators.py:96 , what are the values of data, axis_dim, levels, npart

And for src/parcels/interpolators/_xinterpolators.py:75 what is the value of selection_dict?

Here are the outputs

DATA:
<xarray.DataArray 'U' (time: 361, mockZ: 1, N: 700, M: 100)> Size: 101MB
[25270000 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 3kB 2026-06-30T23:00:00 ... 2026-07-03T11:...
  * N        (N) float64 6kB 800.0 801.0 802.0 ... 1.497e+03 1.498e+03 1.499e+03
  * M        (M) float64 800B 100.0 101.0 102.0 103.0 ... 197.0 198.0 199.0
    lat      (N, M) float64 560kB nan nan nan nan nan ... nan nan nan nan nan
    lon      (N, M) float64 560kB nan nan nan nan nan ... nan nan nan nan nan
Dimensions without coordinates: mockZ
Attributes:
    long_name:     eastward_water_velocity_surface
    units:         m/s
    grid_mapping:  crs
AXIS DIM:
{'Y': 'N', 'X': 'M'}
LEVELS:
{'T': (array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0]),), 'Z': (array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0]),), 'Y': (array([549, 550, 548, 545, 542, 556,  -3,  -3, 552,  -3, 563, 566, 570,
       574,  -3, 569, 572, 576,  -3, 583, 574, 579, 582, 586,  -3],
      dtype=int32), array([550, 551, 549, 546, 543, 557,   0,   0, 553,   0, 564, 567, 571,
       575,   0, 570, 573, 577,   0, 584, 575, 580, 583, 587,   0],
      dtype=int32)), 'X': (array([56, 67, 80, 88, 94, 44, -3, -3, 93, -3, 35, 40, 45, 49, -3, 28, 33,
       38, -3, 45, 21, 27, 32, 36, -3], dtype=int32), array([57, 68, 81, 89, 95, 45,  0,  0, 94,  0, 36, 41, 46, 50,  0, 29, 34,
       39,  0, 46, 22, 28, 33, 37,  0], dtype=int32))}
NPART:
25
SELECTION DICT:
{'time': <xarray.DataArray (points: 100)> Size: 800B
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
Dimensions without coordinates: points, 'N': <xarray.DataArray (points: 100)> Size: 400B
array([549, 550, 548, 545, 542, 556,  -3,  -3, 552,  -3, 563, 566, 570,
       574,  -3, 569, 572, 576,  -3, 583, 574, 579, 582, 586,  -3, 549,
       550, 548, 545, 542, 556,  -3,  -3, 552,  -3, 563, 566, 570, 574,
        -3, 569, 572, 576,  -3, 583, 574, 579, 582, 586,  -3, 550, 551,
       549, 546, 543, 557,   0,   0, 553,   0, 564, 567, 571, 575,   0,
       570, 573, 577,   0, 584, 575, 580, 583, 587,   0, 550, 551, 549,
       546, 543, 557,   0,   0, 553,   0, 564, 567, 571, 575,   0, 570,
       573, 577,   0, 584, 575, 580, 583, 587,   0], dtype=int32)
Dimensions without coordinates: points, 'M': <xarray.DataArray (points: 100)> Size: 400B
array([56, 67, 80, 88, 94, 44, -3, -3, 93, -3, 35, 40, 45, 49, -3, 28, 33,
       38, -3, 45, 21, 27, 32, 36, -3, 57, 68, 81, 89, 95, 45,  0,  0, 94,
        0, 36, 41, 46, 50,  0, 29, 34, 39,  0, 46, 22, 28, 33, 37,  0, 56,
       67, 80, 88, 94, 44, -3, -3, 93, -3, 35, 40, 45, 49, -3, 28, 33, 38,
       -3, 45, 21, 27, 32, 36, -3, 57, 68, 81, 89, 95, 45,  0,  0, 94,  0,
       36, 41, 46, 50,  0, 29, 34, 39,  0, 46, 22, 28, 33, 37,  0],
      dtype=int32)
Dimensions without coordinates: points}

Note that the -3 values are expected, as some particles are outside of the domain. But that shouldn't be the issue?

Also not that this full code works when using to_windowed_arrays() instead of to_cache_chunked_arrays()