WeatherRegime Calculation: GFS and ERA RegridDataPlane, PcpCombine, and WeatherRegime python code

model_applications/ s2s_mid_lat/ UserScript_fcstGFS_obsERA_WeatherRegime.py

Scientific Objective

To perform a weather regime analysis using 500 mb height data. There are 2 pre- processing steps, RegridDataPlane and PcpCombine, and 4 steps in the weather regime analysis, elbow, EOFs, K means, and the Time frequency. The elbow and K means steps begin with K means clustering. Elbow then computes the sum of squared distances for clusters 1 - 14 and draws a straight line from the sum of squared distance for the clusters. This helps determine the optimal cluster number by examining the largest difference between the curve and the straight line. The EOFs step is optional. It computes an empirical orthogonal function analysis. The K means step uses clustering to compute the frequency of occurrence and anomalies for each cluster to give the most common weather regimes. Then, the time frequency computes the frequency of each weather regime over a user specified time frame. Finally, stat_analysis can be run to compute an categorical analysis of the weather regime classification or an anomaly correlation of the time frequency data.

Datasets

  • Forecast dataset: GFS Forecast 500 mb height.

  • Observation dataset: ERA Reanlaysis 500 mb height.

External Dependencies

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

* numpy
* netCDF4
* datetime
* pylab
* scipy
* sklearn
* eofs

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 weather regime driver script which runs the steps the user lists in STEPS_OBS. The possible steps are regridding, time averaging, creating a list of input files for the weather regime calculation, computing the elbow (ELBOW), plotting the elbow (PLOTELBOW), computing EOFs (EOF), plotting EOFs (PLOTEOF), computing K means (KMEANS), plotting the K means (PLOTKMEANS), computing a time frequency of weather regimes (TIMEFREQ) and plotting the time frequency (PLOTFREQ). All variables are set up in the UserScript .conf file. The pre- processing steps and stat_analysis are listed in the process list, and are formatted as follows:

PROCESS_LIST = RegridDataPlane(regrid_obs), PcpCombine(daily_mean_obs), UserScript(script_wr)

The other steps are listed in the [user_env_vars] section of the UserScript .conf file in the following format: OBS_STEPS = ELBOW+PLOTELBOW+EOF+PLOTEOF+KMEANS+PLOTKMEANS+TIMEFREQ+PLOTFREQ FCST_STEPS = ELBOW+PLOTELBOW+EOF+PLOTEOF+KMEANS+PLOTKMEANS+TIMEFREQ+PLOTFREQ

METplus Workflow

The weather regime python code is run for each time for the forecast and observations data. This example loops by valid time. This version is set to only process the weather regime steps (ELBOW, PLOTELBOW, EOF, PLOTEOF, KMEANS, PLOTKMEANS, TIMEFREQ, PLOTFREQ) and stat_analysis, omitting the regridding, time averaging, and creating the file list 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_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime.py. The file UserScript_fcstGFS_obsERA_WeatherRegime.conf runs the python program and sets the variables for all steps of the Weather Regime use case including data paths.

[config]

# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime.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 pre-processing:
# PROCESS_LIST = RegridDataPlane(regrid_obs), PcpCombine(daily_mean_obs), UserScript(script_wr), StatAnalysis(sanal_wrclass), StatAnalysis(sanal_wrfreq)
# Weather Regime Analysis and stat_analysis:

PROCESS_LIST = UserScript(script_wr), StatAnalysis(sanal_wrclass), StatAnalysis(sanal_wrfreq)


###
# 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 = 1979120100
VALID_BEG = 2000120100
VALID_END = 2017022800
VALID_INCREMENT = 86400

LEAD_SEQ = 0

# Only Process DJF
SKIP_TIMES = "%m:begin_end_incr(3,11,1)", "%m%d:1201,0229"


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

# Regridding Pre-Processing Step
[regrid_obs]

VALID_BEG = 2000120200
VALID_END = 2017022818
VALID_INCREMENT = 21600

# REGRID_DATA_PLANE (Pre Processing Step 1), currently turned off
# Run regrid_data_plane on forecast data
OBS_REGRID_DATA_PLANE_RUN = True

