Start by making an output directory for Ensemble-Stat and changing directories:

The behavior of Ensemble-Stat is controlled by the contents of the configuration file passed to it on the command line. The default Ensemble-Stat configuration file may be found in the data/config/EnsembleStatConfig_default file. The configurations used by the test script may be found in the scripts/config/EnsembleStatConfig* files. Prior to modifying the configuration file, users are advised to make a copy of the default:

Open up the EnsembleStatConfig_tutorial file for editing with your preferred text editor.

The configurable items for Ensemble-Stat are broken out into two sections. The first section specifies how the ensemble should be processed to derive summary fields, such as the ensemble mean and spread. The second section specifies how the ensemble should be verified directly, such as the computation of rank histograms and spread/skill. The configurable items include specifications for the following:

• Section 1: Ensemble Processing (ens dictionary)
  • The ensemble fields to be summarized at the specified vertical level or accumulation interval.
  • The threshold values to be applied in computing ensemble relative frequencies (e.g. the percent of ensemble members exceeding some threshold at each point).
  • Thresholds to specify how many of the ensemble members must actually be present with valid data.
• Section 2: Verification (fcst and obs dictionaries)
  • The forecast and observation fields to be verified at the specified vertical level or accumulation interval.
  • The matching time window for point observations.
  • The type of point observations to be matched to the forecasts.
  • The areas over which to aggregate statistics - as predefined grids or configurable lat/lon polylines.
  • The interpolation or smoothing methods to be used.

You may find a complete description of the configurable items in section 9.3.2 of the MET User's Guide. Please take some time to review them.

For this tutorial, we'll configure Ensemble-Stat to summarize and verify 24-hour accumulated precipitation. While we'll run Ensemble-Stat on a single field, please note that it may be configured to operate on multiple fields. The ensemble we're verifying consists of 6 members defined over the west coast of the United States. Edit the EnsembleStatConfig_tutorial file as follows:

• In the ens dictionary, set
  
     field = [
       {
         name       = "APCP";
         level      = [ "A24" ];
         cat_thresh = [ >0, >=5.0, >=10.0 ];
       }
     ];

  To read 24-hour accumulated precipitation from the input GRIB files and compute ensemble relative frequencies for the thresolds listed.

• In the fcst dictionary, set
  
     field = [
       {
         name     = "APCP";
         level      = [ "A24" ];
       }
     ];

  To also verify the 24-hour accumulated precipitation fields.

• In the fcst dictionary, set:

     message_type = [ "ADPSFC" ];

  To verify against surface observations.

   

• In the mask dictionary, set:
  
     grid = [ "FULL" ]; 

  To accumulate statistics over the full model domain.

• In the mask dictionary, set
  
     poly = [ "MET_BASE/poly/NWC.poly",
              "MET_BASE/poly/SWC.poly" ];

  To also verify over the northwest coast (NWC) and southwest coast (SWC) subregions.

• Set:
  
  output_flag = {
     ecnt  = BOTH;
     rhist = BOTH;
     phist = BOTH;
     orank = BOTH;
     ssvar = BOTH;
     relp  = BOTH;
  }

  To compute continuous ensemble statistics (ECNT), ranked histogram (RHIST), probability integral transform histogram (PHIST), observation ranks (ORANK), spread-skill variance (SSVAR), and relative position (RELP).

Save and close this file.

Next, run Ensemble-Stat on the command line using the following command, using wildcards to list the 6 input ensemble member files:

 The command above uses the ASCII2NC output of that was generated by make test when MET was compiled.

Ensemble-Stat is now performing the tasks we requested in the configuration file. Note that we've passed the input ensemble data directly on the command line by specifying the number of ensemble members (6) followed by their names using wildcards. We've also specified one gridded StageIV analysis field (-grid_obs) and one file containing point rain gauge observations (-point_obs) to be used in computing rank histograms. This tool should run pretty quickly.

When Ensemble-Stat is finished, it will have created 9 output files in the current directory: 6 ASCII statistics files (.stat, _ecnt.txt, _rhist.txt, _phist.txt, _orank.txt, _ssvar.txt , and _relp.txt ) , a NetCDF ensemble file (_ens.nc), and a NetCDF matched pairs file (_orank.nc).

The _ens.nc output from Ensemble-Stat is a NetCDF file containing the derived ensemble fields, one or more ASCII files containing statistics summarizing the verification performed, and a NetCDF file containing the gridded matched pairs.

All of the line types are written to the file ending in .stat. The Ensemble-Stat tool currently writes six output line types, ECNT, RHIST, PHIST, RELP, SSVAR, and ORANK.

1. The ECNT line type contains contains continuous ensemble statistics such as spread and skill. Ensemble-Stat uses assumed observation errors to compute both perturbed and unperturbed versions of these statistics. Statistics to which observation error have been applied can be found in columns which include the _OERR (for observation error) suffix.
2. The RHIST line type contains counts for a ranked histogram. This ranks each observation value relative to ensemble member values. Ideally, observation values would fall equally across all available ranks, yiedling a flat rank histogram. In practice, ensembles are often under-(U shape) or over-(inverted U shape) dispersive. In the event of ties, ranks are randomly assigned.
3. The PHIST line type contains counts for a probability integral transform histogram. This scales the observation ranks to a range of values between 0 and 1 and allows ensembles of different size to be compared. Similiary, when ensemble members drop out, RHIST lines cannot be aggregated together but PHIST lines can.
4. The RELP line is the relative position, which indicates how often each ensemble member's value was closest to the observation's value. In the event of ties, credit is divided equally among the tied members.
5. The ORANK line type is similar to the matched pair (MPR) output of Point-Stat. For each point observation value, one ORANK line is written out containing the observation value, it's rank, and the corresponding ensemble values for that point. When verifying against a griddedanalysis, the ranks can be written to the NetCDF output file.
6. The SSVAR line contains binned spread/skill information. For each observation location, the ensemble variance is computed at that point. Those variance values are binned based on the ens_ssvar_bin_size configuration setting. The skill is deteremined by comparing the ensemble mean value to the observation value. One SSVAR line is written for each bin summarizing the all the observation/ensemble mean pairs that it contains.

The STAT file contains all the ASCII output while the _ecnt.txt, _rhist.txt, _phist.txt, _orank.txt, _ssvar.txt, and _relp.txtfiles contain the same data but sorted by line type. Since so much data can be written for the ORANK line type, we recommend disabling the output of the optional text file using the output_flag parameter in the configuration file.

 Since the lines of data in these ASCII file are so long, we strongly recommend configuring your text editor to NOT use dynamic word wrapping. The files will be much easier to read that way.

Open up the ensemble_stat_20100101_120000V_rhist.txt RHIST file using the text editor of your choice and note the following:

