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,158 @@
+"""
+QBO Phase plots and QBO Index: UserScript, Stat-Analysis
+==========================================================================
+
+model_applications/
+s2s_stratosphere/
+UserScript_fcstGFS_obsERA_StratosphereQBO.py
+
+"""
+
+##############################################################################
+# Scientific Objective
+# --------------------
+#
+# This use case performs several differenc calculations and creates multiple plots
+# to evaluate QBO.  First, the zonal and meridional means on U from -10 S to 10N 
+# latitude are computed.  Then, an EOF analysis is performed on the zonal and meridional
+# mean data and and two Phase diagrams of QBO are created.  The two phase diagrams are 
+# phase circuits which gives one plot for each start date listed in the configuration file 
+# for model and observations, and another phase plot for the whole time period.  Additionally
+# the zonal and meridional mean at 30 and 50mb are output as time series, and the bias (ME) 
+# is computed.  The 30 and 50mb U time series are also plotted.
+#
+
+##############################################################################
+# 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 runs the UserScript wrapper tool to run a user provided script,
+# in this case, stratosphere_qbo_driver.py to compute the zonal and meridional
+# mean, compute anomalies, create EOFs, create Phase diagrams, output U at 30 and
+# 50mb, and create time series of 30 and 50mb winds.  StatAnalysis is then run on 
+# the 30 and 50mb winds to compute the bias (ME).
+#
+
+##############################################################################
+# METplus Workflow
+# ----------------
+#
+# This use case does not loop but plots the entire time period of data.  The 
+# following tools are run once: UserScript, StatAnalysis
+# 
+
+##############################################################################
+# 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
+# with the -c option, i.e. -c 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.
+#
+
+##############################################################################
+# Running METplus
+# ---------------
+#
+# This use case can be run two ways:
+#
+# 1) Passing in UserScript_fcstGFS_obsERA_StratosphereQBO.conf, 
+# then a user-specific system configuration file::
+#
+#        run_metplus.py -c /path/to/METplus/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf -c /path/to/user_system.conf
+#
+# 2) Modifying the configurations in parm/metplus_config, then passing in UserScript_fcstGFS_obsERA_StratosphereQBO.conf:
+#
+#        run_metplus.py -c /path/to/METplus/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf
+#
+# The former method is recommended. Whether you add them to a user-specific configuration file or modify the metplus_config files, the following variables must be set correctly:
+#
+# * **INPUT_BASE** - Path to directory where sample data tarballs are unpacked (See Datasets section to obtain tarballs). This is not required to run METplus, but it is required to run the examples in parm/use_cases
+# * **OUTPUT_BASE** - Path where METplus output will be written. This must be in a location where you have write permissions
+# * **MET_INSTALL_DIR** - Path to location where MET is installed locally
+#
+#  and for the [exe] section, you will need to define the location of NON-MET executables.
+#  No executables are required for performing this use case.
+#
+# Example User Configuration File::
+#
+#   [dir]
+#   INPUT_BASE = /path/to/sample/input/data
+#   OUTPUT_BASE = /path/to/output/dir
+#   MET_INSTALL_DIR = /path/to/met-X.Y
+#
+#
+
+##############################################################################
+# 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, ERA_GFS_QBO_circuits.png, ERA5_QBO_PhaseSpace.png,
+# ERA_GFS_timeseries_30mb_u_201710_201802.png, and ERA_GFS_timeseries_30mb_u_201710_201802.png.
+# These graphics will be output to the path specified as OUTPUT_DIR/plots, using 
+# PLOT_PHASE_CIRCUITS_OUTPUT_NAME, PLOT_PHASE_SPACE_OUTPUT_NAME, PLOT_TIME_SERIES_OUTPUT_NAME_30,
+# and PLOT_TIME_SERIES_OUTPUT_NAME_50 in the [user_env_vars] section.  Additionally matched pair
+# .stat files will be output to OUTPUT_DIR/mpr, and the computed bias will be output to 
+# OUTPUT_DIR/StatAnalysis.
+#
+
+##############################################################################
+# 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)