UserScript: Make ERA RMM plots from calculated MJO indices

model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM.py

Scientific Objective

The Madden-Julian Oscillation (MJO) is the largest element of intraseasonal variability in the tropics and is characterized by eastward moving regions of enhanced and suppressed rainfall. These phases are typically grouped into numbers 1 - 8 based on the geographic location of the enhanced and suppressed rainfall. The MJO affects global weather including summer monsoons, tropical cyclone development, and sudden stratospheric warming events, and has teleconnections to mid latitude weather systems.

This use case uses anomalies of outgoing longwave radiation (OLR), 850 hPa wind (U850), and 200 hPa wind (U200) to compute the Real-time Multivariate MJO Index (RMM). In contrast to OMI, which is a convective index of MJO, RMM is a dynamical index. The code for computing RMM came from Maria Gehne at the NOAA Physical Sciences Laboratory (PSL).

Version Added

METplus version 4.1.0

Datasets

Forecast: None

Observation: ERA Reanlaysis Outgoing Longwave Radiation, 850 hPa wind and 200 hPa wind, 2000 - 2002

Climatology:

EOFs: EOF patterns for OLR, U850, and U200 from Matthew Wheeler

Location: All of the input data required for this use case can be found in a sample data tarball. Each use case category will have one or more sample data tarballs. It is only necessary to download the tarball with the use case’s dataset and not the entire collection of sample data. Click here to access the METplus releases page and download sample data for the appropriate release: https://github.com/dtcenter/METplus/releases This tarball should be unpacked into the directory that you will set the value of INPUT_BASE. See Running METplus section for more information.

METplus Components

This use case calls UserScript 5 times and Regrid-Data-Plane three times. The UserScript calls run the harmonic pre-processing on all three variables, create a filelist, and run the RMM calculation. There are four optional pre-processing steps.

This use case requires METcalcpy, METplotpy, and METdataio to run. The METcalcpy scripts accessed include the following:

• metcalcpy/contributed/rmm_omi/compute_mjo_indices.py

The METplotpy scripts accessed include the following:

• metplotpy/contributed/mjo_rmm_omi/plot_mjo_indioces.py

The METdataio scripts accessed include the following:

• METreadnc/util/read_netcdf.py

METplus Workflow

Beginning time (VALID_BEG): 01-01-2000

End time (VALID_END): 12-30-2002

Increment between beginning and end times (INIT_INCREMENT): 1 day

Sequence of forecast leads to process (LEAD_SEQ): 0 hour

The UserScript calls do not loop, but are each run once. The first call creates a list of the mean daily annual data files for OLR, U850, and U200. It is done separately since the mean daily annual files span across all years whereas the RMM calculation can proceed on a different time frame. The second, third, and fourth calls to UserScript run the pre-processing to compute anomalies OLR, U850, and U200 using a harmonic analysis program in python. Then, RegridDataPlane is run for all valid times for the 3 variables. This step cuts the grid to only include -15 to 15 latitude. Finally, the last (fifth) call to UserScript runs the RMM calculation once on the observations.

There are four optional pre-processing steps loop by valid time with different timing settings needed for the different steps. These steps are turned off due to data size and processing time. Two of the steps are calls to Pcp-Combine to compute the mean daily annual data for OLR, wind (U850 and U200). The other two steps also call Pcp-Combine but these compute daily means for OLR and wind. These omitted steps can be turned back on by using the PROCESS_LIST that is commented out in the file:

PROCESS_LIST = PcpCombine(mean_daily_annual_cycle_obs_wind), PcpCombine(mean_daily_annual_cycle_obs_olr), PcpCombine(daily_mean_obs_wind), PcpCombine(daily_mean_obs_olr), UserScript(create_mda_filelist), UserScript(harmonic_anomalies_olr), UserScript(harmonic_anomalies_u850), UserScript(harmonic_anomalies_u200), RegridDataPlane(regrid_obs_olr), RegridDataPlane(regrid_obs_u850), RegridDataPlane(regrid_obs_u200), UserScript(script_rmm)

Settings for the optional pre-processing steps can be found in the respective sections of the configuration, mean_daily_annual_cycle_obs_wind, mean_daily_annual_cycle_obs_olr, daily_mean_obs_wind, and daily_mean_obs_olr. Data is not provided in the tarball to run these steps, but the configurations is provided for reference on how to set up these calculations.

The Phase diagram and timeseries plots are created over a different time frame, 01-01-2002 to 12-30-2002 for both plots.

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_mjo/UserScript_obsERA_obsOnly_RMM.conf.

[config]

# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM.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
###

# All steps, including creating daily means and mean daily annual cycle
#PROCESS_LIST = PcpCombine(mean_daily_annual_cycle_obs_wind), PcpCombine(mean_daily_annual_cycle_obs_olr), PcpCombine(daily_mean_obs_wind), PcpCombine(daily_mean_obs_olr), UserScript(create_mda_filelist), UserScript(harmonic_anomalies_olr), UserScript(harmonic_anomalies_u850), UserScript(harmonic_anomalies_u200), RegridDataPlane(regrid_obs_olr), RegridDataPlane(regrid_obs_u850), RegridDataPlane(regrid_obs_u200), UserScript(script_rmm)
# Computing anomalies, regridding, and RMM Analysis script