• There are 6 lines in this output file resulting from using 3 verification regions in the VX_MASK column (FULL, NWC, and SWC) and two observations datasets in the OBTYPE column (ADPSFC point observations and gridded observations).
• Each line contains columns for the observations ranks (RANK_#) and a handful of ensemble statistics (CRPS, CRPSS, IGN, and SPREAD).
• There is output for 7 ranks - since we verified a 6-member ensemble, there are 7 possible ranks the observation values could attain.

Close this file, and open up the ensemble_stat_20100101_120000V_phist.txt PHIST file, and note the following:

• There are 5 lines in this output file resulting from using 3 verification regions (FULL, NWC, and SWC) and two observations datasets (ADPSFC point observations and gridded observations), where the ADPSFC point observations for the SWC region were all zeros for which the probability integral transform is not defined.
• Each line contains columns for the BIN_SIZE and counts for each bin. The bin size is set in the configuration file using the ens_phist_bin_size field. In this case, it was set to .05, therefore creating 20 bins (1/ens_phist_bin_size).

Close this file, and open up the ensemble_stat_20100101_120000V_orank.txt ORANK file, and note the following:

• This file contains 1866 lines, 1 line for each observation value falling inside each verification region (VX_MASK).
• Each line contains 44 columns... header information, the observation location and value, it's rank, and the 6 values for the ensemble members at that point.
• When there are ties, Ensemble-Stat randomly assigns a rank from all the possible choices. This can be seen in the SWC masking region where all of the observed values are 0 and the ensemble forecasts are 0 as well. Ensemble-Stat randomly assigns a rank between 1 and 7.

Close this file, and use the ncview utility to view the NetCDF ensemble fields file:

This file contains variables for the following:

1. Ensemble Mean
2. Ensemble Standard Deviation
3. Ensemble Mean minus 1 Standard Deviation
4. Ensemble Mean plus 1 Standard Deviation
5. Ensemble Minimum
6. Ensemble Maximum
7. Ensemble Range
8. Ensemble Valid Data Count
9. Ensemble Relative Frequency (for 3 thresholds)

The output of any of these summary fields may be disabled using the output_flag parameter in the configuration file.

Use the ncview utility to view the NetCDF gridded observation rank file:

This file is only created when you've verified using gridded observations and have requested its output using the output_flagparameter in the configuration file. Click through the variables in this file. Note that for each of the three verification areas (FULL, NWC, and SWC) this file contains 4 variables:

1. The gridded observation value
2. The observation rank
3. The probability integral transform
4. The ensemble valid data count

In ncview, the random assignment of tied ranks is evident in areas of zero precipitation.

Close this file.

Feel free to explore using this dataset. Some options to try are:

• Try setting skip_const = TRUE; in the config file to discard points where all ensemble members and the observaiton are tied (i.e. zero precip).
• Try setting obs_thresh = [ >0.01 ]; in the config file to only consider points where the observation meets this threshold. How does this differ from the using skip_const?
• Use wgrib to inventory the input files and add additional entries to the ens.field list. Can you process 10-meter U and V wind?

Stat-Analysis Tool: Configure

Start by making an output directory for Stat-Analysis and changing directories:

The behavior of Stat-Analysis is controlled by the contents of the configuration file or the job command passed to it on the command line. The default Stat-Analysis configuration may be found in the data/config/StatAnalysisConfig_default file. Let's start with a configuration file that's packaged with the met-8.0 test scripts:

Open up the STATAnalysisConfig_tutorial file for editing with your preferred text editor.

You will see that most options are left blank, so the tool will use whatever it finds or whatever is specified in the command or job line.  If you go down to the jobs[]section you will see a list of the jobs run for the test scripts. Remove those existing jobs and add the following 2 analysis jobs:

The first job listed above will select out only the contingency table count lines (CTC) where the threshold applied is >273.0 over the FULL masking region. This should result in 2 lines, one for pressure levels P850-500 and one for pressure P1050-850. So this job will be aggregating contingency table counts across vertical levels.

The second job listed above will perform the same aggregation as the first. However, it'll dump out the corresponding contingency table statistics derived from the aggregated counts.

Stat-Analysis Tool: Run on Point-Stat output

Now, run Stat-Analysis on the command line using the following command:

The output for these two jobs are printed to the screen. Try redirecting their output to a file by adding the -out command line argument:

The output was written to aggr_ctc_lines.out. We'll look at this file in the next section.

Next, try running the first job again, but entirely on the command line without a configuration file:

Note that we had to put double quotes (") around the forecast theshold string for this to work.

Next, run the same command but add the -dump_row command line option. This will redirect all of the STAT lines used by the job to a file. Also, add the -out_stat command line option. This will write a full STAT output file, including the 22 header columns:

Open up the file aggr_ctc_job.stat to see the 2 STAT lines used by this job.

Open up the file aggr_ctc_job_out.stat to see the 1 output STAT line. Notice that the FCST_LEV and OBS_LEV columns contain the input strings concatenated together.

Try re-running this job using -set_hdr FCST_LEV P1050-500 and -set_hdr OBS_LEV P1050-500. How does that affect the output?

The use of the -dump_row option is highly recommended to ensure that your analysis jobs run on the exact set of data that you intended. It's easy to make mistakes here!

Stat-Analysis Tool: Output

On the previous page, we generated the output file aggr_ctc_lines.out by using the -out command line argument. Open that file using the text editor of your choice, and be sure to turn word-wrapping off.

This file contains the output for the two jobs we ran through the configuration file. The output for each job consists of 3 lines as follows:

1. The JOB_LIST line contains the job filtering parameters applied for this job.
2. The COL_NAME line contains the column names for the data to follow in the next line.
3. The third line consists of the line type generated (CTC and CTS in this case) followed by the values computed for that line type.

Next, try running the Stat-Analysis tool on the output file ../point_stat/point_stat_run2_360000L_20070331_120000V.stat. Start by running the following job:

This job should aggregate 2 CTC lines for 2-meter temperature across the EAST and WEST regions. Next, try creating your own Stat-Analysis command line jobs to do the following:

1. Do the same aggregation as above but for the 5x5 interpolation output (i.e. 25 points instead of 1 point).
2. Do the aggregation listed in (1) but compute the corresponding contingency table statistics (CTS) line. Hint: you will need to change the job type to aggregate_stat and specify the desired -out_line_type.
   How do the scores change when you increase the number of interpolation points? Did you expect this?
3. Aggregate the scalar partial sums lines (SL1L2) for 2-meter temperature across the EAST and WEST masking regions.
   How does aggregating the East and West domains affect the output?
4. Do the aggregation listed in (3) but compute the corresponding continuous statistics (CNT) line. Hint: use the aggregate_stat job type.
5. Run an aggregate_stat job directly on the matched pair data (MPR lines), and use the -out_line_type command line argument to select the type of output to be generated. You'll likely have to supply additional command line arguments depending on what computation you request.

    

Now answer this question about this Stat-Analysis output:

1. How do the scores compare to the original (separated by level) scores? What information is gained by aggregating the statistics?

When doing the exercises above, don't forget to use the -dump_row command line option to verify that you're running the job over the STAT lines you intended.

If you get stuck on any of these exercises, you may refer to the exercise answers. We will return to the Stat-Analysis tool in the future practical sessions.

1. Job Number 1:
    
   stat_analysis \
   -lookin ../point_stat/point_stat_run2_360000L_20070331_120000V.stat -v 2 \
   -job aggregate -fcst_var TMP -fcst_lev Z2 -vx_mask EAST -vx_mask WEST -interp_pnts 25 -fcst_thresh ">278.0" \
   -line_type CTC \
   -dump_row job1_ps.stat
2. Job Number 2:
    
   stat_analysis \
   -lookin ../point_stat/point_stat_run2_360000L_20070331_120000V.stat -v 2 \
   -job aggregate_stat -fcst_var TMP -fcst_lev Z2 -vx_mask EAST -vx_mask WEST -interp_pnts 25 -fcst_thresh ">278.0" \
   -line_type CTC -out_line_type CTS \
   -dump_row job2_ps.stat
3. Job Number 3:
    
   stat_analysis \
   -lookin ../point_stat/point_stat_run2_360000L_20070331_120000V.stat -v 2 \
   -job aggregate -fcst_var TMP -fcst_lev Z2 -vx_mask EAST -vx_mask WEST -interp_pnts 25 \
   -line_type SL1L2 \
   -dump_row job3_ps.stat
4. Job Number 4:
    
   stat_analysis \
   -lookin ../point_stat/point_stat_run2_360000L_20070331_120000V.stat -v 2 \
   -job aggregate_stat -fcst_var TMP -fcst_lev Z2 -vx_mask EAST -vx_mask WEST -interp_pnts 25 \
   -line_type SL1L2 -out_line_type CNT \
   -dump_row job4_ps.stat
5. This MPR job recomputes contingency table statistics for 2-meter temperature over G212 using a new threshold of ">=285":
    
   stat_analysis \
   -lookin ../point_stat/point_stat_run2_360000L_20070331_120000V.stat -v 2 \
   -job aggregate_stat -fcst_var TMP -fcst_lev Z2 -vx_mask G212 -interp_pnts 25 \
   -line_type MPR -out_line_type CTS \
   -out_fcst_thresh ge285 -out_obs_thresh ge285 \
   -dump_row job5_ps.stat

Answers to Exercises from Session 3

These are the answers to the exercises from the previous page. Feel free to ask a MET representative if you have any questions!

ANSWER 3.1: accum_3hr - Build a 3 Hour Accumulation Instead of 6

Instructions: Modify the METplus configuration files to build a 3 hour accumulation instead of a 6 hour accumulation from forecast data using Pcp-Combine in the HREF MEAN vs. MRMS QPE example. Then compare 3 hour accumulations in the forecast and observation data with grid_stat.

Answer: In the user_config/qpf-ensemble.accum_3hr.conf file, change the following variables in the [config] section: 

Change:

To:

Change:

To:

ANSWER 3.2: input_1hr - Force Pcp-Combine to only use 1 hour accumulation files

Instructions: Modify the METplus configuration files to force Pcp-Combine to use six 1 hour accumulation files instead of one 6 hour accumulation file of observation data in the PHPT vs. StageIV GRIB example.

Answer: In the user_config/qpf-prob.input_1hr.conf file, add the following variable to the [config] section: