Skip to content

Use datetime64 natively for particle time throughout codebase (and drop TimeConverter and time_origin) #1876

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
erikvansebille wants to merge 1 commit into from
Closed

Use datetime64 natively for particle time throughout codebase (and drop TimeConverter and time_origin) #1876

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 9 additions & 24 deletions parcels/field.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@
)
from parcels.tools._helpers import default_repr, field_repr, timedelta_to_float
from parcels.tools.converters import (
TimeConverter,
UnitConverter,
unitconverters_map,
)
Expand Down Expand Up @@ -137,8 +136,8 @@ class Field:
A numpy array containing the timestamps for each of the files in filenames, for loading
from netCDF files only. Default is None if the netCDF dimensions dictionary includes time.
grid : parcels.grid.Grid
:class:`parcels.grid.Grid` object containing all the lon, lat depth, time
mesh and time_origin information. Can be constructed from any of the Grid objects
:class:`parcels.grid.Grid` object containing all the lon, lat depth, time and mesh information.
Can be constructed from any of the Grid objects
fieldtype : str
Type of Field to be used for UnitConverter (either 'U', 'V', 'Kh_zonal', 'Kh_meridional' or None)
transpose : bool
Expand All @@ -149,8 +148,6 @@ class Field:
Maximum allowed value on the field. Data above this value are set to zero
cast_data_dtype : str
Cast Field data to dtype. Supported dtypes are "float32" (np.float32 (default)) and "float64 (np.float64).
time_origin : parcels.tools.converters.TimeConverter
Time origin of the time axis (only if grid is None)
interp_method : str
Method for interpolation. Options are 'linear' (default), 'nearest',
'linear_invdist_land_tracer', 'cgrid_velocity', 'cgrid_tracer' and 'bgrid_velocity'
Expand Down Expand Up @@ -195,7 +192,6 @@ def __init__(
vmin: float | None = None,
vmax: float | None = None,
cast_data_dtype: type[np.float32] | type[np.float64] | Literal["float32", "float64"] = "float32",
time_origin: TimeConverter | None = None,
interp_method: InterpMethod = "linear",
allow_time_extrapolation: bool | None = None,
time_periodic: TimePeriodic = False,
Expand All @@ -221,12 +217,7 @@ def __init__(
)
self._grid = grid
else:
if (time is not None) and isinstance(time[0], np.datetime64):
time_origin = TimeConverter(time[0])
time = np.array([time_origin.reltime(t) for t in time])
else:
time_origin = TimeConverter(0)
self._grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
self._grid = Grid.create_grid(lon, lat, depth, time, mesh=mesh)
self.igrid = -1
self.fieldtype = self.name if fieldtype is None else fieldtype
self.to_write = to_write
Expand Down Expand Up @@ -439,15 +430,13 @@ def _collect_timeslices(
dataFiles = np.concatenate(dataFiles).ravel()
if time.size == 1 and time[0] is None:
time[0] = 0
time_origin = TimeConverter(time[0])
time = time_origin.reltime(time)

if not np.all((time[1:] - time[:-1]) > 0):
id_not_ordered = np.where(time[1:] < time[:-1])[0][0]
raise AssertionError(
f"Please make sure your netCDF files are ordered in time. First pair of non-ordered files: {dataFiles[id_not_ordered]}, {dataFiles[id_not_ordered + 1]}"
)
return time, time_origin, timeslices, dataFiles
return time, timeslices, dataFiles

@classmethod
def from_netcdf(
Expand Down Expand Up @@ -650,14 +639,14 @@ def from_netcdf(
# Concatenate time variable to determine overall dimension
# across multiple files
if "time" in dimensions or timestamps is not None:
time, time_origin, timeslices, dataFiles = cls._collect_timeslices(
time, timeslices, dataFiles = cls._collect_timeslices(
timestamps, data_filenames, _grid_fb_class, dimensions, indices, netcdf_engine
)
grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
grid = Grid.create_grid(lon, lat, depth, time, mesh=mesh)
grid.timeslices = timeslices
kwargs["dataFiles"] = dataFiles
else: # e.g. for the CROCO CS_w field, see https://github.com/OceanParcels/Parcels/issues/1831
grid = Grid.create_grid(lon, lat, depth, np.array([0.0]), time_origin=TimeConverter(0.0), mesh=mesh)
grid = Grid.create_grid(lon, lat, depth, np.array([0.0]), mesh=mesh)
grid.timeslices = [[0]]
data_filenames = [data_filenames[0]]
elif grid is not None and ("dataFiles" not in kwargs or kwargs["dataFiles"] is None):
Expand Down Expand Up @@ -805,10 +794,7 @@ def from_xarray(
lon = da[dimensions["lon"]].values
lat = da[dimensions["lat"]].values

time_origin = TimeConverter(time[0])
time = time_origin.reltime(time) # type: ignore[assignment]

grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
grid = Grid.create_grid(lon, lat, depth, time, mesh=mesh)
kwargs["time_periodic"] = time_periodic
return cls(
name,
Expand Down Expand Up @@ -1273,7 +1259,7 @@ def write(self, filename, varname=None):
else:
raise NotImplementedError("Field.write only implemented for RectilinearZGrid and CurvilinearZGrid")

attrs = {"units": "seconds since " + str(self.grid.time_origin)} if self.grid.time_origin.calendar else {}
attrs = {} # TODO v4: fix units here
time_counter = xr.DataArray(self.grid.time, dims=["time_counter"], attrs=attrs)
vardata = xr.DataArray(
self.data.reshape((self.grid.tdim, self.grid.zdim, self.grid.ydim, self.grid.xdim)),
Expand Down Expand Up @@ -1340,7 +1326,6 @@ def computeTimeChunk(self, data, tindex):
)
filebuffer.__enter__()
time_data = filebuffer.time
time_data = g.time_origin.reltime(time_data)
filebuffer.ti = (time_data <= g.time[tindex]).argmin() - 1
if self.netcdf_engine != "xarray":
filebuffer.name = filebuffer.parse_name(self.filebuffername)
Expand Down
24 changes: 3 additions & 21 deletions parcels/fieldset.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
from parcels.gridset import GridSet
from parcels.particlefile import ParticleFile
from parcels.tools._helpers import fieldset_repr
from parcels.tools.converters import TimeConverter, convert_xarray_time_units
from parcels.tools.converters import convert_xarray_time_units
from parcels.tools.loggers import logger
from parcels.tools.statuscodes import TimeExtrapolationError
from parcels.tools.warnings import FieldSetWarning
Expand Down Expand Up @@ -43,8 +43,6 @@ def __init__(self, U: Field | NestedField | None, V: Field | NestedField | None,
self._particlefile: ParticleFile | None = None
if U:
self.add_field(U, "U")
# see #1663 for type-ignore reason
self.time_origin = self.U.grid.time_origin if isinstance(self.U, Field) else self.U[0].grid.time_origin # type: ignore
if V:
self.add_field(V, "V")

Expand Down Expand Up @@ -143,14 +141,8 @@ def from_data(
lon = dims["lon"]
lat = dims["lat"]
depth = np.zeros(1, dtype=np.float32) if "depth" not in dims else dims["depth"]
time = np.zeros(1, dtype=np.float64) if "time" not in dims else dims["time"]
time = np.array(time)
if isinstance(time[0], np.datetime64):
time_origin = TimeConverter(time[0])
time = np.array([time_origin.reltime(t) for t in time])
else:
time_origin = kwargs.pop("time_origin", TimeConverter(0))
grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
time = np.zeros(1, dtype=np.datetime64) if "time" not in dims else np.array(dims["time"])
grid = Grid.create_grid(lon, lat, depth, time, mesh=mesh)
if "creation_log" not in kwargs.keys():
kwargs["creation_log"] = "from_data"

Expand Down Expand Up @@ -298,15 +290,6 @@ def check_velocityfields(U, V, W):

for g in self.gridset.grids:
g._check_zonal_periodic()
if len(g.time) == 1:
continue
assert isinstance(
g.time_origin.time_origin, type(self.time_origin.time_origin)
), "time origins of different grids must be have the same type"
g.time = g.time + self.time_origin.reltime(g.time_origin)
if g.defer_load:
g.time_full = g.time_full + self.time_origin.reltime(g.time_origin)
g._time_origin = self.time_origin
self._add_UVfield()

for f in self.get_fields():
Expand Down Expand Up @@ -1506,7 +1489,6 @@ def computeTimeChunk(self, time=0.0, dt=1):
----------
time :
Time around which the FieldSet chunks are to be loaded.
Time is provided as a double, relatively to Fieldset.time_origin.
Default is 0.
dt :
time step of the integration scheme, needed to set the direction of time chunk loading.
Expand Down
53 changes: 14 additions & 39 deletions parcels/grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
import numpy.typing as npt

from parcels._typing import Mesh, UpdateStatus, assert_valid_mesh
from parcels.tools.converters import Geographic, GeographicPolar, TimeConverter, UnitConverter
from parcels.tools.converters import Geographic, GeographicPolar, UnitConverter
from parcels.tools.warnings import FieldSetWarning

__all__ = [
Expand Down Expand Up @@ -46,7 +46,6 @@ def __init__(
lon: npt.NDArray,
lat: npt.NDArray,
time: npt.NDArray | None,
time_origin: TimeConverter | None,
mesh: Mesh,
):
self._ti = -1
Expand All @@ -62,18 +61,11 @@ def __init__(
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.time_full = self.time # needed for deferred_loaded Fields
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._cstruct = None
Expand All @@ -98,7 +90,7 @@ def __repr__(self):
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})"
f"mesh={self.mesh!r})"
)

@property
Expand Down Expand Up @@ -132,10 +124,6 @@ def meridional_halo(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
Expand All @@ -158,7 +146,6 @@ def create_grid(
lat: npt.ArrayLike,
depth,
time,
time_origin,
mesh: Mesh,
**kwargs,
):
Expand All @@ -170,14 +157,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, **kwargs)
return RectilinearZGrid(lon, lat, depth, time, mesh=mesh, **kwargs)
else:
return RectilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return RectilinearSGrid(lon, lat, depth, time, mesh=mesh, **kwargs)
else:
if depth is None or len(depth.shape) <= 1:
return CurvilinearZGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return CurvilinearZGrid(lon, lat, depth, time, mesh=mesh, **kwargs)
else:
return CurvilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return CurvilinearSGrid(lon, lat, depth, time, mesh=mesh, **kwargs)

@property
def ctypes_struct(self):
Expand Down Expand Up @@ -378,14 +365,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)
self.tdim = self.time.size

if self.ydim > 1 and self.lat[-1] < self.lat[0]:
Expand Down Expand Up @@ -465,8 +452,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:
Expand All @@ -476,8 +461,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"

Expand Down Expand Up @@ -514,8 +499,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:
Expand All @@ -531,10 +514,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
Expand Down Expand Up @@ -576,7 +558,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"
Expand All @@ -587,7 +568,7 @@ def __init__(

lon = lon.squeeze()
lat = lat.squeeze()
super().__init__(lon, lat, time, time_origin, mesh)
super().__init__(lon, lat, time, mesh)
self.tdim = self.time.size

@property
Expand Down Expand Up @@ -630,8 +611,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:
Expand All @@ -647,10 +626,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"

Expand Down Expand Up @@ -686,8 +664,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:
Expand All @@ -703,10 +679,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
Expand Down
4 changes: 1 addition & 3 deletions parcels/gridset.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,7 @@ def add_grid(self, field):
existing_grid = True
break
sameGrid = True
if grid.time_origin != g.time_origin:
continue
for attr in ["lon", "lat", "depth", "time"]:
for attr in ["lon", "lat", "depth"]: # HACK removed time because of np.datetime64 support
gattr = getattr(g, attr)
gridattr = getattr(grid, attr)
if gattr.shape != gridattr.shape or not np.allclose(gattr, gridattr):
Expand Down
Loading