PROCESS_LIST = UserScript(create_mda_filelist), UserScript(harmonic_anomalies_olr), UserScript(harmonic_anomalies_u850),  UserScript(harmonic_anomalies_u200), RegridDataPlane(regrid_obs_olr), RegridDataPlane(regrid_obs_u850), RegridDataPlane(regrid_obs_u200), UserScript(script_rmm)


###
# 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 = VALID
VALID_TIME_FMT = %Y%m%d%H
VALID_BEG = 2000010100
VALID_END = 2002123000
VALID_INCREMENT = 86400

LEAD_SEQ = 0

# variables referenced in other sections

# Run the obs for these cases
OBS_RUN = True
FCST_RUN = False


###
# RegridDataPlane Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#regriddataplane
###

REGRID_DATA_PLANE_VERIF_GRID = latlon 144 13 -15 0 2.5 2.5
REGRID_DATA_PLANE_METHOD = NEAREST
REGRID_DATA_PLANE_WIDTH = 1


###
# PCPCombine(mean_daily_annual_cycle_obs_wind) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pcpcombine
###

# Configurations for creating U200 and U850 mean daily annual cycle obs
# Mean daily annual cycle anomalies are computed for 1979 - 2001
[mean_daily_annual_cycle_obs_wind]

LOOP_BY = VALID
VALID_TIME_FMT = %Y%m%d%H
VALID_BEG = 2012010100
VALID_END = 2012123100
VALID_INCREMENT = 86400

OBS_PCP_COMBINE_RUN = {OBS_RUN}

OBS_PCP_COMBINE_METHOD = USER_DEFINED

OBS_PCP_COMBINE_COMMAND = -derive mean {OBS_PCP_COMBINE_INPUT_DIR}/{OBS_PCP_COMBINE_INPUT_TEMPLATE} -field 'name="U_P850_mean"; level="(*,*)"; set_attr_valid = "{valid?fmt=%Y%m%d_%H%M%S}";' -field 'name="U_P200_mean"; level="(*,*)"; set_attr_valid = "{valid?fmt=%Y%m%d_%H%M%S}";' -name U_P850_mean,U_P200_mean

OBS_PCP_COMBINE_INPUT_DIR = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean
OBS_PCP_COMBINE_INPUT_TEMPLATE = ERA_wind_daily_mean_*{valid?fmt=%m%d}.nc

OBS_PCP_COMBINE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/mean_daily_annual_cycle
OBS_PCP_COMBINE_OUTPUT_TEMPLATE = ERA_wind_daily_annual_{valid?fmt=%m%d}.nc


###
# PCPCombine(mean_daily_annual_cycle_obs_olr) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pcpcombine
###

# Configurations for creating OLR mean daily annual cycle obs
# Mean daily annual cycle anomalies are computed for 1979 - 2001
[mean_daily_annual_cycle_obs_olr]

LOOP_BY = VALID
VALID_TIME_FMT = %Y%m%d%H
VALID_BEG = 2012010100
VALID_END = 2012123100
VALID_INCREMENT = 86400

OBS_PCP_COMBINE_RUN = {OBS_RUN}

OBS_PCP_COMBINE_METHOD = USER_DEFINED

OBS_PCP_COMBINE_COMMAND = -derive mean {OBS_PCP_COMBINE_INPUT_DIR}/{OBS_PCP_COMBINE_INPUT_TEMPLATE} -field 'name="olr"; level="(*,*)";'

OBS_PCP_COMBINE_INPUT_DIR = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean
OBS_PCP_COMBINE_INPUT_TEMPLATE = ERA_OLR_daily_mean_*{valid?fmt=%m%d}.nc

OBS_PCP_COMBINE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/mean_daily_annual_cycle
OBS_PCP_COMBINE_OUTPUT_TEMPLATE = ERA_OLR_daily_annual_{valid?fmt=%m%d}.nc


###
# PCPCombine(daily_mean_obs_wind) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pcpcombine
###

# Configurations for creating U200 and U850 daily mean obs
[daily_mean_obs_wind]

LOOP_BY = VALID
VALID_TIME_FMT = %Y%m%d%H
VALID_BEG = 1979010100
VALID_END = 2002123100
VALID_INCREMENT = 86400

OBS_PCP_COMBINE_RUN = {OBS_RUN}

OBS_PCP_COMBINE_METHOD = USER_DEFINED

OBS_PCP_COMBINE_COMMAND = -derive mean {OBS_PCP_COMBINE_INPUT_DIR}/{OBS_PCP_COMBINE_INPUT_TEMPLATE} -field 'name="U"; level="P850"; set_attr_valid = "{valid?fmt=%Y%m%d_%H%M%S}";' -field 'name="U"; level="P200"; set_attr_valid = "{valid?fmt=%Y%m%d_%H%M%S}";'

OBS_PCP_COMBINE_INPUT_DIR = /gpfs/fs1/collections/rda/data/ds627.0/ei.oper.an.pl
OBS_PCP_COMBINE_INPUT_TEMPLATE = {valid?fmt=%Y%m}/ei.oper.an.pl.regn128uv.{valid?fmt=%Y%m%d}*

OBS_PCP_COMBINE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean
OBS_PCP_COMBINE_OUTPUT_TEMPLATE = ERA_wind_daily_mean_{valid?fmt=%Y%m%d}.nc


###
# PCPCombine(daily_mean_obs_olr) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pcpcombine
###

# Configurations for creating mean daily annual cycle obs OLR
[daily_mean_obs_olr]

LOOP_BY = VALID
VALID_TIME_FMT = %Y%m%d%H
VALID_BEG = 1979010100
VALID_END = 2002123100
VALID_INCREMENT = 86400

OBS_PCP_COMBINE_RUN = {OBS_RUN}

OBS_PCP_COMBINE_METHOD = USER_DEFINED

OBS_PCP_COMBINE_COMMAND = -add {OBS_PCP_COMBINE_INPUT_DIR}/{OBS_PCP_COMBINE_INPUT_TEMPLATE} -field 'name="olr"; level="({valid?fmt=%Y%m%d_%H%M%S},*,*)"; file_type=NETCDF_NCCF;'

OBS_PCP_COMBINE_INPUT_DIR = /glade/u/home/kalb/MJO
OBS_PCP_COMBINE_INPUT_TEMPLATE = olr.1x.7920.nc

OBS_PCP_COMBINE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean
OBS_PCP_COMBINE_OUTPUT_TEMPLATE = ERA_OLR_daily_mean_{valid?fmt=%Y%m%d}.nc


###
# UserScript(create_mda_filelist) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript
###

# Creating a file list of the mean daily annual cycle files
# This is run separately since it has different start/end times
[create_mda_filelist]

# Find the files for each lead time
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD

# Valid Begin and End Times for the CBL File Climatology
VALID_BEG = 2012010100
VALID_END = 2012123100
VALID_INCREMENT = 86400
LEAD_SEQ = 0

