''' Name: plot_lead_by_level.py Contact(s): Mallory Row Abstract: Reads average files from plot_time_series.py to make lead-pressue plots History Log: Third version Usage: Called by make_plots_wrapper.py Parameters: None Input Files: Text files Output Files: .png images Condition codes: 0 for success, 1 for failure ''' import os import numpy as np import pandas as pd import itertools import warnings import logging import datetime import re import sys import matplotlib matplotlib.use('agg') import matplotlib.pyplot as plt import matplotlib.dates as md import matplotlib.gridspec as gridspec import plot_util as plot_util from plot_util import get_lead_avg_file # add metplus directory to path so the wrappers and utilities can be found sys.path.insert(0, os.path.abspath(os.environ['HOMEMETplus'])) from metplus.util import do_string_sub #### EMC-verif_global title creation script import plot_title as plot_title # Plot Settings plt.rcParams['font.weight'] = 'bold' plt.rcParams['axes.titleweight'] = 'bold' plt.rcParams['axes.titlesize'] = 16 plt.rcParams['axes.titlepad'] = 5 plt.rcParams['axes.labelweight'] = 'bold' plt.rcParams['axes.labelsize'] = 14 plt.rcParams['axes.labelpad'] = 10 plt.rcParams['axes.formatter.useoffset'] = False plt.rcParams['xtick.labelsize'] = 14 plt.rcParams['xtick.major.pad'] = 5 plt.rcParams['ytick.major.pad'] = 5 plt.rcParams['ytick.labelsize'] = 14 plt.rcParams['figure.subplot.left'] = 0.1 plt.rcParams['figure.subplot.right'] = 0.95 plt.rcParams['figure.titleweight'] = 'bold' plt.rcParams['figure.titlesize'] = 16 if float(matplotlib.__version__[0:3]) >= 3.3: plt.rcParams['date.epoch'] = '0000-12-31T00:00:00' title_loc = 'center' cmap_bias_original = plt.cm.PiYG_r colors_bias = cmap_bias_original( np.append(np.linspace(0,0.3,10), np.linspace(0.7,1,10)) ) cmap_bias = matplotlib.colors.LinearSegmentedColormap.from_list( 'cmap_bias', colors_bias ) cmap = plt.cm.BuPu noaa_logo_img_array = matplotlib.image.imread( os.path.join(os.environ['USHverif_global'], 'plotting_scripts', 'noaa.png') ) noaa_logo_alpha = 0.5 nws_logo_img_array = matplotlib.image.imread( os.path.join(os.environ['USHverif_global'], 'plotting_scripts', 'nws.png') ) nws_logo_alpha = 0.5 # Read environment variables set in make_plots_wrapper.py verif_case = os.environ['VERIF_CASE'] verif_type = os.environ['VERIF_TYPE'] date_type = os.environ['DATE_TYPE'] valid_beg = os.environ['VALID_BEG'] valid_end = os.environ['VALID_END'] init_beg = os.environ['INIT_BEG'] init_end = os.environ['INIT_END'] fcst_valid_hour = os.environ['FCST_VALID_HOUR'] fcst_init_hour = os.environ['FCST_INIT_HOUR'] obs_valid_hour = os.environ['OBS_VALID_HOUR'] obs_init_hour = os.environ['OBS_INIT_HOUR'] fcst_lead_list = [os.environ['FCST_LEAD'].split(', ')] fcst_var_name = os.environ['FCST_VAR'] fcst_var_units = os.environ['FCST_UNITS'] fcst_var_level_list = [os.environ['FCST_LEVEL'].split(', ')] fcst_var_thresh_list = os.environ['FCST_THRESH'].split(', ') obs_var_name = os.environ['OBS_VAR'] obs_var_units = os.environ['OBS_UNITS'] obs_var_level_list = [os.environ['OBS_LEVEL'].split(', ')] obs_var_thresh_list = os.environ['OBS_THRESH'].split(', ') interp_pnts = os.environ['INTERP_PNTS'] vx_mask = os.environ['VX_MASK'] alpha = os.environ['ALPHA'] desc = os.environ['DESC'] obs_lead = os.environ['OBS_LEAD'] cov_thresh = os.environ['COV_THRESH'] stats_list = os.environ['STATS'].split(', ') model_list = os.environ['MODEL'].split(', ') model_obtype_list = os.environ['MODEL_OBTYPE'].split(', ') model_reference_name_list = os.environ['MODEL_REFERENCE_NAME'].split(', ') dump_row_filename_template = os.environ['DUMP_ROW_FILENAME'] average_method = os.environ['AVERAGE_METHOD'] ci_method = os.environ['CI_METHOD'] verif_grid = os.environ['VERIF_GRID'] event_equalization = os.environ['EVENT_EQUALIZATION'] met_version = os.environ['MET_VERSION'] input_base_dir = os.environ['INPUT_BASE_DIR'] output_base_dir = os.environ['OUTPUT_BASE_DIR'] log_metplus = os.environ['LOG_METPLUS'] log_level = os.environ['LOG_LEVEL'] if len(fcst_lead_list[0]) < 2: logger.warning("Must provide more than one forecast lead to " "plot lead average") sys.exit(0) #### EMC-verif_global environment variables var_name = os.environ['var_name'] fcst_var_extra = (os.environ['fcst_var_options'].replace(' ', '') \ .replace('=','').replace(';','').replace('"','') \ .replace("'",'').replace(',','-').replace('_','')) obs_var_extra = (os.environ['obs_var_options'].replace(' ', '') \ .replace('=','').replace(';','').replace('"','') \ .replace("'",'').replace(',','-').replace('_','')) interp_mthd = os.environ['interp'] if date_type == 'VALID': start_date = valid_beg end_date = valid_end elif date_type == 'INIT': start_date = init_beg end_date = init_end img_quality = os.environ['img_quality'] # General set up and settings # Image Quality if img_quality == 'low': plt.rcParams['savefig.dpi'] = 50 elif img_quality == 'medium': plt.rcParams['savefig.dpi'] = 75 # Logging logger = logging.getLogger(log_metplus) logger.setLevel(log_level) formatter = logging.Formatter( '%(asctime)s.%(msecs)03d (%(filename)s:%(lineno)d) %(levelname)s: ' +'%(message)s', '%m/%d %H:%M:%S' ) file_handler = logging.FileHandler(log_metplus, mode='a') file_handler.setFormatter(formatter) logger.addHandler(file_handler) # Output output_data_dir = os.path.join(output_base_dir, 'data') #output_imgs_dir = os.path.join(output_base_dir, 'imgs') #### EMC-verif_global image directory output_imgs_dir = os.path.join(output_base_dir, 'images') # Model info model_info_list = list( zip(model_list, model_reference_name_list, model_obtype_list, ) ) nmodels = len(model_info_list) # Plot info if fcst_var_name == 'O3MR': plot_info_list = list( itertools.product(*[fcst_lead_list, ['all'], fcst_var_thresh_list]) ) else: plot_info_list = list( itertools.product(*[fcst_lead_list, ['all', 'trop', 'lowertrop', 'uppertrop', 'strat'], fcst_var_thresh_list]) ) # Level info fcst_var_level_all_list = [] fcst_var_level_trop_list = [] fcst_var_level_lower_trop_list = [] fcst_var_level_upper_trop_list = [] fcst_var_level_strat_list = [] for fcst_var_level in fcst_var_level_list[0]: fcst_var_level_all_list.append(fcst_var_level) if int(fcst_var_level[1:]) > 100: fcst_var_level_trop_list.append(fcst_var_level) elif int(fcst_var_level[1:]) < 100: fcst_var_level_strat_list.append(fcst_var_level) elif int(fcst_var_level[1:]) == 100: fcst_var_level_trop_list.append(fcst_var_level) fcst_var_level_strat_list.append(fcst_var_level) if int(fcst_var_level[1:]) >= 500: fcst_var_level_lower_trop_list.append(fcst_var_level) if int(fcst_var_level[1:]) <= 500 and int(fcst_var_level[1:]) >= 100: fcst_var_level_upper_trop_list.append(fcst_var_level) # Date and time infomation and build title for plot date_beg = os.environ[date_type+'_BEG'] date_end = os.environ[date_type+'_END'] valid_init_dict = { 'fcst_valid_hour_beg': fcst_valid_hour.split(', ')[0], 'fcst_valid_hour_end': fcst_valid_hour.split(', ')[-1], 'fcst_init_hour_beg': fcst_init_hour.split(', ')[0], 'fcst_init_hour_end': fcst_init_hour.split(', ')[-1], 'obs_valid_hour_beg': obs_valid_hour.split(', ')[0], 'obs_valid_hour_end': obs_valid_hour.split(', ')[-1], 'obs_init_hour_beg': obs_init_hour.split(', ')[0], 'obs_init_hour_end': obs_init_hour.split(', ')[-1], 'valid_hour_beg': '', 'valid_hour_end': '', 'init_hour_beg': '', 'init_hour_end': '' } valid_init_type_list = [ 'valid_hour_beg', 'valid_hour_end', 'init_hour_beg', 'init_hour_end' ] for vitype in valid_init_type_list: if (valid_init_dict['fcst_'+vitype] != '' and valid_init_dict['obs_'+vitype] == ''): valid_init_dict[vitype] = valid_init_dict['fcst_'+vitype] elif (valid_init_dict['obs_'+vitype] != '' and valid_init_dict['fcst_'+vitype] == ''): valid_init_dict[vitype] = valid_init_dict['obs_'+vitype] if valid_init_dict['fcst_'+vitype] == '': if 'beg' in vitype: valid_init_dict['fcst_'+vitype] = '000000' elif 'end' in vitype: valid_init_dict['fcst_'+vitype] = '235959' if valid_init_dict['obs_'+vitype] == '': if 'beg' in vitype: valid_init_dict['obs_'+vitype] = '000000' elif 'end' in vitype: valid_init_dict['obs_'+vitype] = '235959' if valid_init_dict['fcst_'+vitype] == valid_init_dict['obs_'+vitype]: valid_init_dict[vitype] = valid_init_dict['fcst_'+vitype] # Start looping to make plots for plot_info in plot_info_list: fcst_leads = plot_info[0] fcst_lead_timedeltas = np.full_like(fcst_leads, np.nan, dtype=float) for fcst_lead in fcst_leads: fcst_lead_idx = fcst_leads.index(fcst_lead) fcst_lead_timedelta = datetime.timedelta( hours=int(fcst_lead[:-4]), minutes=int(fcst_lead[-4:-2]), seconds=int(fcst_lead[-2:]) ).total_seconds() fcst_lead_timedeltas[fcst_lead_idx] = float(fcst_lead_timedelta) fcst_lead_timedeltas_str = [] for tdelta in fcst_lead_timedeltas: h = int(tdelta/3600) m = int((tdelta-(h*3600))/60) s = int(tdelta-(h*3600)-(m*60)) if h < 10: tdelta_str = f"{h:01d}" elif h < 100: tdelta_str = f"{h:02d}" else: tdelta_str = f"{h:03d}" if m != 0: tdelta_str+=f":{m:02d}" if s != 0: tdelta_str+=f":{s:02d}" fcst_lead_timedeltas_str.append(tdelta_str) if plot_info[1] == 'all': fcst_var_levels = fcst_var_level_all_list elif plot_info[1] == 'trop': fcst_var_levels = fcst_var_level_trop_list elif plot_info[1] == 'lowertrop': fcst_var_levels = fcst_var_level_lower_trop_list elif plot_info[1] == 'uppertrop': fcst_var_levels = fcst_var_level_upper_trop_list elif plot_info[1] == 'strat': fcst_var_levels = fcst_var_level_strat_list logger.info("Working on levels: "+plot_info[1]+" " +' '.join(fcst_var_levels)) obs_var_levels = fcst_var_levels fcst_var_thresh = plot_info[2] obs_var_thresh = obs_var_thresh_list[ fcst_var_thresh_list.index(fcst_var_thresh) ] fcst_var_thresh_symbol, fcst_var_thresh_letter = plot_util.format_thresh( fcst_var_thresh ) obs_var_thresh_symbol, obs_var_thresh_letter = plot_util.format_thresh( obs_var_thresh ) logger.info("Working on forecast lead averages" +" for forecast variable "+fcst_var_name +" "+fcst_var_thresh) # Set up base name for file naming convention for lead average files, # and output data and images base_name = date_type.lower()+date_beg+'to'+date_end if (valid_init_dict['valid_hour_beg'] != '' and valid_init_dict['valid_hour_end'] != ''): base_name+=( '_valid'+valid_init_dict['valid_hour_beg'][0:4] +'to'+valid_init_dict['valid_hour_end'][0:4]+'Z' ) else: base_name+=( '_fcst_valid'+valid_init_dict['fcst_valid_hour_beg'][0:4] +'to'+valid_init_dict['fcst_valid_hour_end'][0:4]+'Z' +'_obs_valid'+valid_init_dict['obs_valid_hour_beg'][0:4] +'to'+valid_init_dict['obs_valid_hour_end'][0:4]+'Z' ) if (valid_init_dict['init_hour_beg'] != '' and valid_init_dict['init_hour_end'] != ''): base_name+=( '_init'+valid_init_dict['init_hour_beg'][0:4] +'to'+valid_init_dict['init_hour_end'][0:4]+'Z' ) else: base_name+=( '_fcst_init'+valid_init_dict['fcst_init_hour_beg'][0:4] +'to'+valid_init_dict['fcst_init_hour_end'][0:4]+'Z' +'_obs_init'+valid_init_dict['obs_init_hour_beg'][0:4] +'to'+valid_init_dict['obs_init_hour_end']+'Z' ) base_name+=( '_fcst_lead_avgs' +'_fcst'+fcst_var_name+'FCSTLEVELHOLDER' +fcst_var_thresh_letter.replace(',', '_')+interp_mthd +'_obs'+obs_var_name+'OBSLEVELHOLDER' +obs_var_thresh_letter.replace(',', '_')+interp_mthd +'_vxmask'+vx_mask ) if desc != '': base_name+='_desc'+desc if obs_lead != '': base_name+='_obs_lead'+obs_lead if interp_pnts != '': base_name+='_interp_pnts'+interp_pnts if cov_thresh != '': cov_thresh_symbol, cov_thresh_letter = plot_util.format_thresh( cov_thresh ) base_name+='_cov_thresh'+cov_thresh_letter.replace(',', '_') if alpha != '': base_name+='_alpha'+alpha # Build date by forecst level grid for plotting fcst_var_levels_int = np.empty(len(fcst_var_levels), dtype=int) for vl in range(len(fcst_var_levels)): fcst_var_levels_int[vl] = fcst_var_levels[vl][1:] xmesh, ymesh = np.meshgrid(fcst_lead_timedeltas, fcst_var_levels_int) for stat in stats_list: logger.debug("Working on "+stat) if stat in ['bias', 'rmse', 'rmse_md', 'rmse_pv', 'fbar_obar']: cmap_diff_original = plt.cm.bwr else: cmap_diff_original = plt.cm.bwr_r colors_diff = cmap_diff_original( np.append(np.linspace(0,0.425,10), np.linspace(0.575,1,10)) ) cmap_diff = matplotlib.colors.LinearSegmentedColormap.from_list( 'cmap_diff', colors_diff ) stat_plot_name = plot_util.get_stat_plot_name(logger, stat) if (stat == 'fbar_obar' or stat == 'orate_frate' or stat == 'baser_frate'): nsubplots = nmodels + 1 else: nsubplots = nmodels if nsubplots == 1: x_figsize, y_figsize = 14, 7 row, col = 1, 1 hspace, wspace = 0, 0 bottom, top = 0.175, 0.825 suptitle_y_loc = 0.92125 noaa_logo_x_scale, noaa_logo_y_scale = 0.1, 0.865 nws_logo_x_scale, nws_logo_y_scale = 0.9, 0.865 cbar_bottom = 0.06 cbar_height = 0.02 elif nsubplots == 2: x_figsize, y_figsize = 14, 7 row, col = 1, 2 hspace, wspace = 0, 0.1 bottom, top = 0.175, 0.825 suptitle_y_loc = 0.92125 noaa_logo_x_scale, noaa_logo_y_scale = 0.1, 0.865 nws_logo_x_scale, nws_logo_y_scale = 0.9, 0.865 cbar_bottom = 0.06 cbar_height = 0.02 elif nsubplots > 2 and nsubplots <= 4: x_figsize, y_figsize = 14, 14 row, col = 2, 2 hspace, wspace = 0.15, 0.1 bottom, top = 0.125, 0.9 suptitle_y_loc = 0.9605 noaa_logo_x_scale, noaa_logo_y_scale = 0.1, 0.9325 nws_logo_x_scale, nws_logo_y_scale = 0.9, 0.9325 cbar_bottom = 0.03 cbar_height = 0.02 elif nsubplots > 4 and nsubplots <= 6: x_figsize, y_figsize = 14, 14 row, col = 3, 2 hspace, wspace = 0.15, 0.1 bottom, top = 0.125, 0.9 suptitle_y_loc = 0.9605 noaa_logo_x_scale, noaa_logo_y_scale = 0.1, 0.9325 nws_logo_x_scale, nws_logo_y_scale = 0.9, 0.9325 cbar_bottom = 0.03 cbar_height = 0.02 elif nsubplots > 6 and nsubplots <= 8: x_figsize, y_figsize = 14, 14 row, col = 4, 2 hspace, wspace = 0.175, 0.1 bottom, top = 0.125, 0.9 suptitle_y_loc = 0.9605 noaa_logo_x_scale, noaa_logo_y_scale = 0.1, 0.9325 nws_logo_x_scale, nws_logo_y_scale = 0.9, 0.9325 cbar_bottom = 0.03 cbar_height = 0.02 elif nsubplots > 8 and nsubplots <= 10: x_figsize, y_figsize = 14, 14 row, col = 5, 2 hspace, wspace = 0.225, 0.1 bottom, top = 0.125, 0.9 suptitle_y_loc = 0.9605 noaa_logo_x_scale, noaa_logo_y_scale = 0.1, 0.9325 nws_logo_x_scale, nws_logo_y_scale = 0.9, 0.9325 cbar_bottom = 0.03 cbar_height = 0.02 else: logger.error("Too many subplots selected, max. is 10") exit(1) suptitle_x_loc = (plt.rcParams['figure.subplot.left'] +plt.rcParams['figure.subplot.right'])/2. fig = plt.figure(figsize=(x_figsize, y_figsize)) gs = gridspec.GridSpec( row, col, bottom = bottom, top = top, hspace = hspace, wspace = wspace, ) noaa_logo_xpixel_loc = ( x_figsize * plt.rcParams['figure.dpi'] * noaa_logo_x_scale ) noaa_logo_ypixel_loc = ( y_figsize * plt.rcParams['figure.dpi'] * noaa_logo_y_scale ) nws_logo_xpixel_loc = ( x_figsize * plt.rcParams['figure.dpi'] * nws_logo_x_scale ) nws_logo_ypixel_loc = ( y_figsize * plt.rcParams['figure.dpi'] * nws_logo_y_scale ) #### EMC_verif-global set up subplots subplot_num = 0 while subplot_num < nsubplots: ax = plt.subplot(gs[subplot_num]) ax.grid(True) if len(fcst_lead_timedeltas) >= 25: ax.set_xticks(fcst_lead_timedeltas[::4]) ax.set_xticklabels(fcst_lead_timedeltas_str[::4]) elif len(fcst_lead_timedeltas) >= 15: ax.set_xticks(fcst_lead_timedeltas[::2]) ax.set_xticklabels(fcst_lead_timedeltas_str[::2]) else: ax.set_xticks(fcst_lead_timedeltas) ax.set_xticklabels(fcst_lead_timedeltas_str) ax.set_xlim([fcst_lead_timedeltas[0], fcst_lead_timedeltas[-1]]) if ax.is_last_row() \ or (nsubplots % 2 != 0 and subplot_num == nsubplots -1): ax.set_xlabel('Forecast Hour') else: plt.setp(ax.get_xticklabels(), visible=False) ax.set_yscale('log') ax.invert_yaxis() ax.minorticks_off() ax.set_yticks(fcst_var_levels_int) ax.set_yticklabels(fcst_var_levels_int) ax.set_ylim([fcst_var_levels_int[0],fcst_var_levels_int[-1]]) if ax.is_first_col() \ or (nsubplots % 2 != 0 and subplot_num == nsubplots -1): ax.set_ylabel('Pressure Level (hPa)') else: plt.setp(ax.get_yticklabels(), visible=False) subplot_num+=1 if (stat == 'fbar_obar' or stat == 'orate_frate' or stat == 'baser_frate'): avg_file_cols = ['LEADS', 'VALS', 'OBS_VALS'] else: avg_file_cols = ['LEADS', 'VALS'] avg_cols_to_array = avg_file_cols[1:] # Reading in model lead average files produced from # plot_time_series.py obs_plotted = False get_clevels = True make_colorbar = False logger.info("Reading in model data") for model_info in model_info_list: model_num = model_info_list.index(model_info) + 1 model_idx = model_info_list.index(model_info) model_name = model_info[0] model_plot_name = model_info[1] model_obtype = model_info[2] model_avg_data = np.empty( [len(avg_cols_to_array), len(fcst_var_levels), len(fcst_leads)] ) model_avg_data.fill(np.nan) for vl in range(len(fcst_var_levels)): fcst_var_level = fcst_var_levels[vl] obs_var_level = obs_var_levels[vl] model_stat_template = dump_row_filename_template string_sub_dict = { 'model': model_name, 'model_reference': model_plot_name, 'obtype': model_obtype, 'fcst_lead': fcst_lead, 'fcst_level': fcst_var_level, 'obs_level': obs_var_level, 'fcst_thresh': fcst_var_thresh, 'obs_thresh': obs_var_thresh, } model_stat_file = do_string_sub(model_stat_template, **string_sub_dict) lead_avg_file = get_lead_avg_file(stat, model_stat_file, fcst_lead, output_base_dir) if os.path.exists(lead_avg_file): nrow = sum(1 for line in open(lead_avg_file)) if nrow == 0: logger.error("Model "+str(model_num)+" " +model_name+" with plot name " +model_plot_name+" file: " +lead_avg_file+" empty") sys.exit(1) else: logger.debug("Model "+str(model_num)+" " +model_name+" with plot name " +model_plot_name+" file: " +lead_avg_file+" exists") model_avg_file_data = pd.read_csv( lead_avg_file, sep=' ', header=None, names=avg_file_cols, dtype=str ) model_avg_file_data_leads = ( model_avg_file_data.loc[:]['LEADS'].tolist() ) for fcst_lead in fcst_leads: fcst_lead_idx = fcst_leads.index(fcst_lead) if fcst_lead in model_avg_file_data_leads: model_fcst_lead_idx = ( model_avg_file_data_leads.index( fcst_lead ) ) for col in avg_cols_to_array: col_idx = avg_cols_to_array.index(col) model_avg_file_data_col = ( model_avg_file_data.loc[:][col].tolist() ) if (model_avg_file_data_col\ [model_fcst_lead_idx] != '--'): model_avg_data[col_idx, vl, fcst_lead_idx] = ( float(model_avg_file_data_col \ [model_fcst_lead_idx]) ) else: logger.error("Model "+str(model_num)+" " +model_name+" with plot name " +model_plot_name+" file: " +lead_avg_file+" does not exist") sys.exit(1) model_avg_data = np.ma.masked_invalid(model_avg_data) if (stat == 'fbar_obar' or stat == 'orate_frate' or stat == 'baser_frate'): if not obs_plotted: ax = plt.subplot(gs[0]) ax.set_title('obs', loc='left') if not model_avg_data[1,:,:].mask.all(): logger.debug("Plotting observations") ax = plt.subplot(gs[0]) obs_avg_data = model_avg_data[1,:,:] ax.set_title('obs', loc='left') CF1 = ax.contourf(xmesh, ymesh, obs_avg_data, cmap=cmap, extend='both') C1 = ax.contour(xmesh, ymesh, obs_avg_data, levels=CF1.levels, colors='k', linewidths=1.0) C1_labels_list = [] for level in C1.levels: if str(level).split('.')[1] == '0': C1_labels_list.append( str(int(level)) ) else: C1_labels_list.append( str(round(level,3)).rstrip('0') ) fmt = {} for lev, label in zip(C1.levels, C1_labels_list): fmt[lev] = label ax.clabel(C1, C1.levels, fmt=fmt, inline=True, fontsize=12.5) obs_plotted = True ax = plt.subplot(gs[model_num]) ax.set_title(model_plot_name+'-obs', loc='left') model_obs_diff = (model_avg_data[0,:,:] - model_avg_data[1,:,:]) if not model_obs_diff.mask.all(): logger.debug("Plotting model "+str(model_num)+" " +model_name+" - obs " +"with name on plot "+model_plot_name+" " +"- obs") if get_clevels: clevels_diff = plot_util.get_clevels(model_obs_diff, 1) CF2 = ax.contourf(xmesh, ymesh, model_obs_diff, levels=clevels_diff, cmap=cmap_diff, extend='both') get_clevels = False make_colorbar = True colorbar_CF = CF2 colorbar_CF_ticks = CF2.levels colorbar_label = 'Difference' else: CF = ax.contourf(xmesh, ymesh, model_obs_diff, levels=CF2.levels, cmap=cmap_diff, extend='both') elif stat == 'bias' or stat == 'fbias': ax = plt.subplot(gs[model_idx]) ax.set_title(model_plot_name, loc='left') if not model_avg_data[0,:,:].mask.all(): logger.debug("Plotting model "+str(model_num) +" "+model_name+" with name on plot " +model_plot_name) if get_clevels: if np.max(np.abs(model_avg_data[0,:,:])) > 100: bias_spacing = 2.25 elif np.max(np.abs(model_avg_data[0,:,:])) > 10: bias_spacing = 2 else: bias_spacing = 1.75 clevels_bias = plot_util.get_clevels( model_avg_data[0,:,:], bias_spacing ) CF1 = ax.contourf(xmesh, ymesh, model_avg_data[0,:,:], levels=clevels_bias, cmap=cmap_bias, extend='both') C1 = ax.contour(xmesh, ymesh, model_avg_data[0,:,:], levels=CF1.levels, colors='k', linewidths=1.0) C1_labels_list = [] for level in C1.levels: if str(level).split('.')[1] == '0': C1_labels_list.append( str(int(level)) ) else: C1_labels_list.append( str(round(level,3)).rstrip('0') ) fmt = {} for lev, label in zip(C1.levels, C1_labels_list): fmt[lev] = label ax.clabel(C1, C1.levels, fmt=fmt, inline=True, fontsize=12.5) get_clevels = False make_colorbar = True colorbar_CF = CF1 colorbar_CF_ticks = CF1.levels colorbar_label = 'Bias' else: CF = ax.contourf(xmesh, ymesh, model_avg_data[0,:,:], levels=CF1.levels, cmap=cmap_bias, extend='both') C = ax.contour(xmesh, ymesh, model_avg_data[0,:,:], levels=CF1.levels, colors='k', linewidths=1.0) C_labels_list = [] for level in C.levels: if str(level).split('.')[1] == '0': C_labels_list.append( str(int(level)) ) else: C_labels_list.append( str(round(level,3)).rstrip('0') ) fmt = {} for lev, label in zip(C.levels, C_labels_list): fmt[lev] = label ax.clabel(C, C.levels, fmt=fmt, inline=True, fontsize=12.5) else: ax = plt.subplot(gs[model_idx]) if model_num == 1: model1_name = model_name model1_plot_name = model_plot_name model1_avg_data = model_avg_data[0,:,:] ax.set_title(model_plot_name, loc='left') if not model1_avg_data.mask.all(): logger.debug("Plotting model "+str(model_num)+" " +model_name+" with name on plot " +model_plot_name) if stat in ['acc']: levels = np.array( [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 1] ) CF1 = ax.contourf(xmesh, ymesh, model_avg_data[0,:,:], levels=levels, cmap=cmap, extend='both') elif stat in ['rmse', 'rmse_md', 'rmse_pv']: cmax = np.nanmax(model_avg_data[0,:,:]) steps = 12 dx = 1.0 / (steps-1) if cmax > 100: spacing = 2.25 elif cmax > 10: spacing = 2 else: spacing = 1.75 levels = np.array( [0+(i*dx)**spacing*cmax for i in range(steps)], dtype=float ) CF1 = ax.contourf(xmesh, ymesh, model_avg_data[0,:,:], levels=levels, cmap=cmap, extend='both') else: CF1 = ax.contourf(xmesh, ymesh, model_avg_data[0,:,:], cmap=cmap, extend='both') C1 = ax.contour(xmesh, ymesh, model_avg_data[0,:,:], levels=CF1.levels, colors='k', linewidths=1.0) C1_labels_list = [] for level in C1.levels: if str(level).split('.')[1] == '0': C1_labels_list.append( str(int(level)) ) else: C1_labels_list.append( str(round(level,3)).rstrip('0') ) fmt = {} for lev, label in zip(C1.levels, C1_labels_list): fmt[lev] = label ax.clabel(C1, C1.levels, fmt=fmt, inline=True, fontsize=12.5) else: ax.set_title(model_plot_name+'-'+model1_plot_name, loc='left') model_model1_diff = (model_avg_data[0,:,:] - model1_avg_data) if not model_model1_diff.mask.all(): logger.debug("Plotting model "+str(model_num)+" " +model_name+" - model 1 "+model1_name+" " +"with name on plot "+model_plot_name+" " +"- "+model1_plot_name) if get_clevels: if stat in ['acc']: clevels_diff = np.array( [-0.5, -0.4, -0.3, -0.2, -0.1, -0.05, 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5] ) else: clevels_diff = plot_util.get_clevels( model_model1_diff, 1.25 ) CF2 = ax.contourf(xmesh, ymesh, model_model1_diff, levels=clevels_diff, cmap=cmap_diff, extend='both') get_clevels = False make_colorbar = True colorbar_CF = CF2 colorbar_CF_ticks = CF2.levels colorbar_label = 'Difference' else: CF = ax.contourf(xmesh, ymesh, model_model1_diff, levels=CF2.levels, cmap=cmap_diff, extend='both') #### EMC-verif_global build formal plot title if verif_grid == vx_mask: grid_vx_mask = verif_grid else: grid_vx_mask = verif_grid+vx_mask if verif_type in ['sfc', 'conus_sfc']: var_info_title = plot_title.get_var_info_title( fcst_var_name, 'all', fcst_var_extra, fcst_var_thresh ) else: var_info_title = plot_title.get_var_info_title( var_name, 'all', fcst_var_extra, fcst_var_thresh ) vx_mask_title = plot_title.get_vx_mask_title(vx_mask) date_info_title = plot_title.get_date_info_title( date_type, fcst_valid_hour.split(', '), fcst_init_hour.split(', '), datetime.datetime.strptime( start_date, '%Y%m%d' ).strftime('%d%b%Y'), datetime.datetime.strptime( end_date, '%Y%m%d' ).strftime('%d%b%Y'), verif_case ) full_title = ( stat_plot_name+'\n' +var_info_title+', '+vx_mask_title+'\n' +date_info_title ) fig.suptitle(full_title, x = suptitle_x_loc, y = suptitle_y_loc, horizontalalignment = title_loc, verticalalignment = title_loc) noaa_img = fig.figimage(noaa_logo_img_array, noaa_logo_xpixel_loc, noaa_logo_ypixel_loc, zorder=1, alpha=noaa_logo_alpha) nws_img = fig.figimage(nws_logo_img_array, nws_logo_xpixel_loc, nws_logo_ypixel_loc, zorder=1, alpha=nws_logo_alpha) if img_quality in ['low', 'medium']: noaa_img.set_visible(False) nws_img.set_visible(False) plt.subplots_adjust( left = noaa_img.get_extent()[1] \ /(plt.rcParams['figure.dpi']*x_figsize), right = nws_img.get_extent()[0] \ /(plt.rcParams['figure.dpi']*x_figsize) ) #### EMC-verif_global add colorbar cbar_left = ( noaa_img.get_extent()[1]/(plt.rcParams['figure.dpi']*x_figsize) ) cbar_width = ( nws_img.get_extent()[0]/(plt.rcParams['figure.dpi']*x_figsize) - noaa_img.get_extent()[1]/(plt.rcParams['figure.dpi']*x_figsize) ) if make_colorbar: cax = fig.add_axes( [cbar_left, cbar_bottom, cbar_width, cbar_height] ) cbar = fig.colorbar(colorbar_CF, cax = cax, orientation = 'horizontal', ticks = colorbar_CF_ticks) cbar.ax.set_xlabel(colorbar_label, labelpad = 0) cbar.ax.xaxis.set_tick_params(pad=0) cbar_tick_labels_list = [] for tick in cbar.get_ticks(): if str(tick).split('.')[1] == '0': cbar_tick_labels_list.append( str(int(tick)) ) else: cbar_tick_labels_list.append( str(round(tick,3)).rstrip('0') ) cbar.ax.set_xticklabels(cbar_tick_labels_list) #### EMC-verif_global build savefig name savefig_name = os.path.join(output_imgs_dir, stat) if date_type == 'VALID': if verif_case == 'grid2obs': savefig_name = ( savefig_name+'_init' +fcst_init_hour.split(', ')[0][0:2]+'Z' ) else: savefig_name = ( savefig_name+'_valid' +fcst_valid_hour.split(', ')[0][0:2]+'Z' ) elif date_type == 'INIT': if verif_case == 'grid2obs': savefig_name = ( savefig_name+'_valid' +fcst_valid_hour.split(', ')[0][0:2]+'Z' ) else: savefig_name = ( savefig_name+'_init' +fcst_init_hour.split(', ')[0][0:2]+'Z' ) if verif_case == 'grid2grid' and verif_type == 'anom': savefig_name = savefig_name+'_'+var_name+'_'+plot_info[1] else: savefig_name = savefig_name+'_'+fcst_var_name+'_'+plot_info[1] if verif_case == 'precip': savefig_name = savefig_name+'_'+fcst_var_thresh savefig_name = savefig_name+'_fhrmean_'+grid_vx_mask+'.png' logger.info("Saving image as "+savefig_name) plt.savefig(savefig_name) plt.clf() plt.close('all')