Merge pull request NCAR#483 from dustinswales/feature/beef_up_plotting

Updates to CI and plotting

.github/workflows/ci_run_scm_rts.yml

@@ -188,7 +188,7 @@ jobs:
     - name: Download SCM RT baselines
       run: |
         cd ${dir_bl}
-        wget ftp://ftp.rap.ucar.edu:/pub/ccpp-scm/rt-baselines-${{matrix.build-type}}.zip
+        wget https://dtcenter.ucar.edu/ccpp/users/rt/rt-baselines-${{matrix.build-type}}.zip
         unzip rt-baselines-${{matrix.build-type}}.zip
 
     - name: Compare SCM RT output to baselines

test/cmp_rt2bl.py

@@ -14,22 +14,20 @@
 
 #
 parser = argparse.ArgumentParser()
-parser.add_argument('-drt', '--dir_rt',  help='Directory containing SCM RT output')
-parser.add_argument('-dbl', '--dir_bl',  help='Directory containing SCM RT baselines')
-parser.add_argument('-plt', '--do_plot', help='If true, create plots of SCM RT output', action='store_true')
+parser.add_argument('-drt',  '--dir_rt',  help='Directory containing SCM RT output',      required=True)
+parser.add_argument('-dbl',  '--dir_bl',  help='Directory containing SCM RT baselines',   required=True)
 
 #
 def parse_args():
     args    = parser.parse_args()
     dir_rt  = args.dir_rt 
     dir_bl  = args.dir_bl
-    do_plot = args.do_plot
-    return (dir_rt, dir_bl, do_plot)
+    return (dir_rt, dir_bl)
 
 #
 def main():
     #
-    (dir_rt, dir_bl, do_plot) = parse_args()
+    (dir_rt, dir_bl) = parse_args()
 
     #
     error_count = 0
@@ -48,7 +46,7 @@ def main():
 
             # Create plots between RTs and baselines (only if differences exist)
             if (result != 0):
-                plot_files = plot_results(file_rt, file_bl, do_plot)
+                plot_files = plot_results(file_bl, file_rt)
 
                 # Setup output directories for plots.
                 result = os.system("mkdir -p scm_rt_out/"+run["case"]+"/"+run["suite"])

test/cmp_scmout.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+
+##############################################################################
+#
+# This script compares SCM RT output to baselines.
+#
+##############################################################################
+import os
+import sys
+from os.path import exists
+import argparse
+from plot_scm_out import plot_results
+
+#
+parser = argparse.ArgumentParser()
+parser.add_argument('-fbl',  '--file_bl', help='File containing SCM RT baselines', required=True)
+parser.add_argument('-frt',  '--file_rt', help='File containing SCM RT output')
+
+#
+def parse_args():
+    args    = parser.parse_args()
+    file_rt = args.file_rt
+    file_bl = args.file_bl
+    return (file_bl, file_rt)
+
+#
+def main():
+    #
+    (file_bl, file_rt) = parse_args()
+
+    plot_files = plot_results(file_bl, file_rt)
+
+    # Put plots in local directory (scm_plots/)
+    result = os.system("mkdir -p scm_plots/")
+    com = "mv"
+    for plot_file in plot_files:
+        com = com + " " + plot_file
+    # end for
+    result = os.system(com+" scm_plots/")
+
+    # Housekeeping
+    if file_rt is not None:
+        result = os.system('tar -cvf scm_out_abs.tar scm_plots/*.png')
+        result = os.system('tar -cvf /scratch1/data_untrusted/Dustin.Swales/scm_out_abs.tar scm_plots/*.png')
+    else:
+        result = os.system('tar -cvf scm_out_diff.tar scm_plots/*.png')
+        result = os.system('tar -cvf /scratch1/data_untrusted/Dustin.Swales/scm_out_diff.tar scm_plots/*.png')
+    # end if
+    result = os.system('rm -rf scm_plots/')
+#
+if __name__ == '__main__':
+    main()

test/plot_scm_out.py

@@ -1,56 +1,73 @@
 #!/usr/bin/env python
 
-##############################################################################
+###############################################################################
 #
-# This script compares SCM RT output to baselines.
+# This script compares SCM output from two simulations.
+# Defined here as Baseline (BL) and Regression Test (RT).
 #
-##############################################################################
+###############################################################################
 from netCDF4 import Dataset
 import numpy as np
 import matplotlib.pyplot as plt
 #
-def plot_results(file_BL, file_RT, plot_RT_only):
-    # List of SCM output fields to plot (This is everything)
-    vars2plot = ["time_inst","time_diag","time_swrad","time_lwrad","time_rad","pres",    \
-                 "pres_i","sigma","sigma_i","pres_s","qv","T","u","v","ql","qi","qc",    \
-                 "qv_force_tend","T_force_tend","u_force_tend","v_force_tend","w_ls",    \
-                 "u_g","v_g","dT_dt_rad_forc","h_advec_thil","h_advec_qt",               \
-                 "v_advec_thil","v_advec_qt","T_s","lhf","shf","tprcp_inst",             \
-                 "tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u","tau_v","upd_mf",  \
-                 "dwn_mf","det_mf","sfc_up_lw_land","sfc_up_lw_ice","sfc_up_lw_water",   \
-                 "sfc_up_sw_dir_nir","sfc_up_sw_dif_nir","sfc_up_sw_dir_vis",            \
-                 "sfc_up_sw_dif_vis","sfc_dwn_sw_dir_nir","sfc_dwn_sw_dif_nir",          \
-                 "sfc_dwn_sw_dir_vis","sfc_dwn_sw_dif_vis","mp_prcp_inst",               \
-                 "dcnv_prcp_inst","scnv_prcp_inst","pwat","dT_dt_lwrad","dT_dt_swrad",   \
-                 "dT_dt_pbl","dT_dt_deepconv","dT_dt_shalconv","dT_dt_micro",            \
-                 "dT_dt_ogwd","dT_dt_cgwd","dT_dt_phys","dT_dt_nonphys","dq_dt_pbl",     \
-                 "dq_dt_deepconv","dq_dt_shalconv","dq_dt_micro","dq_dt_phys",           \
-                 "dq_dt_nonphys","doz_dt_pbl","doz_dt_prodloss","doz_dt_oz","doz_dt_T",  \
-                 "doz_dt_ovhd","doz_dt_phys","doz_dt_nonphys","du_dt_pbl","du_dt_ogwd",  \
-                 "du_dt_deepconv","du_dt_cgwd","du_dt_shalconv","du_dt_phys",            \
-                 "du_dt_nonphys","dv_dt_pbl","dv_dt_ogwd","dv_dt_deepconv","dv_dt_cgwd", \
-                 "dv_dt_shalconv","dv_dt_phys","dv_dt_nonphys","sfc_dwn_sw","sfc_up_sw", \
-                 "sfc_net_sw","sfc_dwn_lw","gflux","u10m","v10m","hpbl","tprcp_accum",   \
-                 "ice_accum","snow_accum","graupel_accum","conv_prcp_accum",             \
-                 "tprcp_rate_accum","ice_rate_accum","snow_rate_accum",                  \
-                 "graupel_rate_accum","conv_prcp_rate_accum","max_cloud_fraction",       \
-                 "toa_total_albedo","vert_int_lwp_mp","vert_int_iwp_mp",                 \
-                 "vert_int_lwp_cf","vert_int_iwp_cf"]
-    # Use subset of available SCM output for plots.
-    vars2plot = ["qv","T","u","v","ql","qi","qc","sfc_dwn_sw","sfc_up_sw",               \
-                 "sfc_net_sw","sfc_dwn_lw", "u10m","v10m","hpbl","gflux",                \
-                 "qv_force_tend","T_force_tend","u_force_tend","v_force_tend","w_ls",    \
-                 "u_g","v_g","dT_dt_rad_forc","h_advec_thil","h_advec_qt",               \
-                 "v_advec_thil","v_advec_qt","T_s","lhf","shf","tprcp_inst",             \
-                 "tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u","tau_v"]
+def plot_results(file_bl, file_rt=None, vars2plt=None):
+#def plot_results(file_bl, file_rt, plot_bl, plot_rt, plot_all, debug):
+    # List of SCM output fields to plot
+    vars2plot_ALL = \
+        ["pres", "pres_i","sigma","sigma_i","pres_s","qv","T","u","v","ql",   \
+         "qi","qc","qv_force_tend","T_force_tend","u_force_tend",             \
+         "v_force_tend","w_ls","u_g","v_g","dT_dt_rad_forc","h_advec_thil",   \
+         "h_advec_qt", "v_advec_thil","v_advec_qt","T_s","lhf","shf",         \
+         "tprcp_inst","tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u",   \
+         "tau_v","upd_mf","dwn_mf","det_mf","sfc_up_lw_land","sfc_up_lw_ice", \
+         "sfc_up_lw_water","sfc_up_sw_dir_nir","sfc_up_sw_dif_nir",           \
+         "sfc_up_sw_dir_vis","sfc_up_sw_dif_vis","sfc_dwn_sw_dir_nir",        \
+         "sfc_dwn_sw_dif_nir","sfc_dwn_sw_dir_vis","sfc_dwn_sw_dif_vis",      \
+         "mp_prcp_inst","dcnv_prcp_inst","scnv_prcp_inst","pwat",             \
+         "dT_dt_lwrad","dT_dt_swrad","dT_dt_pbl","dT_dt_deepconv",            \
+         "dT_dt_shalconv","dT_dt_micro","dT_dt_ogwd","dT_dt_cgwd",            \
+         "dT_dt_phys","dT_dt_nonphys","dq_dt_pbl","dq_dt_deepconv",           \
+         "dq_dt_shalconv","dq_dt_micro","dq_dt_phys","dq_dt_nonphys",         \
+         "doz_dt_pbl","doz_dt_prodloss","doz_dt_oz","doz_dt_T","doz_dt_ovhd", \
+         "doz_dt_phys","doz_dt_nonphys","du_dt_pbl","du_dt_ogwd","dv_dt_cgwd",\
+         "du_dt_deepconv","du_dt_cgwd","du_dt_shalconv","du_dt_phys",         \
+         "du_dt_nonphys","dv_dt_pbl","dv_dt_ogwd","dv_dt_deepconv",           \
+         "dv_dt_shalconv","dv_dt_phys","dv_dt_nonphys","sfc_dwn_sw",          \
+         "sfc_up_sw","sfc_net_sw","sfc_dwn_lw","gflux","u10m","v10m","hpbl",  \
+         "tprcp_accum","ice_accum","snow_accum","graupel_accum",              \
+         "conv_prcp_accum","tprcp_rate_accum","ice_rate_accum",               \
+         "snow_rate_accum","graupel_rate_accum","conv_prcp_rate_accum",       \
+         "max_cloud_fraction","toa_total_albedo","vert_int_lwp_mp",           \
+         "vert_int_iwp_mp","vert_int_lwp_cf","vert_int_iwp_cf"]
+
+    # Smaller subset of SCM outputs to plot.
+    vars2plot_SUB = \
+        ["qv","T","u","v","ql","qi","qc","sfc_dwn_sw","sfc_up_sw",            \
+         "sfc_net_sw","sfc_dwn_lw", "u10m","v10m","hpbl","gflux",             \
+         "qv_force_tend","T_force_tend","u_force_tend","v_force_tend","w_ls", \
+         "u_g","v_g","dT_dt_rad_forc","h_advec_thil","h_advec_qt",            \
+         "v_advec_thil","v_advec_qt","T_s","lhf","shf","tprcp_inst",          \
+         "tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u","tau_v"]
+    vars2plot_DEBUG = \
+        ["qv","u10m"]
 
-    # Open SCM datasets
-    SCM_BL = Dataset(file_BL)
-    SCM_RT = Dataset(file_RT)
+    # Which fields to plot? (default is subset of full fields)
+    if vars2plt is None:
+        vars2plot = vars2plot_SUB
+    else:
+        vars2plot = vars2plt
+    # end if
 
-    plot_diff = True
-    if plot_RT_only:
-        plot_diff = False
+    # Only plot differences if two files provided
+    plot_diff = False
+    if file_rt is not None:
+        plot_diff = True
+    # end if
+
+    # Open SCM datasets
+    SCM_BL = Dataset(file_bl)
+    if file_rt is not None:
+        SCM_RT = Dataset(file_rt)
     # end if
 
     plot_files = []
@@ -62,7 +79,13 @@ def plot_results(file_BL, file_RT, plot_RT_only):
             timeD = SCM_BL[var].dimensions[0]
             timeD = timeD[0:len(timeD)-4]
             x1    = SCM_BL[timeD][:].squeeze()/3600. #seconds -> hours
-            x2    = SCM_RT[timeD][:].squeeze()/3600. #seconds - >hours
+            if file_rt is not None:
+                x2    = SCM_RT[timeD][:].squeeze()/3600. #seconds - >hours
+                # If temporal dimensions disagree, con't compute deltas from experiments, turn off difference plots.
+                if (len(x1) != len(x2)):
+                    plot_diff = False
+                # end if
+            # end if
             # Is this a 2D (time, x) variable? (Really 1D since x=1 in SCM)
             is2D  = False
             if (len(SCM_BL[var].dimensions)==2):
@@ -72,90 +95,88 @@ def plot_results(file_BL, file_RT, plot_RT_only):
             if (len(SCM_BL[var].shape) != 3):
                 if (is2D):
                     y1 = SCM_BL[var][:,0].squeeze()
-                    y2 = SCM_RT[var][:,0].squeeze()
+                    if file_rt is not None:
+                        y2 = SCM_RT[var][:,0].squeeze()
+                    # end if
                 else:
                     y1 = SCM_BL[var][:]
-                    y2 = SCM_RT[var][:]
+                    if file_rt is not None:
+                        y2 = SCM_RT[var][:]
+                    # end if
                 # endif
 
                 # Make figure
                 if (np.size(x1) > 1):
                     fig = plt.figure(figsize=(13,10))
-                    # Baselines and SCM RTs on same plot
-                    if plot_diff:
-                        plt.subplot(2,1,1)
-                    # end if
+
+                    # Baselines and RTs on same plot
+                    if plot_diff: plt.subplot(2,1,1)
                     plt.title(SCM_BL[var].description)
-                    if plot_RT_only:
-                        plt.plot(x1, y1,  color='blue')
-                    else:
-                        plt.plot(x1, y1,  color='blue')
-                        plt.plot(x2, y2,  color='black')
-                    # end if
+                    plt.plot(x1, y1,  color='blue')
+                    if plot_diff: plt.plot(x2, y2,  color='black')
                     plt.ylabel('('+SCM_BL[var].units+')')
                     plt.xlabel('(hours)')
 
-                    # Difference (Baseline-SCMRT)
+                    # Difference (Baseline-MRT)
                     if plot_diff:
                         plt.subplot(2,1,2)
                         plt.title("Difference (blue - black)")
                         plt.plot(x1, y1 - y2,  color='red')
                         plt.plot(x1, np.zeros(len(y1)), color='grey',linestyle='dashed')
-                        plt.ylabel('('+SCM_RT[var].units+')')
+                        plt.ylabel('('+SCM_BL[var].units+')')
                         plt.xlabel('(hours)')
-                    #
+                    # Save figure
                     fileOUT = 'scm.' + var +'.png'
                     plt.savefig(fileOUT)
-                    #
                     plot_files.append(fileOUT)
             # three-dimensional variables
             elif len(SCM_BL[var].shape) == 3:
                 z1 = np.transpose(SCM_BL[var][:,:,0]).squeeze()
-                z2 = np.transpose(SCM_RT[var][:,:,0]).squeeze()
+                if file_rt is not None:
+                    z2 = np.transpose(SCM_RT[var][:,:,0]).squeeze()
+                # end if
 
                 # vertical axis
                 y1 = SCM_BL["pres"][0,:].squeeze()*0.01
-                y2 = SCM_RT["pres"][0,:].squeeze()*0.01
+                if file_rt is not None:
+                    y2 = SCM_RT["pres"][0,:].squeeze()*0.01
+                # end if
                 nlev = SCM_BL[var][:,:,0].shape[1]
                 # Layer (nlev) quantities are the default, so check for case where we have an
                 # interface (nlev+1) variable to plot.
                 if (SCM_BL[var][:,:,0].shape[1] > len(y1)):
                     y1 = SCM_BL["pres_i"][0,:].squeeze()*0.01
-                    y2 = SCM_RT["pres_i"][0,:].squeeze()*0.01
+                    if file_rt is not None:
+                        y2 = SCM_RT["pres_i"][0,:].squeeze()*0.01
+                    # end if
                 # endif
 
-                # Comppute differences and determine valid plot range(s).
-                dz = z1-z2
-
                 # Finally, make figure.
                 if (np.size(x1) > 1):
                     fig = plt.figure(figsize=(13,10))
-                    if plot_RT_only:
-                        plt.contourf(x2, y2, z2, 20, cmap='YlGnBu')
-                        plt.ylim(1000,200)
-                        plt.ylabel('(Pa)')
-                        plt.xlabel('(hours)')
-                        cbr = plt.colorbar()
-                        cbr.set_label('('+SCM_RT[var].units+')')
-                    else:
-                        # Baselines
-                        plt.subplot(3,1,1)
-                        plt.title(SCM_BL[var].description, fontsize=12)
-                        plt.contourf(x1, y1, z1, 20, cmap='YlGnBu')
-                        plt.ylim(1000,200)
-                        plt.ylabel('(Pa)')
-                        cbr = plt.colorbar()
-                        cbr.set_label('('+SCM_RT[var].units+')')
+                    if file_rt is not None: plt.subplot(3,1,1)
+                    plt.contourf(x1, y1, z1, 20, cmap='YlGnBu')
+                    plt.ylim(1000,200)
+                    plt.xlim(0,np.max(x1))
+                    plt.ylabel('(Pa)')
+                    plt.xlabel('(hours)')
+                    cbr = plt.colorbar()
+                    cbr.set_label('('+SCM_BL[var].units+')')
+                    if file_rt is not None:
                         # SCM RTs
                         plt.subplot(3,1,2)
                         plt.contourf(x2, y2, z2, 20, cmap='YlGnBu')
                         plt.ylim(1000,200)
+                        plt.xlim(0,np.max(x1))
                         plt.ylabel('(Pa)')
                         plt.xlabel('(hours)')
                         cbr = plt.colorbar()
                         cbr.set_label('('+SCM_RT[var].units+')')
-                        # Only plot differences if they exist.
-                        if np.count_nonzero(dz) > 0:
+                    # end if
+                    # Only plot differences if requested, and only if they are non-zero.
+                    if plot_diff:
+                        dz = z1-z2
+                        if (np.count_nonzero(dz) > 0):
                             plt.subplot(3,1,3)
                             plt.title("Difference (top - middle)", fontsize=8)
                             plt.contourf(x2, y2, dz, 20, cmap='bwr')
@@ -165,9 +186,11 @@ def plot_results(file_BL, file_RT, plot_RT_only):
                             cbr = plt.colorbar()
                             cbr.set_label('('+SCM_RT[var].units+')')
                         # end if (no differences exist)
+                    # end if     (plot differences)
+                    # Save figure
                     fileOUT = 'scm.' + var +'.png'
+                    plt.show()
                     plt.savefig(fileOUT)
-                    #
                     plot_files.append(fileOUT)
                 # end if (Have enought pts to plot?)
             # end if     (fields exist?)