Repository: OceanParcels/Parcels
Pull request #1891 (2025-02-27 12:15:59)

Removing particle ti
This PR removed the `particle.ti` `Variable` altogether, as it turned out not to be used anywhere in the code. 

It also cleans up the `search_indices` method by 
1. Not setting to default `ti=-1` anymore (as not necessary, and confusing)
2. Reordering arguments list to align with `(time, depth, lat, lon)` order throughout Parcels

This PR supersedes #1881 

- [x] Chose the correct base branch (`main` for v3 changes, `v4-dev` for v4 changes)



Patch:
diff --git a/parcels/_index_search.py b/parcels/_index_search.py
index 1b7e7788..a52baf56 100644
--- a/parcels/_index_search.py
+++ b/parcels/_index_search.py
@@ -137,7 +137,7 @@ def search_indices_vertical_s(
 
 
 def _search_indices_rectilinear(
-    field: Field, time: float, z: float, y: float, x: float, ti=-1, particle=None, search2D=False
+    field: Field, time: float, z: float, y: float, x: float, ti: int, particle=None, search2D=False
 ):
     grid = field.grid
 
@@ -223,7 +223,7 @@ def _search_indices_rectilinear(
     return (zeta, eta, xsi, zi, yi, xi)
 
 
-def _search_indices_curvilinear(field: Field, time, z, y, x, ti=-1, particle=None, search2D=False):
+def _search_indices_curvilinear(field: Field, time, z, y, x, ti, particle=None, search2D=False):
     if particle:
         zi, yi, xi = field.unravel_index(particle.ei)
     else:
diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py
index da62bf89..d35d3063 100644
--- a/parcels/application_kernels/advection.py
+++ b/parcels/application_kernels/advection.py
@@ -185,7 +185,7 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
         ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])
 
     zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices(
-        -1, particle.depth, particle.lat, particle.lon, particle=particle
+        time, particle.depth, particle.lat, particle.lon, ti, particle=particle
     )
     if withW:
         if abs(xsi - 1) < tol:
