Note

Go to the end to download the full example code.

UserScript: Compute Cross Spectra and Make a Plot

model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra.py

Scientific Objective

This use case calls the METcalcpy cross spectra function and then the METplotpy space time plot to compute cross-spectra and create a sample cross spectra diagram using sample data.

The space time plot and cross spectra calculations were created by Maria Gehne at the Physical Sciences Labratory in NOAA.

Version Added

METplus version 5.1

Datasets

Forecast: Unified Forecast System (UFS) Prototype 7

Observation: ERAI

Climatology: None

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 runs the UserScript wrapper tool to run a user provided script, in this case, cross_spectra.py and cross_spectra_plot.py.

METplus Workflow

Beginning time (VALID_BEG): 2014

End time (VALID_END): None

Increment between beginning and end times (VALID_INCREMENT): None

Sequence of forecast leads to process (LEAD_SEQ): None

This use case computes spectra and plots for the entire time period of data. The use case loops over two processes, computing and plotting the cross-spectra. The user is able to edit the process list to turn off the computation part if only plotting is desired, or vice versa.

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/UserScript_fcstS2S_obsERAI_CrossSpectra.conf

[config]

PROCESS_LIST = UserScript(comp_spectra), UserScript(plot_spectra)

# Note: time looping is not used in this use case
LOOP_BY = REALTIME
VALID_TIME_FMT = %Y
VALID_BEG = 2014

USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE

[user_env_vars]
# Make output base avabilable to the script
COMP_SPECTRA_SCRIPT_OUTPUT_DIR = {OUTPUT_BASE}/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/output

# YAML configuration file for the cross spectra calculation
COMP_SPECTRA_YAML_CONFIG_NAME = {METPLUS_BASE}/parm/use_cases/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/spectra_comp.yaml

# Input files for the cross spectra calculation
COMP_SPECTRA_INPUT_FILE_NAMES = {INPUT_BASE}/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/data/precip.erai.sfc.1p0.2x.2014-2016.nc,{INPUT_BASE}/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/data/prate_avg_ufs_p7_2014040100.nc,{INPUT_BASE}/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/data/u850_ufs_p7_2014040100.nc,{INPUT_BASE}/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/data/u200_ufs_p7_2014040100.nc

PLOT_SPECTRA_INPUT_FILE_NAMES = {OUTPUT_BASE}/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/output/SpaceTimeSpectra_ufs_p7_P_D850_symm_4spd.nc,{OUTPUT_BASE}/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/output/SpaceTimeSpectra_ufs_p7_P_D200_symm_4spd.nc,{OUTPUT_BASE}/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/output/SpaceTimeSpectra_ufs_p7_P_D200_symm_4spd.nc

PLOT_SPECTRA_YAML_CONFIG_NAME = {METPLUS_BASE}/parm/use_cases/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/spectra_plot.yaml

PLOT_SPECTRA_OUTPUT_DIR = {OUTPUT_BASE}/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/plots/


[comp_spectra]
# Settings for computing the cross-spectra
USER_SCRIPT_COMMAND = python {METCALCPY_BASE}/metcalcpy/contributed/spacetime/cross_spectra.py

LOG_FILE = "cross_spectra.log"
LOG_LEVEL = "DEBUG"

METCALCPY_BASE = {METPLUS_BASE}/../METcalcpy


[plot_spectra]
# settings for plotting the cross-spectra

USER_SCRIPT_COMMAND = python {METPLUS_BASE}/parm/use_cases/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/cross_spectra_plot.py

LOG_FILE = "cross_spectra_plot.log"
LOG_LEVEL = "INFO"

MET Configuration

There are no MET tools used in this use case.

Python Embedding

There is no Python embedding in this use case.

User Scripting

This use case uses a Python script to perform plotting.

parm/use_cases/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra/cross_spectra_plot.py
#!/usr/bin/env python3

"""
This is an example script for plotting cross spectral components. The script reads in output files computed
by the example_cross_spectra.py script and uses the plotly plotting routines in spacetime_plot.py to generate
a panel plot of coherence spectra.
"""
import numpy as np
import os
import sys
import xarray as xr

import metplotpy.contributed.spacetime_plot.spacetime_plot as stp
import metcalcpy.util.read_env_vars_in_config as readconfig


# Read in the YAML config file
# user can use their own, if none specified at the command line,
# use the "default" example YAML config file, spectra_plot_coh2.py
# Using a custom YAML reader so we can use environment variables
plot_config_file = os.getenv("PLOT_SPECTRA_YAML_CONFIG_NAME","spectra_plot.yaml")

config_dict = readconfig.parse_config(plot_config_file)

# Retrieve settings from config file
#pathdata is now set in the METplus conf file
#pathdata = config_dict['pathdata'][0]
plotpath = config_dict['plotpath'][0]
print("Output path ",plotpath)
model = config_dict['model']
var2 = config_dict['var2']
var3 = config_dict['var3']

