Note
Go to the end to download the full example code.
GridStat: Use binary observation field to verify percentile forecast
model_applications/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE.conf
Scientific Objective
Evaluation of a Probability of Exceedence (POE) field presents several difficulties. Some of these include a fitting verification statistic to report on, choosing a meaningful percentile field, and more. This use case was the culmination of attempting to verify a POE field for extreme temperature (defined as the 85th percentile) in METplus. In order to provide a streamlined process that didn’t require vast reworkings of the MET tools, the observation field was converted to binary: 0s indicating a non-85th percentile temperature was observed, and a 1 indicating the opposite. Those observations are compared to the chosen forecast percentile and the HSS_EC becomes the main statistical focus, as the new hss_ec_value feature allowed the use case to more closely replicate in-house verificaiton that already existed. A final note that because the POE forecast file is a non-standard netCDF, Python Embedding was used to extract the desired field
Version Added
METplus version 5.1
Datasets
Forecast: NOAA Global Ensemble Forecast System (GEFS) 85th percentile of maximum temperature
Observations: NOAA Climate Prediction Center Climate Assessment Data Base (CADB) converted into a binary field relative to the 85th percentile
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 calls a Python script to extract the user-defined percentile forecast. METplus then verifies it against a binary observation field in GridStat and returns the requested line type outputs.
METplus Workflow
Beginning time (INIT_BEG): 20220522
End time (INIT_END): 20220522
Increment between beginning and end times (INIT_INCREMENT): 12H
Sequence of forecast leads to process (LEAD_SEQ): 8d
The following boundary time is used for the entire script:
There is only one time processed for the use case.
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/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE.conf
[config]
PROCESS_LIST = GridStat
LOOP_BY = INIT
INIT_TIME_FMT = %Y%m%d
INIT_BEG=20220522
INIT_END=20220522
INIT_INCREMENT = 12H
LEAD_SEQ = 8d
###
# File I/O
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#directory-and-filename-template-info
###
FCST_GRID_STAT_INPUT_TEMPLATE = PYTHON_NUMPY
OBS_GRID_STAT_INPUT_DIR = {INPUT_BASE}/model_applications/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE
OBS_GRID_STAT_INPUT_TEMPLATE = tmax_cats_{valid?fmt=%Y%m%d}.nc
GRID_STAT_CLIMO_MEAN_INPUT_DIR =
GRID_STAT_CLIMO_MEAN_INPUT_TEMPLATE =
GRID_STAT_CLIMO_STDEV_INPUT_DIR =
GRID_STAT_CLIMO_STDEV_INPUT_TEMPLATE =
GRID_STAT_OUTPUT_DIR = {OUTPUT_BASE}/model_applications/POE_tmax
GRID_STAT_OUTPUT_TEMPLATE =
###
# Field Info
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#field-info
###
MODEL = GEFS
OBTYPE = Obs
GRID_STAT_ONCE_PER_FIELD = False
FCST_IS_PROB = false
#FCST_GRID_STAT_PROB_THRESH = ==0.1
FCST_VAR1_NAME = {PARM_BASE}/use_cases/model_applications/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE/Tmax_fcst_embedded.py {INPUT_BASE}/model_applications/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE/gefs-00z_rfcst-cal_tmax_{init?fmt=%Y%m%d}_day8.nc:85
FCST_VAR1_THRESH = >=0.2
OBS_VAR1_NAME = tmax
OBS_VAR1_LEVELS = "(*,*)"
OBS_VAR1_THRESH = gt0.5
###
# GridStat Settings (optional)
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#gridstat
###
#LOG_GRID_STAT_VERBOSITY = 2
GRID_STAT_CONFIG_FILE = {PARM_BASE}/met_config/GridStatConfig_wrapped
#FCST_GRID_STAT_FILE_TYPE =
OBS_GRID_STAT_FILE_TYPE = NETCDF_NCCF
GRID_STAT_REGRID_TO_GRID = FCST
GRID_STAT_DESC = NA
FCST_GRID_STAT_FILE_WINDOW_BEGIN = 0
FCST_GRID_STAT_FILE_WINDOW_END = 0
OBS_GRID_STAT_FILE_WINDOW_BEGIN = 0
OBS_GRID_STAT_FILE_WINDOW_END = 0
GRID_STAT_NEIGHBORHOOD_WIDTH = 1
GRID_STAT_NEIGHBORHOOD_SHAPE = SQUARE
GRID_STAT_NEIGHBORHOOD_COV_THRESH = >=0.5
GRID_STAT_OUTPUT_PREFIX =
#GRID_STAT_OUTPUT_FLAG_FHO = NONE
GRID_STAT_OUTPUT_FLAG_CTC = BOTH
GRID_STAT_OUTPUT_FLAG_CTS = BOTH
#GRID_STAT_OUTPUT_FLAG_MCTC = NONE
#GRID_STAT_OUTPUT_FLAG_MCTS = NONE
#GRID_STAT_OUTPUT_FLAG_CNT = NONE
#GRID_STAT_OUTPUT_FLAG_SL1L2 = NONE
#GRID_STAT_OUTPUT_FLAG_SAL1L2 = NONE
#GRID_STAT_OUTPUT_FLAG_VL1L2 = NONE
#GRID_STAT_OUTPUT_FLAG_VAL1L2 = NONE
#GRID_STAT_OUTPUT_FLAG_VCNT = NONE
#GRID_STAT_OUTPUT_FLAG_PCT = NONE
#GRID_STAT_OUTPUT_FLAG_PSTD = NONE
#GRID_STAT_OUTPUT_FLAG_PJC = NONE
#GRID_STAT_OUTPUT_FLAG_PRC = NONE
#GRID_STAT_OUTPUT_FLAG_ECLV = BOTH
#GRID_STAT_OUTPUT_FLAG_NBRCTC = NONE
#GRID_STAT_OUTPUT_FLAG_NBRCTS = NONE
#GRID_STAT_OUTPUT_FLAG_NBRCNT = NONE
#GRID_STAT_OUTPUT_FLAG_GRAD = BOTH
#GRID_STAT_OUTPUT_FLAG_DMAP = NONE
GRID_STAT_NC_PAIRS_FLAG_LATLON = FALSE
GRID_STAT_NC_PAIRS_FLAG_RAW = FALSE
GRID_STAT_NC_PAIRS_FLAG_DIFF = FALSE
GRID_STAT_NC_PAIRS_FLAG_CLIMO = FALSE
#GRID_STAT_NC_PAIRS_FLAG_CLIMO_CDP = FALSE
#GRID_STAT_NC_PAIRS_FLAG_WEIGHT = FALSE
#GRID_STAT_NC_PAIRS_FLAG_NBRHD = FALSE
#GRID_STAT_NC_PAIRS_FLAG_FOURIER = FALSE
#GRID_STAT_NC_PAIRS_FLAG_GRADIENT = FALSE
#GRID_STAT_NC_PAIRS_FLAG_DISTANCE_MAP = FALSE
GRID_STAT_NC_PAIRS_FLAG_APPLY_MASK = FALSE
GRID_STAT_HSS_EC_VALUE = 0.15
#GRID_STAT_MASK_GRID =
GRID_STAT_MASK_POLY = {MET_INSTALL_DIR}/share/met/poly/CONUS.poly, {MET_INSTALL_DIR}/share/met/poly/EAST.poly, {MET_INSTALL_DIR}/share/met/poly/WEST.poly, {MET_INSTALL_DIR}/share/met/poly/NPL.poly, {MET_INSTALL_DIR}/share/met/poly/NEC.poly
MET Configuration
METplus sets environment variables based on user settings in the METplus configuration file. See How METplus controls MET config file settings for more details.
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 currently not supported by METplus you’d like to control, please refer to: Overriding Unsupported MET config file settings
GridStatConfig_wrapped
////////////////////////////////////////////////////////////////////////////////
//
// Grid-Stat configuration file.
//
// For additional information, see the MET_BASE/config/README file.
//
////////////////////////////////////////////////////////////////////////////////
//
// Output model name to be written
//
// model =
${METPLUS_MODEL}
//
// Output description to be written
// May be set separately in each "obs.field" entry
//
// desc =
${METPLUS_DESC}
//
// Output observation type to be written
//
// obtype =
${METPLUS_OBTYPE}
////////////////////////////////////////////////////////////////////////////////
//
// Verification grid
//
// regrid = {
${METPLUS_REGRID_DICT}
////////////////////////////////////////////////////////////////////////////////
//censor_thresh =
${METPLUS_CENSOR_THRESH}
//censor_val =
${METPLUS_CENSOR_VAL}
//cat_thresh =
${METPLUS_CAT_THRESH}
cnt_thresh = [ NA ];
cnt_logic = UNION;
wind_thresh = [ NA ];
wind_logic = UNION;
eclv_points = 0.05;
//nc_pairs_var_name =
${METPLUS_NC_PAIRS_VAR_NAME}
nc_pairs_var_suffix = "";
//hss_ec_value =
${METPLUS_HSS_EC_VALUE}
rank_corr_flag = FALSE;
//
// Forecast and observation fields to be verified
//
fcst = {
${METPLUS_FCST_FILE_TYPE}
${METPLUS_FCST_FIELD}
${METPLUS_FCST_CLIMO_MEAN_DICT}
${METPLUS_FCST_CLIMO_STDEV_DICT}
}
obs = {
${METPLUS_OBS_FILE_TYPE}
${METPLUS_OBS_FIELD}
${METPLUS_OBS_CLIMO_MEAN_DICT}
${METPLUS_OBS_CLIMO_STDEV_DICT}
}
////////////////////////////////////////////////////////////////////////////////
//
// Climatology mean data
//
//climo_mean = {
${METPLUS_CLIMO_MEAN_DICT}
//climo_stdev = {
${METPLUS_CLIMO_STDEV_DICT}
//
// May be set separately in each "obs.field" entry
//
//climo_cdf = {
${METPLUS_CLIMO_CDF_DICT}
////////////////////////////////////////////////////////////////////////////////
//
// Verification masking regions
//
// mask = {
${METPLUS_MASK_DICT}
////////////////////////////////////////////////////////////////////////////////
//
// Confidence interval settings
//
ci_alpha = [ 0.05 ];
boot = {
interval = PCTILE;
rep_prop = 1.0;
n_rep = 0;
rng = "mt19937";
seed = "";
}
////////////////////////////////////////////////////////////////////////////////
//
// Data smoothing methods
//
//interp = {
${METPLUS_INTERP_DICT}
////////////////////////////////////////////////////////////////////////////////
//
// Neighborhood methods
//
nbrhd = {
field = BOTH;
// shape =
${METPLUS_NBRHD_SHAPE}
// width =
${METPLUS_NBRHD_WIDTH}
// cov_thresh =
${METPLUS_NBRHD_COV_THRESH}
vld_thresh = 1.0;
}
////////////////////////////////////////////////////////////////////////////////
//
// Fourier decomposition
// May be set separately in each "obs.field" entry
//
//fourier = {
${METPLUS_FOURIER_DICT}
////////////////////////////////////////////////////////////////////////////////
//
// Gradient statistics
// May be set separately in each "obs.field" entry
//
//gradient = {
${METPLUS_GRADIENT_DICT}
////////////////////////////////////////////////////////////////////////////////
//
// Distance Map statistics
// May be set separately in each "obs.field" entry
//
//distance_map = {
${METPLUS_DISTANCE_MAP_DICT}
////////////////////////////////////////////////////////////////////////////////
// Threshold for SEEPS p1 (Probability of being dry)
//seeps_p1_thresh =
${METPLUS_SEEPS_P1_THRESH}
////////////////////////////////////////////////////////////////////////////////
//
// Statistical output types
//
//output_flag = {
${METPLUS_OUTPUT_FLAG_DICT}
//
// NetCDF matched pairs output file
// May be set separately in each "obs.field" entry
//
// nc_pairs_flag = {
${METPLUS_NC_PAIRS_FLAG_DICT}
////////////////////////////////////////////////////////////////////////////////
//ugrid_dataset =
${METPLUS_UGRID_DATASET}
//ugrid_max_distance_km =
${METPLUS_UGRID_MAX_DISTANCE_KM}
//ugrid_coordinates_file =
${METPLUS_UGRID_COORDINATES_FILE}
////////////////////////////////////////////////////////////////////////////////
//grid_weight_flag =
${METPLUS_GRID_WEIGHT_FLAG}
tmp_dir = "${MET_TMP_DIR}";
// output_prefix =
${METPLUS_OUTPUT_PREFIX}
////////////////////////////////////////////////////////////////////////////////
${METPLUS_TIME_OFFSET_WARNING}
${METPLUS_MET_CONFIG_OVERRIDES}
Python Embedding
This use case calls a Python script to parse the user-requested percentile from the forecast dataset. This is controlled in the forecast VAR1 variable setting and is provided in:
parm/use_cases/model_applications/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE/Tmax_fcst_embedded.py
import sys
import datetime as dt
import numpy as np
from netCDF4 import Dataset
try:
#user will input name of the file, as well as a percentile they're interested in
#in future iteration, this may need to change to multiple percentiles (a la list style)
print("1")
input_file,ptile = sys.argv[1].split(':')
f = Dataset(input_file, 'r')
print("2")
v = f['poe']
val_time = f.valid_date_range[1]
val_time = dt.datetime.strptime(val_time,"%Y%m%d")
ini_time = str(input_file.split('_')[-2])
ini_time = dt.datetime.strptime(ini_time,"%Y%m%d")
print("3")
lead, rem = divmod((val_time - ini_time).total_seconds(), 3600)
ptile_ind = np.where(f['ptile'][:] == int(ptile))[0][0]
print("4")
lat = np.float64(f.variables['latitude'][:])
lon = np.float64(f.variables['longitude'][:])
#var = np.float64(v[0,ptile_ind,:,:],fill_value=-9999.)
var = np.float64(v[0,ptile_ind,:,:])
print(np.amax(var),np.amin(var))
met_data = var.copy()
except NameError:
print("Can't find input file")
sys.exit(1)
#ADDED
#for i in range(len(met_data)):
# for j in range(len(met_data[i])):
# if j <=2 or j >=358:
# print("edge of ", met_data[i,j],"at", lat[i],lon[j])
# if lat[i] >=42.0 and lat[i] <= 46.0:
# if lon[j] >= 235.0 and lon[j] <= 239.0:
# print("found",met_data[i,j]," at ",lat[i],lon[j],i,j)
attrs = {
'valid': str(val_time.strftime("%Y%m%d"))+'_000000',
'init': str(ini_time.strftime("%Y%m%d"))+'_000000',
'name': 'poe_P'+str(ptile),
'long_name': v.long_name,
'lead': str(int(lead)),
'accum': '00',
'level': 'SURFACE',
'units': 'PERCENTILES',
'grid': {
'name': 'Global 1 degree',
'type': 'LatLon',
'lat_ll': -90.0,
'lon_ll': 0.0,
'delta_lat': 1.0,
'delta_lon': 1.0,
'Nlon': f.dimensions['longitude'].size,
'Nlat': f.dimensions['latitude'].size,
}
}
#print output for user to show successful run
print("Input file: " + repr(input_file.split('/')[-1]))
print("Attributes:\t"+ repr(attrs))
f.close()
For more information on the basic requirements to utilize Python Embedding in METplus, please refer to the MET User’s Guide section on Python embedding.
User Scripting
User Scripting is not used in this use case.
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/medium_range/GridStat_fcstGEFS_obsCADB_BinaryObsPOE.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 the use case will be found in model_applications/POE_tmax (relative to OUTPUT_BASE). The following files should exist:
grid_stat_1920000L_20220530_000000V_ctc.txt
grid_stat_1920000L_20220530_000000V_cts.txt
grid_stat_1920000L_20220530_000000V.stat
Keywords
Note
GridStatUseCase
PythonEmbeddingFileUseCase
MediumRangeAppUseCase
NETCDFFileUseCase
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.