'''
Name: plot_date_by_level.py
Contact(s): Mallory Row
Abstract: Reads filtered files from stat_analysis_wrapper
run_all_times to make date-pressure plots
History Log: Third version
Usage: Called by make_plots_wrapper.py
Parameters: None
Input Files: MET .stat 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
# 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'
nticks = 2
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_mthd = os.environ['INTERP_MTHD']
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']
#### 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('_',''))
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]
# MET .stat file formatting
stat_file_base_columns = plot_util.get_stat_file_base_columns(met_version)
nbase_columns = len(stat_file_base_columns)
# Start looping to make plots
for plot_info in plot_info_list:
fcst_lead = plot_info[0]
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 "+fcst_lead+" "
+"and forecast variable "+fcst_var_name+" "
+fcst_var_thresh)
# Set up base name for file naming convention for MET .stat 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'+fcst_lead
+'_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
# Set up expected date in MET .stat file and date plot information
plot_time_dates, expected_stat_file_dates = plot_util.get_date_arrays(
date_type, date_beg, date_end,
fcst_valid_hour, fcst_init_hour,
obs_valid_hour, obs_init_hour,
fcst_lead
)
total_dates = len(plot_time_dates)
#### EMC-verif_global date tick interval
if len(plot_time_dates) < nticks:
date_tick_intvl = 1
else:
date_tick_intvl = int(len(plot_time_dates)/nticks)
# 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(plot_time_dates, fcst_var_levels_int)
# Reading in model .stat files from stat_analysis
logger.info("Reading in model data")
for model_info in model_info_list:
model_num = model_info_list.index(model_info) + 1
model_name = model_info[0]
model_plot_name = model_info[1]
model_obtype = model_info[2]
for vl in range(len(fcst_var_levels)):
fcst_var_level = fcst_var_levels[vl]
obs_var_level = obs_var_levels[vl]
model_level_now_data_index = pd.MultiIndex.from_product(
[
[model_plot_name], [fcst_var_level],
expected_stat_file_dates
],
names=['model_plot_name', 'levels', 'dates']
)
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)
if os.path.exists(model_stat_file):
nrow = sum(1 for line in open(model_stat_file))
if nrow == 0:
logger.warning("Model "+str(model_num)+" "+model_name+" "
+"with plot name "+model_plot_name+" "
+"file: "+model_stat_file+" empty")
model_level_now_data = (
pd.DataFrame(np.nan,
index=model_level_now_data_index,
columns=[ 'TOTAL' ])
)
else:
logger.debug("Model "+str(model_num)+" "+model_name+" "
+"with plot name "+model_plot_name+" "
+"file: "+model_stat_file+" exists")
model_level_now_stat_file_data = pd.read_csv(
model_stat_file, sep=" ", skiprows=1,
skipinitialspace=True, header=None
)
model_level_now_stat_file_data.rename(
columns=dict(zip(
model_level_now_stat_file_data \
.columns[:nbase_columns],
stat_file_base_columns
)), inplace=True
)
line_type = model_level_now_stat_file_data['LINE_TYPE'][0]
stat_file_line_type_columns = (
plot_util.get_stat_file_line_type_columns(logger,
met_version,
line_type)
)
model_level_now_stat_file_data.rename(
columns=dict(zip(
model_level_now_stat_file_data \
.columns[nbase_columns:],
stat_file_line_type_columns
)), inplace=True
)
model_level_now_stat_file_data_fcstvaliddates = (
model_level_now_stat_file_data.loc[:] \
['FCST_VALID_BEG'].values
)
model_level_now_data = (
pd.DataFrame(np.nan, index=model_level_now_data_index,
columns=stat_file_line_type_columns)
)
#model_level_now_stat_file_data.fillna(
# {'FCST_UNITS':'NA', 'OBS_UNITS':'NA', 'VX_MASK':'NA'},
# inplace=True
#)
for expected_date in expected_stat_file_dates:
if expected_date in \
model_level_now_stat_file_data_fcstvaliddates:
matching_date_idx = (
model_level_now_stat_file_data_fcstvaliddates \
.tolist().index(expected_date)
)
model_level_now_stat_file_data_idx = (
model_level_now_stat_file_data \
.loc[matching_date_idx][:]
)
for col in stat_file_line_type_columns:
#### EMC-verif_global changes for PRMSL,PRES/Z0
#### SPFH, O3MR
if fcst_var_name == 'PRMSL' \
or \
(fcst_var_name == 'PRES' \
and fcst_var_level == 'Z0'):
if col in ['FBAR', 'OBAR']:
scale = 1/100.
elif col in ['FFBAR', 'FOBAR', 'OOBAR']:
scale = 1/(100.*100.)
else:
scale = 1
elif fcst_var_name == 'SPFH':
if col in ['FBAR', 'OBAR']:
scale = 1000
elif col in ['FFBAR', 'FOBAR', 'OOBAR']:
scale = 1000*1000
else:
scale = 1
elif fcst_var_name == 'O3MR':
if col in ['FBAR', 'OBAR']:
scale = 1e6
elif col in ['FFBAR', 'FOBAR', 'OOBAR']:
scale = 1e6*1e6
else:
scale = 1
else:
scale = 1
model_level_now_data.loc[
(model_plot_name,
fcst_var_level,
expected_date)
][col] = (
model_level_now_stat_file_data_idx \
.loc[:][col] * scale
)
else:
logger.warning("Model "+str(model_num)+" "+model_name+" "
+"with plot name "+model_plot_name+" "
+"file: "+model_stat_file+" does not exist")
model_level_now_data = (
pd.DataFrame(np.nan,
index=model_level_now_data_index,
columns=[ 'TOTAL' ])
)
if vl > 0:
model_now_data = pd.concat(
[model_now_data, model_level_now_data], sort=True
)
else:
model_now_data = model_level_now_data
if model_num > 1:
model_data = pd.concat([model_data, model_now_data], sort=True)
else:
model_data = model_now_data
# Calculate statistics and plots
logger.info("Calculating and plotting statistics")
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_values, stat_values_array, stat_plot_name = (
plot_util.calculate_stat(logger, model_data, stat)
)
if event_equalization == "True":
logger.debug("Doing event equalization")
for l in range(len(stat_values_array[:,0,0])):
for vl in range(len(fcst_var_levels)):
stat_values_array[l,:,vl,:] = (
np.ma.mask_cols(stat_values_array[l,:,vl,:])
)
np.ma.set_fill_value(stat_values_array, np.nan)
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)
ax.set_xlim([plot_time_dates[0],plot_time_dates[-1]])
day_interval = int(len(plot_time_dates)/nticks)
ax.set_xticks(plot_time_dates[::day_interval])
ax.xaxis.set_major_formatter(md.DateFormatter('%d%b%Y'))
ax.get_xticklabels()[0].set_ha('left')
ax.get_xticklabels()[-1].set_ha('right')
if len(plot_time_dates) > 60:
ax.xaxis.set_minor_locator(md.MonthLocator())
else:
ax.xaxis.set_minor_locator(md.DayLocator())
if ax.is_last_row() \
or (nmodels % 2 != 0 and model_num == nmodels -1):
ax.set_xlabel(date_type.title()+" Date")
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
obs_plotted = False
get_clevels = True
make_colorbar = False
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_stat_values_array = stat_values_array[0,model_idx,:,:]
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 stat_values_array[1,model_idx,:,:].mask.all():
logger.debug("Plotting observations from "+model_name)
obs_stat_values_array = (
stat_values_array[1,model_idx,:,:]
)
CF1 = ax.contourf(xmesh, ymesh, obs_stat_values_array,
cmap=cmap, extend='both')
C1 = ax.contour(xmesh, ymesh, obs_stat_values_array,
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_stat_values_array
- stat_values_array[1,model_idx,:,:]
)
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_stat_values_array.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_stat_values_array)) > 100:
bias_spacing = 2.25
elif np.max(np.abs(model_stat_values_array)) > 10:
bias_spacing = 2
else:
bias_spacing = 1.75
clevels_bias = plot_util.get_clevels(
model_stat_values_array, bias_spacing
)
CF1 = ax.contourf(xmesh, ymesh,
model_stat_values_array,
levels=clevels_bias,
cmap=cmap_bias,
extend='both')
C1 = ax.contour(xmesh, ymesh, model_stat_values_array,
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_stat_values_array,
levels=CF1.levels,
cmap=cmap_bias,
extend='both')
C = ax.contour(xmesh, ymesh, model_stat_values_array,
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:
ax.set_title(model_plot_name, loc='left')
model1_name = model_name
model1_plot_name = model_plot_name
model1_stat_values_array = model_stat_values_array
if not model_stat_values_array.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_stat_values_array,
levels=levels, cmap=cmap,
extend='both')
elif stat in ['rmse', 'rmse_md', 'rmse_pv']:
cmax = np.nanmax(model_stat_values_array)
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_stat_values_array,
levels=levels, cmap=cmap,
extend='both')
else:
CF1 = ax.contourf(xmesh, ymesh,
model_stat_values_array,
cmap=cmap,
extend='both')
C1 = ax.contour(xmesh, ymesh, model_stat_values_array,
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_stat_values_array - model1_stat_values_array
)
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(', '),
str(datetime.date.fromordinal(int(
plot_time_dates[0])
).strftime('%d%b%Y')),
str(datetime.date.fromordinal(int(
plot_time_dates[-1])
).strftime('%d%b%Y')),
verif_case
)
forecast_lead_title = plot_title.get_lead_title(fcst_lead[:-4])
full_title = (
stat_plot_name+'\n'
+var_info_title+', '+vx_mask_title+'\n'
+date_info_title+', '+forecast_lead_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+'_fhr'+fcst_lead[:-4]
+'_'+grid_vx_mask+'.png')
logger.info("Saving image as "+savefig_name)
plt.savefig(savefig_name)
plt.clf()
plt.close('all')