# plot layout parameters
flim = 0.5  # maximum frequency in cpd for plotting
nWavePlt = 20  # maximum wavenumber for plotting
contourmin = 0.1  # contour minimum
contourmax = 0.8  # contour maximum
contourspace = 0.1  # contour spacing
N = [1, 2]  # wave modes for plotting
source = ""
spd = 4

symmetry = "symm"      #("symm", "asymm", "latband")
#filenames = os.environ.get("INPUT_FILE_NAMES","ERAI_TRMM_P_symn,ERAI_P_D850_symn,ERAI_P_D200_symn").split(",")
#vars1 = ['ERAI P', 'ERAI P', 'ERAI P']
#vars2 = ['TRMM', 'ERAI D850', 'ERAI D200']
filenames = os.environ.get("PLOT_SPECTRA_INPUT_FILE_NAMES","ERAI_P_D850_symn,ERAI_P_D200_symn,ERAI_P_D200_symn").split(",")
vars1 = [model+' P',model+' P',model+' P']
vars2 = [model+' '+var2,model+' '+var3,model+' '+var3]
nplot = len(vars1)

npanel =3
for pp in np.arange(0, nplot, 1):

    # read data from file
    var1 = vars1[pp]
    var2 = vars2[pp]
    print("Filename ",filenames[pp])
    fin = xr.open_dataset(filenames[pp])
    STC = fin['STC'][:, :, :]
    wnum = fin['wnum']
    freq = fin['freq']

    #ifreq = np.where((freq[:] >= 0) & (freq[:] <= flim))
    #iwave = np.where(abs(wnum[:]) <= nWavePlt)

    STC[:, freq[:] == 0, :] = 0.
    STC = STC.sel(wnum=slice(-nWavePlt, nWavePlt))
    STC = STC.sel(freq=slice(0, flim))
    coh2 = np.squeeze(STC[4, :, :])
    phs1 = np.squeeze(STC[6, :, :])
    phs2 = np.squeeze(STC[7, :, :])
    phs1.where(coh2 <= contourmin, drop=True)
    phs2.where(coh2 <= contourmin, drop=True)
    pow1 = np.squeeze(STC[0, :, :])
    pow2 = np.squeeze(STC[1, :, :])
    pow1.where(pow1 <= 0, drop=True)
    pow2.where(pow2 <= 0, drop=True)

    if pp == 0:
        ifreq = np.where((freq[:] >= 0) & (freq[:] <= flim))
        iwave = np.where(abs(wnum[:]) <= nWavePlt)
        Coh2 = np.full([npanel, len(freq[ifreq]), len(wnum[iwave])], np.nan)
        Phs1 = np.full([npanel, len(freq[ifreq]), len(wnum[iwave])], np.nan)
        Phs2 = np.full([npanel, len(freq[ifreq]), len(wnum[iwave])], np.nan)
        Pow1 = np.full([npanel, len(freq[ifreq]), len(wnum[iwave])], np.nan)
        Pow2 = np.full([npanel, len(freq[ifreq]), len(wnum[iwave])], np.nan)
        k = wnum[iwave]
        w = freq[ifreq]

    Coh2[pp, :, :] = coh2
    Phs1[pp, :, :] = phs1
    Phs2[pp, :, :] = phs2
    Pow1[pp, :, :] = np.log10(pow1)
    Pow2[pp, :, :] = np.log10(pow2)

    phstmp = Phs1
    phstmp = np.square(Phs1) + np.square(Phs2)
    phstmp = np.where(phstmp == 0, np.nan, phstmp)
    scl_one = np.sqrt(1 / phstmp)
    Phs1 = scl_one * Phs1
    Phs2 = scl_one * Phs2

# create output directory if it does not exist
if not os.path.exists(plotpath):
    print(f"Creating output directory: {plotpath}")
    os.makedirs(plotpath)

# plot coherence
stp.plot_coherence(Coh2, Phs1, Phs2, symmetry, source, vars1, vars2, plotpath, flim, 20, contourmin, contourmax,
                   contourspace, npanel, N)

# check if output file exists since plotting function
# doesn't return an error code on failure
expected_file = os.path.join(plotpath,
                             'SpaceTimeCoherence_.png')
if not os.path.exists(expected_file):
    print(f"ERROR: Could not create output file: {expected_file}")
    sys.exit(1)

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 -c /path/to/METplus/parm/use_cases/model_applications/s2s/UserScript_fcstS2S_obsERAI_CrossSpectra.conf -c /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.

Keywords

Note

  • UserScriptUseCase

  • S2SAppUseCase

  • METcalcpyUseCase

  • METplotpyUseCase

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

Gallery generated by Sphinx-Gallery