Cleanign up of search_indices method

Not setting to default ti=-1 anymore (as not necessary, and confusing)

Also reordering arguments list to align with (time, depth, lat, lon) order throughout Parcels

Author: erikvansebille

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:

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:

parcels/field.py

@@ -877,15 +877,15 @@ def cell_areas(self):
         """
         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 @@ def _interpolator2D(self, ti, z, y, x, particle=None):
                 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)
@@ -931,13 +931,13 @@ def temporal_interpolate_fullfield(self, ti, time):
             f1 = self.data[ti + 1, :]
             return f0 + (f1 - f0) * ((time - t0) / (t1 - t0))
 
-    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
@@ -1002,16 +1002,16 @@ def eval(self, time, z, y, x, particle=None, applyConversion=True):
         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))
         else:
             # 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)