From 848788da652659e772590ab9846736a0b499bc91 Mon Sep 17 00:00:00 2001 From: Joe Date: Tue, 6 May 2025 09:38:55 -0400 Subject: [PATCH 1/8] Remove allow_time_extrapolation attribute on Field When a xarray or uxarray dataarray does not have a time dimension, one is added to the Field.data dataarray. When a field.data datarray has a time dimension that is length 1, we assume that extrapolation is allowed meaning that we just keep using the same single time-slice for that field. --- parcels/_index_search.py | 5 +++-- parcels/field.py | 13 +++++++------ parcels/fieldset.py | 1 - parcels/particleset.py | 6 +++--- 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/parcels/_index_search.py b/parcels/_index_search.py index 8af53bee27..ae8258da79 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,9 @@ 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 (len(field.data.time)) > 1 and (time < field.data.time[0] or time > field.data.time[-1]): _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 613a88ac82..a5323f2212 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -149,7 +149,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( @@ -204,10 +203,12 @@ 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 + # Check if time is in self.data.dims + if "time" not in self.data.dims: + # Add time dimension of length 1 if not present + # While the choice of actual date is arbitrary, it should be a datetime object + # Here, we use the current time + self.data = self.data.expand_dims({"time": [datetime.now()]}) def __repr__(self): return field_repr(self) @@ -361,7 +362,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 diff --git a/parcels/fieldset.py b/parcels/fieldset.py index 26543e9e8c..2775e2cafa 100644 --- a/parcels/fieldset.py +++ b/parcels/fieldset.py @@ -151,7 +151,6 @@ def add_constant_field(self, name: str, value, mesh: Mesh = "flat"): name, data, interp_method=None, # TODO : Need to define an interpolation method for constants - allow_time_extrapolation=True, ) ) diff --git a/parcels/particleset.py b/parcels/particleset.py index a2dd9f7a75..6750d245c2 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -161,7 +161,7 @@ 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): + if isinstance(fieldset.U, Field) and (len(fieldset.U.data.time) > 1): _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.U.grid.time) if lonlatdepth_dtype is None: @@ -1130,13 +1130,13 @@ def _warn_particle_times_outside_fieldset_time_bounds(release_times: np.ndarray, 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 before the fieldset's first time and the fields are not constant in time.", 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, ) From d2417f7b13d07ec7e2aeeb24767ec13dae76dc36 Mon Sep 17 00:00:00 2001 From: Joe Schoonover <11430768+fluidnumerics-joe@users.noreply.github.com> Date: Tue, 6 May 2025 11:53:10 -0600 Subject: [PATCH 2/8] Update parcels/particleset.py Co-authored-by: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> --- parcels/particleset.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 6750d245c2..47e896c14d 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -1130,7 +1130,7 @@ def _warn_particle_times_outside_fieldset_time_bounds(release_times: np.ndarray, 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 the fields are not constant in time.", + "Some particles are set to be released outside the FieldSet's executable time domain.", ParticleSetWarning, stacklevel=2, ) From 14ef097f8c45a7d50fbdba25b3b7d61806e0288b Mon Sep 17 00:00:00 2001 From: Joe Date: Wed, 7 May 2025 07:47:52 -0400 Subject: [PATCH 3/8] Switch to using time_interval for conditional checks on current time --- parcels/_index_search.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/_index_search.py b/parcels/_index_search.py index ae8258da79..9716857fb0 100644 --- a/parcels/_index_search.py +++ b/parcels/_index_search.py @@ -37,7 +37,7 @@ def _search_time_index(field: Field, time: datetime): Note that we normalize to either the first or the last index if the sampled value is outside the time value range. """ - if (len(field.data.time)) > 1 and (time < field.data.time[0] or time > field.data.time[-1]): + if field.time_interval is not None and time in field.time_interval: _raise_time_extrapolation_error(time, field=None) time_index = field.data.time <= time From 2afafa0740efdc42e19ac66575b898f1a60ac721 Mon Sep 17 00:00:00 2001 From: Joe Date: Wed, 7 May 2025 07:55:31 -0400 Subject: [PATCH 4/8] Change to using fieldset.time_interval for time out of bounds check --- parcels/particleset.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index f9e1b88b04..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 (len(fieldset.U.data.time) > 1): - _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,7 +1127,7 @@ 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( From 80eba8860e7ce7d76918068c92cde64121308ce9 Mon Sep 17 00:00:00 2001 From: Joe Date: Wed, 7 May 2025 07:58:14 -0400 Subject: [PATCH 5/8] When time is not provided, use hardset date The beginning of time when not provided is the date of the first parcels commit --- parcels/field.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/parcels/field.py b/parcels/field.py index 87edb80b6a..02e8acaace 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -206,8 +206,7 @@ def __init__( if "time" not in self.data.dims: # Add time dimension of length 1 if not present # While the choice of actual date is arbitrary, it should be a datetime object - # Here, we use the current time - self.data = self.data.expand_dims({"time": [datetime.now()]}) + self.data = self.data.expand_dims({"time": [datetime.strptime("2015-09-27", "%Y-%m-%d")]}) def __repr__(self): return field_repr(self) From 9697eca73094e83872c01a78debc82aaff839a19 Mon Sep 17 00:00:00 2001 From: Joe Date: Wed, 7 May 2025 09:03:32 -0400 Subject: [PATCH 6/8] Require time dimension in field.data.dims --- parcels/_index_search.py | 5 ++++- parcels/field.py | 6 ++---- parcels/fieldset.py | 3 ++- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/parcels/_index_search.py b/parcels/_index_search.py index 9716857fb0..d8548d2642 100644 --- a/parcels/_index_search.py +++ b/parcels/_index_search.py @@ -37,7 +37,10 @@ def _search_time_index(field: Field, time: datetime): Note that we normalize to either the first or the last index if the sampled value is outside the time value range. """ - if field.time_interval is not None and time in field.time_interval: + 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 diff --git a/parcels/field.py b/parcels/field.py index 02e8acaace..63dfbf9514 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -204,9 +204,7 @@ def __init__( # Check if time is in self.data.dims if "time" not in self.data.dims: - # Add time dimension of length 1 if not present - # While the choice of actual date is arbitrary, it should be a datetime object - self.data = self.data.expand_dims({"time": [datetime.strptime("2015-09-27", "%Y-%m-%d")]}) + raise ValueError("Field is missing a 'time' dimension. ") def __repr__(self): return field_repr(self) @@ -666,7 +664,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 ba76526794..8cf73b659d 100644 --- a/parcels/fieldset.py +++ b/parcels/fieldset.py @@ -154,7 +154,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"}), @@ -162,6 +162,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) From 699d4776be4a3b95a4f0d1a00858c248ccf515cf Mon Sep 17 00:00:00 2001 From: Joe Schoonover <11430768+fluidnumerics-joe@users.noreply.github.com> Date: Wed, 7 May 2025 11:54:02 -0600 Subject: [PATCH 7/8] Update parcels/field.py Co-authored-by: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> --- parcels/field.py | 1 - 1 file changed, 1 deletion(-) diff --git a/parcels/field.py b/parcels/field.py index bc79c507e1..0484f30e90 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -208,7 +208,6 @@ def __init__( else: raise ValueError("Unsupported mesh type in data array attributes. Choose either: 'spherical' or 'flat'") - # Check if time is in self.data.dims if "time" not in self.data.dims: raise ValueError("Field is missing a 'time' dimension. ") From 1d44b691fc8ab417b111ab58f8cfaa7297382743 Mon Sep 17 00:00:00 2001 From: Joe Date: Tue, 13 May 2025 09:18:47 -0400 Subject: [PATCH 8/8] Gaurd against the case when all fields have length 1 time dimension --- parcels/fieldset.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/parcels/fieldset.py b/parcels/fieldset.py index 873fedae0a..9f9393d81d 100644 --- a/parcels/fieldset.py +++ b/parcels/fieldset.py @@ -254,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: