Repository: OceanParcels/parcels (python)

Issue:
Implement FieldSet.from_xarray_dataset()
Following the suggestion by @rabernat in #462, it would be good to implement a simple way to convert xarray `Datasets` to Parcels `FieldSets`.

Pull request description:
Implementing Fieldset.from_xarray_dataset() method for integration with xarray Datasets
This PR implements a new method, `FieldSet.from_ds()` to run Parcels directly on xarray `Dataset` objects. Following the suggestion by @rabernat, this fixes #467 

An example script would be
```python
from parcels import FieldSet, AdvectionRK4, ParticleSet, JITParticle
from parcels import plotTrajectoriesFile
import numpy as np
import xarray as xr
import gcsfs
from datetime import timedelta as delta


def periodicBC(particle, fieldset, time, dt):
    if particle.lon < 0:
        particle.lon += 360
    elif particle.lon > 360:
        particle.lon -= 360


# create xarray Dataset of AVISO geostrophic velocities using gcsfs
fname = "pangeo-data/dataset-duacs-rep-global-merged-allsat-phy-l4-v3-alt"
gcsmap = gcsfs.mapping.GCSMap(fname)
ds = xr.open_zarr(gcsmap)

# convert Dataset to Parcels FieldSet
variables = {'U': 'ugos', 'V': 'vgos'}
dimensions = {'lon': 'longitude', 'lat': 'latitude', 'time': 'time'}
fset = FieldSet.from_xarray_dataset(ds, variables, dimensions)
fset.add_periodic_halo(zonal=True, meridional=False, halosize=2)

# Create ParticleSet and output file
lons, lats = np.meshgrid(np.arange(60, 80, 2), np.arange(-50, -30, 2))
pset = ParticleSet(fset, JITParticle, lon=lons, lat=lats)
ofile = pset.ParticleFile('aviso_particles.nc', outputdt=delta(days=2))

# Advect ParticleSet with RK4 and periodic Boundary conditions
pset.execute(AdvectionRK4+pset.Kernel(periodicBC), runtime=delta(days=30),
             dt=delta(hours=1), output_file=ofile)

plotTrajectoriesFile('aviso_particles.nc')
```

The key statement here is `fset = FieldSet.from_xarray_dataset(ds, variables, dimensions)`. Note that users will still need to define what the `U` and `V` Fields are named in the Dataset, and that also the dimensions need to be specified.

Patch:
diff --git a/environment_py2_linux.yml b/environment_py2_linux.yml
index c51502b1..a0882889 100644
--- a/environment_py2_linux.yml
+++ b/environment_py2_linux.yml
@@ -23,6 +23,7 @@ dependencies:
   - scipy>=0.16.0
   - six >=1.10.0
   - xarray>=0.10.8
+  - dask
   - cftime
   - pip:
       - pytest>=2.7.0
diff --git a/environment_py2_osx.yml b/environment_py2_osx.yml
index 8a0a25b3..84a57995 100644
--- a/environment_py2_osx.yml
+++ b/environment_py2_osx.yml
@@ -23,6 +23,7 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.10.8
+  - dask
   - cftime
   - pip:
       - pytest>=2.7.0
diff --git a/environment_py3_linux.yml b/environment_py3_linux.yml
index ac9a13c5..eecb2c43 100644
--- a/environment_py3_linux.yml
+++ b/environment_py3_linux.yml
@@ -23,6 +23,7 @@ dependencies:
   - six >=1.10.0
   - xarray>=0.10.8
   - cftime
+  - dask
   - pip:
       - pytest>=2.7.0
       - nbval
diff --git a/environment_py3_osx.yml b/environment_py3_osx.yml
index 41760b19..393488f1 100644
--- a/environment_py3_osx.yml
+++ b/environment_py3_osx.yml
@@ -22,6 +22,7 @@ dependencies:
   - scipy>=0.16.0
   - six>=1.10.0
   - xarray>=0.10.8
+  - dask
   - cftime
   - pip:
       - pytest>=2.7.0
diff --git a/parcels/examples/example_globcurrent.py b/parcels/examples/example_globcurrent.py
index 2d15e0a0..d511ceda 100644
--- a/parcels/examples/example_globcurrent.py
+++ b/parcels/examples/example_globcurrent.py
@@ -4,29 +4,35 @@ from os import path
 from glob import glob
 import numpy as np
 import pytest
+import xarray as xr
 
 
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 
 
-def set_globcurrent_fieldset(filename=None, indices=None, full_load=False):
+def set_globcurrent_fieldset(filename=None, indices=None, full_load=False, use_xarray=False):
     if filename is None:
         filename = path.join(path.dirname(__file__), 'GlobCurrent_example_data',
                              '20*-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc')
     variables = {'U': 'eastward_eulerian_current_velocity', 'V': 'northward_eulerian_current_velocity'}
     dimensions = {'lat': 'lat', 'lon': 'lon', 'time': 'time'}
-    return FieldSet.from_netcdf(filename, variables, dimensions, indices, full_load=full_load)
+    if use_xarray:
+        ds = xr.open_mfdataset(filename)
+        return FieldSet.from_xarray_dataset(ds, variables, dimensions, indices, full_load=full_load)
+    else:
+        return FieldSet.from_netcdf(filename, variables, dimensions, indices, full_load=full_load)
 
 
-def test_globcurrent_fieldset():
-    fieldset = set_globcurrent_fieldset()
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_globcurrent_fieldset(use_xarray):
+    fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
     assert(fieldset.U.lon.size == 81)
     assert(fieldset.U.lat.size == 41)
     assert(fieldset.V.lon.size == 81)
     assert(fieldset.V.lat.size == 41)
 
     indices = {'lon': [5], 'lat': range(20, 30)}
-    fieldsetsub = set_globcurrent_fieldset(indices=indices)
+    fieldsetsub = set_globcurrent_fieldset(indices=indices, use_xarray=use_xarray)
     assert np.allclose(fieldsetsub.U.lon, fieldset.U.lon[indices['lon']])
     assert np.allclose(fieldsetsub.U.lat, fieldset.U.lat[indices['lat']])
     assert np.allclose(fieldsetsub.V.lon, fieldset.V.lon[indices['lon']])
@@ -34,34 +40,32 @@ def test_globcurrent_fieldset():
 
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
-@pytest.mark.parametrize('dt, substart, subend, lonstart, latstart, irange', [
-    (3600., 0, 3, 25, -35, range(3, 9, 1)),
-    (-3600., 8, 10, 20, -39, range(7, 2, -1))
-])
-def test_globcurrent_fieldset_advancetime(mode, dt, substart, subend, lonstart, latstart, irange):
+@pytest.mark.parametrize('dt, lonstart, latstart', [(3600., 25, -35), (-3600., 20, -39)])
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_globcurrent_fieldset_advancetime(mode, dt, lonstart, latstart, use_xarray):
     basepath = path.join(path.dirname(__file__), 'GlobCurrent_example_data',
                          '20*-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc')
     files = sorted(glob(str(basepath)))
 
-    fieldsetsub = set_globcurrent_fieldset(files[0:10])
+    fieldsetsub = set_globcurrent_fieldset(files[0:10], use_xarray=use_xarray)
     psetsub = ParticleSet.from_list(fieldset=fieldsetsub, pclass=ptype[mode], lon=[lonstart], lat=[latstart])
 
-    fieldsetall = set_globcurrent_fieldset(files[0:10], full_load=True)
+    fieldsetall = set_globcurrent_fieldset(files[0:10], full_load=True, use_xarray=use_xarray)
     psetall = ParticleSet.from_list(fieldset=fieldsetall, pclass=ptype[mode], lon=[lonstart], lat=[latstart])
     if dt < 0:
         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)
-        psetall.execute(AdvectionRK4, runtime=delta(days=1), dt=dt)
+    psetsub.execute(AdvectionRK4, runtime=delta(days=7), dt=dt)
+    psetall.execute(AdvectionRK4, runtime=delta(days=7), dt=dt)
 
     assert abs(psetsub[0].lon - psetall[0].lon) < 1e-4
 
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_globcurrent_particles(mode):
-    fieldset = set_globcurrent_fieldset()
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_globcurrent_particles(mode, use_xarray):
+    fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
 
     lonstart = [25]
     latstart = [-35]
@@ -74,6 +78,22 @@ def test_globcurrent_particles(mode):
     assert(abs(pset[0].lat - -35.3) < 1)
 
 
+@pytest.mark.parametrize('dt', [-300, 300])
+def test_globcurrent_xarray_vs_netcdf(dt):
+    fieldsetNetcdf = set_globcurrent_fieldset(use_xarray=False)
+    fieldsetxarray = set_globcurrent_fieldset(use_xarray=True)
+    lonstart, latstart, runtime = (25, -35, delta(days=7))
+
+    psetN = ParticleSet(fieldsetNetcdf, pclass=JITParticle, lon=lonstart, lat=latstart)
+    psetN.execute(AdvectionRK4, runtime=runtime, dt=dt)
+
+    psetX = ParticleSet(fieldsetxarray, pclass=JITParticle, lon=lonstart, lat=latstart)
+    psetX.execute(AdvectionRK4, runtime=runtime, dt=dt)
+
+    assert np.allclose(psetN[0].lon, psetX[0].lon)
+    assert np.allclose(psetN[0].lat, psetX[0].lat)
+
+
 def test__particles_init_time():
     fieldset = set_globcurrent_fieldset()
 
@@ -92,8 +112,9 @@ def test__particles_init_time():
 
 @pytest.mark.xfail(reason="Time extrapolation error expected to be thrown", strict=True)
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_globcurrent_time_extrapolation_error(mode):
-    fieldset = set_globcurrent_fieldset()
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_globcurrent_time_extrapolation_error(mode, use_xarray):
+    fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
 
     pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[25], lat=[-35],
                        time=fieldset.U.time[0]-delta(days=1).total_seconds())
@@ -103,8 +124,9 @@ def test_globcurrent_time_extrapolation_error(mode):
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 @pytest.mark.parametrize('dt', [-300, 300])
-def test_globcurrent_variable_fromfield(mode, dt):
-    fieldset = set_globcurrent_fieldset()
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_globcurrent_variable_fromfield(mode, dt, use_xarray):
+    fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
 
     class MyParticle(ptype[mode]):
         sample_var = Variable('sample_var', initial=fieldset.U)
@@ -115,8 +137,9 @@ def test_globcurrent_variable_fromfield(mode, dt):
 
 
 @pytest.mark.parametrize('full_load', [True, False])
-def test_globcurrent_deferred_fieldset_gradient(full_load):
-    fieldset = set_globcurrent_fieldset(full_load=full_load)
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_globcurrent_deferred_fieldset_gradient(full_load, use_xarray):
+    fieldset = set_globcurrent_fieldset(full_load=full_load, use_xarray=use_xarray)
     (dU_dx, dU_dy) = fieldset.U.gradient()
     fieldset.add_field(dU_dy)
 
diff --git a/parcels/examples/example_ofam.py b/parcels/examples/example_ofam.py
index f7c318cd..5fdf63b9 100644
--- a/parcels/examples/example_ofam.py
+++ b/parcels/examples/example_ofam.py
@@ -1,30 +1,54 @@
 from parcels import FieldSet, ParticleSet, ScipyParticle, JITParticle, AdvectionRK4
 from datetime import timedelta as delta
 import pytest
+import xarray as xr
+import numpy as np
 from os import path
 
 
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
 
 
-def set_ofam_fieldset(full_load=False):
+def set_ofam_fieldset(full_load=False, use_xarray=False):
     filenames = {'U': path.join(path.dirname(__file__), 'OFAM_example_data', 'OFAM_simple_U.nc'),
                  'V': path.join(path.dirname(__file__), 'OFAM_example_data', 'OFAM_simple_V.nc')}
     variables = {'U': 'u', 'V': 'v'}
     dimensions = {'lat': 'yu_ocean', 'lon': 'xu_ocean', 'depth': 'st_ocean',
                   'time': 'Time'}
-    return FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True, full_load=full_load)
+    if use_xarray:
+        ds = xr.open_mfdataset([filenames['U'], filenames['V']])
+        return FieldSet.from_xarray_dataset(ds, variables, dimensions, allow_time_extrapolation=True, full_load=full_load)
+    else:
+        return FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True, full_load=full_load)
 
 
-def test_ofam_fieldset_fillvalues():
-    fieldset = set_ofam_fieldset(full_load=True)
+@pytest.mark.parametrize('use_xarray', [True, False])
+def test_ofam_fieldset_fillvalues(use_xarray):
+    fieldset = set_ofam_fieldset(full_load=True, use_xarray=use_xarray)
     # V.data[0, 0, 150] is a landpoint, that makes NetCDF4 generate a masked array, instead of an ndarray
     assert(fieldset.V.data[0, 0, 150] == 0)
 
 
+@pytest.mark.parametrize('dt', [delta(minutes=-5), delta(minutes=5)])
+def test_ofam_xarray_vs_netcdf(dt):
+    fieldsetNetcdf = set_ofam_fieldset(use_xarray=False)
+    fieldsetxarray = set_ofam_fieldset(use_xarray=True)
+    lonstart, latstart, runtime = (180, 10, delta(days=7))
+
+    psetN = ParticleSet(fieldsetNetcdf, pclass=JITParticle, lon=lonstart, lat=latstart)
+    psetN.execute(AdvectionRK4, runtime=runtime, dt=dt)
+
+    psetX = ParticleSet(fieldsetxarray, pclass=JITParticle, lon=lonstart, lat=latstart)
+    psetX.execute(AdvectionRK4, runtime=runtime, dt=dt)
+
+    assert np.allclose(psetN[0].lon, psetX[0].lon)
+    assert np.allclose(psetN[0].lat, psetX[0].lat)
+
+
+@pytest.mark.parametrize('use_xarray', [True, False])
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
-def test_ofam_particles(mode):
-    fieldset = set_ofam_fieldset()
+def test_ofam_particles(mode, use_xarray):
+    fieldset = set_ofam_fieldset(use_xarray=use_xarray)
 
     lonstart = [180]
     latstart = [10]
diff --git a/parcels/field.py b/parcels/field.py
index 22003f49..73dfd838 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -155,30 +155,35 @@ class Field(object):
         :param netcdf_engine: engine to use for netcdf reading in xarray. Default is 'netcdf',
                but in cases where this doesn't work, setting netcdf_engine='scipy' could help
         """
-
-        if not isinstance(filenames, Iterable) or isinstance(filenames, str):
-            filenames = [filenames]
-
-        data_filenames = filenames['data'] if type(filenames) is dict else filenames
-        if type(filenames) == dict:
-            for k in filenames.keys():
-                assert k in ['lon', 'lat', 'depth', 'data'], \
-                    'filename dimension keys must be lon, lat, depth or data'
-            assert len(filenames['lon']) == 1
-            if filenames['lon'] != filenames['lat']:
-                raise NotImplementedError('longitude and latitude dimensions are currently processed together from one single file')
-            lonlat_filename = filenames['lon'][0]
-            if 'depth' in dimensions:
-                if len(filenames['depth']) != 1:
-                    raise NotImplementedError('Vertically adaptive meshes not implemented for from_netcdf()')
-                depth_filename = filenames['depth'][0]
-
+        if isinstance(variable, xr.core.dataarray.DataArray):
+            lonlat_filename = variable
+            depth_filename = variable
+            data_filenames = variable
+            netcdf_engine = 'xarray'
         else:
-            lonlat_filename = filenames[0]
-            depth_filename = filenames[0]
+            if not isinstance(filenames, Iterable) or isinstance(filenames, str):
+                filenames = [filenames]
+
+            data_filenames = filenames['data'] if type(filenames) is dict else filenames
+            if type(filenames) == dict:
+                for k in filenames.keys():
+                    assert k in ['lon', 'lat', 'depth', 'data'], \
+                        'filename dimension keys must be lon, lat, depth or data'
+                assert len(filenames['lon']) == 1
+                if filenames['lon'] != filenames['lat']:
+                    raise NotImplementedError('longitude and latitude dimensions are currently processed together from one single file')
+                lonlat_filename = filenames['lon'][0]
+                if 'depth' in dimensions:
+                    if len(filenames['depth']) != 1:
+                        raise NotImplementedError('Vertically adaptive meshes not implemented for from_netcdf()')
+                    depth_filename = filenames['depth'][0]
+            else:
+                lonlat_filename = filenames[0]
+                depth_filename = filenames[0]
+            netcdf_engine = kwargs.pop('netcdf_engine', 'netcdf4')
 
         indices = {} if indices is None else indices.copy()
-        netcdf_engine = kwargs.pop('netcdf_engine', 'netcdf4')
+
         with NetcdfFileBuffer(lonlat_filename, dimensions, indices, netcdf_engine) as filebuffer:
             lon, lat = filebuffer.read_lonlat
             indices = filebuffer.indices
@@ -199,16 +204,22 @@ class Field(object):
         if grid is None:
             # Concatenate time variable to determine overall dimension
             # across multiple files
-            timeslices = []
-            dataFiles = []
-            for fname in data_filenames:
-                with NetcdfFileBuffer(fname, dimensions, indices, netcdf_engine) as filebuffer:
-                    ftime = filebuffer.time
-                    timeslices.append(ftime)
-                    dataFiles.append([fname] * len(ftime))
-            timeslices = np.array(timeslices)
-            time = np.concatenate(timeslices)
-            dataFiles = np.concatenate(np.array(dataFiles))
+            if netcdf_engine == 'xarray':
+                with NetcdfFileBuffer(data_filenames, dimensions, indices, netcdf_engine) as filebuffer:
+                    time = filebuffer.time
+                    timeslices = time if isinstance(time, (list, np.ndarray)) else [time]
+                    dataFiles = data_filenames if isinstance(data_filenames, (list, np.ndarray)) else [data_filenames] * len(time)
+            else:
+                timeslices = []
+                dataFiles = []
+                for fname in data_filenames:
+                    with NetcdfFileBuffer(fname, dimensions, indices, netcdf_engine) as filebuffer:
+                        ftime = filebuffer.time
+                        timeslices.append(ftime)
+                        dataFiles.append([fname] * len(ftime))
+                timeslices = np.array(timeslices)
+                time = np.concatenate(timeslices)
+                dataFiles = np.concatenate(np.array(dataFiles))
             if time.size == 1 and time[0] is None:
                 time[0] = 0
             time_origin = TimeConverter(time[0])
@@ -239,17 +250,20 @@ class Field(object):
                 with NetcdfFileBuffer(fname, dimensions, indices, netcdf_engine) as filebuffer:
                     # If Field.from_netcdf is called directly, it may not have a 'data' dimension
                     # In that case, assume that 'name' is the data dimension
-                    filebuffer.name = filebuffer.parse_name(dimensions, variable)
+                    if netcdf_engine == 'xarray':
+                        tslice = [tslice]
+                    else:
+                        filebuffer.name = filebuffer.parse_name(dimensions, variable)
 
-                    if len(filebuffer.dataset[filebuffer.name].shape) == 2:
-                        data[ti:ti+len(tslice), 0, :, :] = filebuffer.data[:, :]
-                    elif len(filebuffer.dataset[filebuffer.name].shape) == 3:
+                    if len(filebuffer.data.shape) == 2:
+                        data[ti:ti+len(tslice), 0, :, :] = filebuffer.data
+                    elif len(filebuffer.data.shape) == 3:
                         if len(filebuffer.indices['depth']) > 1:
-                            data[ti:ti+len(tslice), :, :, :] = filebuffer.data[:, :, :]
+                            data[ti:ti+len(tslice), :, :, :] = filebuffer.data
                         else:
-                            data[ti:ti+len(tslice), 0, :, :] = filebuffer.data[:, :, :]
+                            data[ti:ti+len(tslice), 0, :, :] = filebuffer.data
                     else:
-                        data[ti:ti+len(tslice), :, :, :] = filebuffer.data[:, :, :, :]
+                        data[ti:ti+len(tslice), :, :, :] = filebuffer.data
                 ti += len(tslice)
         else:
             grid.defer_load = True
@@ -264,6 +278,7 @@ class Field(object):
         kwargs['time_periodic'] = time_periodic
         kwargs['netcdf_engine'] = netcdf_engine
 
+        variable = kwargs['var_name'] if netcdf_engine == 'xarray' else variable
         return cls(variable, data, grid=grid,
                    allow_time_extrapolation=allow_time_extrapolation, **kwargs)
 
@@ -866,21 +881,20 @@ class Field(object):
     def computeTimeChunk(self, data, tindex):
         g = self.grid
         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
             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[:, :]
-            elif len(filebuffer.dataset[filebuffer.name].shape) == 3:
+            filebuffer.ti = (time_data <= g.time[tindex]).argmin() - 1
+            if self.netcdf_engine != 'xarray':
+                filebuffer.name = filebuffer.parse_name(self.dimensions, self.name)
+            if len(filebuffer.data.shape) == 2:
+                data[tindex, 0, :, :] = filebuffer.data
+            elif len(filebuffer.data.shape) == 3:
                 if g.zdim > 1:
-                    data[tindex, :, :, :] = filebuffer.data[:, :, :]
+                    data[tindex, :, :, :] = filebuffer.data
                 else:
-                    data[tindex, 0, :, :] = filebuffer.data[ti, :, :]
+                    data[tindex, 0, :, :] = filebuffer.data
             else:
-                data[tindex, :, :, :] = filebuffer.data[ti, :, :, :]
-
-        return data
+                data[tindex, :, :, :] = filebuffer.data
 
     def __add__(self, field):
         if isinstance(self, Field) and isinstance(field, Field):
@@ -1139,12 +1153,16 @@ class NetcdfFileBuffer(object):
 
     def __init__(self, filename, dimensions, indices, netcdf_engine):
         self.filename = filename
-        self.dimensions = dimensions  # Dict with dimension keyes for file data
+        self.dimensions = dimensions  # Dict with dimension keys for file data
         self.indices = indices
         self.dataset = None
         self.netcdf_engine = netcdf_engine
+        self.ti = None
 
     def __enter__(self):
+        if self.netcdf_engine == 'xarray':
+            self.dataset = self.filename
+            return self
         try:
             self.dataset = xr.open_dataset(str(self.filename), decode_cf=True, engine=self.netcdf_engine)
             self.dataset['decoded'] = True
@@ -1159,7 +1177,10 @@ class NetcdfFileBuffer(object):
         return self
 
     def __exit__(self, type, value, traceback):
-        self.dataset.close()
+        if self.netcdf_engine == 'xarray':
+            pass
+        else:
+            self.dataset.close()
 
     def parse_name(self, dimensions, variable):
         name = dimensions['data'] if 'data' in dimensions else variable
@@ -1228,16 +1249,20 @@ class NetcdfFileBuffer(object):
 
     @property
     def data(self):
-        data = self.dataset[self.name]
+        if self.netcdf_engine == 'xarray':
+            data = self.dataset
+        else:
+            data = self.dataset[self.name]
+        ti = range(data.shape[0]) if self.ti is None else self.ti
         if len(data.shape) == 2:
             data = data[self.indices['lat'], self.indices['lon']]
         elif len(data.shape) == 3:
             if len(self.indices['depth']) > 1:
                 data = data[self.indices['depth'], self.indices['lat'], self.indices['lon']]
             else:
-                data = data[:, self.indices['lat'], self.indices['lon']]
+                data = data[ti, self.indices['lat'], self.indices['lon']]
         else:
-            data = data[:, self.indices['depth'], self.indices['lat'], self.indices['lon']]
+            data = data[ti, self.indices['depth'], self.indices['lat'], self.indices['lon']]
 
         if np.ma.is_masked(data):  # convert masked array to ndarray
             data = np.ma.filled(data, np.nan)
@@ -1247,7 +1272,7 @@ class NetcdfFileBuffer(object):
     def time(self):
         try:
             time_da = self.dataset[self.dimensions['time']]
-            if self.dataset['decoded'] and 'Unit' not in time_da.attrs:
+            if self.netcdf_engine != 'xarray' and (self.dataset['decoded'] and 'Unit' not in time_da.attrs):
                 time = np.array([time_da]) if len(time_da.shape) == 0 else np.array(time_da)
             else:
                 if 'units' not in time_da.attrs and 'Unit' in time_da.attrs:
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index 6232ce73..012cd2df 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -331,6 +331,51 @@ class FieldSet(object):
                                dimensions=dimensions, allow_time_extrapolation=allow_time_extrapolation,
                                time_periodic=time_periodic, full_load=full_load, **kwargs)
 
+    @classmethod
+    def from_xarray_dataset(cls, ds, variables, dimensions, indices=None, mesh='spherical', allow_time_extrapolation=None,
+                            time_periodic=False, full_load=False, **kwargs):
+        """Initialises FieldSet data from xarray Datasets.
+
+        :param ds: xarray Dataset
+        :param dimensions: Dictionary mapping data dimensions (lon,
+               lat, depth, time, data) to dimensions in the xarray Dataset.
+               Note that dimensions can also be a dictionary of dictionaries if
+               dimension names are different for each variable
+               (e.g. dimensions['U'], dimensions['V'], etc).
+        :param indices: Optional dictionary of indices for each dimension
+               to read from file(s), to allow for reading of subset of data.
+               Default is to read the full extent of each dimension.
+        :param mesh: 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.
+        :param allow_time_extrapolation: boolean whether to allow for extrapolation
+               (i.e. beyond the last available time snapshot)
+               Default is False if dimensions includes time, else True
+        :param time_periodic: boolean whether to loop periodically over the time component of the FieldSet
+               This flag overrides the allow_time_interpolation and sets it to False
+        :param full_load: boolean whether to fully load the data or only pre-load them. (default: False)
+               It is advised not to fully load the data, since in that case Parcels deals with
+               a better memory management during particle set execution.
+               full_load is however sometimes necessary for plotting the fields.
+        """
+
+        fields = {}
+        for var, name in variables.items():
+
+            # Use dimensions[var] and indices[var] if either of them is a dict of dicts
+            dims = dimensions[var] if var in dimensions else dimensions
+            inds = indices[var] if (indices and var in indices) else indices
+
+            fields[var] = Field.from_netcdf(None, ds[name], dimensions=dims, indices=inds, grid=None, mesh=mesh,
+                                            allow_time_extrapolation=allow_time_extrapolation, var_name=var,
+                                            time_periodic=time_periodic, full_load=full_load, **kwargs)
+        u = fields.pop('U', None)
+        v = fields.pop('V', None)
+        return cls(u, v, fields=fields)
+
     @property
     def fields(self):
         """Returns a list of all the :class:`parcels.field.Field` objects
@@ -446,7 +491,7 @@ class FieldSet(object):
                 data = np.empty((g.tdim, g.zdim, g.ydim-2*g.meridional_halo, g.xdim-2*g.zonal_halo), dtype=np.float32)
                 f.loaded_time_indices = range(3)
                 for tind in f.loaded_time_indices:
-                    data = f.computeTimeChunk(data, tind)
+                    f.computeTimeChunk(data, tind)
                 f.data = f.reshape(data)
             elif g.update_status == 'updated':
                 data = np.empty((g.tdim, g.zdim, g.ydim-2*g.meridional_halo, g.xdim-2*g.zonal_halo), dtype=np.float32)
@@ -456,7 +501,7 @@ class FieldSet(object):
                 else:
                     f.data[1:, :] = f.data[:2, :]
                     f.loaded_time_indices = [0]
-                data = f.computeTimeChunk(data, f.loaded_time_indices[0])
+                f.computeTimeChunk(data, f.loaded_time_indices[0])
                 f.data[f.loaded_time_indices[0], :] = f.reshape(data)[f.loaded_time_indices[0], :]
             else:
                 f.loaded_time_indices = []


Test patch:
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index d4f897f2..bebcda16 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -4,6 +4,7 @@ from parcels.tools.converters import TimeConverter
 from datetime import timedelta as delta
 import datetime
 import numpy as np
+import xarray as xr
 import math
 import pytest
 from os import path
@@ -365,3 +366,31 @@ def test_fieldset_defer_loading_function(zdim, scale_fac, tmpdir, filename='test
     pset.execute(DoNothing, dt=3600)
     assert np.allclose(fieldset.U.data, scale_fac*(zdim-1.)/zdim)
     assert np.allclose(dFdx.data, 0)
+
+
+@pytest.mark.parametrize('maxlatind', [3, pytest.param(2, marks=pytest.mark.xfail(strict=True))])
+def test_fieldset_from_xarray(maxlatind):
+    def generate_dataset(xdim, ydim, zdim=1, tdim=1):
+        lon = np.linspace(0., 12, xdim, dtype=np.float32)
+        lat = np.linspace(0., 12, ydim, dtype=np.float32)
+        depth = np.linspace(0., 20., zdim, dtype=np.float32)
+        time = np.linspace(0., 10, tdim, dtype=np.float64)
+        Uxr = np.ones((tdim, zdim, ydim, xdim), dtype=np.float32)
+        Vxr = np.ones((tdim, zdim, ydim, xdim), dtype=np.float32)
+        for t in range(Uxr.shape[0]):
+            Uxr[t, :, :, :] = t/10.
+        coords = {'lat': lat, 'lon': lon, 'depth': depth, 'time': time}
+        dims = ('time', 'depth', 'lat', 'lon')
+        return xr.Dataset({'Uxr': xr.DataArray(Uxr, coords=coords, dims=dims),
+                           'Vxr': xr.DataArray(Vxr, coords=coords, dims=dims)})
+
+    ds = generate_dataset(3, 3, 2, 10)
+    variables = {'U': 'Uxr', 'V': 'Vxr'}
+    dimensions = {'lat': 'lat', 'lon': 'lon', 'depth': 'depth', 'time': 'time'}
+    indices = {'lat': range(0, maxlatind)}
+    fieldset = FieldSet.from_xarray_dataset(ds, variables, dimensions, indices, mesh='flat')
+
+    pset = ParticleSet(fieldset, JITParticle, 0, 0)
+
+    pset.execute(AdvectionRK4, dt=1)
+    assert pset[0].lon == 4.5 and pset[0].lat == 10