OBS_DATA_PLANE_ONCE_PER_FIELD = False

OBS_REGRID_DATA_PLANE_VAR1_INPUT_FIELD_NAME = Z
OBS_REGRID_DATA_PLANE_VAR1_INPUT_LEVEL = P500
OBS_REGRID_DATA_PLANE_VAR1_OUTPUT_FIELD_NAME = Z500

# Mask to use for regridding
# A 1 degree latitude/longitude grid running 24 to 54 degrees latitude
# and 230 to 300 degrees longitude
REGRID_DATA_PLANE_VERIF_GRID = latlon 71 31 54 230 -1.0 1.0

# Method to run regrid_data_plane, not setting this will default to NEAREST
REGRID_DATA_PLANE_METHOD = BILIN

# Regridding width used in regrid_data_plane, not setting this will default to 1
REGRID_DATA_PLANE_WIDTH = 2

OBS_REGRID_DATA_PLANE_INPUT_DIR = {INPUT_BASE}/model_applications/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/ERA/OrigData
OBS_REGRID_DATA_PLANE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/ERA/Regrid

OBS_REGRID_DATA_PLANE_INPUT_TEMPLATE = {valid?fmt=%Y%m}/ei.oper.an.pl.regn128sc.{valid?fmt=%Y%m%d%H}
OBS_REGRID_DATA_PLANE_OUTPUT_TEMPLATE = {valid?fmt=%Y%m%d}/Z500_6hourly_{init?fmt=%Y%m%d%H}_NH.nc


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

# Daily Mean Pre-Processing Step
[daily_mean_obs]

VALID_BEG = 2000120218
VALID_END = 2017022818

OBS_PCP_COMBINE_RUN = True


OBS_PCP_COMBINE_METHOD = DERIVE
OBS_PCP_COMBINE_STAT_LIST = MEAN

OBS_PCP_COMBINE_INPUT_ACCUMS = 6
OBS_PCP_COMBINE_INPUT_NAMES = Z500
OBS_PCP_COMBINE_INPUT_LEVELS = "(*,*)"
OBS_PCP_COMBINE_INPUT_OPTIONS = convert(x) = x / 9.81; set_attr_valid = "{valid?fmt=%Y%m%d_%H%M%S?shift=-64800}";

# Convert height and derive mean over 24 hours
OBS_PCP_COMBINE_OUTPUT_ACCUM = 24
OBS_PCP_COMBINE_DERIVE_LOOKBACK = 24

# Name output variable Z500
OBS_PCP_COMBINE_OUTPUT_NAME = Z500

OBS_PCP_COMBINE_INPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/ERA/Regrid
OBS_PCP_COMBINE_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/ERA/Daily

OBS_PCP_COMBINE_INPUT_TEMPLATE = {valid?fmt=%Y%m%d}/Z500_6hourly_{valid?fmt=%Y%m%d%H}_NH.nc
OBS_PCP_COMBINE_OUTPUT_TEMPLATE = Z500_daily_{valid?fmt=%Y%m%d?shift=-64800}_NH.nc


# Variables for the Weather Regime code
[user_env_vars]
# Steps to Run
FCST_STEPS = ELBOW+PLOTELBOW+EOF+PLOTEOF+KMEANS+PLOTKMEANS+TIMEFREQ+PLOTFREQ
OBS_STEPS = ELBOW+PLOTELBOW+EOF+PLOTEOF+KMEANS+PLOTKMEANS+TIMEFREQ+PLOTFREQ

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

# Number of Seasons and Days per season that should be available
# The code will fill missing data, but requires the same number of days per
# season for each year.  You may need to omit leap days if February is part of
# the processing
NUM_SEASONS = 17
DAYS_PER_SEASON = 89

# Variable for the Z500 data
OBS_WR_VAR = Z500
FCST_WR_VAR = Z500_P500

# Weather Regime Number
OBS_WR_NUMBER = 6
FCST_WR_NUMBER = {OBS_WR_NUMBER}

# Number of clusters
OBS_NUM_CLUSTERS = 20
FCST_NUM_CLUSTERS = {OBS_NUM_CLUSTERS}

# Number of principal components
OBS_NUM_PCS = 10
FCST_NUM_PCS = {OBS_NUM_PCS}

# Time (in timesteps) over which to compute weather regime frequencies
# i.e. if your data time step is days and you want to average over 7
# days, input 7
# Optional, only needed if you want to compute frequencies
OBS_WR_FREQ = 7
FCST_WR_FREQ = {OBS_WR_FREQ}

# These variables control reordering the forecast weather regime to match the
# observations if their orders are different
# REORDER_FCST_MANUAL will use the order in FCST_ORDER, whereas REORDER_FCST will 
# use a pattern correlation to reorder
# It is recommended to set REORDER_FCST_MANUAL to False if this is the first time running the
# case
REORDER_FCST = True
REORDER_FCST_MANUAL = False
#Order to use if REORDER_FCST_MANUAL = True; will be ignored if REORER_FCST_MANUAL = False
FCST_ORDER = 1,3,4,2,5,6

# Type, name and directory of Output File for weather regime classification
# Type options are text or netcdf
OBS_WR_OUTPUT_FILE_TYPE = text
OBS_WR_OUTPUT_FILE = obs_weather_regime_class
FCST_WR_OUTPUT_FILE_TYPE = text
FCST_WR_OUTPUT_FILE = fcst_weather_regime_class
WR_OUTPUT_FILE_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime

# Directory to send output plots
WR_PLOT_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/plots/

# Elbow Plot Title and output file name
OBS_ELBOW_PLOT_TITLE = ERA Elbow Method For Optimal k
OBS_ELBOW_PLOT_OUTPUT_NAME = obs_elbow
FCST_ELBOW_PLOT_TITLE = GFS Elbow Method For Optimal k
FCST_ELBOW_PLOT_OUTPUT_NAME = fcst_elbow


# EOF plot output name and contour levels
OBS_EOF_PLOT_OUTPUT_NAME = obs_eof
FCST_EOF_PLOT_OUTPUT_NAME = fcst_eof
EOF_PLOT_LEVELS = -50, -45, -40, -35, -30, -25, -20, -15, -10, -5,  0,  5, 10, 15, 20, 25, 30, 35, 40, 45, 50

# K means Plot Output Name and contour levels
OBS_KMEANS_PLOT_OUTPUT_NAME = obs_kmeans
FCST_KMEANS_PLOT_OUTPUT_NAME = fcst_kmeans
KMEANS_PLOT_LEVELS = -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80

# Frequency Plot title and output file name
OBS_FREQ_PLOT_TITLE = ERA Seasonal Cycle of WR Days/Week (2000-2017)
OBS_FREQ_PLOT_OUTPUT_NAME = obs_freq
FCST_FREQ_PLOT_TITLE = GFS Seasonal Cycle of WR Days/Week (2000-2017)
FCST_FREQ_PLOT_OUTPUT_NAME = fcst_freq

# MPR file information
MASK_NAME = FULL
WR_MPR_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/mpr


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

# Run the Weather Regime Script
[script_wr]

LEAD_SEQ = 24

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_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/ERA/Daily/Z500_daily_{valid?fmt=%Y%m%d}_NH.nc,{INPUT_BASE}/model_applications/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/GFS/Daily/Z500_{init?fmt=%Y%m%d}_{lead?fmt=%HHH}_NH.nc

# Name of the file containing the listing of input files
# The options are OBS_INPUT for observations or FCST_INPUT for forecast
# Or, set OBS_INPUT, FCST_INPUT if doing both and make sure the USER_SCRIPT_INPUT_TEMPLATE is ordered:
# observation_template, forecast_template
USER_SCRIPT_INPUT_TEMPLATE_LABELS = OBS_INPUT, FCST_INPUT

# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/WeatherRegime_driver.py


###
# StatAnalysis(sanal_wrclass) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#statanalysis
###

[sanal_wrclass]

VALID_TIME_FMT = %Y%m%d
VALID_BEG = 20001202
VALID_END = 20170228

STAT_ANALYSIS_CONFIG_FILE = {PARM_BASE}/met_config/STATAnalysisConfig_wrapped

MODEL1 = GFS
MODEL1_OBTYPE = ADPUPA

