diff --git a/environment_py3_linux.yml b/environment_py3_linux.yml
index bbc9be01d..6beee1212 100644
--- a/environment_py3_linux.yml
+++ b/environment_py3_linux.yml
@@ -28,4 +28,5 @@ dependencies:
   - nbval
   - scikit-learn
   - pykdtree
+  - cartopy
   - zarr
diff --git a/environment_py3_win.yml b/environment_py3_win.yml
index f4ecf30fd..93e4c9698 100644
--- a/environment_py3_win.yml
+++ b/environment_py3_win.yml
@@ -25,4 +25,5 @@ dependencies:
   - pytest
   - nbval
   - pykdtree
+  - cartopy
   - zarr
diff --git a/parcels/plotting.py b/parcels/plotting.py
index dc68d7106..603657c3f 100644
--- a/parcels/plotting.py
+++ b/parcels/plotting.py
@@ -2,6 +2,7 @@
 from datetime import timedelta as delta
 
 import numpy as np
+import dask.array as da
 import copy
 
 from parcels.field import Field
@@ -101,7 +102,7 @@ def plotparticles(particles, with_particles=True, show_time=None, field=None, do
 
 
 def plotfield(field, show_time=None, domain=None, depth_level=0, projection='PlateCarree', land=True,
-              vmin=None, vmax=None, savefile=None, **kwargs):
+              vmin=None, vmax=None, savefile=None, use3D=False, **kwargs):
     """Function to plot a Parcels Field
 
     :param show_time: Time in seconds from start after which to show the Field
@@ -113,6 +114,7 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla
     :param vmax: maximum colour scale (only in single-plot mode)
     :param savefile: Name of a file to save the plot to
     :param animation: Boolean whether result is a single plot, or an animation
+    :param use3D: tells if requested projection is 2D or 3D
     """
 
     if type(field) is VectorField:
@@ -130,7 +132,7 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla
         logger.warning('Field.show() does not always correctly determine the domain for curvilinear grids. '
                        'Use plotting with caution and perhaps use domain argument as in the NEMO 3D tutorial')
 
-    plt, fig, ax, cartopy = create_parcelsfig_axis(spherical, land, projection=projection, cartopy_features=kwargs.pop('cartopy_features', []))
+    plt, fig, ax, cartopy = create_parcelsfig_axis(spherical, land, projection=projection, cartopy_features=kwargs.pop('cartopy_features', []), use3D=use3D)
     if plt is None:
         return None, None, None, None  # creating axes was not possible
 
@@ -165,6 +167,8 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla
                 data[i] = np.squeeze(fld.data)[depth_level, latS:latN, lonW:lonE]
             else:
                 data[i] = np.squeeze(fld.data)[latS:latN, lonW:lonE]
+        if isinstance(data[i], da.Array):
+            data[i] = np.array(data[i])
 
     if plottype == 'vector':
         if field[0].interp_method == 'cgrid_velocity':
@@ -264,14 +268,18 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection='Pla
     return plt, fig, ax, cartopy
 
 
-def create_parcelsfig_axis(spherical, land=True, projection='PlateCarree', central_longitude=0, cartopy_features=[]):
+def create_parcelsfig_axis(spherical, land=True, projection='PlateCarree', central_longitude=0, cartopy_features=[], use3D=False):
     try:
         import matplotlib.pyplot as plt
     except:
         logger.info("Visualisation is not possible. Matplotlib not found.")
         return None, None, None, None  # creating axes was not possible
 
-    if spherical and projection:
+    if use3D:
+        cartopy = None
+        fig, ax = plt.subplots(1, 1, subplot_kw={'projection': '3d'})
+        ax.grid()
+    elif spherical and projection:
         try:
             import cartopy
         except:
diff --git a/parcels/scripts/plottrajectoriesfile.py b/parcels/scripts/plottrajectoriesfile.py
index f215c43db..825e59d79 100644
--- a/parcels/scripts/plottrajectoriesfile.py
+++ b/parcels/scripts/plottrajectoriesfile.py
@@ -53,13 +53,13 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
 
     if tracerfile is not None and mode != 'hist2d':
         tracerfld = Field.from_netcdf(tracerfile, tracerfield, {'lon': tracerlon, 'lat': tracerlat})
-        plt, fig, ax, cartopy = plotfield(tracerfld)
+        plt, fig, ax, cartopy = plotfield(tracerfld, use3D=(mode == '3d'))
         if plt is None:
             return  # creating axes was not possible
         titlestr = ' and ' + tracerfield
     else:
         spherical = False if mode == '3d' or mesh == 'flat' else True
