Fix field plotting

environment_py3_linux.yml

@@ -28,4 +28,5 @@ dependencies:
   - nbval
   - scikit-learn
   - pykdtree
+  - cartopy
   - zarr

environment_py3_win.yml

@@ -25,4 +25,5 @@ dependencies:
   - pytest
   - nbval
   - pykdtree
+  - cartopy
   - zarr

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,21 +268,26 @@ 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:
             logger.info("Visualisation of field with geographic coordinates is not possible. Cartopy not found.")
             return None, None, None, None  # creating axes was not possible
 
         projection = cartopy.crs.PlateCarree(central_longitude) if projection == 'PlateCarree' else projection
+        # projection = '3d' if use3D else projection
         fig, ax = plt.subplots(1, 1, subplot_kw={'projection': projection})
         try:  # gridlines not supported for all projections
             if isinstance(projection, cartopy.crs.PlateCarree) and central_longitude != 0:

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 = ''
@@ -71,7 +71,8 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
     if mode == '3d':
         from mpl_toolkits.mplot3d import Axes3D  # noqa
         plt.clf()  # clear the figure
-        ax = fig.gca(projection='3d')
+        # ax = fig.gca(projection='3d')
+        ax = fig.gca()
         for p in range(len(lon)):
             ax.plot(lon[p, :], lat[p, :], z[p, :], '.-')
         ax.set_xlabel('Longitude')
@@ -116,7 +117,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))

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)

tests/test_scripts.py

@@ -21,15 +21,15 @@ 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)
 
     fp = dir.join("DelayParticle.nc")
     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)