Note
Go to the end to download the full example code.
UserScript: Calculate Weather Regimes for ERA
model_applications/s2s_mid_lat/UserScript_obsERA_obsOnly_WeatherRegime.py
Scientific Objective
Weather regimes are defined as atmospheric patterns with high likelihoods of recurrence. They can modulate many atmospheric phenomena including tornado and severe weather occurrence. This use case computes the top 6 most frequent weather regimes using K-means clustering for the ERA observations. It also computes the frequency of occurrence of each weather regime over a 7 day period. The code for computing weather regimes comes from Douglas Miller.
Miller, D. E., Wang, Z., Trapp, R. J., &Harnos, D. S., 2020: Hybrid prediction of weekly tornado activity out to Week 3: Utilizing weather regimes. Geophysical Research Letters, 47, https://doi.org/10.1029/2020GL087253.
Version Added
METplus version 4.0.0
Datasets
Forecast: None
Observation: ERA Reanlaysis 500 mb height for DJF 1979 - 2017
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 calles UserScript once and Stat-Analysis twice. There are two optional pre-processing steps, Regrid-Data-Plane and PCP-Combine. Additionally, METcalcpy and METplotpy are required to run this use case. The METcalcpy scripts accessed include the following:
metcalcpy/contributed/blocking_weather_regime/WeatherRegime.py
metcalcpy/contributed/blocking_weather_regime/Blocking_WeatherRegime_util.py
The METplopty scripts accessed include the following:
metplotpy/contributed/weather_regime/plot_weather_regime.py
METplus Workflow
Beginning time (VALID_BEG): 12-01-2000
End time (VALID_END): 02-28-2017
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 that runs the weather regime driver script is run once over the entire time period. The weather regime driver script performs the weather regime calculation for the observations. This calculation is divided up into steps, which the user can select by setting STEPS_OBS in the [user_env_vars] section of the configuration. More information on the steps and how the calculation proceeds is given in the User Scripting section below.
The two optional pre-processing steps loop by valid time when they are turned on, with different timing settings needed for the different steps. These steps are turned off due to data size and processing time. The first optional step calls Regrid-Data-Plane to regrid the data to a 1 degree latitude/longitude grid. The second calls PCP-Combine to compute daily means. These omitted steps can be turned back on by using the PROCESS_LIST that is commented out:
PROCESS_LIST = RegridDataPlane(regrid_obs), PcpCombine(daily_mean_obs), UserScript(script_wr)
Settings for the optional pre-processing steps can be found in the respective sections of the configuration, regrid_obs and daily_mean_obs. Data is not provided in the tarball to run these steps, but the configurations are provided for reference on how to set up these calculations.
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_obsERA_obsOnly_WeatherREgime.conf
[config]
# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/s2s_mid_lat/UserScript_obsERA_obsOnly_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)
# Weather Regime Analysis only:
PROCESS_LIST = UserScript(script_wr)
###
# 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_END = 2017022800
VALID_INCREMENT = 86400
LEAD_SEQ = 0
# Only Process DJF
SKIP_TIMES = "%m:begin_end_incr(3,11,1)", "%m%d:0229"
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE_PER_LEAD
###
# RegridDataPlane(regrid_obs) Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#regriddataplane
###
# Regridding Pre-Processing Step
[regrid_obs]
VALID_BEG = 1979120100
VALID_END = 2017022818
VALID_INCREMENT = 21600
# REGRID_DATA_PLANE (Pre Processing Step 1), currently turned off
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
REGRID_DATA_PLANE_METHOD = BILIN
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 = 1979120118
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
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
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 = 38
DAYS_PER_SEASON = 90
# Variable for the Z500 data
OBS_WR_VAR = Z500
# Weather Regime Number
OBS_WR_NUMBER = 6
# Number of clusters
OBS_NUM_CLUSTERS = 20
# Number of principal components
OBS_NUM_PCS = 10
# 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
# 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
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
# EOF plot output name and contour levels
OBS_EOF_PLOT_OUTPUT_NAME = obs_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
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 (1979-2017)
OBS_FREQ_PLOT_OUTPUT_NAME = obs_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]
# 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
# 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
# Command to run the user script with input configuration file
USER_SCRIPT_COMMAND = {METPLUS_BASE}/parm/use_cases/model_applications/s2s_mid_lat/common/WeatherRegime_driver.py
MET Configuration
This case does not use MET configuration files.
Python Embedding
This use case does not use python embedding.
User Scripting
This use case runs WeatherRegime_driver.py. This driver script runs the selected steps of the weather regime calculation that are specified in OBS_STEPS in the [user_env_vars] section of the UserScript .conf file. All steps are run for this use case, as specified in the following format:
OBS_STEPS = ELBOW+PLOTELBOW+EOF+PLOTEOF+KMEANS+PLOTKMEANS+TIMEFREQ+PLOTFREQ
The possible steps are computing the elbow or optimal number of clusters (ELBOW), plotting the elbow (PLOTELBOW), computing EOFs (EOF), plotting EOFs (PLOTEOF), computing the weather regimes using K means clustering (KMEANS), plotting the weather regimes (PLOTKMEANS), computing a user specified time frequency of weather regimes (TIMEFREQ) and plotting the time frequency (PLOTFREQ). The TIMEFREQ and PLOTFREQ steps require that the KMEANS step be run first, while all other steps can be run individally. Input variables to the WeatherRegime driver are set and described in the [user_env_vars] section of the configuration file.
Elbow computes the optimal number of clusters using 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 computes empirical orthogonal functions. These EOFs are used to reconstruct the height field, with this reconstructed data used in the K means calculation. If EOFs are not compted, the original height field is used in the K means calculation. 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.
parm/use_cases/model_applications/s2s_mid_lat/common/WeatherRegime_driver.py
#!/usr/bin/env python3
import os
import warnings
import numpy as np
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 ValueError('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_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 ValueError('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_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 "ELBOW" not in steps_list_obs:
raise ValueError('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 "ELBOW" not in steps_list_fcst:
raise ValueError('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 "EOF" not in steps_list_obs:
raise ValueError('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 "EOF" not in steps_list_fcst:
raise ValueError('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 "KMEANS" not in steps_list_obs:
raise ValueError('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 "KMEANS" not in steps_list_fcst:
raise ValueError('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 "KMEANS" not in steps_list_obs:
raise ValueError('Must run observed Kmeans before running frequencies.')
wrfreq_obs,dlen_obs,_ = steps_obs.compute_wr_freq(wrc_obs)
if "TIMEFREQ" in steps_list_fcst:
if "KMEANS" not in steps_list_fcst:
raise ValueError('Must run forecast Kmeans before running frequencies.')
wrfreq_fcst,dlen_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')
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 "TIMEFREQ" not in steps_list_obs:
raise ValueError('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 "TIMEFREQ" not in steps_list_fcst:
raise ValueError('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
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_mid_lat/UserScript_obsERA_obsOnly_WeatherRegime.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_mid_lat/WeatherRegime and will contain output for the steps requested. The output includes 4 plots in the plots directory:
* obs_elbow.png
* obs_eof.png
* obs_kmeans.png
* obs_freq.png
The output also includes a daily classification of weather regimes as a text file:
* obs_weather_regime_class.txt
If the pre-processing steps are turned on, the output will include the regridded data and daily averaged files.
Keywords
Note
RegridDataPlaneToolUseCase
PCPCombineToolUseCase
StatAnalysisToolUseCase
S2SAppUseCase
S2SMidLatAppUseCase
UserScriptUseCase
NetCDFFileUseCase
GRIB2FileUseCase
METcalcpyUseCase
METplotpyUseCase
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.