UserScript: Make ERA OMI plot from calculated MJO indices

model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI.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 outgoing longwave radiation (OLR) to compute the OLR based MJO Index (OMI), which is a convective index of the MJO. OMI is computed for the ERA observations and then displayed on a phase diagram to evaluate the model reprentation of this important oscillation. The code for computing OMI came from Maria Gehne at PSL.

Version Added

METplus version 4.1.0

Datasets

Forecast: None

Observation: ERA Reanlaysis Outgoing Longwave Radiation, 1979 - 2012

Climatology: None

EOFs: Observed OMI EOF1 and EOF2 patterns from the PSL Website (https://psl.noaa.gov/mjo/mjoindex/)

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 twice, first to create a list of EOF files and then to run the OMI calculation. In addition, there is one optional pre-processing steps, using Regrid-Data-Plane.

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 calculation (VALID_BEG): 01-01-1979

End time calculation (VALID_END): 12-30-2012

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

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

This use case does not loop, but the UserScript to create and EOF filelist is run once and the OMI driver script is run once. The EOF filelist is done separately since the EOF files are needed for each day of the year while the OMI calculation is on a separate time frame. The optional pre-processing step loops by valid time. The optional pre-processing step uses Regrid-Data-Plane to cut the observation grid to only include -20 to 20 latitude. This omitted step can be turned back on by using the PROCESS_LIST that is commented out in the file:

PROCESS_LIST = RegridDataPlane(regrid_obs_olr), UserScript(create_eof_filelist), UserScript(script_omi)

Settings for the optional pre-processing step can be found in the regrid_obs_olr section of the configuration. Data is not provided in the tarball to run this steps, but the configuration is provided for reference on how to set up this step.

The Phase diagram plot is created over a different time frame than the calculation, 10-01-2012 to 03-30-2012.

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_OMI.conf.

# OMI UserScript wrapper
[config]


# All steps, including pre-processing:
#PROCESS_LIST = RegridDataPlane(regrid_obs_olr), UserScript(create_eof_filelist), UserScript(script_omi)
# Finding EOF files and OMI Analysis script for the observations
PROCESS_LIST = UserScript(create_eof_filelist), UserScript(script_omi)


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

LEAD_SEQ = 0

# variables referenced in other sections

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

OBS_OLR_INPUT_DIR = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI/ERA
OBS_OLR_INPUT_TEMPLATE = OLR_{valid?fmt=%Y%m%d}.nc


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

REGRID_DATA_PLANE_VERIF_GRID = latlon 144 17 -20 0 2.5 2.5

REGRID_DATA_PLANE_METHOD = NEAREST

REGRID_DATA_PLANE_WIDTH = 1


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

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

OBS_REGRID_DATA_PLANE_RUN = {OBS_RUN}

OBS_DATA_PLANE_ONCE_PER_FIELD = False

OBS_REGRID_DATA_PLANE_VAR1_NAME = olr
OBS_REGRID_DATA_PLANE_VAR1_LEVELS = "({valid?fmt=%Y%m%d_%H%M%S},*,*)"
OBS_REGRID_DATA_PLANE_VAR1_OPTIONS = file_type=NETCDF_NCCF; censor_thresh=eq-999.0; censor_val=-9999.0;

OBS_REGRID_DATA_PLANE_VAR1_OUTPUT_FIELD_NAME = olr

OBS_REGRID_DATA_PLANE_INPUT_DIR = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI
OBS_REGRID_DATA_PLANE_OUTPUT_DIR = {OBS_OLR_INPUT_DIR}

OBS_REGRID_DATA_PLANE_INPUT_TEMPLATE = olr.1x.7920.nc
OBS_REGRID_DATA_PLANE_OUTPUT_TEMPLATE = {OBS_OLR_INPUT_TEMPLATE}


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

# Create the EOF filelists
[create_eof_filelist]

# Find the files for each time to create the time list
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE

# Valid Begin and End Times for the EOF files
VALID_BEG = 2012010100
VALID_END = 2012123100

# Find the EOF files for each time
# Filename templates for EOF1 and EOF2
USER_SCRIPT_INPUT_TEMPLATE = {INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI/EOF/eof1/eof{valid?fmt=%j}.txt,{INPUT_BASE}/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI/EOF/eof2/eof{valid?fmt=%j}.txt

# Name of the file containing the listing of input files
# The options are EOF1_INPUT and EOF2_INPUT
# *** Make sure the order is the same as the order of templates listed in USER_SCRIPT_INPUT_TEMPLATE
USER_SCRIPT_INPUT_TEMPLATE_LABELS = EOF1_INPUT, EOF2_INPUT

# 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 EOF1 and EOF2 Input


# Configurations for the OMI 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

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

# Phase Plot start date, end date, output name, and format
PHASE_PLOT_TIME_BEG = 2012010100
PHASE_PLOT_TIME_END = 2012033000
PHASE_PLOT_TIME_FMT = {VALID_TIME_FMT}
OBS_PHASE_PLOT_OUTPUT_NAME = obs_OMI_comp_phase
OBS_PHASE_PLOT_OUTPUT_FORMAT = png                                   


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

# Configurations for UserScript: Run the RMM Analysis driver
[script_omi]
#  Run the script once per lead time
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD

## Template of OLR filenames to input to the user-script
USER_SCRIPT_INPUT_TEMPLATE = {OBS_OLR_INPUT_DIR}/{OBS_OLR_INPUT_TEMPLATE}

## Name of the file containing the listing of OLR input files
## The options are OBS_OLR_INPUT and  FCST_OLR_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

# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mjo/common/OMI_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

Python Scripting

This use case runs the OMI driver which computes OMI and creates a phase diagram. Inputs to the OMI driver include netCDF files formatted in MET’s netCDF version. In addition, a text file containing the listing of these input netCDF files is required, as well as text file listings of the EOF1 and EOF2 files. These text files can be generated using the USER_SCRIPT_INPUT_TEMPLATES in the [create_eof_filelist] and [script_omi] sections.

For the OMI calculation, the OLR data are then projected onto Empirical Orthogonal Function (EOF) data that is computed for each day of the year, latitude, and longitude. 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. The OMI driver script orchestrates the calculation of the MJO indices and the generation of a phase diagram OMI plot. Variables for the OMI calculation are set in the [user_env_vars] section of the .conf file.

parm/use_cases/model_applications/s2s_mjo/common/OMI_driver.py

#!/usr/bin/env python3
"""
Driver Script to Compute OMI from input OLR data. Data is averaged from 20S-20N
"""

import numpy as np
import xarray as xr
import pandas as pd
import datetime
import os
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_omi_eofs(eof1_files, eof2_files):
    """
    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 3D DataArrays
    """

    # observed EOFs from NOAA PSL are saved in individual text files for each doy
    # horizontal resolution of EOFs is 2.5 degree
    EOF1 = xr.DataArray(np.empty([366,17,144]),dims=['doy','lat','lon'],
    coords={'doy':np.arange(1,367,1), 'lat':np.arange(-20,22.5,2.5), 'lon':np.arange(0,360,2.5)})
    EOF2 = xr.DataArray(np.empty([366,17,144]),dims=['doy','lat','lon'],
    coords={'doy':np.arange(1,367,1), 'lat':np.arange(-20,22.5,2.5), 'lon':np.arange(0,360,2.5)})
    nlat = len(EOF1['lat'])
    nlon = len(EOF1['lon'])

    for doy in range(len(eof1_files)):
        tmp1 = pd.read_csv(eof1_files[doy], header=None, delim_whitespace=True, names=['eof1'])
        tmp2 = pd.read_csv(eof2_files[doy], header=None, delim_whitespace=True, names=['eof2'])
        eof1 = xr.DataArray(np.reshape(tmp1.eof1.values,(nlat, nlon)),dims=['lat','lon'])
        eof2 = xr.DataArray(np.reshape(tmp2.eof2.values,(nlat, nlon)),dims=['lat','lon'])
        EOF1[doy,:,:] = eof1.values
        EOF2[doy,:,:] = eof2.values

    return EOF1, EOF2


def run_omi_steps(inlabel, olr_filetxt, spd, eof1, eof2, oplot_dir):

    # Read the listing of EOF 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:]

    # Read in the netCDF data from a list of files

    netcdf_reader = read_netcdf.ReadNetCDF()
    ds_orig = netcdf_reader.read_into_xarray(olr_input_files)

    # Add some needed attributes
    ds_list = []
    time = []
    for din in ds_orig:
        ctime =  datetime.datetime.strptime(din['olr'].valid_time,'%Y%m%d_%H%M%S')
        time.append(ctime.strftime('%Y-%m-%d'))
        din = din.assign_coords(time=ctime)
        din = din.expand_dims("time")
        ds_list.append(din)
    time = np.array(time,dtype='datetime64[D]')

    everything = xr.concat(ds_list,"time")
    olr = everything['olr']
    print(olr.min(), olr.max())

    # project OLR onto EOFs
    PC1, PC2 = cmi.omi(olr, time, spd, eof1, eof2)

    # 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_plot = PC1.sel(time=slice(phase_plot_start_time,phase_plot_end_time))

    # 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+'_OMI_comp_phase'))
    print(phase_plot_name)
    phase_plot_format = os.environ.get(inlabel+'_PHASE_PLOT_OUTPUT_FORMAT','png')

    # plot the PC phase diagram
    pmi.phase_diagram('OMI',PC1,PC2,np.array(pc1_plot['time'].dt.strftime("%Y-%m-%d").values),
        np.array(pc1_plot['time.month'].values),np.array(pc1_plot['time.day'].values),
        phase_plot_name,phase_plot_format)


def main():

    # Get Obs and Forecast OLR file listing
    obs_olr_filetxt = os.environ.get('METPLUS_FILELIST_OBS_OLR_INPUT','')
    fcst_olr_filetxt = os.environ.get('METPLUS_FILELIST_FCST_OLR_INPUT','')

    # Read in EOF filenames
    eof1_filetxt = os.environ['METPLUS_FILELIST_EOF1_INPUT']
    eof2_filetxt = os.environ['METPLUS_FILELIST_EOF2_INPUT']

    # Read the listing of EOF files
    with open(eof1_filetxt) as ef1:
        eof1_input_files = ef1.read().splitlines()
    if (eof1_input_files[0] == 'file_list'):
        eof1_input_files = eof1_input_files[1:]
    with open(eof2_filetxt) as ef2:
        eof2_input_files = ef2.read().splitlines()
    if (eof2_input_files[0] == 'file_list'):
        eof2_input_files = eof2_input_files[1:]

    # Read in the EOFs
    EOF1, EOF2 = read_omi_eofs(eof1_input_files, eof2_input_files)

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

    # Check for an output plot directory in the configs.  Create one if it does not exist
    oplot_dir = os.environ.get('OMI_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)

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

    # Run the steps to compute OMM
    # Observations
    if run_obs_omi == 'true':
        run_omi_steps('OBS', obs_olr_filetxt, spd, EOF1, EOF2, oplot_dir)

    # Forecast
    if run_fcst_omi == 'true':
        run_omi_steps('FCST', fcst_olr_filetxt, spd, EOF1, EOF2, oplot_dir)

    # nothing selected
    if (run_obs_omi == 'false') and (run_fcst_omi == '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()

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_OMI.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_fcstGFS_obsERA_OMI/plots. A Phase diagram plots will be generated:

* obs_OMI_comp_phase.png

If the pre-processing steps are turned on, the output will include the regridded data and daily averaged files.

Keywords

Note

• S2SAppUseCase

• S2SMJOAppUseCase

• RegridDataPlaneUseCase

• PCPCombineUseCase

• UserScriptUseCase

• METdataioUseCase

• METcalcpyUseCase

• METplotpyUseCase

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