Merge pull request #1999 from OceanParcels/remove-time-extrapolation-option

Remove allow_time_extrapolation attribute on Field

Author: fluidnumericsJoe

parcels/_index_search.py

@@ -25,7 +25,7 @@
     # from .grid import Grid
 
 
-def _search_time_index(field: Field, time: datetime, allow_time_extrapolation=True):
+def _search_time_index(field: Field, time: datetime):
     """Find and return the index and relative coordinate in the time array associated with a given time.
 
     Parameters
@@ -37,8 +37,12 @@ def _search_time_index(field: Field, time: datetime, allow_time_extrapolation=Tr
     Note that we normalize to either the first or the last index
     if the sampled value is outside the time value range.
     """
-    if not allow_time_extrapolation and (time < field.data.time[0] or time > field.data.time[-1]):
+    if field.time_interval is None:
+        return 0
+
+    if time in field.time_interval:
         _raise_time_extrapolation_error(time, field=None)
+
     time_index = field.data.time <= time
 
     if time_index.all():

parcels/field.py

@@ -148,7 +148,6 @@ def __init__(
         grid: ux.Grid | Grid,
         mesh_type: Mesh = "flat",
         interp_method: Callable | None = None,
-        allow_time_extrapolation: bool | None = None,
     ):
         if not isinstance(data, (ux.UxDataArray, xr.DataArray)):
             raise ValueError(
@@ -209,10 +208,8 @@ def __init__(
         else:
             raise ValueError("Unsupported mesh type in data array attributes. Choose either: 'spherical' or 'flat'")
 
-        if allow_time_extrapolation is None:
-            self.allow_time_extrapolation = True if len(getattr(self.data, "time", [])) == 1 else False
-        else:
-            self.allow_time_extrapolation = allow_time_extrapolation
+        if "time" not in self.data.dims:
+            raise ValueError("Field is missing a 'time' dimension. ")
 
     def __repr__(self):
         return field_repr(self)
@@ -366,7 +363,7 @@ def _search_indices_structured(self, z, y, x, ei=None, search2D=False):
         return (zeta, eta, xsi, zi, yi, xi)
 
     def _search_indices(self, time: datetime, z, y, x, ei=None, search2D=False):
-        tau, ti = _search_time_index(self, time, self.allow_time_extrapolation)
+        tau, ti = _search_time_index(self, time)
 
         if ei is None:
             _ei = None
@@ -679,7 +676,7 @@ def _assert_compatible_combination(data: xr.DataArray | ux.UxDataArray, grid: ux
 
 
 def get_time_interval(data: xr.DataArray | ux.UxDataArray) -> TimeInterval | None:
-    if "time" not in data.dims:
+    if len(data.time) == 1:
         return None
 
     return TimeInterval(data.time.values[0], data.time.values[-1])

parcels/fieldset.py

@@ -164,14 +164,15 @@ def add_constant_field(self, name: str, value, mesh: Mesh = "flat"):
         """
         da = xr.DataArray(
             data=np.full((1, 1, 1, 1), value),
-            dims=["T", "ZG", "YG", "XG"],
+            dims=["time", "ZG", "YG", "XG"],
             coords={
                 "ZG": (["ZG"], np.arange(1), {"axis": "Z"}),
                 "YG": (["YG"], np.arange(1), {"axis": "Y"}),
                 "XG": (["XG"], np.arange(1), {"axis": "X"}),
                 "lon": (["XG"], np.arange(1), {"axis": "X"}),
                 "lat": (["YG"], np.arange(1), {"axis": "Y"}),
                 "depth": (["ZG"], np.arange(1), {"axis": "Z"}),
+                "time": (["time"], np.arange(1), {"axis": "T"}),
             },
         )
         grid = Grid(da)
@@ -253,6 +254,10 @@ class CalendarError(Exception):  # TODO: Move to a parcels errors module
 
 def assert_compatible_calendars(fields: Iterable[Field]):
     time_intervals = [f.time_interval for f in fields if f.time_interval is not None]
+
+    if len(time_intervals) == 0:  # All time intervals are none
+        return
+
     reference_datetime_object = time_intervals[0].left
 
     for field in fields:

parcels/particleset.py

@@ -11,9 +11,9 @@
 from tqdm import tqdm
 
 from parcels._compat import MPI
+from parcels._core.utils.time import TimeInterval
 from parcels._reprs import particleset_repr
 from parcels.application_kernels.advection import AdvectionRK4
-from parcels.field import Field
 from parcels.grid import GridType
 from parcels.interaction.interactionkernel import InteractionKernel
 from parcels.interaction.neighborsearch import (
@@ -162,8 +162,8 @@ def ArrayClass_init(self, *args, **kwargs):
             raise NotImplementedError("If fieldset.time_origin is not a date, time of a particle must be a double")
         time = np.array([self.time_origin.reltime(t) if _convert_to_reltime(t) else t for t in time])
         assert lon.size == time.size, "time and positions (lon, lat, depth) do not have the same lengths."
-        if isinstance(fieldset.U, Field) and (not fieldset.U.allow_time_extrapolation):
-            _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.U.grid.time)
+        if fieldset.time_interval:
+            _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.time_interval)
 
         if lonlatdepth_dtype is None:
             lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U)
@@ -1127,17 +1127,17 @@ def _warn_outputdt_release_desync(outputdt: float, starttime: float, release_tim
         )
 
 
-def _warn_particle_times_outside_fieldset_time_bounds(release_times: np.ndarray, time: np.ndarray):
+def _warn_particle_times_outside_fieldset_time_bounds(release_times: np.ndarray, time: np.ndarray | TimeInterval):
     if np.any(release_times):
         if np.any(release_times < time[0]):
             warnings.warn(
-                "Some particles are set to be released before the fieldset's first time and allow_time_extrapolation is set to False.",
+                "Some particles are set to be released outside the FieldSet's executable time domain.",
                 ParticleSetWarning,
                 stacklevel=2,
             )
         if np.any(release_times > time[-1]):
             warnings.warn(
-                "Some particles are set to be released after the fieldset's last time and allow_time_extrapolation is set to False.",
+                "Some particles are set to be released after the fieldset's last time and the fields are not constant in time.",
                 ParticleSetWarning,
                 stacklevel=2,
             )