Adding ti index to the ravelled particle index too

.readthedocs.yaml

@@ -3,7 +3,9 @@ build:
   os: ubuntu-22.04
   tools:
     python: mambaforge-22.9
-
+  jobs:
+    post_create_environment:
+      - pip install -e .
 sphinx:
   configuration: docs/conf.py
 

docs/examples/example_globcurrent.py

@@ -13,7 +13,6 @@ def set_globcurrent_fieldset(
     indices=None,
     deferred_load=True,
     use_xarray=False,
-    time_periodic=False,
     timestamps=None,
 ):
     if filename is None:
@@ -32,7 +31,9 @@ def set_globcurrent_fieldset(
     if use_xarray:
         ds = xr.open_mfdataset(filename, combine="by_coords")
         return parcels.FieldSet.from_xarray_dataset(
-            ds, variables, dimensions, time_periodic=time_periodic
+            ds,
+            variables,
+            dimensions,
         )
     else:
         return parcels.FieldSet.from_netcdf(
@@ -41,7 +42,6 @@ def set_globcurrent_fieldset(
             dimensions,
             indices,
             deferred_load=deferred_load,
-            time_periodic=time_periodic,
             timestamps=timestamps,
         )
 
@@ -118,6 +118,8 @@ def test_globcurrent_particles(use_xarray):
     assert abs(pset[0].lat - -35.3) < 1
 
 
+@pytest.mark.v4remove
+@pytest.mark.xfail(reason="time_periodic removed in v4")
 @pytest.mark.parametrize("rundays", [300, 900])
 def test_globcurrent_time_periodic(rundays):
     sample_var = []

docs/examples/tutorial_particle_field_interaction.ipynb

@@ -162,7 +162,7 @@
     "    deltaC = fieldset.a * fieldset.C[particle] - fieldset.b * particle.c\n",
     "\n",
     "    ei = fieldset.C.unravel_index(particle.ei)\n",
-    "    xi, yi = ei[2], ei[1]\n",
+    "    yi, xi = ei[2], ei[3]\n",
     "\n",
     "    if abs(particle.lon - fieldset.C.grid.lon[xi + 1]) < abs(\n",
     "        particle.lon - fieldset.C.grid.lon[xi]\n",

docs/examples/tutorial_timestamps.ipynb

@@ -155,14 +155,6 @@
     "    )"
    ]
   },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Note, by the way, that adding the `time_periodic` argument to `Field.from_netcdf()` will also mean that the climatology can be cycled for multiple years.\n"
-   ]
-  },
   {
    "attachments": {},
    "cell_type": "markdown",

docs/v4/TODO.md

@@ -5,3 +5,5 @@ List of tasks that are important to do before the release of version 4 (but can'
 - [ ] Make migration guide for v3 to v4
 - [ ] Just prior to release: Update conda feedstock recipe dependencies (remove cgen and compiler dependencies). Make sure that recipe is up-to-date.
 - [ ] Revamp the oceanparcels.org landing page, and perhaps also consider new logo/branding?
+- [ ] Look into xarray and whether users can create periodic datasets without increasing the size of the original dataset (i.e., no compromise alternative to `time_periodic` param in v3). Update docs accordingly.
+- [ ] Rerun all the tutorials so that their output is in line with new v4 print statements etc

parcels/_index_search.py

@@ -218,14 +218,14 @@
         _raise_field_sampling_error(z, y, x)
 
     if particle:
-        particle.ei[field.igrid] = field.ravel_index(zi, yi, xi)
+        particle.ei[field.igrid] = field.ravel_index(ti, zi, yi, xi)
 
     return (zeta, eta, xsi, zi, yi, xi)
 
 
 def _search_indices_curvilinear(field: Field, time, z, y, x, ti=-1, particle=None, search2D=False):
     if particle:
-        zi, yi, xi = field.unravel_index(particle.ei)
+        ti, zi, yi, xi = field.unravel_index(particle.ei)
     else:
         xi = int(field.grid.xdim / 2) - 1
         yi = int(field.grid.ydim / 2) - 1
@@ -307,7 +307,7 @@
         _raise_field_sampling_error(z, y, x)
 
     if particle:
-        particle.ei[field.igrid] = field.ravel_index(zi, yi, xi)
+        particle.ei[field.igrid] = field.ravel_index(ti, zi, yi, xi)
 
     return (zeta, eta, xsi, zi, yi, xi)
 

parcels/_typing.py

@@ -6,7 +6,6 @@
 
 """
 
-import datetime
 import os
 from collections.abc import Callable
 from typing import Any, Literal, get_args
@@ -32,7 +31,6 @@
 ChunkMode = Literal["auto", "specific", "failsafe"]  # corresponds with `chunk_mode`
 GridIndexingType = Literal["pop", "mom5", "mitgcm", "nemo", "croco"]  # corresponds with `gridindexingtype`
 UpdateStatus = Literal["not_updated", "first_updated", "updated"]  # corresponds with `_update_status`
-TimePeriodic = float | datetime.timedelta | Literal[False]  # corresponds with `time_periodic`
 NetcdfEngine = Literal["netcdf4", "xarray", "scipy"]
 
 

parcels/application_kernels/advection.py

@@ -210,7 +210,7 @@
                 yi += 1
                 eta = 0
 
-    particle.ei[:] = fieldset.U.ravel_index(zi, yi, xi)
+    particle.ei[:] = fieldset.U.ravel_index(ti, zi, yi, xi)
 
     grid = fieldset.U.grid
     if grid._gtype < 2:

parcels/field.py

@@ -22,12 +22,11 @@
     GridIndexingType,
     InterpMethod,
     Mesh,
-    TimePeriodic,
     VectorType,
     assert_valid_gridindexingtype,
     assert_valid_interp_method,
 )
-from parcels.tools._helpers import default_repr, field_repr, timedelta_to_float
+from parcels.tools._helpers import default_repr, field_repr
 from parcels.tools.converters import (
     TimeConverter,
     UnitConverter,
@@ -155,10 +154,6 @@
     allow_time_extrapolation : bool
         boolean whether to allow for extrapolation in time
         (i.e. beyond the last available time snapshot)
-    time_periodic : bool, float or datetime.timedelta
-        To loop periodically over the time component of the Field. It is set to either False or the length of the period (either float in seconds or datetime.timedelta object).
-        The last value of the time series can be provided (which is the same as the initial one) or not (Default: False)
-        This flag overrides the allow_time_extrapolation and sets it to False
     chunkdims_name_map : str, optional
         Gives a name map to the FieldFileBuffer that declared a mapping between chunksize name, NetCDF dimension and Parcels dimension;
         required only if currently incompatible OCM field is loaded and chunking is used by 'chunksize' (which is the default)
@@ -174,7 +169,6 @@
     """
 
     allow_time_extrapolation: bool
-    time_periodic: TimePeriodic
     _cast_data_dtype: type[np.float32] | type[np.float64]
 
     def __init__(
@@ -196,7 +190,6 @@
         time_origin: TimeConverter | None = None,
         interp_method: InterpMethod = "linear",
         allow_time_extrapolation: bool | None = None,
-        time_periodic: TimePeriodic = False,
         gridindexingtype: GridIndexingType = "nemo",
         to_write: bool = False,
         **kwargs,
@@ -261,28 +254,6 @@
         else:
             self.allow_time_extrapolation = allow_time_extrapolation
 
-        self.time_periodic = time_periodic
-        if self.time_periodic is not False and self.allow_time_extrapolation:
-            warnings.warn(
-                "allow_time_extrapolation and time_periodic cannot be used together. allow_time_extrapolation is set to False",
-                FieldSetWarning,
-                stacklevel=2,
-            )
-            self.allow_time_extrapolation = False
-        if self.time_periodic is True:
-            raise ValueError(
-                "Unsupported time_periodic=True. time_periodic must now be either False or the length of the period (either float in seconds or datetime.timedelta object."
-            )
-        if self.time_periodic is not False:
-            self.time_periodic = timedelta_to_float(self.time_periodic)
-
-            if not np.isclose(self.grid.time[-1] - self.grid.time[0], self.time_periodic):
-                if self.grid.time[-1] - self.grid.time[0] > self.time_periodic:
-                    raise ValueError("Time series provided is longer than the time_periodic parameter")
-                self.grid._add_last_periodic_data_timestep = True
-                self.grid.time = np.append(self.grid.time, self.grid.time[0] + self.time_periodic)
-                self.grid.time_full = self.grid.time
-
         self.vmin = vmin
         self.vmax = vmax
 
@@ -312,16 +283,11 @@
             if self.vmax is not None:
                 self.data[self.data > self.vmax] = 0.0
 
-            if self.grid._add_last_periodic_data_timestep:
-                self.data = lib.concatenate((self.data, self.data[:1, :]), axis=0)
-
         self._scaling_factor = None
 
         self._dimensions = kwargs.pop("dimensions", None)
         self.indices = kwargs.pop("indices", None)
         self._dataFiles = kwargs.pop("dataFiles", None)
-        if self.grid._add_last_periodic_data_timestep and self._dataFiles is not None:
-            self._dataFiles = np.append(self._dataFiles, self._dataFiles[0])
         self._field_fb_class = kwargs.pop("FieldFileBuffer", None)
         self._netcdf_engine = kwargs.pop("netcdf_engine", "netcdf4")
         self._creation_log = kwargs.pop("creation_log", "")
@@ -457,7 +423,6 @@
         mesh: Mesh = "spherical",
         timestamps=None,
         allow_time_extrapolation: bool | None = None,
-        time_periodic: TimePeriodic = False,
         deferred_load: bool = True,
         **kwargs,
     ) -> "Field":
@@ -491,9 +456,6 @@
             boolean whether to allow for extrapolation in time
             (i.e. beyond the last available time snapshot)
             Default is False if dimensions includes time, else True
-        time_periodic : bool, float or datetime.timedelta
-            boolean whether to loop periodically over the time component of the FieldSet
-            This flag overrides the allow_time_extrapolation and sets it to False (Default value = False)
         deferred_load : bool
             boolean whether to only pre-load data (in deferred mode) or
             fully load them (default: True). It is advised to deferred load the data, since in
@@ -743,7 +705,6 @@
 
         kwargs["dimensions"] = dimensions.copy()
         kwargs["indices"] = indices
-        kwargs["time_periodic"] = time_periodic
         kwargs["netcdf_engine"] = netcdf_engine
 
         return cls(
@@ -764,7 +725,6 @@
         dimensions,
         mesh: Mesh = "spherical",
         allow_time_extrapolation: bool | None = None,
-        time_periodic: TimePeriodic = False,
         **kwargs,
     ):
         """Create field from xarray Variable.
@@ -788,9 +748,6 @@
             boolean whether to allow for extrapolation in time
             (i.e. beyond the last available time snapshot)
             Default is False if dimensions includes time, else True
-        time_periodic : bool, float or datetime.timedelta
-            boolean whether to loop periodically over the time component of the FieldSet
-            This flag overrides the allow_time_extrapolation and sets it to False (Default value = False)
         **kwargs :
             Keyword arguments passed to the :class:`Field` constructor.
         """
@@ -806,7 +763,6 @@
         time = time_origin.reltime(time)  # type: ignore[assignment]
 
         grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
-        kwargs["time_periodic"] = time_periodic
         return cls(
             name,
             data,
@@ -1001,24 +957,10 @@
         Note that we normalize to either the first or the last index
         if the sampled value is outside the time value range.
         """
-        if (
-            not self.time_periodic
-            and not self.allow_time_extrapolation
-            and (time < self.grid.time[0] or time > self.grid.time[-1])
-        ):
+        if not self.allow_time_extrapolation and (time < self.grid.time[0] or time > self.grid.time[-1]):
             raise TimeExtrapolationError(time, field=self)
         time_index = self.grid.time <= time
-        if self.time_periodic:
-            if time_index.all() or np.logical_not(time_index).all():
-                periods = int(
-                    math.floor((time - self.grid.time_full[0]) / (self.grid.time_full[-1] - self.grid.time_full[0]))
-                )
-                self.grid.periods = periods
-                time -= periods * (self.grid.time_full[-1] - self.grid.time_full[0])
-                time_index = self.grid.time <= time
-                ti = time_index.argmin() - 1 if time_index.any() else 0
-                return (ti, periods)
-            return (time_index.argmin() - 1 if time_index.any() else 0, 0)
+
         if time_index.all():
             # If given time > last known field time, use
             # the last field frame without interpolation
@@ -1313,11 +1255,13 @@
         self.filebuffers[tindex] = filebuffer
         return data
 
-    def ravel_index(self, zi, yi, xi):
+    def ravel_index(self, ti, zi, yi, xi):
         """Return the flat index of the given grid points.
 
         Parameters
         ----------
+        ti : int
+            time index
         zi : int
             z index
         yi : int
@@ -1330,10 +1274,10 @@
         int
             flat index
         """
-        return xi + self.grid.xdim * (yi + self.grid.ydim * zi)
+        return xi + self.grid.xdim * (yi + self.grid.ydim * (zi + self.grid.tdim * ti))
 
     def unravel_index(self, ei):
-        """Return the zi, yi, xi indices for a given flat index.
+        """Return the ti, zi, yi, xi indices for a given flat index.
 
         Parameters
         ----------
@@ -1342,6 +1286,8 @@
 
         Returns
         -------
+        ti : int
+            The time index.
         zi : int
             The z index.
         yi : int
@@ -1349,12 +1295,14 @@
         xi : int
             The x index.
         """
         _ei = ei[self.igrid]
+        ti = _ei // (self.grid.xdim * self.grid.ydim * self.grid.zdim)
+        _ei = abs(_ei) % (self.grid.xdim * self.grid.ydim * self.grid.zdim)
         zi = _ei // (self.grid.xdim * self.grid.ydim)
         _ei = _ei % (self.grid.xdim * self.grid.ydim)
         yi = _ei // self.grid.xdim
         xi = _ei % self.grid.xdim
-        return zi, yi, xi
+        return ti, zi, yi, xi
 
 
 class VectorField: