From 5fb29d4d2f7ae990883eb3abe735c3e914d902bc Mon Sep 17 00:00:00 2001 From: Dustin Swales Date: Tue, 21 May 2024 15:44:36 +0000 Subject: [PATCH] Added plotting step to SCM RT Baseline comparison in CI --- test/cmp_rt2bl.py | 177 ++++++++++++++++++++++++++++++++++++++++++ test/rt_test_cases.py | 2 +- 2 files changed, 178 insertions(+), 1 deletion(-) diff --git a/test/cmp_rt2bl.py b/test/cmp_rt2bl.py index 683e511db..02eedaeeb 100755 --- a/test/cmp_rt2bl.py +++ b/test/cmp_rt2bl.py @@ -10,6 +10,9 @@ from rt_test_cases import run_list from os.path import exists import argparse +from netCDF4 import Dataset +import numpy as np +import matplotlib.pyplot as plt # parser = argparse.ArgumentParser() @@ -17,6 +20,7 @@ parser.add_argument('-drt', '--dir_rt', help='Directory containing SCM RT output') parser.add_argument('-dbl', '--dir_bl', help='Directory containing SCM RT baselines') +# def parse_args(): args = parser.parse_args() build_type = args.build_type @@ -25,6 +29,155 @@ def parse_args(): return (build_type,dir_rt,dir_bl) +# +def plot_results(file_BL,file_RT): + # List of SCM output fields to plot (This 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"] + + # Open SCM datasets + SCM_BL = Dataset(file_BL) + SCM_RT = Dataset(file_RT) + + plot_files = [] + for var in SCM_BL.variables.keys(): + if (var in vars2plot): + # Handle temporal axis. + # There are 4 different dimensions in the SCM output, identified by the suffix "_dim". + # Here the suffix is stripped and used to identify the temporal dimenesion (index 0 in netcdf file) + 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 + # Is this a 2D (time, x) variable? (Really 1D since x=1 in SCM) + is2D = False + if (len(SCM_BL[var].dimensions)==2): + is2D = True + # endif + # one/two-dimensional variables + if (len(SCM_BL[var].shape) != 3): + print("Creating 1D field",var) + if (is2D): + y1 = SCM_BL[var][:,0].squeeze() + y2 = SCM_RT[var][:,0].squeeze() + else: + y1 = SCM_BL[var][:] + y2 = SCM_RT[var][:] + # endif + plt_range = [np.min(SCM_BL[var]),np.max(SCM_BL[var])] + plt_ranged = [-1*abs(np.max(SCM_BL[var])),abs(np.max(SCM_BL[var]))] + + # Make figure + fig = plt.figure(figsize=(13,10)) + # Baselines and SCM RTs on same plot + plt.subplot(2,1,1) + plt.title(SCM_BL[var].description) + plt.plot(x1, y1, color='blue') + plt.plot(x2, y2, color='black') + plt.ylim(plt_range) + plt.ylabel('('+SCM_BL[var].units+')') + plt.xlabel('(hours)') + # Difference (Baseline-SCMRT) + plt.subplot(2,1,2) + plt.title("Difference (blue - black)") + plt.plot(x1, y1 - y2, color='red') + plt.plot(x1, np.zeros(len(x1)), color='grey',linestyle='dashed') + plt.ylim(plt_ranged) + plt.ylabel('('+SCM_RT[var].units+')') + plt.xlabel('(hours)') + # + 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() + + # vertical axis + y1 = SCM_BL["pres"][0,:].squeeze()*0.01 + y2 = SCM_RT["pres"][0,:].squeeze()*0.01 + 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 + # endif + + # Comppute differences and determine valid plot range(s). + dz = z1-z2 + if np.count_nonzero(dz) > 0: + clev = np.arange(np.min(z1),np.max(z1),(np.max(z1)-np.min(z1))*0.05) + clevd = np.arange(np.min(dz),np.max(dz),(np.max(dz)-np.min(dz))*0.05) + else: + clev = 0 + clevd = 0 + # end if + + # Finally, make figure. + fig = plt.figure(figsize=(13,10)) + # Baselines + plt.subplot(3,1,1) + plt.title(SCM_BL[var].description, fontsize=12) + plt.contourf(x1, y1, z1, clev, cmap='YlGnBu') + plt.ylim(1000,200) + plt.ylabel('(Pa)') + cbr = plt.colorbar() + cbr.set_label('('+SCM_RT[var].units+')') + # SCM RTs + plt.subplot(3,1,2) + plt.contourf(x2, y2, z2, clev, cmap='YlGnBu') + plt.ylim(1000,200) + 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: + plt.subplot(3,1,3) + plt.title("Difference (top - middle)", fontsize=8) + plt.contourf(x2, y2, dz, clevd, cmap='bwr') + plt.ylim(1000,200) + plt.ylabel('(Pa)') + plt.xlabel('(hours)') + cbr = plt.colorbar() + cbr.set_label('('+SCM_RT[var].units+')') + # end if (no differences exist) + fileOUT = 'scm.' + var +'.png' + plt.savefig(fileOUT) + # + plot_files.append(fileOUT) + # end if (fields exist?) + # end if (field requested?) + # end for (fields in file) + + return(plot_files) + # def main(): # @@ -43,12 +196,35 @@ def main(): error_count = error_count + 1 else: print("Output for "+run["case"]+"_"+run["suite"]+ " is IDENTICAL to baseline") + # end if + + # Create plots between RTs and baselines + plot_files = plot_results(file_rt,file_bl) + + # Setup output directories for plots. + result = os.system("mkdir -p scm_rt_out/"+run["case"]+"/"+run["suite"]) + + # Archive plots. + com = "mv" + for plot_file in plot_files: + com = com + " " + plot_file + # end if + com = com + " scm_rt_out/" + run["case"] + "/" + run["suite"] + result = os.system(com) else: if not exists(file_rt): print("Output for "+run["case"]+"_"+run["suite"]+ " is MISSING from output") + # end if if not exists(file_bl): print("Output for "+run["case"]+"_"+run["suite"]+ " is MISSING from baseline") + # end if error_count = error_count + 1 + # end if + # end for + + # Bundle all plots into tarball for artifact in github action. + com = 'tar -cvf scm_rt_plots.tar scm_rt_out/*' + result = os.system(com) # if error_count == 0: @@ -56,6 +232,7 @@ def main(): else: print("ALL TESTS PASSED, BUT OUTPUT DIFFERS FROM BASELINE.") #1/0 + # end if # if __name__ == '__main__': diff --git a/test/rt_test_cases.py b/test/rt_test_cases.py index 7c52f88b5..0a2bd76d0 100644 --- a/test/rt_test_cases.py +++ b/test/rt_test_cases.py @@ -21,7 +21,7 @@ {"case": "astex", "suite": "SCM_GFS_v17_HR3_RRTMGP"}, \ {"case": "astex", "suite": "SCM_GFS_v16"}, \ {"case": "astex", "suite": "SCM_WoFS_v0"}, \ - {"case": "astex", "suite": "SCM_RRFS_v1"}, \ + {"case": "astex", "suite": "SCM_RRFS_v1"}, \ {"case": "LASSO_2016051812", "suite": "SCM_GFS_v17_HR3"}, \ {"case": "LASSO_2016051812", "suite": "SCM_GFS_v17_HR3_RRTMGP"}, \ {"case": "LASSO_2016051812", "suite": "SCM_GFS_v16"}, \