Merge branch 'master' into fix-mitgcm-time-dimension

Author: erikvansebille

environment_py2_linux.yml

@@ -23,6 +23,8 @@ dependencies:
   - scipy>=0.16.0
   - six >=1.10.0
   - xarray>=0.5.1
+  - netcdftime
+  - cftime
   - pip:
       - pytest>=2.7.0
       - nbval

environment_py2_osx.yml

@@ -23,6 +23,8 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.5.1
+  - netcdftime
+  - cftime
   - pip:
       - pytest>=2.7.0
       - nbval

environment_py2_win.yml

@@ -22,6 +22,8 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.5.1
+  - netcdftime
+  - cftime
   - pip:
       - pytest>=2.7.0
       - nbval

environment_py3_linux.yml

@@ -22,6 +22,8 @@ dependencies:
   - scipy>=0.16.0
   - six >=1.10.0
   - xarray>=0.5.1
+  - netcdftime
+  - cftime
   - pip:
       - pytest>=2.7.0
       - nbval

environment_py3_osx.yml

@@ -22,6 +22,8 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.5.1
+  - netcdftime
+  - cftime
   - pip:
       - pytest>=2.7.0
       - nbval

environment_py3_win.yml

@@ -21,6 +21,8 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.5.1
+  - netcdftime
+  - cftime
   - pip:
       - pytest>=2.7.0
       - nbval

parcels/examples/example_globcurrent.py

@@ -49,8 +49,8 @@ def test_globcurrent_fieldset_advancetime(mode, dt, substart, subend, lonstart,
     fieldsetall = set_globcurrent_fieldset(files[0:10], full_load=True)
     psetall = ParticleSet.from_list(fieldset=fieldsetall, pclass=ptype[mode], lon=[lonstart], lat=[latstart])
     if dt < 0:
-        psetsub[0].time = fieldsetsub.U.time[-1]
-        psetall[0].time = fieldsetall.U.time[-1]
+        psetsub[0].time = fieldsetsub.U.grid.time[-1]
+        psetall[0].time = fieldsetall.U.grid.time[-1]
 
     for i in irange:
         psetsub.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)
@@ -108,7 +108,7 @@ def test_globcurrent_variable_fromfield(mode, dt):
 
     class MyParticle(ptype[mode]):
         sample_var = Variable('sample_var', initial=fieldset.U)
-    time = fieldset.U.time[0] if dt > 0 else fieldset.U.time[-1]
+    time = fieldset.U.grid.time[0] if dt > 0 else fieldset.U.grid.time[-1]
     pset = ParticleSet(fieldset, pclass=MyParticle, lon=[25], lat=[-35], time=time)
 
     pset.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)

parcels/field.py

@@ -1,5 +1,6 @@
 from parcels.tools.loggers import logger
 from parcels.tools.converters import unitconverters_map, UnitConverter, Geographic, GeographicPolar
+from parcels.tools.converters import TimeConverter
 from parcels.tools.error import FieldSamplingError, TimeExtrapolationError
 from collections import Iterable
 from py import path
@@ -44,7 +45,7 @@ class Field(object):
     :param transpose: Transpose data to required (lon, lat) layout
     :param vmin: Minimum allowed value on the field. Data below this value are set to zero
     :param vmax: Maximum allowed value on the field. Data above this value are set to zero
-    :param time_origin: Time origin (datetime or np.datetime64 object) of the time axis (only if grid is None)
+    :param time_origin: Time origin (TimeConverter object) of the time axis (only if grid is None)
     :param interp_method: Method for interpolation. Either 'linear' or 'nearest'
     :param allow_time_extrapolation: boolean whether to allow for extrapolation in time
            (i.e. beyond the last available time snapshot)
@@ -57,6 +58,7 @@ def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=N
                  interp_method='linear', allow_time_extrapolation=None, time_periodic=False, **kwargs):
         self.name = name
         self.data = data
+        time_origin = TimeConverter(0) if time_origin is None else time_origin
         if grid:
             self.grid = grid
         else:
@@ -68,7 +70,6 @@ def __init__(self, name, data, lon=None, lat=None, depth=None, time=None, grid=N
         self.lon = self.grid.lon
         self.lat = self.grid.lat
         self.depth = self.grid.depth
-        self.time = self.grid.time
         fieldtype = self.name if fieldtype is None else fieldtype
         if self.grid.mesh == 'flat' or (fieldtype not in unitconverters_map.keys()):
             self.units = UnitConverter()
@@ -181,15 +182,12 @@ def from_netcdf(cls, filenames, variable, dimensions, indices=None, grid=None,
             timeslices = np.array(timeslices)
             time = np.concatenate(timeslices)
             dataFiles = np.concatenate(np.array(dataFiles))
-            if isinstance(time[0], np.datetime64):
-                time_origin = time[0]
-                time = (time - time_origin) / np.timedelta64(1, 's')
-            else:
-                time_origin = None
-            assert(np.all((time[1:]-time[:-1]) > 0))
-
             if time.size == 1 and time[0] is None:
                 time[0] = 0
+            time_origin = TimeConverter(time[0])
+            time = time_origin.reltime(time)
+            assert(np.all((time[1:]-time[:-1]) > 0))
+
             if len(lon.shape) == 1:
                 if len(depth.shape) == 1:
                     grid = RectilinearZGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
@@ -815,7 +813,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 else {}
+        attrs = {'units': 'seconds since ' + str(self.grid.time_origin)} if self.grid.time_origin.calendar else {}
         time_counter = xr.DataArray(self.grid.time,
                                     dims=['time_counter'],
                                     attrs=attrs)
@@ -842,9 +840,7 @@ def computeTimeChunk(self, data, tindex):
         with NetcdfFileBuffer(self.dataFiles[g.ti+tindex], self.dimensions, self.indices, self.netcdf_engine) as filebuffer:
             filebuffer.name = filebuffer.parse_name(self.dimensions, self.name)
             time_data = filebuffer.time
-            if isinstance(time_data[0], np.datetime64):
-                assert isinstance(time_data[0], type(g.time_origin)), ('Field %s stores times as dates, but time_origin is not defined ' % self.name)
-                time_data = (time_data - g.time_origin) / np.timedelta64(1, 's')
+            time_data = g.time_origin.reltime(time_data)
             ti = (time_data <= g.time[tindex]).argmin() - 1
             if len(filebuffer.dataset[filebuffer.name].shape) == 2:
                 data[tindex, 0, :, :] = filebuffer.data[:, :]

parcels/fieldset.py

@@ -2,6 +2,7 @@
 from parcels.gridset import GridSet
 from parcels.grid import RectilinearZGrid
 from parcels.tools.loggers import logger
+from parcels.tools.converters import TimeConverter
 import numpy as np
 from os import path
 from glob import glob
@@ -77,7 +78,7 @@ def from_data(cls, data, dimensions, transpose=False, mesh='spherical',
             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']
-            grid = RectilinearZGrid(lon, lat, depth, time, mesh=mesh)
+            grid = RectilinearZGrid(lon, lat, depth, time, time_origin=TimeConverter(), mesh=mesh)
 
             fields[name] = Field(name, datafld, grid=grid, transpose=transpose,
                                  allow_time_extrapolation=allow_time_extrapolation, time_periodic=time_periodic, **kwargs)
@@ -124,11 +125,10 @@ def check_complete(self):
             if len(g.time) == 1:
                 continue
             assert isinstance(g.time_origin, type(self.time_origin)), 'time origins of different grids must be have the same type'
-            if g.time_origin:
-                g.time = g.time + (g.time_origin - self.time_origin) / np.timedelta64(1, 's')
-                if g.defer_load:
-                    g.time_full = g.time_full + (g.time_origin - self.time_origin) / np.timedelta64(1, 's')
-                g.time_origin = self.time_origin
+            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
         if not hasattr(self, 'UV'):
             if isinstance(self.U, SummedField):
                 self.add_vector_field(SummedVectorField('UV', self.U, self.V))
@@ -184,7 +184,8 @@ def from_netcdf(cls, filenames, variables, dimensions, indices=None,
             if not isinstance(paths, list):
                 paths = sorted(glob(str(paths)))
             if len(paths) == 0:
-                raise IOError("FieldSet files not found: %s" % str(paths))
+                notfound_paths = filenames[var] if type(filenames) is dict else filenames
+                raise IOError("FieldSet files not found: %s" % str(notfound_paths))
             for fp in paths:
                 if not path.exists(fp):
                     raise IOError("FieldSet file not found: %s" % str(fp))

parcels/grid.py

@@ -1,8 +1,8 @@
 from parcels.tools.loggers import logger
+from parcels.tools.converters import TimeConverter
 import numpy as np
 from ctypes import Structure, c_int, c_float, c_double, POINTER, cast, c_void_p, pointer
 from enum import IntEnum
-import datetime
 
 __all__ = ['GridCode', 'RectilinearZGrid', 'RectilinearSGrid', 'CurvilinearZGrid', 'CurvilinearSGrid', 'CGrid']
 
@@ -39,11 +39,8 @@ def __init__(self, lon, lat, time, time_origin, mesh):
             assert isinstance(self.time[0], (np.integer, np.floating, float, int)), 'Time vector must be an array of int or floats'
             logger.warning_once("Casting time data to np.float64")
             self.time = self.time.astype(np.float64)
-        self.time_origin = time_origin
-        if self.time_origin:
-            if isinstance(self.time_origin, datetime.datetime):
-                self.time_origin = np.datetime64(self.time_origin)
-            assert isinstance(self.time_origin, np.datetime64), 'If defined, time_origin must be a datetime.datetime or a np.datetime64'
+        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'
         self.mesh = mesh
         self.cstruct = None
         self.cell_edge_sizes = {}
@@ -214,7 +211,7 @@ class RectilinearZGrid(RectilinearGrid):
     :param depth: Vector containing the vertical coordinates of the grid, which are z-coordinates.
            The depth of the different layers is thus constant.
     :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (datetime or np.datetime64 object) of the time axis
+    :param time_origin: Time origin (TimeConverter object) of the time axis
     :param mesh: String indicating the type of mesh coordinates and
            units used during velocity interpolation:
 
@@ -251,7 +248,7 @@ class RectilinearSGrid(RectilinearGrid):
            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].
     :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (datetime or np.datetime64 object) of the time axis
+    :param time_origin: Time origin (TimeConverter object) of the time axis
     :param mesh: String indicating the type of mesh coordinates and
            units used during velocity interpolation:
 
@@ -343,7 +340,7 @@ class CurvilinearZGrid(CurvilinearGrid):
     :param depth: Vector containing the vertical coordinates of the grid, which are z-coordinates.
            The depth of the different layers is thus constant.
     :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (datetime or np.datetime64 object) of the time axis
+    :param time_origin: Time origin (TimeConverter object) of the time axis
     :param mesh: String indicating the type of mesh coordinates and
            units used during velocity interpolation:
 
@@ -379,7 +376,7 @@ class CurvilinearSGrid(CurvilinearGrid):
            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].
     :param time: Vector containing the time coordinates of the grid
-    :param time_origin: Time origin (datetime or np.datetime64 object) of the time axis
+    :param time_origin: Time origin (TimeConverter object) of the time axis
     :param mesh: String indicating the type of mesh coordinates and
            units used during velocity interpolation: