Ability to distribute filtering over dask workers

[angus-g]

Since Parcels' change to MPI-based parallelism, and moving away from OpenMP on the Parcels fork, we lost the ability to parallelise filtering. This PR adds a new metaclass which is responsible for setting up the filtering workflow to distribute over dask workers.

Closes #61.

Issues

The mechanism here is to chunk the data grid, and compute each chunk as a separate task. Because Parcels is highly stateful, we need to create new objects for each task. In particular, we end up recompiling the kernel several times, when it could (and should) be shared between workers.

Dask operates under an IO-limiting assumption, and prefers to transfer "live" chunks of data between workers, rather than reload them from disk. This ends up putting a lot of extra overhead on workers and generally doesn't work very well in this case. Instead, this method relies on being given a function that recreates the input data as a (delayed) xarray dataset, so that workers don't share the input data. This doesn't seem ideal, but it's hard to imagine a way around this at the moment.

[angus-g]

This is a little entangled, but there are some performance tweaks that need to go here to make it feasible to load deferred data from xarray. We also need a patch on Parcels:

diff --git a/parcels/field.py b/parcels/field.py
index 66be9287..f653e6a5 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -436,7 +436,7 @@ class Field(object):
 
     @classmethod
     def from_xarray(cls, da, name, dimensions, mesh='spherical', allow_time_extrapolation=None,
-                    time_periodic=False, **kwargs):
+                    time_periodic=False, deferred_load=False, **kwargs):
         """Create field from xarray Variable
 
         :param da: Xarray DataArray
@@ -467,6 +467,11 @@ class Field(object):
         time = time_origin.reltime(time)
 
         grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
+
+        if deferred_load:
+            grid.defer_load = True
+            grid.ti = -1
+
         return cls(name, data, grid=grid, allow_time_extrapolation=allow_time_extrapolation,
                    interp_method=interp_method, **kwargs)
 
@@ -1091,7 +1096,7 @@ class Field(object):
             for block_id in range(len(self.grid.load_chunk)):
                 if self.grid.load_chunk[block_id] == 1 or self.grid.load_chunk[block_id] > 1 and self.data_chunks[block_id] is None:
                     block = self.get_block(block_id)
-                    self.data_chunks[block_id] = np.array(self.data.blocks[(slice(self.grid.tdim),) + block])
+                    self.data_chunks[block_id] = np.array(self.data.blocks[(slice(self.grid.ti, self.grid.ti + self.grid.tdim),) + block])
                 elif self.grid.load_chunk[block_id] == 0:
                     if isinstance(self.data_chunks, list):
                         self.data_chunks[block_id] = None
@@ -2201,7 +2206,7 @@ class NetcdfFileBuffer(object):
                             self.chunking_finalized = True
                     else:
                         # ==== I think this can be "pass" too ==== #
-                        data = data.rechunk(self.chunk_mapping)
+                        #data = data.rechunk(self.chunk_mapping)
                         self.chunking_finalized = True
             else:
                 da_data = da.from_array(data, chunks=self.field_chunksize)
diff --git a/parcels/fieldset.py b/parcels/fieldset.py
index 251178de..e417057b 100644
--- a/parcels/fieldset.py
+++ b/parcels/fieldset.py
@@ -1009,7 +1009,7 @@ class FieldSet(object):
                                 f.data_chunks[block_id][0] = np.array(f.data.blocks[(slice(3),)+block][0])
         # do user-defined computations on fieldset data
         if self.compute_on_defer:
-            self.compute_on_defer(self)
+            self.compute_on_defer(self, signdt)
 
         # update time varying grid depth
         for f in self.get_fields():