# Template of filenames to input to the user-script
USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/mean_daily_annual_cycle/ERA_OLR_daily_annual_{valid?fmt=%m%d}.nc,{INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/mean_daily_annual_cycle/ERA_wind_daily_annual_{valid?fmt=%m%d}.nc

# Name of the file containing the listing of input files
USER_SCRIPT_INPUT_TEMPLATE_LABELS = input_mean_daily_annual_infiles_olr,input_mean_daily_annual_infiles_wind

# Placeholder command just to build the file list
# This just states that it's building the file list
USER_SCRIPT_COMMAND = echo Populated file list for Mean daily annual cycle Input


###
# UserScript(harmonic_anomalies_olr) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript
###

# Configurations to create anomalies for OLR
[harmonic_anomalies_olr]

# list of strings to loop over for each run time.
# Run the user script once per lead
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD

# Template of filenames to input to the user-script
USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean/ERA_OLR_daily_mean_{valid?fmt=%Y%m%d}.nc

# Name of the file containing the listing of input files
# The options are OBS_OLR_INPUT, OBS_U850_INPUT, OBS_U200_INPUT, FCST_OLR_INPUT, FCST_U850_INPUT, and FCST_U200_INPUT
# *** Make sure the order is the same as the order of templates listed in USER_SCRIPT_INPUT_TEMPLATE
USER_SCRIPT_INPUT_TEMPLATE_LABELS = input_daily_mean_infiles

# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/compute_harmonic_anomalies.py 'METPLUS_FILELIST_INPUT_MEAN_DAILY_ANNUAL_INFILES_OLR' 'olr' 'olr_NA_mean' '{OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Anomaly' 'ERA_OLR_anom'


###
# UserScript(harmonic_anomalies_u850) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript
###

# Configurations to create anomalies for U850
[harmonic_anomalies_u850]

# list of strings to loop over for each run time.
# Run the user script once per lead
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD

# Template of filenames to input to the user-script
USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean/ERA_wind_daily_mean_{valid?fmt=%Y%m%d}.nc

# Name of the file containing the listing of input files
# The options are OBS_OLR_INPUT, OBS_U850_INPUT, OBS_U200_INPUT, FCST_OLR_INPUT, FCST_U850_INPUT, and FCST_U200_INPUT
# *** Make sure the order is the same as the order of templates listed in USER_SCRIPT_INPUT_TEMPLATE
USER_SCRIPT_INPUT_TEMPLATE_LABELS = input_daily_mean_infiles

# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/compute_harmonic_anomalies.py 'METPLUS_FILELIST_INPUT_MEAN_DAILY_ANNUAL_INFILES_WIND' 'U_P850_mean' 'U_P850_mean' '{OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Anomaly' 'ERA_U850_anom'


###
# UserScript(harmonic_anomalies_u200) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript
###

# Configurations to create anomalies for U200
[harmonic_anomalies_u200]

# list of strings to loop over for each run time.
# Run the user script once per lead
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD

# Template of filenames to input to the user-script
USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/daily_mean/ERA_wind_daily_mean_{valid?fmt=%Y%m%d}.nc

# Name of the file containing the listing of input files
# The options are OBS_OLR_INPUT, OBS_U850_INPUT, OBS_U200_INPUT, FCST_OLR_INPUT, FCST_U850_INPUT, and FCST_U200_INPUT
# *** Make sure the order is the same as the order of templates listed in USER_SCRIPT_INPUT_TEMPLATE
USER_SCRIPT_INPUT_TEMPLATE_LABELS = input_daily_mean_infiles

# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/compute_harmonic_anomalies.py 'METPLUS_FILELIST_INPUT_MEAN_DAILY_ANNUAL_INFILES_WIND' 'U_P200_mean' 'U_P200_mean' '{OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Anomaly' 'ERA_U200_anom'


###
# RegridDataPlane(regrid_obs_olr) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#regriddataplane
###

# Configurations for regrid_data_plane: Regrid OLR to -15 to 15 latitude
[regrid_obs_olr]

OBS_REGRID_DATA_PLANE_RUN = {OBS_RUN}

REGRID_DATA_PLANE_ONCE_PER_FIELD = False

OBS_REGRID_DATA_PLANE_VAR1_NAME = olr_anom
OBS_REGRID_DATA_PLANE_VAR1_LEVELS = "(*,*)"
OBS_REGRID_DATA_PLANE_VAR1_OUTPUT_FIELD_NAME = OLR_anom

OBS_REGRID_DATA_PLANE_INPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Anomaly
OBS_REGRID_DATA_PLANE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Regrid

OBS_REGRID_DATA_PLANE_INPUT_TEMPLATE = ERA_OLR_anom_{lead?fmt=%H%M%S}L_{valid?fmt=%Y%m%d}_{valid?fmt=%H%M%S}V.nc
OBS_REGRID_DATA_PLANE_OUTPUT_TEMPLATE = ERA_OLR_{valid?fmt=%Y%m%d}.nc


###
# RegridDataPlane(regrid_obs_u850) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#regriddataplane
###

# Configurations for regrid_data_plane: Regrid u850 to -15 to 15 latitude
[regrid_obs_u850]

OBS_REGRID_DATA_PLANE_RUN = {OBS_RUN}

REGRID_DATA_PLANE_ONCE_PER_FIELD = False

OBS_REGRID_DATA_PLANE_VAR1_NAME = U_P850_mean_anom
OBS_REGRID_DATA_PLANE_VAR1_LEVELS = "(*,*)"
OBS_REGRID_DATA_PLANE_VAR1_OUTPUT_FIELD_NAME = U_P850_anom

OBS_REGRID_DATA_PLANE_INPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Anomaly
OBS_REGRID_DATA_PLANE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Regrid

OBS_REGRID_DATA_PLANE_INPUT_TEMPLATE = ERA_U850_anom_{lead?fmt=%H%M%S}L_{valid?fmt=%Y%m%d}_{valid?fmt=%H%M%S}V.nc
OBS_REGRID_DATA_PLANE_OUTPUT_TEMPLATE = ERA_U850_{valid?fmt=%Y%m%d}.nc


###
# RegridDataPlane(regrid_obs_u200) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#regriddataplane
###

# Configurations for regrid_data_plane: Regrid u200 to -15 to 15 latitude
[regrid_obs_u200]

OBS_REGRID_DATA_PLANE_RUN = {OBS_RUN}

REGRID_DATA_PLANE_ONCE_PER_FIELD = False

OBS_REGRID_DATA_PLANE_VAR1_NAME = U_P200_mean_anom
OBS_REGRID_DATA_PLANE_VAR1_LEVELS = "(*,*)"
OBS_REGRID_DATA_PLANE_VAR1_OUTPUT_FIELD_NAME = U_P200_anom

OBS_REGRID_DATA_PLANE_INPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Anomaly
OBS_REGRID_DATA_PLANE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Regrid

OBS_REGRID_DATA_PLANE_INPUT_TEMPLATE = ERA_U200_anom_{lead?fmt=%H%M%S}L_{valid?fmt=%Y%m%d}_{valid?fmt=%H%M%S}V.nc
OBS_REGRID_DATA_PLANE_OUTPUT_TEMPLATE = ERA_U200_{valid?fmt=%Y%m%d}.nc


# Configurations for the RMM analysis script
[user_env_vars]
# Whether to Run the model or obs
RUN_OBS = {OBS_RUN}
RUN_FCST = {FCST_RUN}

# Make OUTPUT_BASE Available to the script
SCRIPT_OUTPUT_BASE = {OUTPUT_BASE}

# Number of obs per day
OBS_PER_DAY = 1

# Variable names for OLR, U850, U200
OBS_OLR_VAR_NAME = OLR_anom
OBS_U850_VAR_NAME = U_P850_anom
OBS_U200_VAR_NAME = U_P200_anom

# EOF Filename
OLR_EOF_INPUT_TEXTFILE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/EOF/rmm_olr_eofs.txt
U850_EOF_INPUT_TEXTFILE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/EOF/rmm_u850_eofs.txt
U200_EOF_INPUT_TEXTFILE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/EOF/rmm_u200_eofs.txt

# Normalization factors for RMM
RMM_OLR_NORM = 15.11623
RMM_U850_NORM = 1.81355
RMM_U200_NORM = 4.80978
PC1_NORM = 8.618352504159244
PC2_NORM = 8.40736449709697

# Output Directory for the plots
# If not set, it this will default to {OUTPUT_BASE}/plots
RMM_PLOT_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/plots

# EOF plot information
EOF_PLOT_OUTPUT_NAME = RMM_EOFs
EOF_PLOT_OUTPUT_FORMAT = png

# Phase Plot start date, end date, output name, and format
PHASE_PLOT_TIME_BEG = 2002010100
PHASE_PLOT_TIME_END = 2002123000
PHASE_PLOT_TIME_FMT = {VALID_TIME_FMT}
OBS_PHASE_PLOT_OUTPUT_NAME = obs_RMM_comp_phase
OBS_PHASE_PLOT_OUTPUT_FORMAT = png

# Time Series Plot start date, end date, output name, and format
TIMESERIES_PLOT_TIME_BEG = 2002010100
TIMESERIES_PLOT_TIME_END = 2002123000
TIMESERIES_PLOT_TIME_FMT = {VALID_TIME_FMT}
OBS_TIMESERIES_PLOT_OUTPUT_NAME = obs_RMM_time_series
OBS_TIMESERIES_PLOT_OUTPUT_FORMAT = png


###
# UserScript(script_rmm) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#userscript
###

# Configurations for UserScript: Run the RMM Analysis driver
[script_rmm]
# list of strings to loop over for each run time.
# Run the user script once per lead
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD

# Template of filenames to input to the user-script
USER_SCRIPT_INPUT_TEMPLATE = {OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Regrid/ERA_OLR_{valid?fmt=%Y%m%d}.nc,{OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Regrid/ERA_U850_{valid?fmt=%Y%m%d}.nc,{OUTPUT_BASE}/s2s_mjo/UserScript_obsERA_obsOnly_RMM/ERA/Regrid/ERA_U200_{valid?fmt=%Y%m%d}.nc

# Name of the file containing the listing of input files
# The options are OBS_OLR_INPUT, OBS_U850_INPUT, OBS_U200_INPUT, FCST_OLR_INPUT, FCST_U850_INPUT, and FCST_U200_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_OLR_INPUT,OBS_U850_INPUT, OBS_U200_INPUT

# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/RMM_driver.py

MET Configuration

There are no MET configuration files used in this use case.

Python Embedding

This use case does not use Python embedding.

User Scripting

There are two Python scripts used in this use case. The RMM driver script (RMM_driver.py) orchestrates the calculation of the MJO indices and the generation of three RMM plots. The calculation proceeds using OLR, U850, and U200 data between 15N and 15S. The 120 day mean is first removed and the data are normalized by normalization factors (generally the square root of the average variance). The anomalies are projected onto Empirical Orthogonal Function (EOF) data. The OLR is then filtered for 20 - 96 days, and regressed onto the daily EOFs. Finally, it’s normalized and these normalized components are plotted on a phase diagram and timeseries plot.

The anomalies are created using a harmonic analysis for OLR, U850, and U200 with the Python script compute_harmonic_anomalies.py. Input to the harmonic analysis script include a text file containing the list of input files, the daily mean variable name, mean daily average varable name, output directory and output file basename. These are defined in the [harmonic_anomalies_olr], [harmonic_anomalies_u850], and harmonic_anomalies_u200] sections of the configuration file. Additional variables for the RMM calculation and harmonic analysis are set in the [user_env_vars] section of the .conf file.

parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/RMM_driver.py

#!/usr/bin/env python3
"""
Driver Script to Compute RMM index from input U850, U200 and OLR data. Data is averaged from 20S-20N
"""

import numpy as np
import xarray as xr
import pandas as pd
import datetime
import glob
import os
import sys
import warnings

import metcalcpy.contributed.rmm_omi.compute_mjo_indices as cmi
import metplotpy.contributed.mjo_rmm_omi.plot_mjo_indices as pmi
import METreadnc.util.read_netcdf as read_netcdf


def read_rmm_eofs(olrfile, u850file, u200file):
    """
    Read the OMI EOFs from file and into a xarray DataArray.
    :param eofpath: filepath to the location of the eof files
    :return: EOF1 and EOF2 2D DataArrays
    """

    # observed EOFs from BOM Australia are saved in individual text files for each variable
    # horizontal resolution of EOFs is 2.5 degree and longitudes go from 0 - 375.5, column1 is eof1
    # column 2 is eof2 in each file
    EOF1 = xr.DataArray(np.empty([3,144]),dims=['var','lon'],
    coords={'var':['olr','u850','u200'], 'lon':np.arange(0,360,2.5)})
    EOF2 = xr.DataArray(np.empty([3,144]),dims=['var','lon'],
    coords={'var':['olr','u850','u200'], 'lon':np.arange(0,360,2.5)})
    nlon = len(EOF1['lon'])

    tmp = pd.read_csv(olrfile, header=None, sep=r'\s+', names=['eof1','eof2'])
    EOF1[0,:] = tmp.eof1.values
    EOF2[0,:] = tmp.eof2.values
    tmp = pd.read_csv(u850file, header=None, sep=r'\s+', names=['eof1','eof2'])
    EOF1[1,:] = tmp.eof1.values
    EOF2[1,:] = tmp.eof2.values
    tmp = pd.read_csv(u200file, header=None, sep=r'\s+', names=['eof1','eof2'])
    EOF1[2,:] = tmp.eof1.values
    EOF2[2,:] = tmp.eof2.values

    return EOF1, EOF2


def run_rmm_steps(inlabel, spd, EOF1, EOF2, oplot_dir):

    # Get OLR, U850, U200 file listings and variable names
    olr_filetxt = os.environ['METPLUS_FILELIST_'+inlabel+'_OLR_INPUT']
    u850_filetxt = os.environ['METPLUS_FILELIST_'+inlabel+'_U850_INPUT']
    u200_filetxt = os.environ['METPLUS_FILELIST_'+inlabel+'_U200_INPUT']

    olr_var = os.environ[inlabel+'_OLR_VAR_NAME']
    u850_var = os.environ[inlabel+'_U850_VAR_NAME']
    u200_var = os.environ[inlabel+'_U200_VAR_NAME']

    # Read the listing of OLR, U850, U200 files
    with open(olr_filetxt) as ol:
        olr_input_files = ol.read().splitlines()
    if (olr_input_files[0] == 'file_list'):
            olr_input_files = olr_input_files[1:]
    with open(u850_filetxt) as u8:
        u850_input_files = u8.read().splitlines()
    if (u850_input_files[0] == 'file_list'):
            u850_input_files = u850_input_files[1:]
    with open(u200_filetxt) as u2:
        u200_input_files = u2.read().splitlines()
    if (u200_input_files[0] == 'file_list'):
            u200_input_files = u200_input_files[1:]

    # Check the input data to make sure it's not all missing
    olr_allmissing = all(elem == 'missing' for elem in olr_input_files)
    if olr_allmissing:
        raise IOError ('No input OLR files were found, check file paths')
    u850_allmissing = all(elem == 'missing' for elem in u850_input_files)
    if u850_allmissing:
        raise IOError('No input U850 files were found, check file paths')
    u200_allmissing = all(elem == 'missing' for elem in u200_input_files)
    if u200_allmissing:
        raise IOError('No input U200 files were found, check file paths')
    

    # Read OLR, U850, U200 data from file
    netcdf_reader_olr = read_netcdf.ReadNetCDF()
    ds_olr = netcdf_reader_olr.read_into_xarray(olr_input_files)

    netcdf_reader_u850 = read_netcdf.ReadNetCDF()
    ds_u850 = netcdf_reader_u850.read_into_xarray(u850_input_files)

    netcdf_reader_u200 = read_netcdf.ReadNetCDF()
    ds_u200 = netcdf_reader_u200.read_into_xarray(u200_input_files)


    time = []
    for din in range(len(ds_olr)):
        colr = ds_olr[din]
        ctime =  datetime.datetime.strptime(colr[olr_var].valid_time,'%Y%m%d_%H%M%S')
        time.append(ctime.strftime('%Y-%m-%d'))
        colr = colr.assign_coords(time=ctime)
        ds_olr[din] = colr.expand_dims("time")

        cu850 = ds_u850[din]
        cu850 = cu850.assign_coords(time=ctime)
        ds_u850[din] = cu850.expand_dims("time")

        cu200 = ds_u200[din]
        cu200 = cu200.assign_coords(time=ctime)
        ds_u200[din] = cu200.expand_dims("time")

    time = np.array(time,dtype='datetime64[D]')

    everything_olr = xr.concat(ds_olr,"time")
    olr = everything_olr[olr_var]
    olr = olr.mean('lat')
    print(olr.min(), olr.max())

    everything_u850 = xr.concat(ds_u850,"time")
    u850 = everything_u850[u850_var]
    u850 = u850.mean('lat')
    print(u850.min(), u850.max())

    everything_u200 = xr.concat(ds_u200,"time")
    u200 = everything_u200[u200_var]
    u200 = u200.mean('lat')
    print(u200.min(), u200.max())
    

    # Get normalization factors for use 
    rmm_norm = [float(os.environ['RMM_OLR_NORM']),float(os.environ['RMM_U850_NORM']),float(os.environ['RMM_U200_NORM'])]
    pc_norm = [float(os.environ['PC1_NORM']),float(os.environ['PC2_NORM'])]

    # project data onto EOFs
    PC1, PC2 = cmi.rmm(olr, u850, u200, time, spd, EOF1, EOF2, rmm_norm, pc_norm)
    print(PC1.min(), PC1.max())

    # Get times for the PC phase diagram
    plase_plot_time_format = os.environ['PHASE_PLOT_TIME_FMT']
    phase_plot_start_time = datetime.datetime.strptime(os.environ['PHASE_PLOT_TIME_BEG'],plase_plot_time_format)
    phase_plot_end_time = datetime.datetime.strptime(os.environ['PHASE_PLOT_TIME_END'],plase_plot_time_format)
    PC1_pcplot = PC1.sel(time=slice(phase_plot_start_time,phase_plot_end_time))
    PC2_pcplot = PC2.sel(time=slice(phase_plot_start_time,phase_plot_end_time))
    pc_ntim_plot = len(PC1_pcplot)
    PC1_pcplot = PC1_pcplot[0:pc_ntim_plot]
    PC2_pcplot = PC2_pcplot[0:pc_ntim_plot]
    
    # Get the output name and format for the PC plase diagram
    phase_plot_name = os.path.join(oplot_dir,os.environ.get(inlabel+'_PHASE_PLOT_OUTPUT_NAME',inlabel+'_RMM_comp_phase'))
    phase_plot_format = os.environ.get(inlabel+'_PHASE_PLOT_OUTPUT_FORMAT','png')

    # plot the PC phase diagram
    pmi.phase_diagram('RMM',PC1_pcplot,PC2_pcplot,np.array(PC1_pcplot['time'].dt.strftime("%Y-%m-%d").values),
        np.ndarray.tolist(PC1_pcplot['time.month'].values),np.ndarray.tolist(PC1_pcplot['time.day'].values),
        phase_plot_name, phase_plot_format)

    # Get times for the PC time series plot
    timeseries_plot_time_format = os.environ['TIMESERIES_PLOT_TIME_FMT']
    timeseries_plot_start_time = datetime.datetime.strptime(os.environ['TIMESERIES_PLOT_TIME_BEG'],plase_plot_time_format)
    timeseries_plot_end_time = datetime.datetime.strptime(os.environ['TIMESERIES_PLOT_TIME_END'],plase_plot_time_format)
    PC1_tsplot = PC1.sel(time=slice(timeseries_plot_start_time,timeseries_plot_end_time))
    PC2_tsplot = PC2.sel(time=slice(timeseries_plot_start_time,timeseries_plot_end_time))
    ts_ntim_plot = len(PC1_tsplot)
    PC1_tsplot = PC1_tsplot[0:ts_ntim_plot]
    PC2_tsplot = PC2_tsplot[0:ts_ntim_plot]

    # Get the ouitput name and format for the PC Timeseries plot
    timeseries_plot_name = os.path.join(oplot_dir,os.environ.get(inlabel+'_TIMESERIES_PLOT_OUTPUT_NAME',
        inlabel+'_RMM_timeseries'))
    timeseries_plot_format = os.environ.get(inlabel+'_TIMESERIES_PLOT_OUTPUT_FORMAT','png')

    # plot PC time series
    pmi.pc_time_series('RMM',PC1_tsplot,PC2_tsplot,np.array(PC1_tsplot['time'].values),
        np.ndarray.tolist(PC1_tsplot['time.month'].values),np.ndarray.tolist(PC1_tsplot['time.day'].values),
        timeseries_plot_name,timeseries_plot_format)


def main():

    # Get the EOF files and EOF plot variables
    olr_eoffile = os.environ['OLR_EOF_INPUT_TEXTFILE']
    u850_eoffile = os.environ['U850_EOF_INPUT_TEXTFILE']
    u200_eoffile = os.environ['U200_EOF_INPUT_TEXTFILE']

    # Get Number of Obs per day
    spd = os.environ.get('OBS_PER_DAY',1)

    # Check for an output plot directory
    oplot_dir = os.environ.get('RMM_PLOT_OUTPUT_DIR','')
    if not oplot_dir:
        obase = os.environ['SCRIPT_OUTPUT_BASE']
        oplot_dir = os.path.join(obase,'plots')
    if not os.path.exists(oplot_dir):
        os.makedirs(oplot_dir)

    # Read in the EOFs and plot and get plot name and format
    EOF1, EOF2 = read_rmm_eofs(olr_eoffile, u850_eoffile, u200_eoffile)
    eof_plot_name = os.path.join(oplot_dir,os.environ.get('EOF_PLOT_OUTPUT_NAME','RMM_EOFs'))
    eof_plot_format = os.environ.get('EOF_PLOT_OUTPUT_FORMAT','png')
    pmi.plot_rmm_eofs(EOF1, EOF2, eof_plot_name, eof_plot_format)

    #  Determine if doing forecast or obs
    run_obs_rmm = os.environ.get('RUN_OBS', 'False').lower()
    run_fcst_rmm = os.environ.get('RUN_FCST', 'False').lower()

    if (run_obs_rmm == 'true'):
        run_rmm_steps('OBS', spd, EOF1, EOF2, oplot_dir)

    if (run_fcst_rmm == 'true'):
        run_rmm_steps('FCST', spd, EOF1, EOF2, oplot_dir)

    # nothing selected
    if (run_obs_rmm == 'false') and (run_fcst_rmm == 'false'):
        warnings.warn('Forecast and Obs runs not selected, nothing will be calculated')
        warnings.warn('Set RUN_FCST or RUN_OBS in the [user_en_vars] section to generate output')


if __name__ == "__main__":
    main()

parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/compute_harmonic_anomalies.py

#!/usr/bin/env python3
import numpy as np
import xarray as xr
import glob
import os
import sys
import datetime
import METreadnc.util.read_netcdf as read_netcdf

input_mean_daily_annual_infiles_list = os.environ[sys.argv[1]]
dm_var = sys.argv[2]
mda_var = sys.argv[3]
anom_output_dir = sys.argv[4]
anom_output_base = sys.argv[5]
input_daily_mean_infiles_list = os.environ['METPLUS_FILELIST_INPUT_DAILY_MEAN_INFILES']

# Environment variables for script
nobs = int(os.environ.get('OBS_PER_DAY',1))
out_var = dm_var+'_anom'

# Read the listing of files
with open(input_daily_mean_infiles_list) as idm:
    input_daily_mean_infiles = idm.read().splitlines()
if (input_daily_mean_infiles[0] == 'file_list'):
    input_daily_mean_infiles = input_daily_mean_infiles[1:]

with open(input_mean_daily_annual_infiles_list) as imda:
    input_mean_daily_annual_infiles = imda.read().splitlines()
if (input_mean_daily_annual_infiles[0] == 'file_list'):
    input_mean_daily_annual_infiles = input_mean_daily_annual_infiles[1:]


# Read in the data
netcdf_reader = read_netcdf.ReadNetCDF()
dm_orig = netcdf_reader.read_into_xarray(input_daily_mean_infiles)
# Add some needed attributes
dm_list = []
time_dm = []
yr_dm = []
doy_dm = []
for din in dm_orig:
    ctime =  datetime.datetime.strptime(din[dm_var].valid_time,'%Y%m%d_%H%M%S')
    time_dm.append(ctime.strftime('%Y-%m-%d'))
    yr_dm.append(int(ctime.strftime('%Y')))
    doy_dm.append(int(ctime.strftime('%j')))
    din = din.assign_coords(time=ctime)
    din = din.expand_dims("time")
    dm_list.append(din)
time_dm = np.array(time_dm,dtype='datetime64[D]')
yr_dm = np.array(yr_dm)
doy_dm = np.array(doy_dm)
everything = xr.concat(dm_list,"time")
dm_data = np.array(everything[dm_var])

netcdf_reader2 = read_netcdf.ReadNetCDF()
mda_orig = netcdf_reader2.read_into_xarray(input_mean_daily_annual_infiles)
# Add some needed attributes
mda_list = []
time_mda = []
for din in mda_orig:
    ctime =  datetime.datetime.strptime(din[mda_var].valid_time,'%Y%m%d_%H%M%S')
    time_mda.append(ctime.strftime('%Y-%m-%d'))
    din = din.assign_coords(time=ctime)
    din = din.expand_dims("time")
    mda_list.append(din)
time_mda = np.array(time_mda,dtype='datetime64[D]')
everything2 = xr.concat(mda_list,"time")
mda_data = np.array(everything2[mda_var])

# Harmonic Analysis, first step is Forward Fast Fourier Transform
clmfft =  np.fft.rfft(mda_data,axis=0)

smthfft = np.zeros(clmfft.shape,dtype=complex)
for f in np.arange(0,3):
    smthfft[f,:,:] = clmfft[f,:,:]

clmout = np.fft.irfft(smthfft,axis=0)

# Subtract the clmout from the data to create anomalies, each year at a time
yrstrt = yr_dm[0]
yrend = yr_dm[-1]
anom = np.zeros(dm_data.shape)

for y in np.arange(yrstrt,yrend+1,1):
    curyr = np.where(yr_dm == y)
    dd = doy_dm[curyr] - 1
    ndd = len(curyr[0])
    clmshp = [np.arange(dd[0]*nobs,dd[0]*nobs+ndd,1)]
    anom[curyr,:,:] = dm_data[curyr,:,:] - clmout[clmshp,:,:]

# Assign to an xarray and write output
if not os.path.exists(anom_output_dir):
    os.makedirs(anom_output_dir)
for o in np.arange(0,len(dm_orig)):
    dm_orig_cur = dm_orig[o]
    dout = xr.Dataset({out_var: (("lat", "lon"),anom[o,:,:])},
    coords={"lat": dm_orig_cur.coords['lat'], "lon": dm_orig_cur.coords['lon']},
    attrs=dm_orig_cur.attrs)
    dout[out_var].attrs = dm_orig_cur[dm_var].attrs
    dout[out_var].attrs['long_name'] = dm_orig_cur[dm_var].attrs['long_name']+' Anomalies'
    dout[out_var].attrs['name'] = out_var

    # write to a file
    cvtime = datetime.datetime.strptime(dm_orig_cur[dm_var].valid_time,'%Y%m%d_%H%M%S')
    citime = datetime.datetime.strptime(dm_orig_cur[dm_var].init_time,'%Y%m%d_%H%M%S')
    cltime = (cvtime - citime)
    leadmin,leadsec = divmod(cltime.total_seconds(), 60)
    leadhr,leadmin = divmod(leadmin,60)
    lead_str = str(int(leadhr)).zfill(2)+str(int(leadmin)).zfill(2)+str(int(leadsec)).zfill(2)
    dout.to_netcdf(os.path.join(anom_output_dir,anom_output_base+'_'+lead_str+'L_'+cvtime.strftime('%Y%m%d')+'_'+cvtime.strftime('%H%M%S')+'V.nc'))

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_obsERA_obsOnly_RMM.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_mjo/UserScript_obsERA_obsOnly_RMM/plots. The output may include the regridded data and daily averaged files if those steps are turned on. Three output plots will be generated, a phase diagram, time series, and EOF plot:

* obs_RMM_comp_phase.png
* obs_RMM_time_series.png
* RMM_EOFs.png

Output from the Harmonic analysis pre-processing includes multiple files, one for each variable and day in the following format:

* ERA_OLR_anom_000000L_YYYYMMDD_000000V.nc
* ERA_U200_anom_000000L_YYYYMMDD_000000V.nc
* ERA_U850_anom_000000L_YYYYMMDD_000000V.nc

One variabe is output in each of the above files (not including the lat/lon fields). Those variables are:

* olr_anom(lat, lon)
* U_P200_mean_anom(lat, lon)
* U_P850_mean_anom(lat, lon)

Output from the regridding also includes one file for each day and variable with the following format:

* ERA_OLR_YYYYMMDD.nc
* ERA_U200_YYYYMMDD.nc
* ERA_U850_YYYYMMDD.nc

One variabe is output in each of the regriddwed files (not including the lat/lon fields). Those variables are:

* OLR_anom(lat, lon)
* U_P200_anom(lat, lon)
* U_P850_anom(lat, lon)

Keywords

Note

• S2SAppUseCase

• S2SMJOAppUseCase

• UserScriptUseCase

• NetCDFFileUseCase

• RegridDataPlaneUseCase

• PCPCombineUseCase

• METdataioUseCase

• METcalcpyUseCase

• METplotpyUseCase

Navigate to METplus Quick Search for Use Cases to discover other similar use cases.