STAT_ANALYSIS_JOB_NAME = aggregate_stat
STAT_ANALYSIS_JOB_ARGS = -out_line_type MCTS -out_thresh >=1,>=2,>=3,>=4,>=5 -out_stat [out_stat_file]

MODEL_LIST = {MODEL1}
FCST_LEAD_LIST = 24
LINE_TYPE_LIST = MPR

GROUP_LIST_ITEMS = MODEL_LIST
LOOP_LIST_ITEMS = FCST_LEAD_LIST

MODEL1_STAT_ANALYSIS_LOOKIN_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/mpr/WeatherRegime

STAT_ANALYSIS_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime

MODEL1_STAT_ANALYSIS_OUT_STAT_TEMPLATE = {model?fmt=%s}_ERA_WRClass_{lead?fmt=%H%M%S}L_MCTS.stat


###
# StatAnalysis(sanal_wrfreq) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#statanalysis
###

[sanal_wrfreq]

VALID_TIME_FMT = %Y%m%d
VALID_BEG = 20001202
VALID_END = 20170228

STAT_ANALYSIS_CONFIG_FILE = {PARM_BASE}/met_config/STATAnalysisConfig_wrapped

MODEL1 = GFS
MODEL1_OBTYPE = ADPUPA

STAT_ANALYSIS_JOB_NAME = aggregate_stat
STAT_ANALYSIS_JOB_ARGS = -out_line_type CNT -by DESC -out_stat [out_stat_file]

MODEL_LIST = {MODEL1}
FCST_LEAD_LIST = 24
LINE_TYPE_LIST = MPR

GROUP_LIST_ITEMS = MODEL_LIST
LOOP_LIST_ITEMS = FCST_LEAD_LIST

MODEL1_STAT_ANALYSIS_LOOKIN_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/mpr/freq

STAT_ANALYSIS_OUTPUT_DIR = {OUTPUT_BASE}/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime

MODEL1_STAT_ANALYSIS_OUT_STAT_TEMPLATE = {model?fmt=%s}_ERA_WR_freq_{lead?fmt=%H%M%S}L_CNT.stat

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.

See the following files for more information about the environment variables set in this configuration file.

parm/use_cases/met_tool_wrapper/RegridDataPlane/RegridDataPlane.py parm/use_cases/met_tool_wrapper/PCPCombine/PCPCOmbine_derive.py parm/use_cases/met_tool_wrapper/StatAnalysis/StatAnalysis.py

Python Scripts

This use case uses Python scripts to perform the blocking calculation

parm/use_cases/model_applications/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime/WeatherRegime_driver.py: This script calls the requested steps in the blocking analysis for a forecast, observation, or both. The possible steps are computing the elbow, computing EOFs, and computing weather regimes using k means clustering.

metcalcpy/contributed/blocking_weather_regime/WeatherRegime.py: This script runs the requested steps, containing the code for computing the bend in the elbow, computing EOFs, and computing weather regimes using k means clustering. See the METcalcpy Weather Regime Calculation Script for more information.

metcalcpy/contributed/blocking_weather_regime/Blocking_WeatherRegime_util.py: This script contains functions used by both the blocking anwd weather regime analysis, including the code for determining which steps the user wants to run, and finding and reading the input files in the format from the output pre-processing steps. See the METcalcpy Utility script for more information.

#!/usr/bin/env python3
import sys
import os
import numpy as np
import netCDF4
import warnings

from metcalcpy.contributed.blocking_weather_regime.WeatherRegime import WeatherRegimeCalculation
from metcalcpy.contributed.blocking_weather_regime.Blocking_WeatherRegime_util import parse_steps, read_nc_met, write_mpr_file, reorder_fcst_regimes,reorder_fcst_regimes_correlate
from metplotpy.contributed.weather_regime import plot_weather_regime as pwr


