Note

Go to the end to download the full example code

UserScript: Make OMI plot from calculated MJO indices

model_applications/ s2s_mjo/ UserScript_obsERA_obsOnly_OMI.py

Scientific Objective

To use Outgoing Longwave Radiation (OLR) to compute the OLR based MJO Index (OMI). Specifically, OMI is computed using OLR data between 20N and 20S. 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.

Datasets

  • Forecast dataset: None

  • Observation dataset: ERA Reanlaysis Outgoing Longwave Radiation.

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 OMI driver which computes OMI and creates a phase diagram. Inputs to the OMI driver include netCDF files that are in MET’s netCDF version. In addition, a txt 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. Some optional pre-processing steps include using regrid_data_plane to either regrid your data or cut the domain to 20N - 20S.

METplus Workflow

The OMI 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 OMI calculation and creating a text file listing of the EOF files, omitting the creation of daily means for the model and the regridding 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_OMI.conf. The file UserScript_obsERA_obsOnly_OMI/OMI_driver.py runs the python program and UserScript_fcstGFS_obsERA_OMI.conf sets the variables for all steps of the OMI use case.

# 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/UserScript_obsERA_obsOnly_OMI/OMI_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 OMI driver script orchestrates the calculation of the MJO indices and the generation of a phase diagram OMI plot: parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI/OMI_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 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)):
        doystr = str(doy).zfill(3)
        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))
    PC2_plot = PC2.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

This use case is run in the following ways:

  1. Passing in UserScript_obsERA_obsOnly_OMI.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_OMI.conf -c /path/to/user_system.conf

  2. Modifying the configurations in parm/metplus_config, then passing in UserScript_obsERA_obsOnly_OMI.py:

    run_metplus.py -c /path/to/METplus/parm/use_cases/model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI.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_OMI. This may include the regridded data and daily averaged files. In addition, the phase diagram plots will be generated and the output location can be specified as OMI_PLOT_OUTPUT_DIR. If it is not specified, plots will be sent to model_applications/s2s_mjo/UserScript_obsERA_obsOnly_OMI/plots (relative to OUTPUT_BASE).

Keywords

Note

  • S2SAppUseCase

  • S2SMJOAppUseCase

  • RegridDataPlaneUseCase

  • PCPCombineUseCase

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

sphinx_gallery_thumbnail_path = ‘_static/s2s_mjo-UserScript_obsERA_obsOnly_OMI.png’

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

Gallery generated by Sphinx-Gallery