diff --git a/parcels/field.py b/parcels/field.py
index 6467157e..374b4f9e 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -877,15 +877,15 @@ class Field:
         """
         return _calc_cell_areas(self.grid)
 
-    def _search_indices(self, time, z, y, x, ti=-1, particle=None, search2D=False):
+    def _search_indices(self, time, z, y, x, ti, particle=None, search2D=False):
         if self.grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
             return _search_indices_rectilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D)
         else:
             return _search_indices_curvilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D)
 
-    def _interpolator2D(self, ti, z, y, x, particle=None):
+    def _interpolator2D(self, time, z, y, x, ti, particle=None):
         """Impelement 2D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019.."""
-        (_, eta, xsi, _, yi, xi) = self._search_indices(-1, z, y, x, particle=particle)
+        (_, eta, xsi, _, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle)
         ctx = InterpolationContext2D(self.data, eta, xsi, ti, yi, xi)
 
         try:
@@ -899,7 +899,7 @@ class Field:
                 raise RuntimeError(self.interp_method + " is not implemented for 2D grids")
         return f(ctx)
 
-    def _interpolator3D(self, ti, z, y, x, time, particle=None):
+    def _interpolator3D(self, time, z, y, x, ti, particle=None):
         """Impelement 3D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019.."""
         (zeta, eta, xsi, zi, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle)
         ctx = InterpolationContext3D(self.data, zeta, eta, xsi, ti, zi, yi, xi, self.gridindexingtype)
@@ -910,13 +910,13 @@ class Field:
             raise RuntimeError(self.interp_method + " is not implemented for 3D grids")
         return f(ctx)
 
-    def _spatial_interpolation(self, ti, z, y, x, time, particle=None):
-        """Interpolate horizontal field values using a SciPy interpolator."""
+    def _spatial_interpolation(self, time, z, y, x, ti, particle=None):
+        """Interpolate spatial field values."""
         try:
             if self.grid.zdim == 1:
-                val = self._interpolator2D(ti, z, y, x, particle=particle)
+                val = self._interpolator2D(time, z, y, x, ti, particle=particle)
             else:
-                val = self._interpolator3D(ti, z, y, x, time, particle=particle)
+                val = self._interpolator3D(time, z, y, x, ti, particle=particle)
 
             if np.isnan(val):
                 # Detect Out-of-bounds sampling and raise exception
@@ -980,8 +980,8 @@ class Field:
         if self.gridindexingtype == "croco" and self not in [self.fieldset.H, self.fieldset.Zeta]:
             z = _croco_from_z_to_sigma_scipy(self.fieldset, time, z, y, x, particle=particle)
         if ti < self.grid.tdim - 1 and time > self.grid.time[ti]:
-            f0 = self._spatial_interpolation(ti, z, y, x, time, particle=particle)
-            f1 = self._spatial_interpolation(ti + 1, z, y, x, time, particle=particle)
+            f0 = self._spatial_interpolation(time, z, y, x, ti, particle=particle)
+            f1 = self._spatial_interpolation(time, z, y, x, ti + 1, particle=particle)
             t0 = self.grid.time[ti]
             t1 = self.grid.time[ti + 1]
             value = f0 + (f1 - f0) * ((time - t0) / (t1 - t0))
@@ -989,7 +989,7 @@ class Field:
             # Skip temporal interpolation if time is outside
             # of the defined time range or if we have hit an
             # exact value in the time array.
-            value = self._spatial_interpolation(ti, z, y, x, self.grid.time[ti], particle=particle)
+            value = self._spatial_interpolation(self.grid.time[ti], z, y, x, ti, particle=particle)
 
         if applyConversion:
             return self.units.to_target(value, z, y, x)
diff --git a/parcels/particledata.py b/parcels/particledata.py
index 8a9e69c1..2aed1ee6 100644
--- a/parcels/particledata.py
+++ b/parcels/particledata.py
@@ -122,8 +122,8 @@ class ParticleData:
         self._ncount = len(lon)
 
         for v in self.ptype.variables:
-            if v.name in ["ei", "ti"]:
-                self._data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype)
+            if v.name == "ei":
+                self._data[v.name] = np.empty((len(lon), ngrid), dtype=v.dtype)  # TODO len(lon) can be self._ncount?
             else:
                 self._data[v.name] = np.empty(self._ncount, dtype=v.dtype)
 
diff --git a/parcels/particleset.py b/parcels/particleset.py
index 0fa834b9..e6946e2d 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -128,13 +128,11 @@ class ParticleSet:
                 self.ngrids = type(self).ngrids.initial
                 if self.ngrids >= 0:
                     self.ei = np.zeros(self.ngrids, dtype=np.int32)
-                    self.ti = -1 * np.ones(self.ngrids, dtype=np.int32)
                 super(type(self), self).__init__(*args, **kwargs)
 
             array_class_vdict = {
                 "ngrids": Variable("ngrids", dtype=np.int32, to_write=False, initial=-1),
                 "ei": Variable("ei", dtype=np.int32, to_write=False),
-                "ti": Variable("ti", dtype=np.int32, to_write=False, initial=-1),
                 "__init__": ArrayClass_init,
             }
             array_class = type(class_name, (pclass,), array_class_vdict)
@@ -719,7 +717,6 @@ class ParticleSet:
                 v.name
                 not in [
                     "ei",
-                    "ti",
                     "dt",
                     "depth",
                     "id",


Test patch:
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index bab9acae..4ab303b9 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -571,7 +571,8 @@ def test_fieldset_write(tmp_zarrfile):
     def UpdateU(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[particle]
         _, yi, xi = fieldset.U.unravel_index(particle.ei)
-        fieldset.U.data[particle.ti, yi, xi] += 1
+        ti = fieldset.U._time_index(time)
+        fieldset.U.data[ti, yi, xi] += 1
         fieldset.U.grid.time[0] = time
 
     pset = ParticleSet(fieldset, pclass=Particle, lon=5, lat=5)