UserScript: Make RMM plots from calculated MJO indices

model_applications/ s2s_mjo/ UserScript_obsERA_obsOnly_RMM.py

Scientific Objective

To compute the Real-time Multivariate MJO Index (RMM) using Outgoing Longwave Radiation (OLR), 850 hPa wind (U850), and 200 hPa wind (U200). Specifically, RMM is computed using OLR, U850, and U200 data between 15N and 15S. Anomalies of OLR, U850, and U200 are created using a harmonic analysis, 120 day day mean 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.

Datasets

  • Forecast dataset: None

  • Observation dataset: ERA Reanlaysis Outgoing Longwave Radiation, 850 hPa wind and 200 hPa wind

External Dependencies

You will need to use a version of Python 3.6+ that has the following packages installed:

* numpy
* netCDF4
* datetime
* xarray
* matplotlib
* scipy
* pandas

If the version of Python used to compile MET did not have these libraries at the time of compilation, you will need to add these packages or create a new Python environment with these packages.

If this is the case, you will need to set the MET_PYTHON_EXE environment variable to the path of the version of Python you want to use. If you want this version of Python to only apply to this use case, set it in the [user_env_vars] section of a METplus configuration file.:

[user_env_vars] MET_PYTHON_EXE = /path/to/python/with/required/packages/bin/python

METplus Components

This use case runs the RMM driver which computes first computes anomalies of outgoing longwave raidation, 850 hPa wind and 200 hPa wind. Then, it regrids the data to 15S to 15N. Next, RMM is computed and a phase diagram, time series, and EOF plot are created. Inputs to the RMM driver include netCDF files that are in MET’s netCDF version. In addition, a text file containing the listing of these input netCDF files for OLR, u850 and u200 is required. Some optional pre-processing steps include using pcp_combine to compute daily means and the mean daily annual cycle for the data.

METplus Workflow

The RMM driver script python code is run for each lead time on the forecast and observations data. This example loops by valid time for the model pre-processing, and valid time for the other steps. This version is set to only process the creation of anomalies, regridding, and RMM calculation, omitting the caluclation of daily means and the mean daily annucal cycle pre-processing steps. However, the configurations for pre-processing are available for user reference.

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. The file UserScript_obsERA_obsOnly_RMM/RMM_driver.py runs the python program and UserScript_obsERA_obsOnly_RMM.conf sets the variables for all steps of the RMM use case.

[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

METplus sets environment variables based on the values in the METplus configuration file. These variables are referenced in the MET configuration file. 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 not controlled by an environment variable, you can add additional environment variables to be set only within the METplus environment using the [user_env_vars] section of the METplus configuration files. See the ‘User Defined Config’ section on the ‘System Configuration’ page of the METplus User’s Guide for more information.

Python Scripts

The RMM driver script orchestrates the calculation of the MJO indices and the generation of three RMM plots: parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/RMM_driver.py: The harmonic anomalies script creates anomalies of input data using a harmonic analysis: parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/compute_harmonic_anomalies.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, delim_whitespace=True, names=['eof1','eof2'])
    EOF1[0,:] = tmp.eof1.values
    EOF2[0,:] = tmp.eof2.values
    tmp = pd.read_csv(u850file, header=None, delim_whitespace=True, names=['eof1','eof2'])
    EOF1[1,:] = tmp.eof1.values
    EOF2[1,:] = tmp.eof2.values
    tmp = pd.read_csv(u200file, header=None, delim_whitespace=True, 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()
#!/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

This use case is run in the following ways:

  1. Passing in UserScript_obsERA_obsOnly_RMM.conf then a user-specific system configuration file:

    run_metplus.py -c /path/to/METplus/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM.conf -c /path/to/user_system.conf
    
  2. Modifying the configurations in parm/metplus_config, then passing in UserScript_obsERA_obsOnly_RMM.py:

    run_metplus.py -c /path/to/METplus/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM.conf
    

The following variables must be set correctly:

  • INPUT_BASE - Path to directory where sample data tarballs are unpacked (See Datasets section to obtain tarballs). This is not required to run METplus, but it is required to run the examples in parm/use_cases

  • OUTPUT_BASE - Path where METplus output will be written. This must be in a location where you have write permissions

  • MET_INSTALL_DIR - Path to location where MET is installed locally

Example User Configuration File:

[dir]
INPUT_BASE = /path/to/sample/input/data
OUTPUT_BASE = /path/to/output/dir
MET_INSTALL_DIR = /path/to/met-X.Y

Expected Output

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 model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM. This may include the regridded data and daily averaged files. In addition, three plots will be generated, a phase diagram, time series, and EOF plot, and the output location can be specified as RMM_PLOT_OUTPUT_DIR. If it is not specified, plots will be sent to model_applications/s2s_mjo/UserScript_obsERA_obsOnly_RMM/plots (relative to OUTPUT_BASE).

Keywords

Note

  • S2SAppUseCase

  • S2SMJOAppUseCase

  • NetCDFFileUseCase

  • RegridDataPlaneUseCase

  • PCPCombineUseCase

  • METcalcpyUseCase

  • METplotpyUseCase

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

sphinx_gallery_thumbnail_path = ‘_static/s2s_mjo-UserScript_obsERA_obsOnly_RMM.png’

Total running time of the script: (0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery