Update plotting routines in QC script (#355)

* Add an option to plot the ice thickness data on a map

* Update the t-test script to include a new option to plot the mean ice thickness data for both the baseline and test dataset.

* Modify the QC script to always try to plot the data on a map.  Add new optional plot_type argument, which can plot the data using scatter, contour, or pcolor.  The default map (if no plot_type argument is passed) is pcolor.

* Add a brief paragraph about the plotting functionality of the QC script

* Fix typo in documentation

configuration/scripts/tests/QC/cice.t-test.py

@@ -358,48 +358,102 @@ def skill_test(path_a, fname, data_a, data_b, num_files, hemisphere):
         logger.info('Quadratic Skill Test Failed for %s Hemisphere', hemisphere)
         return False
 
-def plot_data(data, lat, lon, units):
-    '''This function plots CICE data and creates a .png file (ice_thickness_map.png).'''
+def plot_data(data, lat, lon, units, case, plot_type):
+    '''This function plots CICE data and creates a .png file.'''
 
     try:
-        logger.info('Creating map of the data (ice_thickness_map.png)')
         # Load the necessary plotting libraries
         import matplotlib.pyplot as plt
         from mpl_toolkits.basemap import Basemap
         from mpl_toolkits.axes_grid1 import make_axes_locatable
     except ImportError:
         logger.warning('Error loading necessary Python modules in plot_data function')
+        return
 
     # Suppress Matplotlib deprecation warnings
     import warnings
     warnings.filterwarnings("ignore", category=UserWarning)
 
     # Create the figure and axis
-    fig = plt.figure(figsize=(12, 8))
-    ax = fig.add_axes([0.05, 0.08, 0.9, 0.9])
+    fig, axes = plt.subplots(nrows=1, ncols=2,figsize=(14, 8))
 
+    # Plot the northern hemisphere data as a scatter plot
     # Create the basemap, and draw boundaries
-    m = Basemap(projection='kav7', lon_0=180., resolution='l')
-    m.drawmapboundary(fill_color='white')
+    plt.sca(axes[0])
+    m = Basemap(projection='npstere', boundinglat=35,lon_0=270, resolution='l')
     m.drawcoastlines()
+    m.fillcontinents()
     m.drawcountries()
 
-    # Plot the data as a scatter plot
-    x, y = m(lon, lat)
-    sc = m.scatter(x, y, c=data, cmap='jet', lw=0)
+    if plot_type == 'scatter':
+        x, y = m(lon,lat)
+        sc = m.scatter(x, y, c=data, cmap='jet', lw=0, s=4)
+    else:
+        # Create new arrays to add 1 additional longitude value to prevent a 
+        # small amount of whitespace around longitude of 0/360 degrees.
+        lon_cyc = np.zeros((lon.shape[0],lon.shape[1]+1))
+        mask = np.zeros((data.shape[0],data.shape[1]+1))
+        lat_cyc = np.zeros((lat.shape[0],lat.shape[1]+1))
+
+        mask[:,0:-1] = data.mask[:,:]
+        mask[:,-1] = data.mask[:,0]
+        lon_cyc[:,0:-1] = lon[:,:]; lon_cyc[:,-1] = lon[:,0]
+        lat_cyc[:,0:-1] = lat[:,:]; lat_cyc[:,-1] = lat[:,0]
+
+        lon1 = np.ma.masked_array(lon_cyc, mask=mask)
+        lat1 = np.ma.masked_array(lat_cyc, mask=mask)
+
+        d = np.zeros((data.shape[0],data.shape[1]+1))
+        d[:,0:-1] = data[:,:]
+        d[:,-1] = data[:,0]
+        d1 = np.ma.masked_array(d,mask=mask)
+
+        x, y = m(lon1.data, lat1.data)
+
+        if plot_type == 'contour':
+            sc = m.contourf(x, y, d1, cmap='jet')
+        else:  # pcolor
+            sc = m.pcolor(x, y, d1, cmap='jet')
+
+    m.drawparallels(np.arange(-90.,120.,15.),labels=[1,0,0,0]) # draw parallels
+    m.drawmeridians(np.arange(0.,420.,30.),labels=[1,1,1,1]) # draw meridians
+
+    # Plot the southern hemisphere data as a scatter plot
+    plt.sca(axes[1])
+    m = Basemap(projection='spstere', boundinglat=-45,lon_0=270, resolution='l')
+    m.drawcoastlines()
+    m.fillcontinents()
+    m.drawcountries()
+
+    if plot_type == 'scatter':
+        x, y = m(lon,lat)
+        sc = m.scatter(x, y, c=data, cmap='jet', lw=0, s=4)
+    else:
+        x, y = m(lon1.data, lat1.data)
 
-    m.drawmeridians(np.arange(0, 360, 60), labels=[0, 0, 0, 1], fontsize=10)
-    m.drawparallels(np.arange(-90, 90, 30), labels=[1, 0, 0, 0], fontsize=10)
+        # Bandaid for a bug in the version of Basemap used during development
+        outside = (x <= m.xmin) | (x >= m.xmax) | (y <= m.ymin) | (y >= m.ymax)
+        tmp = np.ma.masked_where(outside,d1)
 
-    plt.title('CICE Ice Thickness')
+        if plot_type == 'contour':
+            sc = m.contourf(x, y, tmp, cmap='jet')
+        else:  # pcolor
+            sc = m.pcolor(x, y, tmp, cmap='jet')
 
-    # Create the colorbar and add Pass / Fail labels
-    divider = make_axes_locatable(ax)
-    cax = divider.append_axes("bottom", size="5%", pad=0.5)
-    cb = plt.colorbar(sc, cax=cax, orientation="horizontal", format="%.2f")
+    m.drawparallels(np.arange(-90.,120.,15.),labels=[1,0,0,0]) # draw parallels
+    m.drawmeridians(np.arange(0.,420.,30.),labels=[1,1,1,1]) # draw meridians
+
+    plt.suptitle('CICE Mean Ice Thickness\n{}'.format(case), y=0.95)
+
+    # Make some room at the bottom of the figure, and create a colorbar
+    fig.subplots_adjust(bottom=0.2)
+    cbar_ax = fig.add_axes([0.11,0.1,0.8,0.05])
+    cb = plt.colorbar(sc, cax=cbar_ax, orientation="horizontal", format="%.2f")
     cb.set_label(units, x=1.0)
 
-    plt.savefig('ice_thickness_map.png', dpi=300)
+    outfile = 'ice_thickness_{}.png'.format(case.replace('\n- ','_minus_'))
+    logger.info('Creating map of the data ({})'.format(outfile))
+    plt.savefig(outfile, dpi=300, bbox_inches='tight')
 
 def plot_two_stage_failures(data, lat, lon):
     '''This function plots each grid cell and whether or not it Passed or Failed
@@ -428,7 +482,7 @@ def plot_two_stage_failures(data, lat, lon):
         ax = fig.add_axes([0.05, 0.08, 0.9, 0.9])
 
         # Create the basemap, and draw boundaries
-        m = Basemap(projection='kav7', lon_0=180., resolution='l')
+        m = Basemap(projection='moll', lon_0=0., resolution='l')
         m.drawmapboundary(fill_color='white')
         m.drawcoastlines()
         m.drawcountries()
@@ -440,7 +494,7 @@ def plot_two_stage_failures(data, lat, lon):
 
         # Plot the data as a scatter plot
         x, y = m(lon, lat)
-        sc = m.scatter(x, y, c=int_data, cmap=cm, lw=0, vmin=0, vmax=1)
+        sc = m.scatter(x, y, c=int_data, cmap=cm, lw=0, vmin=0, vmax=1, s=4)
 
         m.drawmeridians(np.arange(0, 360, 60), labels=[0, 0, 0, 1], fontsize=10)
         m.drawparallels(np.arange(-90, 90, 30), labels=[1, 0, 0, 0], fontsize=10)
@@ -482,8 +536,11 @@ def main():
                 help='Path to the test history (iceh_inst*) files.  REQUIRED')
     parser.add_argument('-v', '--verbose', dest='verbose', help='Print debug output?', \
                         action='store_true')
+    parser.add_argument('-pt','--plot_type', dest='plot_type', help='Specify type of plot \
+                        to create', choices=['scatter','contour','pcolor'])
 
     parser.set_defaults(verbose=False)
+    parser.set_defaults(plot_type='pcolor')
 
     # If no arguments are provided, print the help message
     if len(sys.argv) == 1:
@@ -528,6 +585,17 @@ def main():
     # If test failed, attempt to create a plot of the failure locations
     if not PASSED:
         plot_two_stage_failures(H1_array, t_lat, t_lon)
+
+        # Create plots of mean ice thickness
+        baseDir = os.path.abspath(args.base_dir).rstrip('history/').rstrip(\
+                                                        'history').split('/')[-1]
+        testDir = os.path.abspath(args.test_dir).rstrip('history/').rstrip( \
+                                                        'history').split('/')[-1]
+        plot_data(np.mean(data_base,axis=0), t_lat, t_lon, 'm', baseDir, args.plot_type)
+        plot_data(np.mean(data_test,axis=0), t_lat, t_lon, 'm', testDir, args.plot_type)
+        plot_data(np.mean(data_base-data_test,axis=0), t_lat, t_lon, 'm', '{}\n- {}'.\
+                  format(baseDir,testDir), args.plot_type)
+
         logger.error('Quality Control Test FAILED')
         sys.exit(-1)
 
@@ -559,6 +627,16 @@ def main():
 
     PASSED_SKILL = PASSED_NH and PASSED_SH
 
+    # Plot the ice thickness data for the base and test cases
+    baseDir = os.path.abspath(args.base_dir).rstrip('history/').rstrip( \
+                                                    'history').split('/')[-1]
+    testDir = os.path.abspath(args.test_dir).rstrip('history/').rstrip( \
+                                                    'history').split('/')[-1]
+    plot_data(np.mean(data_base,axis=0), t_lat, t_lon, 'm', baseDir, args.plot_type)
+    plot_data(np.mean(data_test,axis=0), t_lat, t_lon, 'm', testDir, args.plot_type)
+    plot_data(np.mean(data_base-data_test,axis=0), t_lat, t_lon, 'm', '{}\n- {}'.\
+              format(baseDir,testDir), args.plot_type)
+
     logger.info('')
     if not PASSED_SKILL:
         logger.error('Quality Control Test FAILED')

doc/source/user_guide/ug_testing.rst

@@ -847,7 +847,7 @@ To install the necessary Python packages, the ``pip`` Python utility can be used
   pip install --user matplotlib
 
 To run the compliance test, setup a baseline run with the original baseline model and then 
-a perturbation run based on recent model changes.  Use ``--sets qc`` in both runs in addition
+a perturbation run based on recent model changes.  Use ``--set qc`` in both runs in addition
 to other settings needed.  Then use the QC script to compare history output,
 
 .. code-block:: bash
@@ -871,11 +871,17 @@ Implementation notes: 1) Provide a pass/fail on each of the confidence
 intervals, 2) Facilitate output of a bitmap for each test so that
 locations of failures can be identified.
 
-The cice.t-test.py requires memory to store multiple two-dimensional fields spanning 
+The ``cice.t-test.py`` requires memory to store multiple two-dimensional fields spanning 
 1825 unique timesteps, a total of several GB.  An appropriate resource is needed to 
 run the script.  If the script runs out of memory on an interactive resource, try
 logging into a batch resource or finding a large memory node.
 
+The ``cice.t-test.py`` script will also attempt to generate plots of the mean ice thickness
+for both the baseline and test cases. Additionally, if the 2-stage test fails then the 
+script will attempt to plot a map showing the grid cells that failed the test.  For a 
+full list of options, run ``python cice.t-test.py -h``.
+
+
 
 End-To-End Testing Procedure
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~