-        plt, fig, ax, cartopy = create_parcelsfig_axis(spherical=spherical, central_longitude=central_longitude)
+        plt, fig, ax, cartopy = create_parcelsfig_axis(spherical=spherical, central_longitude=central_longitude, use3D=(mode == '3d'))
         if plt is None:
             return  # creating axes was not possible
         titlestr = ''
@@ -116,7 +116,11 @@ def timestr(plottimes, index):
                 return str(plottimes[index])
 
         if cartopy:
-            scat = ax.scatter(lon[b], lat[b], s=20, color='k', transform=cartopy.crs.Geodetic())
+            scat = None
+            try:
+                scat = ax.scatter(lon[b], lat[b], s=20, color='k', transform=cartopy.crs.Geodetic())
+            except (ValueError,):
+                scat = ax.scatter(lon[b], lat[b], s=20, color='k', transform=cartopy.crs.PlateCarree())
         else:
             scat = ax.scatter(lon[b], lat[b], s=20, color='k')
         ttl = ax.set_title('Particles' + titlestr + ' at time ' + timestr(plottimes, 0))
diff --git a/tests/test_plotting.py b/tests/test_plotting.py
new file mode 100644
index 000000000..5d9812c01
--- /dev/null
+++ b/tests/test_plotting.py
@@ -0,0 +1,54 @@
+from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4, ErrorCode
+from datetime import timedelta as delta
+import numpy as np
+import dask.array as da  # NOQA
+from os import path
+from matplotlib import pyplot as plt  # NOQA
+import pytest
+
+
+def periodicBC(particle, fieldset, time):
+    if particle.lon > 180:
+        particle.lon -= 360
+
+
+def test_field_from_netcdf():
+    data_path = path.join(path.dirname(__file__), 'test_data/')
+
+    filenames = {'U': {'lon': data_path + 'mask_nemo_cross_180lon.nc',
+                       'lat': data_path + 'mask_nemo_cross_180lon.nc',
+                       'data': data_path + 'Uu_eastward_nemo_cross_180lon.nc'},
+                 'V': {'lon': data_path + 'mask_nemo_cross_180lon.nc',
+                       'lat': data_path + 'mask_nemo_cross_180lon.nc',
+                       'data': data_path + 'Vv_eastward_nemo_cross_180lon.nc'}
+                 }
+    variables = {'U': 'U',
+                 'V': 'V'}
+    dimensions = {'lon': 'glamf', 'lat': 'gphif'}
+    return FieldSet.from_netcdf(filenames, variables, dimensions, interp_method='cgrid_velocity', chunksize='auto', allow_time_extrapolation=True)
+
+
+@pytest.fixture(name="fieldset")
+def fieldset_fixture():
+    return test_field_from_netcdf()
+
+
+def test_pset_create_field(fieldset, npart=100):
+    lonp = -180 * np.ones(npart)
+    latp = [i for i in np.linspace(-70, 88, npart)]
+    pset = ParticleSet.from_list(fieldset, JITParticle, lon=lonp, lat=latp)
+    return pset
+
+
+def DeleteParticle(particle, fieldset, time):
+    particle.delete()
+
+
+if __name__ == '__main__':
+    fset = test_field_from_netcdf()
+    print(fset)
+    pset = test_pset_create_field(fset)
+    kernels = pset.Kernel(AdvectionRK4) + periodicBC
+    pset.execute(kernels, dt=delta(hours=1), output_file=None,
+                 recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
+    pset.show(field=fset.U)
diff --git a/tests/test_scripts.py b/tests/test_scripts.py
index 757a104a0..05fe9b276 100644
--- a/tests/test_scripts.py
+++ b/tests/test_scripts.py
@@ -21,7 +21,7 @@ def create_outputfiles(dir, pset_mode):
     npart = 10
     delaytime = delta(hours=1)
     endtime = delta(hours=24)
-    x = 3. * (1. / 1.852 / 60)
+    x = 3. * (1. / 1.852 / 60.)
     y = (fieldset.U.lat[0] + x, fieldset.U.lat[-1] - x)
     lat = np.linspace(y[0], y[1], npart)
 
@@ -29,7 +29,7 @@ def create_outputfiles(dir, pset_mode):
     output_file = pset.ParticleFile(name=fp, outputdt=delaytime)
 
     for t in range(npart):
-        time = 0 if len(pset) == 0 else pset[0].time
+        time = 0. if len(pset) == 0 else pset[0].time
         pset.add(pset_type[pset_mode]['pset'](pclass=JITParticle, lon=x, lat=lat[t], fieldset=fieldset, time=time))
         pset.execute(AdvectionRK4, runtime=delaytime, dt=delta(minutes=5), output_file=output_file)