diff --git a/parcels/_index_search.py b/parcels/_index_search.py index 8af53bee27..d8548d2642 100644 --- a/parcels/_index_search.py +++ b/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(): diff --git a/parcels/field.py b/parcels/field.py index c0af36f4c9..0484f30e90 100644 --- a/parcels/field.py +++ b/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]) diff --git a/parcels/fieldset.py b/parcels/fieldset.py index 71affd6cf2..9f9393d81d 100644 --- a/parcels/fieldset.py +++ b/parcels/fieldset.py @@ -164,7 +164,7 @@ 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"}), @@ -172,6 +172,7 @@ def add_constant_field(self, name: str, value, mesh: Mesh = "flat"): "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: diff --git a/parcels/particleset.py b/parcels/particleset.py index fa09d48451..9fd5f65017 100644 --- a/parcels/particleset.py +++ b/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, )