Feature 2557 stratosphere qbo

.github/parm/use_case_groups.json

@@ -197,12 +197,12 @@
   {
     "category": "s2s_mid_lat",
     "index_list": "0-2",
-    "run": false
+    "run": true
   },
   {
     "category": "s2s_mid_lat",
     "index_list": "3",
-    "run": false
+    "run": true
   },
   {
     "category": "s2s_mjo",

docs/Contributors_Guide/add_use_case.rst

@@ -167,7 +167,8 @@ Use Cases that Exceed GitHub Actions Memory Limit
 
 Use Cases that Exceed GitHub Actions Disk Space
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-- *model_applications/s2s/UserScript_fcstGFS_obsERA_StratospherePolar*
+- *model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar*
+- *model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO*
 
 .. _use_case_documentation:
 

docs/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.py

@@ -0,0 +1,159 @@
+"""
+QBO Phase plots and QBO Index: UserScript, Stat-Analysis
+==========================================================================
+
+model_applications/
+s2s_stratosphere/
+UserScript_fcstGFS_obsERA_StratosphereQBO.py
+
+"""
+
+##############################################################################
+# Scientific Objective
+# --------------------
+#
+# Many common modes of variability in the troposphere have stratospheric teloconnection 
+# pathways.  This use case performs evaluation of the Quasi-biennial Oscillation (QBO), 
+# one of the key players of stratosphic variability, using several different calculations 
+# and plots.  Specifically, phase diagrams can be used to compare the QBO phase progression 
+# between the model and observations.  Additionally, timeseries of U at 30 and 50 mb are also 
+# plotted to compare the speed of propagation of the model versus the observations.  Continuous 
+# statistics (bias, RMSE, etc) are calculated for U at 30 and 50mb, and are also computed
+# separately to evaluate QBO in the easterly phase (U < 0) versus the westerly phase (U > 0).
+#
+
+##############################################################################
+# Datasets
+# --------
+#
+#  * GFS: 0-24 hour forecasts for 10/2017 - 2/2018
+#  * ERA: 30 year climatology for EOFs and 10/2017 - 2/2018  
+#
+
+##############################################################################
+# METplus Components
+# ------------------
+#
+# This use case calls UserScript and StatAnalysis.  The UserScript accesses calculations
+# as part of METcalcpy, METplotpy, and METdataio.  For it to run, the following versions 
+# of those repositories are needed:
+#  * METcalcpy 3.0.0
+#  * METplotpy 3.0.0
+#  * METdataio 2.1
+#
+
+##############################################################################
+# METplus Workflow
+# ----------------
+#
+# This use case does not loop but UserScript and StatAnalysis are each run once. 
+# The UserScript call runs the driver script stratosphere_qbo_driver.py which first 
+# computes zonal and meridional means using directional_means.py in METcalcpy on U from 
+# -10 S to 10N latitude.  Then, an EOF analysis is performed on this zonal and meridional 
+# mean data, and two phase diagrams of QBO are created using the plot_qbo_phase_circuits and 
+# plot_qbo_phase_space functions from stratosphere_plots.py in METplotpy.  Additionally the 
+# zonal and meridional mean at 30 and 50mb are output as time series to matched pair (MPR) 
+# files using write_mpr.py in METcalcpy and are also plotted as timeseries using the 
+# plot_u_timeseries function from stratosphere_plots.py in METplotpy.  Finally StatAnalysis is 
+# run on the 30 and 50mb U mpr files to compute the bias (ME).
+# 
+
+##############################################################################
+# METplus Configuration
+# ---------------------
+#
+# METplus first loads all of the configuration files found in parm/metplus_config,
+# then it loads any configuration files passed to METplus via the command line, i.e
+# parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf
+#
+# .. highlight:: bash
+# .. literalinclude:: ../../../../parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf
+#
+
+#############################################################################
+# MET Configuration
+# ---------------------
+#
+# METplus sets environment variables based on user settings in the METplus configuration file. 
+# See :ref:`How METplus controls MET config file settings<metplus-control-met>` for more details. 
+#
+# **YOU SHOULD NOT SET ANY OF THESE ENVIRONMENT VARIABLES YOURSELF! THEY WILL BE OVERWRITTEN BY METPLUS WHEN IT CALLS THE MET TOOLS!**
+#
+# If there is a setting in the MET configuration file that is currently not supported by METplus you'd like to control, please refer to:
+# :ref:`Overriding Unsupported MET config file settings<met-config-overrides>`
+#
+# **StatAnalysisConfig_wrapped**
+#
+# .. note:: See the :ref:`Stat-Analysis MET Configuration<series-analysis-met-conf>` section of the User's Guide for more information on the environment variables used in the file below:
+#
+# .. highlight:: bash
+# .. literalinclude:: ../../../../parm/met_config/STATAnalysisConfig_wrapped
+#
+
+##############################################################################
+# Python Embedding
+# ----------------
+#
+# This use case does not use python embedding.
+#
+
+##############################################################################
+# Python Scripting
+# ----------------
+#
+# This use case runs the stratospher_qbo_driver.py python script.  The processing
+# performed by the script are detailed in the :ref: `metplus-workflow` section.
+#
+
+##############################################################################
+# Running METplus
+# ---------------
+#
+# Pass the use case configuration file to the run_metplus.py script along with any 
+# user-specific system configuration files if desired:
+#
+#        run_metplus.py /path/to/METplus/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf /path/to/user_system.conf
+#
+# See :ref:`running-metplus` for more information.
+#
+
+##############################################################################
+# Expected Output
+# ---------------
+#
+# A successful run will output the following both to the screen and to the logfile::
+#
+#   INFO: METplus has successfully finished running.
+#
+# There should be 4 graphics output to the plot directory in the location set as OUTPUT_DIR
+# in the [user_env_vars] section of the configuration file:
+#  * ERA_GFS_QBO_circuits.png
+#  * ERA5_QBO_PhaseSpace.png
+#  * ERA_GFS_timeseries_30mb_u_201710_201802.png
+#  * ERA_GFS_timeseries_50mb_u_201710_201802.png
+# The name of the output graphics can be changed using PLOT_PHASE_CIRCUITS_OUTPUT_NAME, 
+# PLOT_PHASE_SPACE_OUTPUT_NAME, PLOT_TIME_SERIES_OUTPUT_NAME_30, and PLOT_TIME_SERIES_OUTPUT_NAME_50 
+# also in the [user_env_vars] section.  Additionally many matched pair .stat files will be output to 
+# OUTPUT_DIR/mpr, and tow computed continuous statistics will be output to OUTPUT_DIR/StatAnalysis:
+#  * GFS_ERA_20171001_20180228_210000L_zonal_wind_byphase_CNT.stat
+#  * GFS_ERA_20171001_20180228_210000L_zonal_wind_CNT.stat
+#
+
+##############################################################################
+# Keywords
+# --------
+#
+# .. note::
+#
+#   * UserScriptUseCase
+#   * S2SAppUseCase
+#   * S2SStratosphereAppUseCase
+#   * StatAnalysisUseCase
+#   * METcalcpyUseCase
+#   * METplotpyUseCase
+#
+#   Navigate to the :ref:`quick-search` page to discover other similar use cases.
+#
+#
+#
+# sphinx_gallery_thumbnail_path = '_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratosphereQBO.png'

internal/tests/use_cases/all_use_cases.txt

@@ -163,6 +163,7 @@ Category: s2s_mjo
 Category: s2s_stratosphere
 0:: UserScript_fcstGFS_obsERA_StratosphereBias:: model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereBias.conf:: weatherregime_env,cartopy,metdataio
 #X::UserScript_fcstGFS_obsERA_StratospherePolar:: model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar.conf:: weatherregime_env,cartopy,metdataio
+#X::UserScript_fcstGFS_obsERA_StratosphereQBO:: model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf:: weatherregime_env,cartopy,metdataio
 
 
 Category: short_range

parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar.conf

@@ -12,6 +12,7 @@
 ###
 
 PROCESS_LIST = UserScript(means), StatAnalysis(sanal_cnt), UserScript(plots_t), UserScript(plots_u)
+#PROCESS_LIST = UserScript(means)
 
 
 SCRUB_STAGING_DIR = False
@@ -33,7 +34,9 @@ VALID_TIME_FMT = %Y%m%d
 VALID_BEG = 20180201
 VALID_END = 20180228
 VALID_INCREMENT = 30d
-LEAD_SEQ = begin_end_incr(0,240,3),begin_end_incr(252,384,12)
+#LEAD_SEQ = begin_end_incr(0,240,3),begin_end_incr(252,384,12)
+LEAD_SEQ = begin_end_incr(0,384,24)
+#LEAD_SEQ = 0
 
 LOOP_ORDER = processes
 
@@ -54,28 +57,41 @@ USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_stratosphere/Us
 # observation_template, forecast_template
 USER_SCRIPT_INPUT_TEMPLATE_LABELS = OBS_INPUT, FCST_INPUT
 
-USER_SCRIPT_COMMAND = {PARM_BASE}/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/polar_t_u_driver.py time {lead?fmt=%HHH}
+USER_SCRIPT_COMMAND = {PARM_BASE}/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/polar_t_u_driver.py {lead?fmt=%HHH}
 
 
 
 [user_env_vars]
 MODEL_NAME = GFS
 
+FCST_T_VAR = T
+FCST_U_VAR = u
+FCST_TIME_VAR = time
+FCST_LAT_DIM = lat
+FCST_LON_DIM = lon
+FCST_PRES_DIM = pres
+
+OBS_T_VAR = T
+OBS_U_VAR = u
+OBS_TIME_VAR = time
+OBS_LAT_DIM = lat
+OBS_LON_DIM = lon
+OBS_PRES_DIM = pres
+
 OUTPUT_DIR = {OUTPUT_BASE}/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar
 
 PLOT_OUTPUT_DIR = {OUTPUT_DIR}/plots
 
-PLOT_INPUT_FILE = {OUTPUT_BASE}/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/SeriesAnalysis/zonal_mean_stats_2018_02.nc
-PLOT_T_BIAS_LEVELS = -6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6
+PLOT_T_BIAS_LEVELS = -7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9
 PLOT_T_BIAS_TITLE = GFS vs ERA Polar Cap Temperature Bias (ME) 02/2018
-PLOT_T_RMSE_LEVELS = 0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6
+PLOT_T_RMSE_LEVELS = 0,1,2,3,4,5,6,7,8,9,10
 PLOT_T_RMSE_TITLE = GFS vs ERA Polar Cap Temperature RMSE 02/2018
 PLOT_T_BIAS_OUTPUT_FILE = ME_2018_02_polar_cap_T.png
 PLOT_T_RMSE_OUTPUT_FILE = RMSE_2018_02_polar_cap_T.png
 
-PLOT_U_BIAS_LEVELS = -6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6
+PLOT_U_BIAS_LEVELS = -4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12
 PLOT_U_BIAS_TITLE = GFS vs ERA Polar Vortex U Bias (ME) 02/2018
-PLOT_U_RMSE_LEVELS = 0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5
+PLOT_U_RMSE_LEVELS = 0,2,4,6,8,10,12,14,16,18,20,22,24
 PLOT_U_RMSE_TITLE = GFS vs ERA Polar Vortex U RMSE 02/2018
 PLOT_U_BIAS_OUTPUT_FILE = ME_2018_02_polar_vortex_U.png
 PLOT_U_RMSE_OUTPUT_FILE = RMSE_2018_02_polar_vortex_U.png

parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/bias_rmse_plot_driver.py

@@ -80,8 +80,6 @@ def main():
     """
     Create plots
     """
-    print(plot_bias_levels)
-    print(plot_rmse_levels)
     plot_polar_bias(pleads,plevels,pme,plot_bias_output_file,plot_bias_title,plot_bias_levels)
     plot_polar_rmse(pleads,plevels,prmse,plot_rmse_output_file,plot_rmse_title,plot_rmse_levels)
 

parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/polar_t_u_driver.py

@@ -17,8 +17,7 @@ def main():
     """
     Read arguments
     """
-    timevar = sys.argv[1]
-    leadvar = sys.argv[2]
+    leadvar = sys.argv[1]
 
     """
     Read METplus filename lists
@@ -40,52 +39,74 @@ def main():
     output_dir = os.environ['OUTPUT_DIR']
     full_output_dir = os.path.join(output_dir,'mpr')
 
+    """
+    Read variable/dimension names
+    """
+    obs_tvar = os.environ.get('OBS_T_VAR','T')
+    obs_uvar = os.environ.get('OBS_U_VAR','u')
+    obs_latvar = os.environ.get('OBS_LAT_VAR','latitude')
+    obs_timevar = os.environ.get('OBS_TIME_VAR','time')
+    obs_latdim = os.environ.get('OBS_LAT_DIM','lat')
+    obs_londim = os.environ.get('OBS_LON_DIM','lon')
+    obs_presdim = os.environ.get('OBS_PRES_DIM','pres')
+
+    fcst_tvar = os.environ.get('FCST_T_VAR','T')
+    fcst_uvar = os.environ.get('FCST_U_VAR','u')
+    fcst_latvar = os.environ.get('FCST_LAT_VAR','latitude')
+    fcst_timevar = os.environ.get('FCST_TIME_VAR','time')
+    fcst_latdim = os.environ.get('FCST_LAT_DIM','lat')
+    fcst_londim = os.environ.get('FCST_LON_DIM','lon')
+    fcst_presdim = os.environ.get('FCST_PRES_DIM','pres')
+
     """
     Read dataset
     """
     file_reader = read_netcdf.ReadNetCDF()
     dsO = file_reader.read_into_xarray(obs_infiles)[0]
-    dsO = dsO.rename({'lat':'latitude',
-                    'lon':'longitude'})
+    dsO = dsO.rename({obs_latdim:'latitude',
+                    obs_londim:'longitude',
+                    obs_presdim:'pres'})
     file_reader2 = read_netcdf.ReadNetCDF()
     dsF = file_reader2.read_into_xarray(fcst_infiles)[0]
-    dsF = dsF.rename({'lat':'latitude',
-                    'lon':'longitude'})
+    dsF = dsF.rename({fcst_latdim:'latitude',
+                    fcst_londim:'longitude',
+                    fcst_presdim:'pres'})
 
     """
     Limit Dataset to 100 - 1 hPa
     """
     dsO = dsO.sel(pres=slice(1,100))
     dsF = dsF.sel(pres=slice(1,100))
 
+    """
+    Assign Latitude Coordinate since it doesn't work
+    """
+    dsO = dsO.assign_coords({obs_latvar:dsO[obs_latvar].values})
+    dsF = dsF.assign_coords({fcst_latvar:dsF[fcst_latvar].values})
 
     """
     Create Polar Cap Temparatures for Forecast and Obs
     """
-    TzmO = directional_means.zonal_mean(dsO.T)
-    TzmO.assign_coords(lat=("latitude",dsO.latitude.values[:,0]))
+    dsO = dsO.assign_coords({obs_latvar:dsO[obs_latvar].values})
+    TzmO = directional_means.zonal_mean(dsO[obs_tvar])
     TO_6090 = directional_means.meridional_mean(TzmO, 60, 90)
-    TzmF = directional_means.zonal_mean(dsF.T)
-    TzmF.assign_coords(lat=("latitude",dsF.latitude.values[:,0]))
+    TzmF = directional_means.zonal_mean(dsF[fcst_tvar])
     TF_6090 = directional_means.meridional_mean(TzmF, 60, 90)
 
     """
     Create Polar Vortex Winds
     """
-    UzmO = directional_means.zonal_mean(dsO.u)
-    UzmO.assign_coords(lat=("latitude",dsO.latitude.values[:,0]))
+    UzmO = directional_means.zonal_mean(dsO[obs_uvar])
     UO_6090 = directional_means.meridional_mean(UzmO, 50, 80)
-    UzmF = directional_means.zonal_mean(dsF.u)
-    UzmF.assign_coords(lat=("latitude",dsF.latitude.values[:,0]))
+    UzmF = directional_means.zonal_mean(dsF[fcst_uvar])
     UF_6090 = directional_means.meridional_mean(UzmF, 50, 80)
 
-
     """
     Add P to the levels since they are pressure levels
     """
-    obs_lvls = ['P'+str(int(op)) for op in dsO.pres.values]
-    obs_lvls2 = [str(int(op)) for op in dsO.pres.values]
-    fcst_lvls = ['P'+str(int(fp)) for fp in dsF.pres.values]
+    obs_lvls = ['P'+str(int(op)) for op in dsO['pres'].values]
+    obs_lvls2 = [str(int(op)) for op in dsO['pres'].values]
+    fcst_lvls = ['P'+str(int(fp)) for fp in dsF['pres'].values]
 
     """
     Write output MPR files
@@ -94,14 +115,14 @@ def main():
     dlength = dlength1*2 
     modname = os.environ.get('MODEL_NAME','GFS')
     maskname = os.environ.get('MASK_NAME','FULL')
-    datetimeindex = dsF.indexes[timevar].to_datetimeindex()
+    datetimeindex = dsF.indexes[fcst_timevar]
     for i in range(len(datetimeindex)):
         valid_str = datetimeindex[i].strftime('%Y%m%d_%H%M%S')
         leadstr = str(int(leadvar)).zfill(2)+'0000'
         outobs = np.concatenate((TO_6090[i,:].values,UO_6090[i,:].values))
         outfcst = np.concatenate((TF_6090[i,:].values,UF_6090[i,:].values))
         write_mpr_file(outfcst,outobs,[0.0]*dlength,[0.0]*dlength,[leadstr]*dlength,[valid_str]*dlength,
-            ['000000']*dlength,[valid_str]*dlength,modname,'NA',['PolarCapT']*dlength1+['PolarVortexU']*dlength1,
+            ['000000']*dlength,[valid_str]*dlength,modname,['NA']*dlength,['PolarCapT']*dlength1+['PolarVortexU']*dlength1,
             ['K']*dlength1+['m/s']*dlength1,fcst_lvls*2,['PolarCapT']*dlength1+['PolarVortexU']*dlength1,
             ['K']*dlength1+['m/s']*dlength1,obs_lvls*2,maskname,obs_lvls2*2,full_output_dir,'polar_cap_T_stat_'+modname)