diff --git a/.github/parm/use_case_groups.json b/.github/parm/use_case_groups.json index ea51000c83..b32abff669 100644 --- a/.github/parm/use_case_groups.json +++ b/.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", diff --git a/docs/Contributors_Guide/add_use_case.rst b/docs/Contributors_Guide/add_use_case.rst index f66345ccbd..10fe4c9a7a 100644 --- a/docs/Contributors_Guide/add_use_case.rst +++ b/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: diff --git a/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratospherePolar.png b/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratospherePolar.png index 253480cd5a..327f899188 100644 Binary files a/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratospherePolar.png and b/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratospherePolar.png differ diff --git a/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratosphereQBO.png b/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratosphereQBO.png new file mode 100644 index 0000000000..e27a1e2a44 Binary files /dev/null and b/docs/_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratosphereQBO.png differ diff --git a/docs/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.py b/docs/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.py new file mode 100644 index 0000000000..79dacccb40 --- /dev/null +++ b/docs/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.py @@ -0,0 +1,164 @@ +""" +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 stratosphereic varability, using several different calculations +# and plots. Specifically, phase diagrams can be used to compare the QBO plhase progression +# between the model and observations. Additionally, timeseries of U at 30 and 50mb are also +# plotted to compare the speed of propagation of QBO in the model versus observations. +# Continuous statistics (ME, RMSE, etc.) are calculated for U at 30 and 50mb, and are also +# computed separately to evaluate QBO in the easterly phase versus the westerly phase. +# + +############################################################################## +# Datasets +# -------- +# +# * GFS: 0-24 hour forecasts for 10/2017 - 2/2018 +# * ERA: 30 year climatology for EOFs and 10/2017 - 2/2018 +# +# Data for this use case is not contained in the sample data tar files due to +# its size. Rather, it is stored as additional data in a separate tar file, +# named additional_data_UserScript_fcstGFS_obsERA_StratosphereQBO.tar.gz at +# https://dtcenter.ucar.edu/dfiles/code/METplus/METplus_Data/v6.0/. +# + +############################################################################## +# METplus Components +# ------------------ +# +# This use case calls UserScript and StatAnalysis. The UserScript accesses calculation +# in METcalcpy, METplotpy, and METdataio. For it to run, the following versions of those +# repositories are needed, METcalcpy 3.0.0, METplotpy 3.0.0, and 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` 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` +# +# **StatAnalysisConfig_wrapped** +# +# .. note:: See the :ref:`Stat-Analysis MET Configuration` 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 METplus Workflow section above. +# + +############################################################################## +# 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' diff --git a/internal/tests/use_cases/all_use_cases.txt b/internal/tests/use_cases/all_use_cases.txt index daef1a124e..9c35982c01 100644 --- a/internal/tests/use_cases/all_use_cases.txt +++ b/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 diff --git a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar.conf b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar.conf index c3cbbfc4ce..d79a37d461 100644 --- a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar.conf +++ b/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 diff --git a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/README b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/README deleted file mode 100644 index de1e1d0abe..0000000000 --- a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/README +++ /dev/null @@ -1,9 +0,0 @@ -These files are a use case to show how to compute meridial and zonal means - -You need METcalcpy and METdatadb in your python path or your conda environment -i.e. -export PYTHONPATH=/METdatadb:/METcalcpy - -The file SSWC_v1.0_varFull_ERAi_d20130106_s20121107_e20130307_c20160701.nc needs to be -on disk somewhere on your computer and referenced correctly in the file -meridial_mean.yaml diff --git a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/bias_rmse_plot_driver.py b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/bias_rmse_plot_driver.py index 451a50aabc..cc9f0b1553 100755 --- a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/bias_rmse_plot_driver.py +++ b/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) diff --git a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/polar_t_u_driver.py b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/polar_t_u_driver.py index 87e2e433c5..774dec04a1 100755 --- a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratospherePolar/polar_t_u_driver.py +++ b/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,17 +39,38 @@ 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 @@ -58,34 +78,35 @@ def main(): 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) diff --git a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf new file mode 100644 index 0000000000..9c496077fc --- /dev/null +++ b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.conf @@ -0,0 +1,169 @@ +[config] + +# Documentation for this use case can be found at +# https://metplus.readthedocs.io/en/latest/generated/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereBias.html + +# For additional information, please see the METplus Users Guide. +# https://metplus.readthedocs.io/en/latest/Users_Guide + +### +# Processes to run +# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#process-list +### +#PROCESS_LIST = UserScript(create_obs_eofs_filelist), UserScript(compute_plot_qbo), StatAnalysis(qbo_bias) +PROCESS_LIST = UserScript(compute_plot_qbo), StatAnalysis(qbo_bias) + + +### +# Time Info +# LOOP_BY options are INIT, VALID, RETRO, and REALTIME +# If set to INIT or RETRO: +# INIT_TIME_FMT, INIT_BEG, INIT_END, and INIT_INCREMENT must also be set +# If set to VALID or REALTIME: +# VALID_TIME_FMT, VALID_BEG, VALID_END, and VALID_INCREMENT must also be set +# LEAD_SEQ is the list of forecast leads to process +# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#timing-control +### + +LOOP_BY = INIT +INIT_TIME_FMT = %Y%m%d%H +INIT_BEG = 2017100100 +INIT_END = 2018022821 +INIT_INCREMENT = 86400 +LEAD_SEQ = begin_end_incr(0,21,3) + + +USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE + +LOOP_ORDER = processes + +### +# User Environment Variables +### +[user_env_vars] +# Variable to determine if Zonal and Meridional means are need to be computed on the data used +# to create EOFS. Set to True to compute zonal and meridional means. If set to False, a file +# should be provided in EOF_DATA_FILE_NAME that specifies the location of a dataset with zonal +# and meridional means. The data should have the coordinates (time, pres), and the variable +# names should be the same as what are set in OBS_TIME_VAR, OBS_LAT_VAR, and OBS_U_VAR below. +# Note, this is only for EOF calculation, and not for the data used to create plots. +COMPUTE_EOF_ZONAL_MERIDIONAL_MEAN = False +# Ignored if the above is set to True +EOF_DATA_FILE_NAME = {INPUT_BASE}/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/erai_zonal_mean_u.nc + +# Option to save out the Zonal/Meridional Mean data that was computed for the EOF calculation +SAVE_EOF_ZONAL_MERIDIONAL_MEAN = False +# Name of the output Zonal/meridional mean file. Ignored if The above is set to False +ZONAL_MERIDIONAL_MEAN_EOF_FILE_NAME = {INPUT_BASE}/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/erai_zonal_mean_u.nc + +# Name of forecast variables for time, latitude, longitude, and U, and +# the dimension names for latitude and longitude +FCST_TIME_VAR = time +FCST_LAT_VAR = latitude +FCST_U_VAR = u +FCST_LON_DIM = lon +FCST_LAT_DIM = lat + +# Name of observation variables for time, latitude, longitude, and U and +# the dimension names for latitude and longitude +OBS_TIME_VAR = time +OBS_LAT_VAR = latitude +OBS_U_VAR = u +OBS_LON_DIM = lon +OBS_LAT_DIM = lat + +# Minimum and maximum pressure levels to use in the QBO calculation +PRES_LEV_MIN = 10 +PRES_LEV_MAX = 100 + +# Minimum and maximim latitude to use in the QBO calculation +LAT_MIN = -10 +LAT_MAX = 10 + +# Name of the model used +MODEL_NAME = GFS + +# Directory where output plots and statistics should be sent +OUTPUT_DIR = {OUTPUT_BASE}/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO + +# Plot Settings +# Settings for the Circut Phase Diagram(s) +# List or single value of start times, each of which will have a phase diagram created +PLOT_START_DATES = 10-1-2017 +# Number of days for each diagram. If the PLOT_START_DATES contains more than one value and +# PLOT_STAT_DATES[0]+PLOT_PERIODS != PLOT_START_DATES[1], the data between +# PLOT_STAT_DATES[0]+PLOT_PERIODS and PLOT_START_DATES[1] will be omitted +PLOT_PERIODS = 151 +# Full path to the output name of the phase circuits plot +PLOT_PHASE_CIRCUTS_OUTPUT_NAME = ERA_GFS_QBO_circuits.png + +# Title and output name of the phase space plot +PLOT_PHASE_SPACE_TITLE = QBO Phase Space from ERA5 (10/2017 - 2/2018) +PLOT_PHASE_SPACE_OUTPUT_NAME = ERA5_QBO_PhaseSpace.png + +# Time Series plot titles and output names for 30mb and 50mb +PLOT_TIME_SERIES_TITLE_30 = ERA5 and GFS Zonal Mean U at 30mb (10/2017 - 2/2018) +PLOT_TIME_SERIES_OUTPUT_NAME_30 = ERA_GFS_timeseries_30mb_u_201710_201802.png +PLOT_TIME_SERIES_TITLE_50 = ERA5 and GFS Zonal Mean U at 50mb (10/2017 - 2/2018) +PLOT_TIME_SERIES_OUTPUT_NAME_50 = ERA_GFS_timeseries_50mb_u_201710_201802.png + + +### +# Creates the file listing for the EOF creation UserScript Settings +# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript +### +[create_obs_eofs_filelist] +INIT_BEG = 1991010100 +INIT_END = 2020123121 +INIT_INCREMENT = 86400 +LEAD_SEQ = 0 + +# Template of filenames to input to the user-script +USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/ERA/{valid?fmt=%Y%m}/ERA_{valid?fmt=%Y%m%d%H}.nc + +# Name of the file containing the listing of input files +# The option is OBS_EOF_INPUT; FCST_EOF_INPUT is not supported +# *** Make sure the order is the same as the order of templates listed in USER_SCRIPT_INPUT_TEMPLATE +USER_SCRIPT_INPUT_TEMPLATE_LABELS = OBS_EOF_INPUT + +USER_SCRIPT_COMMAND = echo Populated file list for obs EOF Input + + +### +# QBO Calculation and plotting UserScript Settings +# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript +### +[compute_plot_qbo] +USER_SCRIPT_RUNTIME_FREQUENCY = RUN_ONCE + +# Template of filenames to input to the user-script +USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/ERA/{valid?fmt=%Y%m}/ERA_{valid?fmt=%Y%m%d%H}.nc, {INPUT_BASE}/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/GFS/GFS_{init?fmt=%Y%m%d%H}_f{lead?fmt=%HHH}h.nc + +# Name of the file containing the listing of input files +# The options are OBS_INPUT and FCST_INPUT +# *** Make sure the order is the same as the order of templates listed in USER_SCRIPT_INPUT_TEMPLATE +USER_SCRIPT_INPUT_TEMPLATE_LABELS = OBS_INPUT, FCST_INPUT + +USER_SCRIPT_COMMAND = {PARM_BASE}/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/stratosphere_qbo_driver.py + + +[qbo_bias] +MODEL1 = GFS +MODEL1_OBTYPE = ADPUPA + +STAT_ANALYSIS_CONFIG_FILE = {PARM_BASE}/met_config/STATAnalysisConfig_wrapped + +STAT_ANALYSIS_JOB1 = -job aggregate_stat -line_type MPR -out_line_type CNT -fcst_var QBO_U -by FCST_LEV -out_stat [out_stat_file]_zonal_wind_CNT.stat +STAT_ANALYSIS_JOB2 = -job aggregate_stat -line_type MPR -out_line_type CNT -fcst_var QBO_U -by FCST_LEV,DESC -out_stat [out_stat_file]_zonal_wind_byphase_CNT.stat + +MODEL_LIST = {MODEL1} +FCST_LEAD_LIST = 21 + +GROUP_LIST_ITEMS = MODEL_LIST +LOOP_LIST_ITEMS = + +MODEL1_STAT_ANALYSIS_LOOKIN_DIR = {OUTPUT_BASE}/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/mpr + +STAT_ANALYSIS_OUTPUT_DIR = {OUTPUT_BASE}/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/StatAnalysis + +MODEL1_STAT_ANALYSIS_OUT_STAT_TEMPLATE = {model?fmt=%s}_ERA_{obs_valid_beg?fmt=%Y%m%d}_{obs_valid_end?fmt=%Y%m%d}_{lead?fmt=%H%M%S}L diff --git a/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/stratosphere_qbo_driver.py b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/stratosphere_qbo_driver.py new file mode 100755 index 0000000000..07c46f8157 --- /dev/null +++ b/parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/stratosphere_qbo_driver.py @@ -0,0 +1,284 @@ +#!/usr/bin/env python3 + +""" +Create QBO + +""" +import os +import datetime +import numpy as np +import xarray as xr +import pandas as pd +from eofs.xarray import Eof + +import metcalcpy.pre_processing.directional_means as directional_means +import METreadnc.util.read_netcdf as read_netcdf +from metplotpy.contributed.stratosphere_diagnostics.stratosphere_plots import plot_qbo_phase_circuits,plot_qbo_phase_space,plot_u_timeseries +from metcalcpy.util.write_mpr import write_mpr_file + + +def process_single_file(infile,latvar,latdim,londim,lat_min,lat_max,pres_min,pres_max): + + file_reader = read_netcdf.ReadNetCDF() + ds = file_reader.read_into_xarray(infile)[0] + ds = ds.assign_coords({latdim:ds[latvar].values}) + #ds = ds.assign_coords({londim:ds[lonvar].values}) + + # Compute zonal mean + Uzm = directional_means.zonal_mean(ds,dimvar=londim) + + # Limit the data to 100 - 10 hPa + Uzm = Uzm.sel(pres=slice(pres_min,pres_max)) + + #Compute Meridional Mean + U_1010 = directional_means.meridional_mean(Uzm, lat_min, lat_max, dimvar=latdim) + + return U_1010 + + +def get_qbo_data(input_files,latvar,latdim,londim,timevar,lat_min,lat_max,pres_min,pres_max): + + # **** This could be edited to read all files at the same time, and it would run faster. + # However I did not have enough memory to run it this way. + # Read the first file + ds = process_single_file(input_files[0],latvar,latdim,londim,lat_min,lat_max,pres_min,pres_max) + # Read the rest of the files and concatenate + for f in input_files[1:]: + # Read in the data file + ds2 = process_single_file(f,latvar,latdim,londim,lat_min,lat_max,pres_min,pres_max) + # Save to array + ds = xr.concat([ds,ds2],dim=timevar) + + return ds + + +def main(): + + """ + Get the input variables + """ + # Read the variable names for lat, lon, and time + obs_latvar = os.environ.get('OBS_LAT_VAR','latitude') + #obs_lonvar = os.environ.get('OBS_LON_VAR','longitude') + obs_latdim = os.environ.get('OBS_LAT_DIM','lat') + obs_londim = os.environ.get('OBS_LON_DIM','lon') + obs_timevar = os.environ.get('OBS_TIME_VAR','time') + obs_uvar = os.environ.get('OBS_U_VAR','u') + + fcst_latvar = os.environ.get('FCST_LAT_VAR','latitude') + #fcst_lonvar = os.environ.get('FCST_LON_VAR','longitude') + fcst_latdim = os.environ.get('FCST_LAT_DIM','lat') + fcst_londim = os.environ.get('FCST_LON_DIM','lon') + fcst_timevar = os.environ.get('FCST_TIME_VAR','time') + fcst_uvar = os.environ.get('FCST_U_VAR','u') + + # Read the lat bounds, default -10 to 10 + lat_min = int(os.environ.get('LAT_MIN','-10')) + lat_max = int(os.environ.get('LAT_MAX','10')) + + # Read the Pressure level bounds, default is 100 to 10 + pres_min = int(os.environ.get('PRES_LEV_MIN','10')) + pres_max = int(os.environ.get('PRES_LEV_MAX','100')) + + # Read output directory + output_dir = os.environ['OUTPUT_DIR'] + + # Read in plotting inits and period + input_plot_inits = os.environ['PLOT_START_DATES'] + plot_period = int(os.environ['PLOT_PERIODS']) + + # Read in plotting outfile names and titles + plot_circuits_outname = os.environ.get('PLOT_PHASE_CIRCUTS_OUTPUT_NAME','QBO_circuits.png') + plot_phasespace_title = os.environ.get('PLOT_PHASE_SPACE_TITLE','QBO Phase Space') + plot_phasespace_outname = os.environ.get('PLOT_PHASE_SPACE_OUTPUT_NAME','QBO_PhaseSpace.png') + plot_timeseries_30_title = os.environ.get('PLOT_TIME_SERIES_TITLE_30','U 30mb') + plot_timeseries_30_outname = os.environ.get('PLOT_TIME_SERIES_OUTPUT_NAME_30','QBO_U_time_series_30.png') + plot_timeseries_50_title = os.environ.get('PLOT_TIME_SERIES_TITLE_50','U 50mb') + plot_timeseries_50_outname = os.environ.get('PLOT_TIME_SERIES_OUTPUT_NAME_50','QBO_U_time_series_50.png') + + + """ + Make output directories if they don't exist + """ + mpr_output_dir = os.path.join(output_dir,'mpr') + if not os.path.exists(mpr_output_dir): + os.makedirs(mpr_output_dir) + plot_output_dir = os.path.join(output_dir,'plots') + if not os.path.exists(plot_output_dir): + os.makedirs(plot_output_dir) + + """ + Read in the data for EOFs + Compute zonal and meridional means if not already computed + """ + compute_zmm = os.environ.get('COMPUTE_EOF_ZONAL_MERIDIONAL_MEAN','True') + if compute_zmm.lower() == 'true': + # Get input files + obs_eof_filetxt = os.environ.get('METPLUS_FILELIST_OBS_EOF_INPUT','') + with open(obs_eof_filetxt) as oef: + obs_input_files_eofs = oef.read().splitlines() + if (obs_input_files_eofs[0] == 'file_list'): + obs_input_files_eofs = obs_input_files_eofs[1:] + # Get the data + dsE = get_qbo_data(obs_input_files_eofs,obs_latvar,obs_latdim,obs_londim,obs_timevar,lat_min,lat_max,pres_min,pres_max) + else: + # Read in the name of the output data file + eof_data_file_name = os.environ['EOF_DATA_FILE_NAME'] + # Load data + dsE = xr.open_dataset(eof_data_file_name) + # Rename time dimension if it's not called time + if obs_timevar != 'time': + dsE.rename({obs_timevar:'time'}) + + + """ + Save the Data for the EOF calculation for future use, if desired + """ + save_zmm = os.environ.get('SAVE_EOF_ZONAL_MERIDIONAL_MEAN','False') + if save_zmm.lower() == 'true': + savefile = os.environ['ZONAL_MERIDIONAL_MEAN_EOF_FILE_NAME'] + dsE.to_netcdf(savefile) + + + """ + Read the other datasets + """ + # Get Obs listing of files for plotting and stats + obs_filetxt = os.environ.get('METPLUS_FILELIST_OBS_INPUT','') + with open(obs_filetxt) as of: + obs_input_files = of.read().splitlines() + if (obs_input_files[0] == 'file_list'): + obs_input_files = obs_input_files[1:] + # Get obs + dsO = get_qbo_data(obs_input_files,obs_latvar,obs_latdim,obs_londim,obs_timevar,lat_min,lat_max,pres_min,pres_max) + if obs_timevar != 'time': + dsO.rename({obs_timevar:'time'}) + + # Get model listing of files for plotting and stats + fcst_filetxt = os.environ.get('METPLUS_FILELIST_FCST_INPUT','') + with open(fcst_filetxt) as ff: + fcst_input_files = ff.read().splitlines() + if (fcst_input_files[0] == 'file_list'): + fcst_input_files = fcst_input_files[1:] + # Forecast + dsF = get_qbo_data(fcst_input_files,fcst_latvar,fcst_latdim,fcst_londim,fcst_timevar,lat_min,lat_max,pres_min,pres_max) + if fcst_timevar != 'time': + dsF.rename({fcst_timevar:'time'}) + + + """ + Compute Anomalies + """ + # Take Daily averages + dsE_daily = dsE.resample(time='1D').mean() + rean_trop_u = dsE_daily[obs_uvar] + + # Take monthly averages + trop_u_monthly = rean_trop_u.resample(time='1MS').mean() + + # Get monthly zonal wind anomalies + trop_u_anom = trop_u_monthly.groupby('time.month') - trop_u_monthly.groupby('time.month').mean('time') + + # get rid of the month coordinate which gets added + # (it can cause an error in the EOF analysis package) + trop_u_anom = trop_u_anom.drop_vars('month') + + + """ + Compute EOFs + """ + # Compute EOFs + solver = Eof(trop_u_anom) + pcs = solver.pcs(npcs=2, pcscaling=1) + eofs = solver.eofs(neofs=2) + + # Daily means from entire dataset + mmdd = rean_trop_u.time.dt.strftime("%m-%d") + mmdd.name = "mmdd" + trop_u_daily_clim = rean_trop_u.groupby(mmdd).mean('time') + + # Daily means for plot time period + dsO = dsO.drop_duplicates(dim='time') + dsO_daily = dsO.resample(time='1D').mean(dim='time') + dsO_daily = dsO_daily[obs_uvar] + mmdd1 = dsO_daily.time.dt.strftime("%m-%d") + mmdd1.name = "mmdd" + + # Daily anomalies for observations time period + rean_u_daily_anom = dsO_daily.groupby(mmdd1) - trop_u_daily_clim + rean_u_daily_anom = rean_u_daily_anom.drop_vars("mmdd") + + # Daily anomalies for forecast plot time period + dsF = dsF.sortby('time') + dsF_daily = dsF.resample(time='1D').mean(dim='time') + max_leads = dsF.resample(time='1D').max(dim='time').lead_time + #min_leads = dsF.resample(time='1D').min(dim='time').lead_time + #mean_leads = dsF_daily.lead_time + dsF_daily = dsF_daily[fcst_uvar] + + utrop_fcst_anoms = dsF_daily - trop_u_daily_clim.sel(mmdd=dsF_daily.time.dt.strftime("%m-%d")) + utrop_fcst_anoms.time.attrs['axis'] = 'T' + utrop_fcst_anoms = utrop_fcst_anoms.drop_vars('mmdd') # get rid of trailing dimension from the anomaly calculation + + # project daily reanalysis zonal winds + rean_qbo_pcs = solver.projectField(rean_u_daily_anom, eofscaling=1, neofs=2) + + # project daily rfcst zonal winds + rfcst_qbo_pcs = solver.projectField(utrop_fcst_anoms, eofscaling=1, neofs=2) + + + """ + EOF Phase Diagram Plots + """ + # Create Circuits plot + plot_inits = pd.DatetimeIndex([input_plot_inits]) + plot_qbo_phase_circuits(plot_inits,plot_period,rean_qbo_pcs,rfcst_qbo_pcs, + os.path.join(plot_output_dir,plot_circuits_outname)) + + # Create Phase Space plot + plot_qbo_phase_space(rean_qbo_pcs,eofs,plot_phasespace_title, + os.path.join(plot_output_dir,plot_phasespace_outname)) + + + """ + Time Series of U at 30 and 50mb Plot + """ + dsO_3050 = dsO_daily.sel(pres=[30,50]) + dsF_3050 = dsF_daily.sel(pres=[30,50]) + plot_u_timeseries(dsO_3050.sel(pres=30)['time'].values,dsO_3050.sel(pres=30).values, + dsF_3050.sel(pres=30)['time'].values,dsF_3050.sel(pres=30).values, + plot_timeseries_30_title,os.path.join(plot_output_dir,plot_timeseries_30_outname)) + plot_u_timeseries(dsO_3050.sel(pres=50)['time'].values,dsO_3050.sel(pres=50).values, + dsF_3050.sel(pres=50)['time'].values,dsF_3050.sel(pres=50).values, + plot_timeseries_50_title,os.path.join(plot_output_dir,plot_timeseries_50_outname)) + + + """ + Write out matched pair files + """ + # 30mb and 50mb wind + dlength = 2 + modname = os.environ.get('MODEL_NAME','GFS') + maskname = 'FULL' + lead_hr = np.floor(max_leads) + lead_min = np.round(np.remainder(max_leads,1) * 60) + lead_sec = np.round(np.remainder(lead_min,1) * 60) + datetimeindex = dsO_3050.indexes[obs_timevar] + for i in range(len(datetimeindex)): + valid_str = datetimeindex[i].strftime('%Y%m%d_%H%M%S') + dsO_cur = dsO_3050.sel(time=datetimeindex[i]) + dsF_cur = dsF_3050.sel(time=datetimeindex[i]) + qbo_phase = np.where(dsO_cur < 0, 'East','West') + obs_lvls = ['P'+str(int(op)) for op in dsO_cur.pres] + obs_lvls2 = [str(int(op)) for op in dsO_cur.pres] + fcst_lvls = ['P'+str(int(fp)) for fp in dsF_cur.pres] + leadstr = str(int(lead_hr[i])).zfill(2)+str(int(lead_min[i])).zfill(2)+str(int(lead_sec[i])).zfill(2) + write_mpr_file(dsF_cur.values,dsO_cur.values,[0.0]*dlength,[0.0]*dlength, + [leadstr]*dlength,[valid_str]*dlength,['000000']*dlength,[valid_str]*dlength,modname,qbo_phase, + ['QBO_U']*dlength,['m/s']*dlength,fcst_lvls,['QBO_U']*dlength,['m/s']*dlength,obs_lvls, + maskname,obs_lvls2,mpr_output_dir,'qbo_u_'+modname) + + + +if __name__ == '__main__': + main()