def main():

    steps_list_fcst,steps_list_obs = parse_steps()

    if not steps_list_obs and not steps_list_fcst:
        warnings.warn('No processing steps requested for either the model or observations,')
        warnings.warn(' nothing will be run')
        warnings.warn('Set FCST_STEPS and/or OBS_STEPS in the [user_env_vars] section to process data')


    ######################################################################
    # Blocking Calculation and Plotting
    ######################################################################
    # Set up the data
    steps_obs = WeatherRegimeCalculation('OBS')
    steps_fcst = WeatherRegimeCalculation('FCST')

    # Check to see if there is a plot directory
    oplot_dir = os.environ.get('WR_PLOT_OUTPUT_DIR','')
    obase = os.environ['SCRIPT_OUTPUT_BASE']
    if not oplot_dir:
        oplot_dir = os.path.join(obase,'plots')
    if not os.path.exists(oplot_dir):
        os.makedirs(oplot_dir)

     # Check to see if there is a mpr output directory
    mpr_outdir = os.environ.get('WR_MPR_OUTPUT_DIR','')
    if not mpr_outdir:
        mpr_outdir = os.path.join(obase,'mpr')

     # Get number of seasons and days per season
    nseasons = int(os.environ['NUM_SEASONS'])
    dseasons = int(os.environ['DAYS_PER_SEASON'])

    # Grab the Daily text files
    obs_wr_filetxt = os.environ.get('METPLUS_FILELIST_OBS_INPUT','')
    fcst_wr_filetxt = os.environ.get('METPLUS_FILELIST_FCST_INPUT','')

    if ("ELBOW" in steps_list_obs) or ("EOF" in steps_list_obs) or ("KMEANS" in steps_list_obs):
        with open(obs_wr_filetxt) as owl:
            obs_infiles = owl.read().splitlines()
        # Remove the first line if it's there
        if (obs_infiles[0] == 'file_list'):
            obs_infiles = obs_infiles[1:]
        if len(obs_infiles) != (nseasons*dseasons):
            raise Exception('Invalid Obs data; each year must contain the same date range to calculate seasonal averages.')
        obs_invar = os.environ.get('OBS_WR_VAR','')
        z500_obs,lats_obs,lons_obs,timedict_obs = read_nc_met(obs_infiles,obs_invar,nseasons,dseasons)
        z500_detrend_obs,z500_detrend_2d_obs = steps_obs.weights_detrend(lats_obs,lons_obs,z500_obs)

    if ("ELBOW" in steps_list_fcst) or ("EOF" in steps_list_fcst) or("KMEANS" in steps_list_fcst):
        with open(fcst_wr_filetxt) as fwl:
            fcst_infiles = fwl.read().splitlines()
        # Remove the first line if it's there
        if (fcst_infiles[0] == 'file_list'):
            fcst_infiles = fcst_infiles[1:]
        if len(fcst_infiles) != (nseasons*dseasons):
            raise Exception('Invalid Obs data; each year must contain the same date range to calculate seasonal averages.')
        fcst_invar = os.environ.get('FCST_WR_VAR','')
        z500_fcst,lats_fcst,lons_fcst,timedict_fcst = read_nc_met(fcst_infiles,fcst_invar,nseasons,dseasons)
        z500_detrend_fcst,z500_detrend_2d_fcst = steps_fcst.weights_detrend(lats_fcst,lons_fcst,z500_fcst)


    if ("ELBOW" in steps_list_obs):
        print('Running Obs Elbow')
        K_obs,d_obs,mi_obs,line_obs,curve_obs = steps_obs.run_elbow(z500_detrend_2d_obs)

    if ("ELBOW" in steps_list_fcst):
        print('Running Forecast Elbow')
        K_fcst,d_fcst,mi_fcst,line_fcst,curve_fcst = steps_fcst.run_elbow(z500_detrend_2d_fcst)

    if ("PLOTELBOW" in steps_list_obs):
        if not ("ELBOW" in steps_list_obs):
            raise Exception('Must run observed Elbow before plotting observed elbow.')
        print('Creating Obs Elbow plot')
        elbow_plot_title = os.environ.get('OBS_ELBOW_PLOT_TITLE','Elbow Method For Optimal k')
        elbow_plot_outname = os.path.join(oplot_dir,os.environ.get('OBS_ELBOW_PLOT_OUTPUT_NAME','obs_elbow'))
        pwr.plot_elbow(K_obs,d_obs,mi_obs,line_obs,curve_obs,elbow_plot_title,elbow_plot_outname)

    if ("PLOTELBOW" in steps_list_fcst):
        if not ("ELBOW" in steps_list_fcst):
            raise Exception('Must run forecast Elbow before plotting forecast elbow.')
        print('Creating Forecast Elbow plot')
        elbow_plot_title = os.environ.get('FCST_ELBOW_PLOT_TITLE','Elbow Method For Optimal k')
        elbow_plot_outname = os.path.join(oplot_dir,os.environ.get('FCST_ELBOW_PLOT_OUTPUT_NAME','fcst_elbow'))
        pwr.plot_elbow(K_fcst,d_fcst,mi_fcst,line_fcst,curve_fcst,elbow_plot_title,elbow_plot_outname)


    if ("EOF" in steps_list_obs):
        print('Running Obs EOF')
        eof_obs,pc_obs,wrnum_obs,variance_fractions_obs = steps_obs.Calc_EOF(z500_obs)
        z500_detrend_2d_obs = steps_obs.reconstruct_heights(eof_obs,pc_obs,z500_detrend_2d_obs.shape)

    if ("EOF" in steps_list_fcst):
        print('Running Forecast EOF')
        eof_fcst,pc_fcst,wrnum_fcst,variance_fractions_fcst = steps_fcst.Calc_EOF(z500_fcst)
        z500_detrend_2d_fcst = steps_fcst.reconstruct_heights(eof_fcst,pc_fcst,z500_detrend_2d_fcst.shape)

    if ("PLOTEOF" in steps_list_obs):
        if not ("EOF" in steps_list_obs):
            raise Exception('Must run observed EOFs before plotting observed EOFs.')
        print('Plotting Obs EOFs')
        pltlvls_str = os.environ['EOF_PLOT_LEVELS'].split(',')
        pltlvls = [float(pp) for pp in pltlvls_str]
        eof_plot_outname = os.path.join(oplot_dir,os.environ.get('OBS_EOF_PLOT_OUTPUT_NAME','obs_eof'))
        pwr.plot_eof(eof_obs,wrnum_obs,variance_fractions_obs,lons_obs,lats_obs,eof_plot_outname,pltlvls)

    if ("PLOTEOF" in steps_list_fcst):
        if not ("EOF" in steps_list_fcst):
            raise Exception('Must run forecast EOFs before plotting forecast EOFs.')
        print('Plotting Forecast EOFs')
        pltlvls_str = os.environ['EOF_PLOT_LEVELS'].split(',')
        pltlvls = [float(pp) for pp in pltlvls_str]
        eof_plot_outname = os.path.join(oplot_dir,os.environ.get('FCST_EOF_PLOT_OUTPUT_NAME','fcst_eof'))
        pwr.plot_eof(eof_fcst,wrnum_fcst,variance_fractions_fcst,lons_fcst,lats_fcst,eof_plot_outname,pltlvls)


    if ("KMEANS" in steps_list_obs):
        print('Running Obs K Means')
        kmeans_obs,wrnum_obs,perc_obs,wrc_obs= steps_obs.run_K_means(z500_detrend_2d_obs,timedict_obs,z500_obs.shape)
        steps_obs.write_K_means_file(timedict_obs,wrc_obs)

    if ("KMEANS" in steps_list_fcst):
        print('Running Forecast K Means')
        kmeans_fcst,wrnum_fcst,perc_fcst,wrc_fcst = steps_fcst.run_K_means(z500_detrend_2d_fcst,timedict_fcst,
            z500_fcst.shape)
        reorder_fcst = os.environ.get('REORDER_FCST','False').lower()
        reorder_fcst_manual = os.environ.get('REORDER_FCST_MANUAL','False').lower()
        if (reorder_fcst == 'true') and ("KMEANS" in steps_list_obs): 
            kmeans_fcst,perc_fcst,wrc_fcst = reorder_fcst_regimes_correlate(kmeans_obs,kmeans_fcst,perc_fcst,wrc_fcst,wrnum_fcst)
        if reorder_fcst_manual == 'true':
            fcst_order_str = os.environ['FCST_ORDER'].split(',')
            fcst_order = [int(fo) for fo in fcst_order_str]
            kmeans_fcst,perc_fcst,wrc_fcst = reorder_fcst_regimes(kmeans_fcst,perc_fcst,wrc_fcst,wrnum_fcst,fcst_order)
        steps_fcst.write_K_means_file(timedict_fcst,wrc_fcst)


        # Write matched pair output for weather regime classification
        modname = os.environ.get('MODEL_NAME','GFS')
        maskname = os.environ.get('MASK_NAME','FULL')
        mpr_full_outdir = os.path.join(mpr_outdir,'WeatherRegime')
        wr_outfile_prefix = os.path.join(mpr_full_outdir,'weather_regime_stat_'+modname)
        wrc_obs_mpr = wrc_obs[:,:,np.newaxis]
        wrc_fcst_mpr = wrc_fcst[:,:,np.newaxis]
        if not os.path.exists(mpr_full_outdir):
            os.makedirs(mpr_full_outdir)
        write_mpr_file(wrc_obs_mpr,wrc_fcst_mpr,[0.0],[0.0],timedict_obs,timedict_fcst,modname,'NA',
            'WeatherRegimeClass','class','Z500','WeatherRegimeClass','class','Z500',maskname,'500',wr_outfile_prefix)

    if ("PLOTKMEANS" in steps_list_obs):
        if not ("KMEANS" in steps_list_obs):
            raise Exception('Must run observed Kmeans before plotting observed Kmeans.')
        print('Plotting Obs K Means')
        pltlvls_str = os.environ['KMEANS_PLOT_LEVELS'].split(',')
        pltlvls = [float(pp) for pp in pltlvls_str]
        kmeans_plot_outname = os.path.join(oplot_dir,os.environ.get('OBS_KMEANS_PLOT_OUTPUT_NAME','obs_kmeans'))
        pwr.plot_K_means(kmeans_obs,wrnum_obs,lons_obs,lats_obs,perc_obs,kmeans_plot_outname,pltlvls)

    if ("PLOTKMEANS" in steps_list_fcst):
        if not ("KMEANS" in steps_list_fcst):
            raise Exception('Must run forecast Kmeans before plotting forecast Kmeans.')
        print('Plotting Forecast K Means')
        pltlvls_str = os.environ['KMEANS_PLOT_LEVELS'].split(',')
        pltlvls = [float(pp) for pp in pltlvls_str]
        kmeans_plot_outname = os.path.join(oplot_dir,os.environ.get('FCST_KMEANS_PLOT_OUTPUT_NAME','fcst_kmeans'))
        pwr.plot_K_means(kmeans_fcst,wrnum_fcst,lons_fcst,lats_fcst,perc_fcst,kmeans_plot_outname,pltlvls)


    if ("TIMEFREQ" in steps_list_obs):
        if not ("KMEANS" in steps_list_obs):
            raise Exception('Must run observed Kmeans before running frequencies.')
        wrfreq_obs,dlen_obs,ts_diff_obs = steps_obs.compute_wr_freq(wrc_obs)

    if ("TIMEFREQ" in steps_list_fcst):
        if not ("KMEANS" in steps_list_fcst):
            raise Exception('Must run forecast Kmeans before running frequencies.')
        wrfreq_fcst,dlen_fcst,ts_diff_fcst = steps_fcst.compute_wr_freq(wrc_fcst)

    if ("TIMEFREQ" in steps_list_obs) and ("TIMEFREQ" in steps_list_fcst):
        # Write matched pair output for frequency of each weather regime
        modname = os.environ.get('MODEL_NAME','GFS')
        maskname = os.environ.get('MASK_NAME','FULL')
        mpr_full_outdir = os.path.join(mpr_outdir,'freq')
        timedict_obs_mpr = {'init':timedict_obs['init'][:,ts_diff_obs-1:],
            'valid':timedict_obs['valid'][:,ts_diff_obs-1:],'lead':timedict_obs['lead'][:,ts_diff_obs-1:]}
        timedict_fcst_mpr = {'init':timedict_fcst['init'][:,ts_diff_fcst-1:],
            'valid':timedict_fcst['valid'][:,ts_diff_fcst-1:],'lead':timedict_fcst['lead'][:,ts_diff_fcst-1:]}
        wrfreq_obs_mpr = wrfreq_obs[:,:,:,np.newaxis]
        wrfreq_fcst_mpr = wrfreq_fcst[:,:,:,np.newaxis]
        if not os.path.exists(mpr_full_outdir):
            os.makedirs(mpr_full_outdir)
        for wrn in np.arange(wrnum_obs):
            wr_outfile_prefix = os.path.join(mpr_full_outdir,'weather_regime'+str(wrn+1).zfill(2)+'_freq_stat_'+modname)
            write_mpr_file(wrfreq_obs_mpr[wrn,:,:,:],wrfreq_fcst_mpr[wrn,:,:,:],[0.0],[0.0],timedict_obs,
                timedict_fcst,modname,str(wrn+1).zfill(2),'WeatherRegimeFreq','percent','Z500','WeatherRegimeFreq',
                'percent','Z500',maskname,'500',wr_outfile_prefix)

    if ("PLOTFREQ" in steps_list_obs):
        if not ("TIMEFREQ" in steps_list_obs):
            raise Exception('Must run observed Frequency calculation before plotting the frequencies.')
        freq_plot_title = os.environ.get('OBS_FREQ_PLOT_TITLE','Seasonal Cycle of WR Days/Week')
        freq_plot_outname = os.path.join(oplot_dir,os.environ.get('OBS_FREQ_PLOT_OUTPUT_NAME','obs_freq'))
        # Compute mean
        wrmean_obs = np.nanmean(wrfreq_obs,axis=1)
        pwr.plot_wr_frequency(wrmean_obs,wrnum_obs,dlen_obs,freq_plot_title,freq_plot_outname)

    if ("PLOTFREQ" in steps_list_fcst):
        if not ("TIMEFREQ" in steps_list_fcst):
            raise Exception('Must run forecast Frequency calculation before plotting the frequencies.')
        freq_plot_title = os.environ.get('FCST_FREQ_PLOT_TITLE','Seasonal Cycle of WR Days/Week')
        freq_plot_outname = os.path.join(oplot_dir,os.environ.get('FCST_FREQ_PLOT_OUTPUT_NAME','fcst_freq'))
        # Compute mean
        wrmean_fcst = np.nanmean(wrfreq_fcst,axis=1)
        pwr.plot_wr_frequency(wrmean_fcst,wrnum_fcst,dlen_fcst,freq_plot_title,freq_plot_outname)


if __name__ == "__main__":
    main()

Running METplus

This use case is run in the following ways:

  1. Passing in UserScript_fcstGFS_obsERA_WeatherRegime.py then a user-specific system configuration file:

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

    run_metplus.py -c /path/to/METplus/parm/use_cases/model_applications/s2s_mid_lat/UserScript_fcstGFS_obsERA_WeatherRegime.py
    

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_mid_lat/WeatherRegime (relative to OUTPUT_BASE) and will contain output for the steps requested. This may include the regridded data, daily averaged files, a text file containing the list of input files, and text files for the weather regime classification and time frequency (if KMEANS and TIMEFREQ are run for both the forecast and observation data). In addition, output elbow, EOF, and Kmeans weather regime plots can be generated. The location of these output plots can be specified as WR_OUTPUT_DIR. If it is not specified, plots will be sent to {OUTPUT_BASE}/plots. The output location for the matched pair files can be specified as WR_MPR_OUTPUT_DIR. If it is not specified, it will be sent to {OUTPUT_BASE}/mpr. The output weather regime text or netCDF file location is set in WR_OUTPUT_FILE_DIR. If this is not specified, the output text/netCDF file will be sent to {OUTPUT_BASE}. The stat_analysis contingency table statistics and anomaly correlation files will be sent to the locations given in STAT_ANALYSIS_OUTPUT_DIR for their respective configuration sections.

Keywords

Note

  • RegridDataPlaneToolUseCase

  • PCPCombineToolUseCase

  • StatAnalysisToolUseCase

  • S2SAppUseCase

  • S2SMidLatAppUseCase

  • NetCDFFileUseCase

  • GRIB2FileUseCase

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

sphinx_gallery_thumbnail_path = ‘_static/s2s_mid_lat-UserScript_fcstGFS_obsERA_WeatherRegime.png’

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

Gallery generated by Sphinx-Gallery