From 4cc5d0a82560be2473d601784e4b6bf40fe282c3 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 3 Jul 2025 15:33:58 +0200 Subject: [PATCH 1/4] Removing time_origin for v4 As time will become a datetime64?timedelta64 object, no need anymore for time_origin and TimeConverter --- parcels/particlefile.py | 9 +-- parcels/particleset.py | 10 +-- parcels/tools/converters.py | 118 --------------------------------- parcels/tools/statuscodes.py | 15 +---- parcels/xgrid.py | 11 --- tests/tools/test_converters.py | 55 --------------- tests/v4/test_xgrid.py | 3 - 7 files changed, 4 insertions(+), 217 deletions(-) delete mode 100644 tests/tools/test_converters.py diff --git a/parcels/particlefile.py b/parcels/particlefile.py index f39dd7b3b8..f8e5aea658 100644 --- a/parcels/particlefile.py +++ b/parcels/particlefile.py @@ -150,10 +150,6 @@ def fname(self): def vars_to_write(self): return self._vars_to_write - @property - def time_origin(self): - return self.particleset.time_origin - def _create_variables_attribute_dict(self): """Creates the dictionary with variable attributes. @@ -173,9 +169,8 @@ def _create_variables_attribute_dict(self): "lat": {"long_name": "", "standard_name": "latitude", "units": "degrees_north", "axis": "Y"}, } - if self.time_origin.calendar is not None: - attrs["time"]["units"] = "seconds since " + str(self.time_origin) - attrs["time"]["calendar"] = _set_calendar(self.time_origin.calendar) + attrs["time"]["units"] = "seconds since " # TODO fix units + attrs["time"]["calendar"] = "None" # TODO fix calendar for vname in self.vars_to_write: if vname not in ["time", "lat", "lon", "depth", "id"]: diff --git a/parcels/particleset.py b/parcels/particleset.py index f000aa4416..36f51b25ac 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -17,7 +17,7 @@ from parcels.particle import Particle from parcels.particledata import ParticleData, ParticleDataIterator from parcels.particlefile import ParticleFile -from parcels.tools.converters import _get_cftime_calendars, convert_to_flat_array +from parcels.tools.converters import convert_to_flat_array from parcels.tools.loggers import logger from parcels.tools.statuscodes import StatusCode from parcels.tools.warnings import ParticleSetWarning @@ -25,13 +25,6 @@ __all__ = ["ParticleSet"] -def _convert_to_reltime(time): - """Check to determine if the value of the time parameter needs to be converted to a relative value (relative to the time_origin).""" - if isinstance(time, np.datetime64) or (hasattr(time, "calendar") and time.calendar in _get_cftime_calendars()): - return True - return False - - class ParticleSet: """Class for storing particle and executing kernel over them. @@ -126,7 +119,6 @@ def __init__( else: raise NotImplementedError("particle time must be a datetime or date object") - 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 fieldset.time_interval: _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.time_interval) diff --git a/parcels/tools/converters.py b/parcels/tools/converters.py index 4b60204e1f..8aa5209b85 100644 --- a/parcels/tools/converters.py +++ b/parcels/tools/converters.py @@ -1,7 +1,6 @@ from __future__ import annotations import inspect -from datetime import timedelta from math import cos, pi import cftime @@ -13,7 +12,6 @@ "GeographicPolar", "GeographicPolarSquare", "GeographicSquare", - "TimeConverter", "UnitConverter", "convert_to_flat_array", "unitconverters_map", @@ -42,122 +40,6 @@ def _get_cftime_calendars() -> list[str]: return [getattr(cftime, cf_datetime)(1990, 1, 1).calendar for cf_datetime in _get_cftime_datetimes()] -class TimeConverter: - """Converter class for dates with different calendars in FieldSets - - Parameters - ---------- - time_origin : float, integer, numpy.datetime64 or cftime.DatetimeNoLeap - time origin of the class. - """ - - def __init__(self, time_origin: float | np.datetime64 | np.timedelta64 | cftime.datetime = 0): - self.time_origin = time_origin - self.calendar: str | None = None - if isinstance(time_origin, np.datetime64): - self.calendar = "np_datetime64" - elif isinstance(time_origin, np.timedelta64): - self.calendar = "np_timedelta64" - elif isinstance(time_origin, cftime.datetime): - self.calendar = time_origin.calendar - - def reltime(self, time: TimeConverter | np.datetime64 | np.timedelta64 | cftime.datetime) -> float | npt.NDArray: - """Method to compute the difference, in seconds, between a time and the time_origin - of the TimeConverter - - Parameters - ---------- - time : - - - Returns - ------- - type - time - self.time_origin - - """ - time = time.time_origin if isinstance(time, TimeConverter) else time - if self.calendar in ["np_datetime64", "np_timedelta64"]: - return (time - self.time_origin) / np.timedelta64(1, "s") # type: ignore - elif self.calendar in _get_cftime_calendars(): - if isinstance(time, (list, np.ndarray)): - try: - return np.array([(t - self.time_origin).total_seconds() for t in time]) # type: ignore - except ValueError: - raise ValueError( - f"Cannot subtract 'time' (a {type(time)} object) from a {self.calendar} calendar.\n" - f"Provide 'time' as a {type(self.time_origin)} object?" - ) - else: - try: - return (time - self.time_origin).total_seconds() # type: ignore - except ValueError: - raise ValueError( - f"Cannot subtract 'time' (a {type(time)} object) from a {self.calendar} calendar.\n" - f"Provide 'time' as a {type(self.time_origin)} object?" - ) - elif self.calendar is None: - return time - self.time_origin # type: ignore - else: - raise RuntimeError(f"Calendar {self.calendar} not implemented in TimeConverter") - - def fulltime(self, time): - """Method to convert a time difference in seconds to a date, based on the time_origin - - Parameters - ---------- - time : - - - Returns - ------- - type - self.time_origin + time - - """ - time = time.time_origin if isinstance(time, TimeConverter) else time - if self.calendar in ["np_datetime64", "np_timedelta64"]: - if isinstance(time, (list, np.ndarray)): - return [self.time_origin + np.timedelta64(int(t), "s") for t in time] - else: - return self.time_origin + np.timedelta64(int(time), "s") - elif self.calendar in _get_cftime_calendars(): - return self.time_origin + timedelta(seconds=time) - elif self.calendar is None: - return self.time_origin + time - else: - raise RuntimeError(f"Calendar {self.calendar} not implemented in TimeConverter") - - def __repr__(self): - return f"{self.time_origin}" - - def __eq__(self, other): - other = other.time_origin if isinstance(other, TimeConverter) else other - return self.time_origin == other - - def __ne__(self, other): - other = other.time_origin if isinstance(other, TimeConverter) else other - if not isinstance(other, type(self.time_origin)): - return True - return self.time_origin != other - - def __gt__(self, other): - other = other.time_origin if isinstance(other, TimeConverter) else other - return self.time_origin > other - - def __lt__(self, other): - other = other.time_origin if isinstance(other, TimeConverter) else other - return self.time_origin < other - - def __ge__(self, other): - other = other.time_origin if isinstance(other, TimeConverter) else other - return self.time_origin >= other - - def __le__(self, other): - other = other.time_origin if isinstance(other, TimeConverter) else other - return self.time_origin <= other - - class UnitConverter: """Interface class for spatial unit conversion during field sampling that performs no conversion.""" diff --git a/parcels/tools/statuscodes.py b/parcels/tools/statuscodes.py index 501741ec2c..94589b8d89 100644 --- a/parcels/tools/statuscodes.py +++ b/parcels/tools/statuscodes.py @@ -69,8 +69,6 @@ class TimeExtrapolationError(RuntimeError): """Utility error class to propagate erroneous time extrapolation sampling.""" def __init__(self, time, field=None): - if field is not None and field.grid.time_origin and time is not None: - time = field.grid.time_origin.fulltime(time) message = ( f"{field.name if field else 'Field'} sampled outside time domain at time {time}." " Try setting allow_time_extrapolation to True." @@ -86,23 +84,12 @@ class KernelError(RuntimeError): """General particle kernel error with optional custom message.""" def __init__(self, particle, fieldset=None, msg=None): - message = ( - f"{particle.state}\n" - f"Particle {particle}\n" - f"Time: {_parse_particletime(particle.time, fieldset)}\n" - f"timestep dt: {particle.dt}\n" - ) + message = f"{particle.state}\nParticle {particle}\nTime: {particle.time}\ntimestep dt: {particle.dt}\n" if msg: message += msg super().__init__(message) -def _parse_particletime(time, fieldset): - if fieldset is not None and fieldset.time_origin: - time = fieldset.time_origin.fulltime(time) - return time - - AllParcelsErrorCodes = { FieldSamplingError: StatusCode.Error, FieldOutOfBoundError: StatusCode.ErrorOutOfBounds, diff --git a/parcels/xgrid.py b/parcels/xgrid.py index 5d536a39ed..d8f65b43b9 100644 --- a/parcels/xgrid.py +++ b/parcels/xgrid.py @@ -8,7 +8,6 @@ from parcels import xgcm from parcels._index_search import _search_indices_curvilinear_2d from parcels.basegrid import BaseGrid -from parcels.tools.converters import TimeConverter _XGRID_AXES_ORDERING = "ZYX" _XGRID_AXES = Literal["X", "Y", "Z"] @@ -130,16 +129,6 @@ def zdim(self): def tdim(self): return get_cell_edge_count_along_dim(self.xgcm_grid.axes.get("T")) - @property - def time_origin(self): - """ - Note - ---- - Included for compatibility with v3 codebase. May be removed in future. - TODO v4: Evaluate - """ - return TimeConverter(self._datetimes[0]) - @property def _z4d(self) -> Literal[0, 1]: """ diff --git a/tests/tools/test_converters.py b/tests/tools/test_converters.py deleted file mode 100644 index 3d35363b1d..0000000000 --- a/tests/tools/test_converters.py +++ /dev/null @@ -1,55 +0,0 @@ -import cftime -import numpy as np -import pytest - -from parcels.tools.converters import TimeConverter, _get_cftime_datetimes - -cf_datetime_classes = [getattr(cftime, c) for c in _get_cftime_datetimes()] -cf_datetime_objects = [c(1990, 1, 1) for c in cf_datetime_classes] - - -@pytest.mark.parametrize( - "cf_datetime", - cf_datetime_objects, -) -def test_TimeConverter_cf(cf_datetime): - assert TimeConverter(cf_datetime).calendar == cf_datetime.calendar - assert TimeConverter(cf_datetime).time_origin == cf_datetime - - -def test_TimeConverter_standard(): - dt = np.datetime64("2001-01-01T12:00") - assert TimeConverter(dt).calendar == "np_datetime64" - assert TimeConverter(dt).time_origin == dt - - dt = np.timedelta64(1, "s") - assert TimeConverter(dt).calendar == "np_timedelta64" - assert TimeConverter(dt).time_origin == dt - - assert TimeConverter(0).calendar is None - assert TimeConverter(0).time_origin == 0 - - -def test_TimeConverter_reltime_one_day(): - ONE_DAY = 24 * 60 * 60 - first_jan = [c(1990, 1, 1) for c in cf_datetime_classes] + [0] - second_jan = [c(1990, 1, 2) for c in cf_datetime_classes] + [ONE_DAY] - - for time_origin, time in zip(first_jan, second_jan, strict=True): - tc = TimeConverter(time_origin) - assert tc.reltime(time) == ONE_DAY - - -@pytest.mark.parametrize( - "x, y", - [ - pytest.param(np.datetime64("2001-01-01T12:00"), 0, id="datetime64 float"), - pytest.param(cftime.DatetimeNoLeap(1990, 1, 1), 0, id="cftime float"), - pytest.param(cftime.DatetimeNoLeap(1990, 1, 1), cftime.DatetimeAllLeap(1991, 1, 1), id="cftime cftime"), - ], -) -def test_TimeConverter_reltime_errors(x, y): - """All of these should raise a ValueError when doing reltime""" - tc = TimeConverter(x) - with pytest.raises((ValueError, TypeError)): - tc.reltime(y) diff --git a/tests/v4/test_xgrid.py b/tests/v4/test_xgrid.py index af84f125dd..ba501346f2 100644 --- a/tests/v4/test_xgrid.py +++ b/tests/v4/test_xgrid.py @@ -22,7 +22,6 @@ GridTestCase(datasets["ds_2d_left"], "ydim", Y), GridTestCase(datasets["ds_2d_left"], "zdim", Z), GridTestCase(datasets["ds_2d_left"], "tdim", T), - GridTestCase(datasets["ds_2d_left"], "time_origin", TimeConverter(datasets["ds_2d_left"].time.values[0])), ] @@ -56,7 +55,6 @@ def test_xgrid_properties_ground_truth(ds, attr, expected): "ydim", "zdim", "tdim", - "time_origin", "_gtype", ], ) @@ -69,7 +67,6 @@ def test_xgrid_against_old(ds, attr): lat=ds.lat.values, depth=ds.depth.values, time=ds.time.values.astype("float64") / 1e9, - time_origin=TimeConverter(ds.time.values[0]), mesh="spherical", ) actual = getattr(grid, attr) From 61dd7963a9b5ed3c781c822a491de9c4c4d8de54 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 7 Jul 2025 16:19:32 +0200 Subject: [PATCH 2/4] More removals of time_origin and TImeConverter --- parcels/grid.py | 51 +++++++++++------------------------------- tests/test_reprs.py | 10 +++------ tests/v4/test_xgrid.py | 3 --- 3 files changed, 16 insertions(+), 48 deletions(-) diff --git a/parcels/grid.py b/parcels/grid.py index ed8c67e295..53ece05f5c 100644 --- a/parcels/grid.py +++ b/parcels/grid.py @@ -4,7 +4,6 @@ import numpy.typing as npt from parcels._typing import Mesh, assert_valid_mesh -from parcels.tools.converters import TimeConverter __all__ = [ "CurvilinearSGrid", @@ -31,7 +30,6 @@ def __init__( lon: npt.NDArray, lat: npt.NDArray, time: npt.NDArray | None, - time_origin: TimeConverter | None, mesh: Mesh, ): lon = np.array(lon) @@ -52,8 +50,6 @@ def __init__( self._lat = lat self.time = time self.tdim = time.size - self._time_origin = TimeConverter() if time_origin is None else time_origin - assert isinstance(self.time_origin, TimeConverter), "time_origin needs to be a TimeConverter object" assert_valid_mesh(mesh) self._mesh = mesh self._zonal_periodic = False @@ -63,11 +59,7 @@ def __init__( def __repr__(self): with np.printoptions(threshold=5, suppress=True, linewidth=120, formatter={"float": "{: 0.2f}".format}): - return ( - f"{type(self).__name__}(" - f"lon={self.lon!r}, lat={self.lat!r}, time={self.time!r}, " - f"time_origin={self.time_origin!r}, mesh={self.mesh!r})" - ) + return f"{type(self).__name__}(lon={self.lon!r}, lat={self.lat!r}, time={self.time!r}, " @property def lon(self): @@ -89,10 +81,6 @@ def mesh(self): def lonlat_minmax(self): return self._lonlat_minmax - @property - def time_origin(self): - return self._time_origin - @property def zonal_periodic(self): return self._zonal_periodic @@ -103,7 +91,6 @@ def create_grid( lat: npt.ArrayLike, depth, time, - time_origin, mesh: Mesh, ): lon = np.array(lon) @@ -114,14 +101,14 @@ def create_grid( if len(lon.shape) <= 1: if depth is None or len(depth.shape) <= 1: - return RectilinearZGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh) + return RectilinearZGrid(lon, lat, depth, time, mesh=mesh) else: - return RectilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh) + return RectilinearSGrid(lon, lat, depth, time, mesh=mesh) else: if depth is None or len(depth.shape) <= 1: - return CurvilinearZGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh) + return CurvilinearZGrid(lon, lat, depth, time, mesh=mesh) else: - return CurvilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh) + return CurvilinearSGrid(lon, lat, depth, time, mesh=mesh) def _check_zonal_periodic(self): if self.zonal_periodic or self.mesh == "flat" or self.lon.size == 1: @@ -167,14 +154,14 @@ class RectilinearGrid(Grid): """ - def __init__(self, lon, lat, time, time_origin, mesh: Mesh): + def __init__(self, lon, lat, time, mesh: Mesh): assert isinstance(lon, np.ndarray) and len(lon.shape) <= 1, "lon is not a numpy vector" assert isinstance(lat, np.ndarray) and len(lat.shape) <= 1, "lat is not a numpy vector" assert isinstance(time, np.ndarray) or not time, "time is not a numpy array" if isinstance(time, np.ndarray): assert len(time.shape) == 1, "time is not a vector" - super().__init__(lon, lat, time, time_origin, mesh) + super().__init__(lon, lat, time, mesh) @property def xdim(self): @@ -199,8 +186,6 @@ class RectilinearZGrid(RectilinearGrid): The depth of the different layers is thus constant. time : Vector containing the time coordinates of the grid - time_origin : parcels.tools.converters.TimeConverter - Time origin of the time axis mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: @@ -210,8 +195,8 @@ class RectilinearZGrid(RectilinearGrid): 2. flat: No conversion, lat/lon are assumed to be in m. """ - def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh: Mesh = "flat"): - super().__init__(lon, lat, time, time_origin, mesh) + def __init__(self, lon, lat, depth=None, time=None, mesh: Mesh = "flat"): + super().__init__(lon, lat, time, mesh) if isinstance(depth, np.ndarray): assert len(depth.shape) <= 1, "depth is not a vector" @@ -247,8 +232,6 @@ class RectilinearSGrid(RectilinearGrid): depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim]. time : Vector containing the time coordinates of the grid - time_origin : parcels.tools.converters.TimeConverter - Time origin of the time axis mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: @@ -264,10 +247,9 @@ def __init__( lat: npt.NDArray, depth: npt.NDArray, time: npt.NDArray | None = None, - time_origin: TimeConverter | None = None, mesh: Mesh = "flat", ): - super().__init__(lon, lat, time, time_origin, mesh) + super().__init__(lon, lat, time, mesh) assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 3D or 4D numpy array" self._gtype = GridType.RectilinearSGrid @@ -306,7 +288,6 @@ def __init__( lon: npt.NDArray, lat: npt.NDArray, time: npt.NDArray | None = None, - time_origin: TimeConverter | None = None, mesh: Mesh = "flat", ): assert isinstance(lon, np.ndarray) and len(lon.squeeze().shape) == 2, "lon is not a 2D numpy array" @@ -317,7 +298,7 @@ def __init__( lon = lon.squeeze() lat = lat.squeeze() - super().__init__(lon, lat, time, time_origin, mesh) + super().__init__(lon, lat, time, mesh) @property def xdim(self): @@ -342,8 +323,6 @@ class CurvilinearZGrid(CurvilinearGrid): The depth of the different layers is thus constant. time : Vector containing the time coordinates of the grid - time_origin : parcels.tools.converters.TimeConverter - Time origin of the time axis mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: @@ -359,10 +338,9 @@ def __init__( lat: npt.NDArray, depth: npt.NDArray | None = None, time: npt.NDArray | None = None, - time_origin: TimeConverter | None = None, mesh: Mesh = "flat", ): - super().__init__(lon, lat, time, time_origin, mesh) + super().__init__(lon, lat, time, mesh) if isinstance(depth, np.ndarray): assert len(depth.shape) == 1, "depth is not a vector" @@ -396,8 +374,6 @@ class CurvilinearSGrid(CurvilinearGrid): depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim]. time : Vector containing the time coordinates of the grid - time_origin : parcels.tools.converters.TimeConverter - Time origin of the time axis mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: @@ -413,10 +389,9 @@ def __init__( lat: npt.NDArray, depth: npt.NDArray, time: npt.NDArray | None = None, - time_origin: TimeConverter | None = None, mesh: Mesh = "flat", ): - super().__init__(lon, lat, time, time_origin, mesh) + super().__init__(lon, lat, time, mesh) assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 4D numpy array" self._gtype = GridType.CurvilinearSGrid diff --git a/tests/test_reprs.py b/tests/test_reprs.py index 8b6fc79f77..18407e583d 100644 --- a/tests/test_reprs.py +++ b/tests/test_reprs.py @@ -5,7 +5,7 @@ import numpy as np import parcels -from parcels import Grid, ParticleFile, TimeConverter, Variable +from parcels import Grid, ParticleFile, Variable from parcels.grid import RectilinearGrid from tests.utils import create_fieldset_unit_mesh, create_simple_pset @@ -58,9 +58,7 @@ def test_particletype_repr(): def test_grid_repr(): """Test arguments are in the repr of a Grid object""" - kwargs = dict( - lon=np.array([1, 2, 3]), lat=np.array([4, 5, 6]), time=None, time_origin=TimeConverter(), mesh="spherical" - ) + kwargs = dict(lon=np.array([1, 2, 3]), lat=np.array([4, 5, 6]), time=None, mesh="spherical") assert_simple_repr(Grid, kwargs) @@ -76,9 +74,7 @@ def test_rectilineargrid_repr(): Mainly to test inherited repr is correct. """ - kwargs = dict( - lon=np.array([1, 2, 3]), lat=np.array([4, 5, 6]), time=None, time_origin=TimeConverter(), mesh="spherical" - ) + kwargs = dict(lon=np.array([1, 2, 3]), lat=np.array([4, 5, 6]), time=None, mesh="spherical") assert_simple_repr(RectilinearGrid, kwargs) diff --git a/tests/v4/test_xgrid.py b/tests/v4/test_xgrid.py index ba501346f2..f5eac0d184 100644 --- a/tests/v4/test_xgrid.py +++ b/tests/v4/test_xgrid.py @@ -8,7 +8,6 @@ from parcels import xgcm from parcels._datasets.structured.generic import T, X, Y, Z, datasets from parcels.grid import Grid as OldGrid -from parcels.tools.converters import TimeConverter from parcels.xgrid import XGrid, _search_1d_array GridTestCase = namedtuple("GridTestCase", ["Grid", "attr", "expected"]) @@ -28,8 +27,6 @@ def assert_equal(actual, expected): if expected is None: assert actual is None - elif isinstance(expected, TimeConverter): - assert actual == expected elif isinstance(expected, np.ndarray): assert actual.shape == expected.shape assert_allclose(actual, expected) From ec6f33fccd9ed8e01fa3433210d791a79d8010c5 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 8 Jul 2025 08:07:37 +0200 Subject: [PATCH 3/4] Removing grid.py from Parcels v4 --- parcels/__init__.py | 1 - parcels/_index_search.py | 3 +- parcels/basegrid.py | 8 + parcels/field.py | 10 +- parcels/grid.py | 423 --------------------------------------- parcels/kernel.py | 2 +- parcels/particleset.py | 2 +- tests/v4/test_xgrid.py | 31 --- 8 files changed, 12 insertions(+), 468 deletions(-) delete mode 100644 parcels/grid.py diff --git a/parcels/__init__.py b/parcels/__init__.py index fe8d4bfb47..c452c9ad96 100644 --- a/parcels/__init__.py +++ b/parcels/__init__.py @@ -5,7 +5,6 @@ from parcels.application_kernels import * from parcels.field import * from parcels.fieldset import * -from parcels.grid import * from parcels.interaction import * from parcels.kernel import * from parcels.particle import * diff --git a/parcels/_index_search.py b/parcels/_index_search.py index 747ee2b186..4638b18486 100644 --- a/parcels/_index_search.py +++ b/parcels/_index_search.py @@ -18,13 +18,12 @@ _raise_time_extrapolation_error, ) -from .grid import GridType +from .basegrid import GridType if TYPE_CHECKING: from parcels.xgrid import XGrid from .field import Field - # from .grid import Grid def _search_time_index(field: Field, time: datetime): diff --git a/parcels/basegrid.py b/parcels/basegrid.py index a7ebcfb42a..e1a866e225 100644 --- a/parcels/basegrid.py +++ b/parcels/basegrid.py @@ -1,12 +1,20 @@ from __future__ import annotations from abc import ABC, abstractmethod +from enum import IntEnum from typing import TYPE_CHECKING if TYPE_CHECKING: import numpy as np +class GridType(IntEnum): + RectilinearZGrid = 0 + RectilinearSGrid = 1 + CurvilinearZGrid = 2 + CurvilinearSGrid = 3 + + class BaseGrid(ABC): @abstractmethod def search(self, z: float, y: float, x: float, ei=None) -> dict[str, tuple[int, float | np.ndarray]]: diff --git a/parcels/field.py b/parcels/field.py index ce82d4c781..473ecbc833 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -4,7 +4,6 @@ import warnings from collections.abc import Callable from datetime import datetime -from enum import IntEnum import numpy as np import uxarray as ux @@ -34,14 +33,7 @@ from ._index_search import _search_time_index -__all__ = ["Field", "GridType", "VectorField"] - - -class GridType(IntEnum): - RectilinearZGrid = 0 - RectilinearSGrid = 1 - CurvilinearZGrid = 2 - CurvilinearSGrid = 3 +__all__ = ["Field", "VectorField"] def _isParticle(key): diff --git a/parcels/grid.py b/parcels/grid.py deleted file mode 100644 index 53ece05f5c..0000000000 --- a/parcels/grid.py +++ /dev/null @@ -1,423 +0,0 @@ -from enum import IntEnum - -import numpy as np -import numpy.typing as npt - -from parcels._typing import Mesh, assert_valid_mesh - -__all__ = [ - "CurvilinearSGrid", - "CurvilinearZGrid", - "Grid", - "GridType", - "RectilinearSGrid", - "RectilinearZGrid", -] - - -class GridType(IntEnum): - RectilinearZGrid = 0 - RectilinearSGrid = 1 - CurvilinearZGrid = 2 - CurvilinearSGrid = 3 - - -class Grid: - """Grid class that defines a (spatial and temporal) grid on which Fields are defined.""" - - def __init__( - self, - lon: npt.NDArray, - lat: npt.NDArray, - time: npt.NDArray | None, - mesh: Mesh, - ): - lon = np.array(lon) - lat = np.array(lat) - time = np.zeros(1, dtype=np.float64) if time is None else time - time = np.array(time) - if not lon.dtype == np.float32: - lon = lon.astype(np.float32) - if not lat.dtype == np.float32: - lat = lat.astype(np.float32) - if not time.dtype == np.float64: - assert isinstance( - time[0], (np.integer, np.floating, float, int) - ), "Time vector must be an array of int or floats" - time = time.astype(np.float64) - - self._lon = lon - self._lat = lat - self.time = time - self.tdim = time.size - assert_valid_mesh(mesh) - self._mesh = mesh - self._zonal_periodic = False - self._lonlat_minmax = np.array( - [np.nanmin(lon), np.nanmax(lon), np.nanmin(lat), np.nanmax(lat)], dtype=np.float32 - ) - - def __repr__(self): - with np.printoptions(threshold=5, suppress=True, linewidth=120, formatter={"float": "{: 0.2f}".format}): - return f"{type(self).__name__}(lon={self.lon!r}, lat={self.lat!r}, time={self.time!r}, " - - @property - def lon(self): - return self._lon - - @property - def lat(self): - return self._lat - - @property - def depth(self): - return self._depth - - @property - def mesh(self): - return self._mesh - - @property - def lonlat_minmax(self): - return self._lonlat_minmax - - @property - def zonal_periodic(self): - return self._zonal_periodic - - @staticmethod - def create_grid( - lon: npt.ArrayLike, - lat: npt.ArrayLike, - depth, - time, - mesh: Mesh, - ): - lon = np.array(lon) - lat = np.array(lat) - - if depth is not None: - depth = np.array(depth) - - if len(lon.shape) <= 1: - if depth is None or len(depth.shape) <= 1: - return RectilinearZGrid(lon, lat, depth, time, mesh=mesh) - else: - return RectilinearSGrid(lon, lat, depth, time, mesh=mesh) - else: - if depth is None or len(depth.shape) <= 1: - return CurvilinearZGrid(lon, lat, depth, time, mesh=mesh) - else: - return CurvilinearSGrid(lon, lat, depth, time, mesh=mesh) - - def _check_zonal_periodic(self): - if self.zonal_periodic or self.mesh == "flat" or self.lon.size == 1: - return - dx = (self.lon[1:] - self.lon[:-1]) if len(self.lon.shape) == 1 else self.lon[0, 1:] - self.lon[0, :-1] - dx = np.where(dx < -180, dx + 360, dx) - dx = np.where(dx > 180, dx - 360, dx) - self._zonal_periodic = sum(dx) > 359.9 - - def _add_Sdepth_periodic_halo(self, zonal, meridional, halosize): - if zonal: - if len(self.depth.shape) == 3: - self._depth = np.concatenate( - (self.depth[:, :, -halosize:], self.depth, self.depth[:, :, 0:halosize]), - axis=len(self.depth.shape) - 1, - ) - assert self.depth.shape[2] == self.xdim, "Third dim must be x." - else: - self._depth = np.concatenate( - (self.depth[:, :, :, -halosize:], self.depth, self.depth[:, :, :, 0:halosize]), - axis=len(self.depth.shape) - 1, - ) - assert self.depth.shape[3] == self.xdim, "Fourth dim must be x." - if meridional: - if len(self.depth.shape) == 3: - self._depth = np.concatenate( - (self.depth[:, -halosize:, :], self.depth, self.depth[:, 0:halosize, :]), - axis=len(self.depth.shape) - 2, - ) - assert self.depth.shape[1] == self.ydim, "Second dim must be y." - else: - self._depth = np.concatenate( - (self.depth[:, :, -halosize:, :], self.depth, self.depth[:, :, 0:halosize, :]), - axis=len(self.depth.shape) - 2, - ) - assert self.depth.shape[2] == self.ydim, "Third dim must be y." - - -class RectilinearGrid(Grid): - """Rectilinear Grid class - - Private base class for RectilinearZGrid and RectilinearSGrid - - """ - - def __init__(self, lon, lat, time, mesh: Mesh): - assert isinstance(lon, np.ndarray) and len(lon.shape) <= 1, "lon is not a numpy vector" - assert isinstance(lat, np.ndarray) and len(lat.shape) <= 1, "lat is not a numpy vector" - assert isinstance(time, np.ndarray) or not time, "time is not a numpy array" - if isinstance(time, np.ndarray): - assert len(time.shape) == 1, "time is not a vector" - - super().__init__(lon, lat, time, mesh) - - @property - def xdim(self): - return self.lon.size - - @property - def ydim(self): - return self.lat.size - - -class RectilinearZGrid(RectilinearGrid): - """Rectilinear Z Grid. - - Parameters - ---------- - lon : - Vector containing the longitude coordinates of the grid - lat : - Vector containing the latitude coordinates of the grid - depth : - Vector containing the vertical coordinates of the grid, which are z-coordinates. - The depth of the different layers is thus constant. - time : - Vector containing the time coordinates of the grid - mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation: - - 1. spherical (default): Lat and lon in degree, with a - correction for zonal velocity U near the poles. - 2. flat: No conversion, lat/lon are assumed to be in m. - """ - - def __init__(self, lon, lat, depth=None, time=None, mesh: Mesh = "flat"): - super().__init__(lon, lat, time, mesh) - if isinstance(depth, np.ndarray): - assert len(depth.shape) <= 1, "depth is not a vector" - - self._gtype = GridType.RectilinearZGrid - self._depth = np.zeros(1, dtype=np.float32) if depth is None else depth - self._depth = np.array(self.depth) - self._z4d = -1 # only used in RectilinearSGrid - if not self.depth.dtype == np.float32: - self._depth = self.depth.astype(np.float32) - - @property - def zdim(self): - return self.depth.size - - -class RectilinearSGrid(RectilinearGrid): - """Rectilinear S Grid. Same horizontal discretisation as a rectilinear z grid, - but with s vertical coordinates - - Parameters - ---------- - lon : - Vector containing the longitude coordinates of the grid - lat : - Vector containing the latitude coordinates of the grid - depth : - 4D (time-evolving) or 3D (time-independent) array containing the vertical coordinates of the grid, - which are s-coordinates. - s-coordinates can be terrain-following (sigma) or iso-density (rho) layers, - or any generalised vertical discretisation. - The depth of each node depends then on the horizontal position (lon, lat), - the number of the layer and the time is depth is a 4D array. - depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim]. - time : - Vector containing the time coordinates of the grid - mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation: - - 1. spherical (default): Lat and lon in degree, with a - correction for zonal velocity U near the poles. - 2. flat: No conversion, lat/lon are assumed to be in m. - """ - - def __init__( - self, - lon: npt.NDArray, - lat: npt.NDArray, - depth: npt.NDArray, - time: npt.NDArray | None = None, - mesh: Mesh = "flat", - ): - super().__init__(lon, lat, time, mesh) - assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 3D or 4D numpy array" - - self._gtype = GridType.RectilinearSGrid - self._depth = depth - self._depth = np.array(self.depth) - self._z4d = 1 if len(self.depth.shape) == 4 else 0 - if self._z4d: - # self.depth.shape[0] is 0 for S grids loaded from netcdf file - assert ( - self.tdim == self.depth.shape[0] or self.depth.shape[0] == 0 - ), "depth dimension has the wrong format. It should be [tdim, zdim, ydim, xdim]" - assert ( - self.xdim == self.depth.shape[-1] or self.depth.shape[-1] == 0 - ), "depth dimension has the wrong format. It should be [tdim, zdim, ydim, xdim]" - assert ( - self.ydim == self.depth.shape[-2] or self.depth.shape[-2] == 0 - ), "depth dimension has the wrong format. It should be [tdim, zdim, ydim, xdim]" - else: - assert ( - self.xdim == self.depth.shape[-1] - ), "depth dimension has the wrong format. It should be [zdim, ydim, xdim]" - assert ( - self.ydim == self.depth.shape[-2] - ), "depth dimension has the wrong format. It should be [zdim, ydim, xdim]" - if not self.depth.dtype == np.float32: - self._depth = self.depth.astype(np.float32) - - @property - def zdim(self): - return self.depth.shape[-3] - - -class CurvilinearGrid(Grid): - def __init__( - self, - lon: npt.NDArray, - lat: npt.NDArray, - time: npt.NDArray | None = None, - mesh: Mesh = "flat", - ): - assert isinstance(lon, np.ndarray) and len(lon.squeeze().shape) == 2, "lon is not a 2D numpy array" - assert isinstance(lat, np.ndarray) and len(lat.squeeze().shape) == 2, "lat is not a 2D numpy array" - assert isinstance(time, np.ndarray) or not time, "time is not a numpy array" - if isinstance(time, np.ndarray): - assert len(time.shape) == 1, "time is not a vector" - - lon = lon.squeeze() - lat = lat.squeeze() - super().__init__(lon, lat, time, mesh) - - @property - def xdim(self): - return self.lon.shape[1] - - @property - def ydim(self): - return self.lon.shape[0] - - -class CurvilinearZGrid(CurvilinearGrid): - """Curvilinear Z Grid. - - Parameters - ---------- - lon : - 2D array containing the longitude coordinates of the grid - lat : - 2D array containing the latitude coordinates of the grid - depth : - Vector containing the vertical coordinates of the grid, which are z-coordinates. - The depth of the different layers is thus constant. - time : - Vector containing the time coordinates of the grid - mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation: - - 1. spherical (default): Lat and lon in degree, with a - correction for zonal velocity U near the poles. - 2. flat: No conversion, lat/lon are assumed to be in m. - """ - - def __init__( - self, - lon: npt.NDArray, - lat: npt.NDArray, - depth: npt.NDArray | None = None, - time: npt.NDArray | None = None, - mesh: Mesh = "flat", - ): - super().__init__(lon, lat, time, mesh) - if isinstance(depth, np.ndarray): - assert len(depth.shape) == 1, "depth is not a vector" - - self._gtype = GridType.CurvilinearZGrid - self._depth = np.zeros(1, dtype=np.float32) if depth is None else depth - self._z4d = -1 # only for SGrid - if not self.depth.dtype == np.float32: - self._depth = self.depth.astype(np.float32) - - @property - def zdim(self): - return self.depth.size - - -class CurvilinearSGrid(CurvilinearGrid): - """Curvilinear S Grid. - - Parameters - ---------- - lon : - 2D array containing the longitude coordinates of the grid - lat : - 2D array containing the latitude coordinates of the grid - depth : - 4D (time-evolving) or 3D (time-independent) array containing the vertical coordinates of the grid, - which are s-coordinates. - s-coordinates can be terrain-following (sigma) or iso-density (rho) layers, - or any generalised vertical discretisation. - The depth of each node depends then on the horizontal position (lon, lat), - the number of the layer and the time is depth is a 4D array. - depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim]. - time : - Vector containing the time coordinates of the grid - mesh : str - String indicating the type of mesh coordinates and - units used during velocity interpolation: - - 1. spherical (default): Lat and lon in degree, with a - correction for zonal velocity U near the poles. - 2. flat: No conversion, lat/lon are assumed to be in m. - """ - - def __init__( - self, - lon: npt.NDArray, - lat: npt.NDArray, - depth: npt.NDArray, - time: npt.NDArray | None = None, - mesh: Mesh = "flat", - ): - super().__init__(lon, lat, time, mesh) - assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 4D numpy array" - - self._gtype = GridType.CurvilinearSGrid - self._depth = depth - self._z4d = 1 if len(self.depth.shape) == 4 else 0 - if self._z4d: - # self.depth.shape[0] is 0 for S grids loaded from netcdf file - assert ( - self.tdim == self.depth.shape[0] or self.depth.shape[0] == 0 - ), "depth dimension has the wrong format. It should be [tdim, zdim, ydim, xdim]" - assert ( - self.xdim == self.depth.shape[-1] or self.depth.shape[-1] == 0 - ), "depth dimension has the wrong format. It should be [tdim, zdim, ydim, xdim]" - assert ( - self.ydim == self.depth.shape[-2] or self.depth.shape[-2] == 0 - ), "depth dimension has the wrong format. It should be [tdim, zdim, ydim, xdim]" - else: - assert ( - self.xdim == self.depth.shape[-1] - ), "depth dimension has the wrong format. It should be [zdim, ydim, xdim]" - assert ( - self.ydim == self.depth.shape[-2] - ), "depth dimension has the wrong format. It should be [zdim, ydim, xdim]" - if not self.depth.dtype == np.float32: - self._depth = self.depth.astype(np.float32) - - @property - def zdim(self): - return self.depth.shape[-3] diff --git a/parcels/kernel.py b/parcels/kernel.py index 287afce4c9..c35ed0314b 100644 --- a/parcels/kernel.py +++ b/parcels/kernel.py @@ -16,7 +16,7 @@ AdvectionRK4_3D_CROCO, AdvectionRK45, ) -from parcels.grid import GridType +from parcels.basegrid import GridType from parcels.tools.statuscodes import ( StatusCode, TimeExtrapolationError, diff --git a/parcels/particleset.py b/parcels/particleset.py index a4b90a3c62..4be041ac39 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -11,7 +11,7 @@ from parcels._core.utils.time import TimeInterval from parcels._reprs import particleset_repr from parcels.application_kernels.advection import AdvectionRK4 -from parcels.grid import GridType +from parcels.basegrid import GridType from parcels.interaction.interactionkernel import InteractionKernel from parcels.kernel import Kernel from parcels.particle import Particle diff --git a/tests/v4/test_xgrid.py b/tests/v4/test_xgrid.py index f5eac0d184..301a16c131 100644 --- a/tests/v4/test_xgrid.py +++ b/tests/v4/test_xgrid.py @@ -7,7 +7,6 @@ from parcels import xgcm from parcels._datasets.structured.generic import T, X, Y, Z, datasets -from parcels.grid import Grid as OldGrid from parcels.xgrid import XGrid, _search_1d_array GridTestCase = namedtuple("GridTestCase", ["Grid", "attr", "expected"]) @@ -41,36 +40,6 @@ def test_xgrid_properties_ground_truth(ds, attr, expected): assert_equal(actual, expected) -@pytest.mark.parametrize( - "attr", - [ - "lon", - "lat", - "depth", - "time", - "xdim", - "ydim", - "zdim", - "tdim", - "_gtype", - ], -) -@pytest.mark.parametrize("ds", [pytest.param(ds, id=key) for key, ds in datasets.items()]) -def test_xgrid_against_old(ds, attr): - grid = XGrid(xgcm.Grid(ds, periodic=False)) - - old_grid = OldGrid.create_grid( - lon=ds.lon.values, - lat=ds.lat.values, - depth=ds.depth.values, - time=ds.time.values.astype("float64") / 1e9, - mesh="spherical", - ) - actual = getattr(grid, attr) - expected = getattr(old_grid, attr) - assert_equal(actual, expected) - - @pytest.mark.parametrize("ds", [pytest.param(ds, id=key) for key, ds in datasets.items()]) def test_xgrid_init_on_generic_datasets(ds): XGrid(xgcm.Grid(ds, periodic=False)) From 231241b4e47c5f5f0c22f203c3862e45be3ffdcc Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 8 Jul 2025 10:35:18 +0200 Subject: [PATCH 4/4] Also removing _get_cftime_calendars and _get_cftime_datetimes --- parcels/tools/converters.py | 13 ------------- tests/test_particlefile.py | 10 ---------- 2 files changed, 23 deletions(-) diff --git a/parcels/tools/converters.py b/parcels/tools/converters.py index 8aa5209b85..17311d2055 100644 --- a/parcels/tools/converters.py +++ b/parcels/tools/converters.py @@ -1,9 +1,7 @@ from __future__ import annotations -import inspect from math import cos, pi -import cftime import numpy as np import numpy.typing as npt @@ -29,17 +27,6 @@ def convert_to_flat_array(var: npt.ArrayLike) -> npt.NDArray: return np.array(var).flatten() -def _get_cftime_datetimes() -> list[str]: - # Is there a more elegant way to parse these from cftime? - cftime_calendars = tuple(x[1].__name__ for x in inspect.getmembers(cftime._cftime, inspect.isclass)) - cftime_datetime_names = [ca for ca in cftime_calendars if "Datetime" in ca] - return cftime_datetime_names - - -def _get_cftime_calendars() -> list[str]: - return [getattr(cftime, cf_datetime)(1990, 1, 1).calendar for cf_datetime in _get_cftime_datetimes()] - - class UnitConverter: """Interface class for spatial unit conversion during field sampling that performs no conversion.""" diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index d8d7c18250..0895774144 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -2,7 +2,6 @@ import tempfile from datetime import timedelta -import cftime import numpy as np import pytest import xarray as xr @@ -17,8 +16,6 @@ ParticleSet, Variable, ) -from parcels.particlefile import _set_calendar -from parcels.tools.converters import _get_cftime_calendars, _get_cftime_datetimes from tests.common_kernels import DoNothing from tests.utils import create_fieldset_zeros_simple @@ -305,13 +302,6 @@ def SampleP(particle, fieldset, time): # pragma: no cover ds.close() -def test_set_calendar(): - for _calendar_name, cf_datetime in zip(_get_cftime_calendars(), _get_cftime_datetimes(), strict=True): - date = getattr(cftime, cf_datetime)(1990, 1, 1) - assert _set_calendar(date.calendar) == date.calendar - assert _set_calendar("np_datetime64") == "standard" - - def test_reset_dt(fieldset, tmp_zarrfile): # Assert that p.dt gets reset when a write_time is not a multiple of dt # for p.dt=0.02 to reach outputdt=0.05 and endtime=0.1, the steps should be [0.2, 0.2, 0.1, 0.2, 0.2, 0.1], resulting in 6 kernel executions