Adding support for depth_level argument in Field.show()

parcels/field.py

@@ -756,22 +756,23 @@ class CField(Structure):
                          pointer(self.grid.ctypes_struct))
         return cstruct
 
-    def show(self, animation=False, show_time=None, domain=None, projection=None, land=True,
+    def show(self, animation=False, show_time=None, domain=None, depth_level=0, projection=None, land=True,
              vmin=None, vmax=None, savefile=None, **kwargs):
         """Method to 'show' a Parcels Field
 
         :param animation: Boolean whether result is a single plot, or an animation
         :param show_time: Time at which to show the Field (only in single-plot mode)
         :param domain: Four-vector (latN, latS, lonE, lonW) defining domain to show
+        :param depth_level: depth level to be plotted (default 0)
         :param projection: type of cartopy projection to use (default PlateCarree)
         :param land: Boolean whether to show land. This is ignored for flat meshes
         :param vmin: minimum colour scale (only in single-plot mode)
         :param vmax: maximum colour scale (only in single-plot mode)
         :param savefile: Name of a file to save the plot to
         """
         from parcels.plotting import plotfield
-        plt, _, _, _ = plotfield(self, animation=animation, show_time=show_time, domain=domain, projection=projection,
-                                 land=land, vmin=vmin, vmax=vmax, savefile=savefile, **kwargs)
+        plt, _, _, _ = plotfield(self, animation=animation, show_time=show_time, domain=domain, depth_level=depth_level,
+                                 projection=projection, land=land, vmin=vmin, vmax=vmax, savefile=savefile, **kwargs)
         if plt:
             plt.show()
 

parcels/plotting.py

@@ -93,12 +93,13 @@ def plotparticles(particles, with_particles=True, show_time=None, field=None, do
         plt.close()
 
 
-def plotfield(field, show_time=None, domain=None, projection=None, land=True,
+def plotfield(field, show_time=None, domain=None, depth_level=0, projection=None, land=True,
               vmin=None, vmax=None, savefile=None, **kwargs):
     """Function to plot a Parcels Field
 
     :param show_time: Time at which to show the Field
     :param domain: Four-vector (latN, latS, lonE, lonW) defining domain to show
+    :param depth_level: depth level to be plotted (default 0)
     :param projection: type of cartopy projection to use (default PlateCarree)
     :param land: Boolean whether to show land. This is ignored for flat meshes
     :param vmin: minimum colour scale (only in single-plot mode)
@@ -144,9 +145,15 @@ def plotfield(field, show_time=None, domain=None, projection=None, land=True,
         if i > 0 and not np.allclose(plotlon[i], plotlon[0]):
             raise RuntimeError('VectorField needs to be on an A-grid for plotting')
         if fld.grid.time.size > 1:
-            data[i] = np.squeeze(fld.temporal_interpolate_fullfield(idx, show_time))[latS:latN, lonW:lonE]
+            if fld.grid.zdim > 1:
+                data[i] = np.squeeze(fld.temporal_interpolate_fullfield(idx, show_time))[depth_level, latS:latN, lonW:lonE]
+            else:
+                data[i] = np.squeeze(fld.temporal_interpolate_fullfield(idx, show_time))[latS:latN, lonW:lonE]
         else:
-            data[i] = np.squeeze(fld.data)[latS:latN, lonW:lonE]
+            if fld.grid.zdim > 1:
+                data[i] = np.squeeze(fld.data)[depth_level, latS:latN, lonW:lonE]
+            else:
+                data[i] = np.squeeze(fld.data)[latS:latN, lonW:lonE]
 
     if plottype is 'vector':
         spd = data[0] ** 2 + data[1] ** 2
@@ -186,10 +193,20 @@ def plotfield(field, show_time=None, domain=None, projection=None, land=True,
 
     timestr = parsetimestr(field[0].grid.time_origin, show_time)
     titlestr = kwargs.pop('titlestr', '')
+    if field[0].grid.zdim > 1:
+        if field[0].grid.gtype in [GridCode.CurvilinearZGrid, GridCode.RectilinearZGrid]:
+            gphrase = 'depth'
+            depth_or_level = field[0].grid.depth[depth_level]
+        else:
+            gphrase = 'level'
+            depth_or_level = depth_level
+        depthstr = ' at %s %g ' % (gphrase, depth_or_level)
+    else:
+        depthstr = ''
     if plottype is 'vector':
-        ax.set_title(titlestr + 'Velocity field' + timestr)
+        ax.set_title(titlestr + 'Velocity field' + depthstr + timestr)
     else:
-        ax.set_title(titlestr + field[0].name + timestr)
+        ax.set_title(titlestr + field[0].name + depthstr + timestr)
 
     if not spherical:
         ax.set_xlabel('Zonal distance [m]')