Note
Go to the end to download the full example code.
UserScript and StatAnalysis: Compute QBO Phase Plots and QBO Index
model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO.py
Scientific Objective
Many common modes of variability in the troposphere have stratospheric teleconnection pathways. This use case performs evaluation of the Quasi-biennial Oscillation (QBO), one of the key players of stratospheric 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 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.
Version Added
METplus version 6.0
Datasets
Forecast: GFS 0-24 hour forecasts for 10/2017 - 2/2018
Observation: ERA 30 year climatology for EOFs and 10/2017 - 2/2018
EOFs: EOFs computed for 10/1991 - 12/2020
Location: 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 and can be downloaded at https://dtcenter.ucar.edu/dfiles/code/METplus/METplus_Data/v6.0/.
METplus Components
This use case calls UserScript once and StatAnalysis once. Additionally, METcalcpy, METplotpy, and METdataio are required to run. The METcalcpy scripts accessed include the following:
metcalcpy/pre_processing/directional_means.py
metcalcpy/util/write_mpr.py
The METplotpy scripts accessed include the following:
metplotpy/contributed/stratosphere_diagnostics/stratosphere_plots.py
The METdataio scripts accessed include the following:
METreadnc/util/read_netcdf.py
METplus Workflow
Beginning time (INIT_BEG): 2018-02-01
End time (INIT_END): 2018-02-28
Increment between beginning and end times (INIT_INCREMENT): 30 days
Sequence of forecast leads to process (LEAD_SEQ): 0, 3, 6, 9, 12, 15, 18, 21 hours
This use case does not loop but UserScript and StatAnalysis are each run once using data from the entire time period. The UserScript call runs stratosphere_qbo_driver.py to compute a QBO index. The output QBO index is run through Stat-Analysis to compute the continuous statistics.
There is an optional step to calculate EOFs rather than providing an input file. This calculation is performed inside qbo_driver.py, but requires building a list of input files to use to create EOFs using a call to UserScript. This step is turned off due to the length of time it needs to run, but can be turned back on using the PROCESS_LIST that is commented out.
PROCESS_LIST = UserScript(create_obs_eofs_filelist), UserScript(compute_plot_qbo), StatAnalysis(qbo_bias)
Other changes needed to calculate EOFs are described in the User Scripting section below. Data is not provided for the EOF calculation.
METcalcpy 3.0.0 or higher, METplotpy 3.0.0 or higher, and METdataio 2.1 or higher are needed for this use case.
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
[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
MET Configuration
METplus sets environment variables based on user settings in the METplus configuration file. See 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: Overriding Unsupported MET config file settings
StatAnalysisConfig_wrapped
////////////////////////////////////////////////////////////////////////////////
//
// STAT-Analysis configuration file.
//
// For additional information, see the MET_BASE/config/README file.
//
////////////////////////////////////////////////////////////////////////////////
//
// Filtering input STAT lines by the contents of each column
//
//model = [
${METPLUS_MODEL}
//desc = [
${METPLUS_DESC}
//fcst_lead = [
${METPLUS_FCST_LEAD}
//obs_lead = [
${METPLUS_OBS_LEAD}
//fcst_valid_beg =
${METPLUS_FCST_VALID_BEG}
//fcst_valid_end =
${METPLUS_FCST_VALID_END}
fcst_valid_inc = [];
fcst_valid_exc = [];
//fcst_valid_hour = [
${METPLUS_FCST_VALID_HOUR}
//obs_valid_beg =
${METPLUS_OBS_VALID_BEG}
//obs_valid_end =
${METPLUS_OBS_VALID_END}
obs_valid_inc = [];
obs_valid_exc = [];
//obs_valid_hour = [
${METPLUS_OBS_VALID_HOUR}
//fcst_init_beg =
${METPLUS_FCST_INIT_BEG}
//fcst_init_end =
${METPLUS_FCST_INIT_END}
fcst_init_inc = [];
fcst_init_exc = [];
//fcst_init_hour = [
${METPLUS_FCST_INIT_HOUR}
//obs_init_beg =
${METPLUS_OBS_INIT_BEG}
//obs_init_end =
${METPLUS_OBS_INIT_END}
obs_init_inc = [];
obs_init_exc = [];
//obs_init_hour = [
${METPLUS_OBS_INIT_HOUR}
//fcst_var = [
${METPLUS_FCST_VAR}
//obs_var = [
${METPLUS_OBS_VAR}
//fcst_units = [
${METPLUS_FCST_UNITS}
//obs_units = [
${METPLUS_OBS_UNITS}
//fcst_lev = [
${METPLUS_FCST_LEVEL}
//obs_lev = [
${METPLUS_OBS_LEVEL}
//obtype = [
${METPLUS_OBTYPE}
//vx_mask = [
${METPLUS_VX_MASK}
//interp_mthd = [
${METPLUS_INTERP_MTHD}
//interp_pnts = [
${METPLUS_INTERP_PNTS}
//fcst_thresh = [
${METPLUS_FCST_THRESH}
//obs_thresh = [
${METPLUS_OBS_THRESH}
//cov_thresh = [
${METPLUS_COV_THRESH}
//alpha = [
${METPLUS_ALPHA}
//line_type = [
${METPLUS_LINE_TYPE}
column = [];
weight = [];
////////////////////////////////////////////////////////////////////////////////
//
// Array of STAT-Analysis jobs to be performed on the filtered data
//
//jobs = [
${METPLUS_JOBS}
////////////////////////////////////////////////////////////////////////////////
//
// Confidence interval settings
//
out_alpha = 0.05;
boot = {
interval = PCTILE;
rep_prop = 1.0;
n_rep = 0;
rng = "mt19937";
seed = "";
}
////////////////////////////////////////////////////////////////////////////////
//
// WMO mean computation logic
//
wmo_sqrt_stats = [ "CNT:FSTDEV", "CNT:OSTDEV", "CNT:ESTDEV",
"CNT:RMSE", "CNT:RMSFA", "CNT:RMSOA",
"VCNT:FS_RMS", "VCNT:OS_RMS", "VCNT:RMSVE",
"VCNT:FSTDEV", "VCNT:OSTDEV" ];
wmo_fisher_stats = [ "CNT:PR_CORR", "CNT:SP_CORR",
"CNT:KT_CORR", "CNT:ANOM_CORR" ];
////////////////////////////////////////////////////////////////////////////////
//hss_ec_value =
${METPLUS_HSS_EC_VALUE}
rank_corr_flag = FALSE;
vif_flag = FALSE;
tmp_dir = "${MET_TMP_DIR}";
//version = "V10.0";
${METPLUS_MET_CONFIG_OVERRIDES}
Python Embedding
This use case does not use python embedding
User Scripting
This use case runs the stratosphere_qbo_driver.py Python script, 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 in MET’s matched pair (MPR) format using write_mpr.py in METcalcpy. They 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).
Variables input to this script are given in the [user_env_vars] section of the configuration file. As mentioned above, the option exists to compute EOFs inside the script. To do this, the COMPUTE_EOF_ZONAL_MERIDIONAL_MEAN variable should be set to true, and an input file list provided by turning on the UserScipt call for [create_eof_filelist].
parm/use_cases/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO/stratosphere_qbo_driver.py
#!/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()
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 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.
Refer to the value set for OUTPUT_BASE to find where the output data was generated. Output for this use case will be found in {OUTPUT_BASE}/model_applications/s2s_stratosphere/UserScript_fcstGFS_obsERA_StratosphereQBO There should be 4 graphics output to the plot directory:
* 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
Additionally matched pair .stat files for each day will be output to the mpr directory. These files have the following format:
* qbo_u_GFS_210000L_YYYYMMDD_000000V.stat
Two computed continuous statistics will be output to the StatAnalysis directory:
* GFS_ERA_20171001_20180228_210000L_zonal_wind_byphase_CNT.stat
* GFS_ERA_20171001_20180228_210000L_zonal_wind_CNT.stat
Keywords
Note
S2SAppUseCase
S2SStratosphereAppUseCase
UserScriptUseCase
StatAnalysisUseCase
METdataioUseCase
METcalcpyUseCase
METplotpyUseCase
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.
sphinx_gallery_thumbnail_path = ‘_static/s2s_stratosphere-UserScript_fcstGFS_obsERA_StratosphereQBO.png’