Deterministic Scalar Statistics

Most meteorological forecasts would be described as non-probabilistic, meaning the forecast value given is provided with no additional information of certainty in that value. When verifying these deterministic forecasts, the only important factor is whether or not the event occurred: the assumed certainty in the forecast is 100%.

Numerous computationally-easy (and very popular) scalar statistics are within reach without too much manipulation of a contingency table’s counts.

Accuracy (Acc)

The scalar attribute of Accuracy can be found as a simple ratio between the forecasts that correctly predicted the event and the total number of occurrences and non-occurrences, n. In equation format,

It’s very easily computed and does address how often a forecast is correctly predicting an event. As most statistical resources will warn you, however, this measure should be used with caution, especially for an event that happens only rarely. The Finley tornado forecasts (1884) is an excellent example of this, with a 96.6% Accuracy for predicting a tornado due to the overwhelming count of correct negatives. Peers were quick to point out that a higher Accuracy (98.2%) could have been achieved with a persistence forecast of No Tornado! See how to use this statistic in METplus!

Probability of Detection (POD)

Probability of Detection (POD), also referred to as the Hit Rate, is a measure of the forecast’s discrimination. Rather than finding a ratio of the forecasts to the entire occurrence and non-occurrence count, it focuses on only the times the forecast predicted an event would occur. This becomes

and can be useful for rare events (tornadoes, 100-year floods, etc.) as it will be penalized (i.e. go toward 0) for too many missed forecasts. See how to use this statistic in METplus! 

Probability of False Detection (POFD)

If there’s a POD, there should also be a probability of false detection, too. The POFD, or false alarm rate, is a discrimination statistic of the frequency of only false alarm forecasts to the count of an event that does not occur.

See how to use this statistic in METplus!

Frequency bias (Bias)

Frequency bias (a measure of, you guessed it, bias!) looks at the count of “yes” forecasts and compares them to the “yes” count of events being observed.

This relative frequency does not hold any information on individual forecasts and their corresponding observation values, so it is not a statistic for understanding Accuracy. See how to use this statistic in METplus!

False Alarm Ratio (FAR)

The False Alarm Ratio (FAR) tackles both the reliability and resolution attributes of a forecast. It looks at the ratio of “yes” forecasts that did not occur to the total number of times a “yes” forecast was made.

This also marks the first statistic covered in this tutorial that has a negative orientation: a FAR of 0 indicates the forecast has high reliability and resolution, while a FAR of 1 indicates little to no reliability or resolution. See how to use this statistic in METplus! 

Critical Success Index (CSI)

The Critical Success Index (CSI), also known as the Threat Score, reflects the accuracy of the forecast. Accuracy pertains to the agreement of individual forecast-observation pairs, so CSI can be calculated with

Note how CSI can be considered a ratio between the times the forecast correctly called for an event and the total times the forecast called for an event or the event was observed. So a CSI of 1 would indicate a highly accurate forecast, while a 0 indicates no accuracy. See how to use this statistic in METplus!

Deterministic Skill Scores

Skill scores can be a more rounded way of speaking to a forecast’s quality. They often combine aspects of the previously-listed scalar statistics and can serve as a starting point for creating your own skill score that is better suited to your forecasts’ properties.

Three of the most popular scores are Heidke Skill Score (HSS), the Hanssen-Kuipers Discriminant (HK), and the Gilbert Skill Score (GSS) and it would be worthwhile to mention them here.

Heidke Skill Score (HSS)

The HSS is a measure of the proportion correct that would be achieved by a “random” forecast, which is located in the denominator of the HSS equation. In this instance, random denotes a forecast that is completely independent of the observation dataset. In practice, the random forecast is often a climatology or persistence forecast. By combining the probability of a correct “yes” forecast (i.e. a hit) with the probability of a correct “no” forecast (i.e. a correct rejection) the resulting equation is

HSS can range from -1 to 1, with a perfect forecast receiving a score of 1. The equation presented above is a compact version which uses a grouped term, C2. The C2 term expands to

This is an important version of HSS, as METplus also calculates a modified HSS, in addition to the more traditional approach, that allows users to control what the C2 term is. See how to use these skill scores in METplus!

Hanssen-Kuipers Discriminant (HK)

HK is known by several names, including the Peirce Skill Score and the true skill statistic. How ever it is referred to, the theory and equation of HK remains the same. This skill score is similar to HSS (ranges from -1 to 1, perfect forecast is 1, etc.), but instead of a comparison to a truly random forecast, HK is formulated for a random forecast that shows no bias in the denominator. This changes the focus of the verification away from the skill of the forecast compared to random chance, and instead shows how well the forecast delineated between observed “yes” events and observed “no” events. That appears in equation form as

which is the expanded version of the POD minus the POFD. Because of its basis in POD, HK can be similarly affected by infrequent events and is suggested as a more useful skill score for frequent events. See how to use this skill score in METplus!

Gilbert Skill Score (GSS)

Finally, GSS looks at the correspondence between a forecasted event and the observed “yes” event. Sometimes called the Equitable Threat Score, the GSS is a good option for those forecasted events where the observed “yes” event is rare. This is due to its omission of correct negatives in the equation, which, with a rare observed “yes” event, would be fairly large. The GSS is given as

And ranges from -1 to 1, with a perfect forecast receiving a score of 1. Similar to HSS, a compact version of GSS is presented using the C1 term. This term expands to

 See how to use this skill score in METplus!

METplus solutions for Deterministic forecast verification

Now that you know a bit more about dichotomous, deterministic forecasts and how to extract information on the scalar attributes through statistics, it’s time to show how you can access those same statistics in METplus!

MET solutions

The MET User’s Guide provides an Appendix that dives into each and every statistical measure that it calculates, as well as the line type it is a part of. Statistics are grouped together by application and category and are available to METplus users in line types. For example, many of the statistics that were discussed above can be found in the Contingency Table Statistics (CTS) line type, which logically groups together statistics based directly on contingency table counts. In fact, MET allows users to access the direct contingency table counts through the aptly named Contingency Table Counts (CTC) line type.

For MET, what line types that are output depend on your selection of the appropriate line type using the output_flag dictionary. Note that certain line types may or may not be available in every tool: for example, both Point-Stat and Grid-Stat produce CTS line types, which allows users access to the various contingency table statistics for both point-based observations and gridded observations. But Ensemble-Stat is the only tool that can generate a Ranked Probability Score (RPS) line type which contains statistics relevant to the analysis of ensemble forecasts. If you don’t see your desired statistic in the line type or tool you’d expect it to be in, be sure to check the Appendix to see if the statistic is available in MET and which line type it’s currently grouped with.

As for the previous statistics that were discussed, here’s a link to the User’s Guide Appendix entry that discusses its use in MET:

• Accuracy
• POD
• POFD
• Bias
• FAR
• CSI
• HSS
• HSS_EC (HSS Expected Correct)
• HK
• GSS

Remember that for categorical statistics, including those that deal with probabilistic datasets, you will need to provide an appropriate threshold that divides the observations and forecasts into two mutually exclusive categories. For more information on the available thresholding options, please review this section of the MET User’s Guide.

METplus Wrapper Solutions

Those same statistics are also available from the METplus wrappers. To better understand how MET configuration options for statistics translate to METplus wrapper configuration options, you can utilize the Statistics and Diagnostics Section of the METplus wrappers User’s Guide, which lists all of the available statistics through the wrappers, including what tools can output what statistics. To access the line type through the tool, find your desired tool in the list of available commands for that tool. Once you do, you’ll see the tool will have several options that contain _OUTPUT_FLAG_. These will exhibit the same behavior and accept the same settings as the line types in MET’s output_flag dictionary, so be sure to review the available settings to get the line type output you want.

MET example of Deterministic forecast verification

Here is an example that demonstrates deterministic forecast verification in MET.

For this example, let’s examine Grid-Stat. Assume we wanted to verify a dichotomous temperature forecast of greater than 86 degrees Fahrenheit. Starting with the general Grid-Stat configuration file, the following would resemble minimum necessary settings/changes for the fcst and obs dictionaries:

We can see that the forecast field name in the forecast input file is named TMP, and is set accordingly in the fcst dictionary. Similarly, the Z0 level is used to grab the lowest (0th) vertical level the TMP variable appears on. Finally, cat_thresh, which controls the categorical threshold that the contingency table will be created with, is set to greater than 86.0. This assumes that the temperature units in the file are in Fahrenheit. The obs dictionary is simply copying the settings from the fcst dictionary, which is a method that can be used if both the forecast and observation input files share the same variable structure (e.g. both inputs use the TMP variable name, in Fahrenheit, with the lowest vertical level being the desired verification level).

Now all that’s necessary would be to adjust the output_flag dictionary settings to have Grid-Stat print out the desired line types:

In this example, we have told MET to output the CTC and CTS line types, which will contain all of the scalar statistics that were discussed in this section. Running this set up would produce one .stat file with the two line types that were selected, CTC and CTS. The CTC line would look something like:

While the stat file full header column contents are discussed in the User’s Guide, the CTC line types are the final 6 columns of the line, beginning after the “CTC” column. The first value is MET’s TOTAL column which is the “total number of matched pairs”. You might better recognize this value as n, the summation of every cell in the contingency table. In fact, the following four columns of the CTC line type are synonymous with the contingency table terms, which have their corresponding MET terms provided in this table for your convenience:

Further descriptions of each of the CTC columns can be found in the MET User’s Guide. Note that the final column of the CTC line type, EC_VALUE, is only relevant to users verifying probabilistic data with the HSS_EC skill score.

The CTS line type is also present in the .stat file and is the second row. It has many more columns than the CTC line, where all of the scalar statistics and skill scores discussed previously are located. Focusing on the first few columns of the example output, you would find:

These columns can be understood by reviewing the MET User’s Guide guidance for CTS line type. After the familiar TOTAL or n column, we find statistics such as Base Rate, forecast mean, Accuracy, plus many more, all with their appropriate lower and upper confidence intervals and the bootstrap confidence intervals. Note that because the bootstrap library’s n_rep variable was kept at its default value of 0, bootstrap methods were not used and appear as NA in the stat file. While all of these statistics could be obtained from the CTC line type values with additional post-processing, the simplicity of having all of them already calculated and ready for additional group statistics or to advise forecast adjustments is one of the many advantages of using the METplus system.

METplus wrapper example of Deterministic forecast verification

To achieve the same success as the previous example but utilizing METplus wrappers instead of MET, very few differences would need to be made. Starting with the standard GridStat configuration file, we would need to set the _VAR1 settings appropriately:

Note how the BOTH option is utilized here (as opposed to individual FCST_ and OBS_ settings) since the forecast and observation datasets utilize the same name and level information. Because the loop/timing information is controlled inside the configuration file for METplus wrappers (as opposed to MET’s non-looping option), that information must also be set accordingly:

Finally, the desired line types need to be selected for output. In the wrappers, that looks like this:

After a successful run of METplus, the same .stat output file that was created in the MET example would be produced here, complete with CTC and CTS line type rows.

Prerequisites: Software

The Requirements section in the Software Installation chapter of the METplus User's Guide lists the software and Python packages that are required to run the METplus wrappers. Note that there is a core set of requirements needed to run the METplus wrappers and additional requirements needed to utilize some of the more advanced features.

Prerequisites: Environment

The following instructions are required so the commands in this tutorial can be copied and run without modification. The steps involve creating a working directory to store files used/generated by the tutorial and configuring a simple script that can be run to set up the shell environment. Once configured correctly, the script can be run upon returning to the tutorial content to easily resume progress.

Pre-Configured Environments

Setting up the Tutorial Environment on Hera (NOAA)

Setting up the Tutorial Environment on Jet (NOAA)

Setting up the Tutorial Environment on Cheyenne (NCAR)

Setting up the Tutorial Environment on Seneca (NCAR)

User Configured Environments

Setting up the Tutorial Environment (bash)

Setting up the Tutorial Environment (csh)

Verify Environment is Set Correctly

Run the Tutorial Setup Script

The tutorial setup script sets the paths for METPLUS_TUTORIAL_DIR, METPLUS_BUILD_BASE, MET_BUILD_BASE, and METPLUS_DATA. It also appends the $PATH environment variable to include the directory where the METplus scripts are located. If necessary, it may also load modules needed for the METplus software to run correctly.

Check Path

Make sure that all of the environment variables are set to the appropriate values and that the path is set up to locate the METplus components.

You should see the usage statement for Point-Stat. The version number listed should correspond to the version listed in MET_BUILD_BASE. If it does not, you will need to either reload the met module, or add ${MET_BUILD_BASE}/bin to your PATH.

$METPLUS_TUTORIAL_DIR

The directory you created to store all of your tutorial files

Example value:

Example contents:

$MET_BUILD_BASE

The directory where MET is installed

Example value:

Example contents

Check contents of MET bin directory

Example contents:

$METPLUS_BUILD_BASE

The directory where METplus is installed

Example value:

Example contents:

$METPLUS_DATA

The directory containing sample input data to use for the tutorial

Example value:

Example contents:

METplus directory structure

A brief description and overview of the METplus/ directory structure can be found in the METplus User's Guide section called METplus Wrappers Directory Structure. The files/directories in ${METPLUS_BUILD_BASE} should match this list.

METplus default configuration file

Look inside the directory ${METPLUS_BUILD_BASE}/parm

Look at the METplus default configuration file:

The METplus default configuration file (defaults.conf) is always read first. Any additional configuration files passed in on the command line are then processed in the order in which they are specified. This allows for each successive conf file the ability to override variables defined in any previously processed conf files. It also allows for defining and setting up conf files from a general (settings used by all use cases, ie. MET install dir) to more specific (Plot type when running track and intensity plotter) structure. The idea is to created a hiearchy of conf files that is easier to maintain, read, and manage. It is important to note, running METplus creates a single configuration file, which can be viewed to understand the result of all the conf file processing.

When METplus is run, the final metplus conf file is generated here:
metplus_runtime.conf:METPLUS_CONF={OUTPUT_BASE}/metplus_final.conf.YYYYMMDDHHmmss

This is a unique final config file for each METplus run. Use this file to see the result of all the conf file processing; this can be very helpful when troubleshooting.

NOTE: The syntax for METplus configuration files MUST include a "[config]" section header with the variable names and values on subsequent lines.

More information about the default configuration variables can be found in the METplus User's Guide section called Default Configuration File.
 

The met_config directory (in ${METPLUS_BUILD_BASE}/parm) contains "wrapped" MET configuration files that are used by calls to the MET applications via the METplus wrappers. The wrappers set environment variables that control settings in the wrapped MET configuration files through these environment variables. See the METplus User's Guide section called How METplus controls MET configuration variables for more information.

METplus Use Cases

The use_cases directory - this is where the use cases you will be running exist. Under the use_cases directory are two directories: met_tool_wrapper and model_applications. The met_tool_wrapper directory contains use cases that run a single METplus wrapper. They are a good starting point to see how the wrapper scripts generate commands that run the MET tools. The model_applications directory contains more complex use cases that often run multiple wrappers and demonstrate real evaluations from users.

MET Tool Wrapper Use Cases

Look at the MET Tool Wrapper Use Cases

The met_tool_wrapper use case files are organized into subdirectories by wrapper, e.g. Example or GridStat. They contain METplus configuration files (ending with .conf). If a MET tool that is called by a use case uses a MET configuration file, the file used for the met_tool_wrapper use cases is found in parm/met_config.

• met_tool_wrapper/Example - directory
• met_tool_wrapper/Example/Example.conf - use case configuration file
• met_tool_wrapper/GridStat - directory
• met_tool_wrapper/GridStat/GridStat.conf - use case configuration file
• parm/met_config/GridStatConfig_wrapped - MET configuration file used in the GridStat.conf use case

Model Application Use Cases

Look at the Model Application use cases

The model_applications use case files are organized in directories by category, e.g. precipitation or data_assimilation. They contain METplus configuration files (ending with .conf). If additional files such as Python Embedding scripts are included with a use case, these files are found in a directory named after the METplus configuration file without the .conf extension.

• model_applications/data_assimilation - directory
• model_applications/data_assimilation/StatAnalysis_fcstHAFS_obsPrepBufr_JEDI_IODA_interface.conf - use case configuration file
• model_applications/data_assimilation/StatAnalysis_fcstHAFS_obsPrepBufr_JEDI_IODA_interface - directory containing supplemental files for the use case
• model_applications/data_assimilation/StatAnalysis_fcstHAFS_obsPrepBufr_JEDI_IODA_interface/read_ioda_mpr.py - script called by the use case

Example Use Case (aka "Hello World" Example - METplus style)

Let's look at the Example use case, Example.conf, under met_tool_wrapper/Example

In this file are variables, denoted in ALL_CAPS. These can be modified as needed. In the METplus system, unmodified config files remain in the directories under ${MET_BUILD_BASE}, while user-modified files are stored in ${METPLUS_TUTORIAL_DIR}/user_config.

No changes are needed in Example.conf. Close it and continue to the next page.

Unless otherwise indicated, all directories are relative to your ${METPLUS_TUTORIAL_DIR} directory.

Modify your Tutorial/User conf files

In this section you will modify the configuration files that will be read for each call to METplus.

The paths in this practical session guide assume:

• You have created a user_config directory in your ${METPLUS_TUTORIAL_DIR} directory
• You have added the shared METplus ush directory to your PATH (done in the Tutorial Setup script)
• You are using the shared installation of MET.

If not, then you need to adjust accordingly.

1. Change to the ${METPLUS_TUTORIAL_DIR} directory.  Try running run_metplus.py. You should see the usage statement output to the screen.

2. Now try to pass in the example.conf configuration file found in your parm directory under use_cases/met_tool_wrapper/Examples tutorial.conf

You should see output like this:

Note it ends with an error message stating that OUTPUT_BASE was not set correctly. You will need to configure the METplus wrappers to be able to run a use case.

Some variables in the system conf are set to '/path/to' and must be overridden to run METplus, such as OUTPUT_BASE in defaults.conf.

3. View the defaults.conf file and notice how OUTPUT_BASE = /path/to . This implies it is REQUIRED to be overridden to a valid path.

4. View the tutorial configuration files in your ${METPLUS_TUTORIAL_DIR} directory.

The INPUT_BASE, OUTPUT_BASE, and MET_INSTALL_DIR variables must all be set to run METplus. Since MET_INSTALL_DIR (and possibly INPUT_BASE) should already be set in the default METplus configuration file (completed on install of METplus), only OUTPUT_BASE is required to run. If INPUT_BASE is not set in the tutorial configuration file, it should be set correctly in the defaults configuration file.

We will test out using these configurations on the next page.

Running METplus

Running METplus involves invoking the python script run_metplus.py followed by a list of configuration files.

Reminder: The default configuration file (defaults.conf) is always read in and processed first before the configuration files passed in on the command line.

If you have configured METplus correctly and call run_metplus.py without passing in any configuration files, it will generate a usage statement to indicate that other config files are required to perform a useful task. It will generate an error statement if something is amiss.

1. Review the Example.conf configuration file - which is METplus' version of a "Hello World" example

2. Call the run_metplus.py script again, this time passing in the Example.conf configuration file and the tutorial.conf configuration file. You should see logs output to the screen.

Note: The environment variable METPLUS_BUILD_BASE determines where to look for paths to use_case files. The METPLUS_TUTORIAL_DIR environment variable determines where to look for user-modified config files.

3. Check the directory specified by the OUTPUT_BASE configuration variable. You should see that files and sub-directories have been created

4. Review the METplus log file to see what was run. Compare the log output to the Example.conf configuration file to see how they correspond to each other. The log file will have today's date in the filename. Since METplus was configured to list today's timestamp in YYYYMMDDHHMMSS format, each run of METplus will generate a separate log file. This is the same behavior as the METplus final configuration file. List all of the log files and view the latest METplus log file:

You will notice that METplus ran for 5 valid times, processing 4 forecast hours for each valid time. For each run time, it ran twice using two different input templates to find files.

5. Now run METplus passing in the Example.conf and tutorial.conf files from the previous run AND an explicit override of the OUTPUT_BASE variable.

6. Check the directory that corresponds to the config.OUTPUT_BASE value that you set in the command line override. You should see that files and sub-directories have been created in the new location. Note that you can reference an environment variable when setting values on the command line.

Remember: Additional conf files and config variable overrides are processed after the default METplus config file (defaults.conf). OUTPUT_BASE was set in tutorial.conf and then overridden with config.OUTPUT_BASE.
Order matters, since each successive conf file and/or explicit variable override will take precedence over any value set for variables defined previously.

Note: The processing order allows for structuring your conf files to contain system/user configurations (settings used for every run) and use case specific configurations (settings only used for a given use case).

Timing Control in METplus

METplus configuration variables that control timing information are described in the Timing Control section of the System Configuration chapter in the METplus User's Guide. The Example wrapper is a good tool to help understand how these settings control what is run by the METplus wrappers.

1. Copy the Example.conf configuration file in your user_config directory, renaming it Example_timing.conf

2. Open the new Example_timing.conf file with an editor.

3. Make the following changes:

3a. Change VALID_BEG and VALID_END:

    to:

3b. Change LEAD_SEQ:

    to:

3c. Change EXAMPLE_CUSTOM_LOOP_LIST:

    to:

4. Call the run_metplus.py script again, this time passing in the Example_timing.conf configuration file and the tutorial.conf configuration file. You should see logs output to the screen.

5. Review the screen output and/or log file output.

6. Open Example_timing.conf again to increase the frequency of output by making the following changes:

6a. Change VALID_END:

    to:

7. Call the run_metplus.py script again and see how the output has changed.

8. Change the timing settings to loop by initialization (or retrospective) time. Open Example_timing.conf again and make the following changes:

8a. Change LOOP_BY:

    to:

8b. Change VALID_TIME_FMT to INIT_TIME_FMT:

    to:

8c. Change VALID_BEG to INIT_BEG:

    to:

8d. Change VALID_END to INIT_END:

    to:

8e. Change VALID_INCREMENT to INIT_INCREMENT:

    to:

9. Call the run_metplus.py script once again and see how the output has changed.

10. Change the value for LOOP_BY from INIT to its nickname RETRO. Open Example_timing.conf again and make the following changes:

10a. Change LOOP_BY:

    to:

11. Call the run_metplus.py script another time.

Filename Template Settings in METplus

METplus configuration variables that control filename templates are described in the Directory and Filename Template Info section of the System Configuration chapter in the METplus User's Guide. The Example wrapper is a good tool to help understand how these settings control what is run by the METplus wrappers.

1. List the sample available in the ${METPLUS_DATA} directory

The output should list:

2. Copy the Example_timing.conf configuration file in your user_config directory, renaming it Example_filename.conf

3. Open the new Example_filename.conf file with an editor.

4. Make the following changes:

4a. Change EXAMPLE_INPUT_DIR:

    to:

4b. Change EXAMPLE_INPUT_TEMPLATE:

    to:

5. Call the run_metplus.py script again, this time passing in the Example_timing.conf configuration file and the tutorial.conf configuration file. You should see logs output to the screen.

6. Review the screen output and/or log file output.

7. Open the new Example_filename.conf file with an editor.

8. Make the following changes:

8a. Change INIT_BEG and INIT_END:

    to:

9. Call the run_metplus.py script again, passing in the modified Example_timing.conf configuration file and the tutorial.conf configuration file. You should see logs output to the screen.

10. Review the screen output and/or log file output.  It appears it is now looking for the right file.

The METplus configuration variable named LOG_TIMESTAMP_TEMPLATE controls to timestamp that is included in the log file names. The default value of LOG_TIMESTAMP_TEMPLATE is %Y%m%d%H%M%S, which will include the year, month, day, hour, minute, and second when the run_metplus.py command was executed. This will create a new log file each time run_metplus.py is called from the command line.

Starting in METplus v5.0.0, the log timestamp is also included by default in the final configuration file that is generated by the METplus wrappers. This helps with debugging because the final config file and its corresponding log file(s) are more easily identified. The path to the final config file is set by METPLUS_CONF.

The values of these settings can be changed by including them in a METplus configuration file that is passed into run_metplus.py. They can also be set directly in the run_metplus.py command using the syntax config.VARIABLE_NAME=VALUE. The following examples will use the latter.

Changing the Log Timestamp Template

In this example we will change the log timestamp template to only include the year, month, and day of the run. This will create a single log file each day. Each call to run_metplus.py will add its log output to the daily file.

1. Run the Example_timing use case from the previous exercise and override the log timestamp template to only include year, month, and day.

2. List the contents of the log directory and notice that there is now a log file that matches the format metplus.log.YYYYMMDD

3. Run the Example_filename use case from the previous exercise and again override the log timestamp template to only include YYYYMMDD.

4. Review the log file and notice that it contains the log output from both runs.

5. Search for the word "Running" to see each command.

Adding Run ID in Log Timestamp Template

The METplus configuration variable RUN_ID was added in METplus v5.0.0. This variable contains an 8 character string that is automatically generated by and is unique to each call to run_metplus.py. This variable can be referenced in other METplus configuration variables. This can be useful in a variety of ways. For example, it can be used to distinguish log files for METplus runs that may have started within the same second.

1. Run the Example_timing use case and override the log timestamp template to include the RUN_ID.

2. List the contents of the log directory and notice that there is now a log file that matches the format metplus.log.YYYYMMDDHHMMSS.XXXXXXXX where XXXXXXXX is a random set of characters.

Since PCP-Combine performs a simple operation and reformatting step, no configuration file is needed.

1. Start by making an output directory for PCP-Combine and changing directories:

2. Now let's run PCP-Combine twice using some sample data that's included with the MET tarball:

The "\" symbols in the commands above are used for ease of reading. They are line continuation markers enabling us to spread a long command line across multiple lines. They should be followed immediately by "Enter". You may copy and paste the command line OR type in the entire line with or without the "\".

Both commands run the sum command which searches the contents of the -pcpdir directory for the data required to create the requested accmululation interval.

In the first command, PCP-Combine summed up 4 3-hourly accumulation forecast files into a single 12-hour accumulation forecast. In the second command, PCP-Combine summed up 12 1-hourly accumulation observation files into a single 12-hour accumulation observation. PCP-Combine performs these tasks very quickly.

We'll use these PCP-Combine output files as input for Grid-Stat. So make sure that these commands have run successfully!

When PCP-Combine is finished, you may view the output NetCDF files it wrote using the ncdump and ncview utilities. Run the following commands to view contents of the NetCDF files:

The ncview windows display plots of the precipitation data in these files. The output of ncdump indicates that the gridded fields are named APCP_12, the GRIB code abbreviation for accumulated precipitation. The accumulation interval is 12 hours for both the forecast (3-hourly * 4 files = 12 hours) and the observation (1-hourly * 12 files = 12 hours).

Note, if ncview is not found when you run it on your system, you may need to load it first.  For example, on hera, you can use this command:

Plot-Data-Plane Tool

The Plot-Data-Plane tool can be run to visualize any gridded data that the MET tools can read. It is a very helpful utility for making sure that MET can read data from your file, orient it correctly, and plot it at the correct spot on the earth. When using new gridded data in MET, it's a great idea to run it through Plot-Data-Plane first:

Ghostview (gv) can take a little while before it displays.  If you don't have gv on your computer, try using display, or any tool that can visualize PostScript files, e.g.:

Another option is to create a PNG file from the PS file, also rotating it to appear the right way:

Next try re-running the command list above, but add the convert(x)=x/25.4; function to the config string (Hint: after the level setting and ; but before the last closing tick) to change units from millimeters to inches. What happened to the values in the colorbar?

Now, try re-running again, but add the censor_thresh=lt1.0; censor_val=0.0; options to the config string to reset any data values less 1.0 to a value of 0.0. How has your plot changed?

The convert(x) and censor_thresh/censor_val options can be used in config strings and MET config files to transform your data in simple ways.

We have run examples of the PCP-Combine -sum command, but the tool also supports the -add, -subtract, and -derive commands. While the -sum command defines a directory to be searched, for -add, -subtract, and -derive we tell PCP-Combine exactly which files to read and what data to process. The following command adds together 3-hourly precipitation from 4 forecast files, just like we did in the previous step with the -sum command:

By default, PCP-Combine looks for accumulated precipitation, and the 03 tells it to look for 3-hourly accumulations. However, that 03 string can be replaced with a configuration string describing the data to be processed, which doesn't have to be accumulated precipation. The configuration string should be enclosed in single quotes. Below, we add together the U and V components of 10-meter wind from the same input file. You would not typically want to do this, but this demonstrates the functionality. We also use the -name command line option to define a descriptive output NetCDF variable name:

While the -add command can be run on one or more input files, the -subtract command requires exactly two. Let's rerun the wind example from above but do a subtraction instead:

Now run Plot-Data-Plane to visualize this output. Use the -plot_range option to specify a the desired plotting range, the -title option to add a title, and the -color_table option to switch from the default color table to one that's good for positive and negative values:

Now view the results:

While the PCP-Combine -add and -subtract commands compute exactly one output field of data, the -derive command can compute multiple output fields in a single run. This command reads data from one or more input files and derives the output fields requested on the command line (sum, min, max, range, mean, stdev, vld_count).

Run the following command to derive several summary metrics for both the 10-meter U and V wind components:

In the above example, we used a wildcard to list multiple input file names.  And we used the -field command line option twice to specify two input fields.  For each input field, PCP-Combine loops over the input files, derives the requested metrics, and writes them to the output NetCDF file.  Run ncview to visualize this output:

This output file contains 8 variables: 2 input fields * 4 metrics. Note the output variable names the tool chose.  You can still override those names using the -name command line argument, but you would have to specify a comma-separated list of 8 names, one for each output variable.

Defining the field string

As you'll see throughout these exercises, the behavior of the MET and METplus tools is controlled using ASCII configuration files, and you will learn more about those options in the coming sessions. The field_string command line argument is actually processed as a miniature configuration file. In fact, that string is written to a temporary file which is then read by MET's configuration file library code.

In general, the name and level entries are required to extract a gridded field of data from a supported input file format. The conventions for specifying them vary based on the input file type:

1. For GRIB1 or GRIB2 inputs, set name as the abbreviation for the desired variable or data type that appears in the GRIB tables and set level to a single letter (A, Z, P, L, or R) to define the level type followed by a number to define the level value. For example 'name = "TMP"; level = "P500";' extracts 500 millibar temperature from a GRIB file.
2. For NetCDF inputs, set name as the NetCDF variable name and level to define how to extract a 2-dimensional slice of gridded data from that variable. For example, 'name = "temperature"; level = "(0,1,*,*)";' extracts a 2-dimensional field of data from the last two dimensions of a NetCDF temperature variable using the first two indicies as constant values.
3. For Python embedding, set name as the python script to be run along with any arguments for that script and do not set level. For example, 'name = "read_my_data.py input.txt";' runs a python script to read data from the specified input file.

Note that all field strings should be enclosed in single quotes, as shown above, so that they are processed on the command line as a single string which may contain embedded whitespace.

Examples for each of these input types are provided in the coming exercises. More details about setting the field string can be found in the Configuration File Overview section of the MET User's Guide. For example, if a field string matches multiple records in an input GRIB file, additional filtering criteria may be specified to further refine them. Additional options exist to explicitly specify the input file_type, define a function to convert the data, or define an operation to censor the data. Each configuration entry should be terminated with a semi-colon (;).

Since the field_string is processed using a temporary file, any syntax errors will produce a parsing error log message similar to the following:

Error messages like this typically mean there is a problem in a configuration string or configuration file being read by MET.

Plot GRIB Data

Start by creating a directory for our Plot-Data-Plane output:

Next, run Plot-Data-Plane to plot 2-meter temperature from a GRIB1 input file:

The default verbosity level of 2 only prints log messages about input and output files.

Next, re-run at verbosity level 4 to see more detailed log message about the grid being read, timing information, and the range of the data values:

It is a great idea to inspect the log messages to sanity check the metadata. Without needing to understand all the details, does the grid definition look reasonable? Are the range of values reasonable for this variable type? Are the timestamps of the data consistent with the input file name?

Next, open the PostScript output file. On some machines, the ghostview utility (or common gv alias) can display PostScript files. On other, the display command works well. On Macs, simply run open. The example below uses gv which is available in the METplus AMI:

Optionally, if the convert utility is in your path, it can be run to change the image file format. Indicate the desired output file format by specifying the suffix (.png shown below).

The plot is created using the default color table (met_default.ctable) and is scaled to the range of valid data (275 to 305). By default, no title is provided and the input file is listed as a sub-title. Does the pattern of the data look reasonable? Does it correspond well to the background map data?

If the plot of the data and the metadata listed in the log messages look reasonable, you can be confident that MET is reading your data well. In addition, the field_string you used to retrieve this data can be used in the configuration strings and configuration files for other MET tools.

While this Plot-Data-Plane validation step is not necessary for every input file, it is very useful when getting started with new input data sources.

Plot NetCDF Data

The NetCDF file format is very flexible and enables the creation of self-describing data files. However, that flexibility makes it impossible to write general purpose software to interpret all NetCDF files. For that reason, MET supports a few types of NetCDF file formats, but does not support all NetCDF files, in general. It can ingest NetCDF files that follow the Climate-Forecast Convention, are created by the WRF-Interp utility, or are created by other MET tools. Additional details can be found in the MET Data I/O chapter of the MET User's Guide.

As described in The Field String, set name to the name of the desired NetCDF variable and level to define how to index into the dimensions of that variable. In the NetCDF level strings, use *,* to indicate the two gridded dimensions. For other, non-gridded dimensions, pick a 0-based integer to specify the value to be used for that dimension. For the time dimension, if present, selecting a 0-based integer does work, however you can also specify a time string in YYYYMMDD[_HH[MMSS]] format. The square braces indicate optional elements of the format. So 19770807, 19770807_12, and 19770807_120000 are all valid time strings. It is often easier to specify a time string directly rather than finding the integer index corresponding to that time string.

Run Plot-Data-Plane to plot quantitative precipitation estimate (QPE) data from a CF-compliant NetCDF file. First run ncdump -h to take a look at the header:

Note the following from the header:

The variable P06M_NONE has 3 dimensions, but the time dimension only has length 1. Therefore, setting level = "(0,*,*)"; should do the trick. Also note that the global Conventions attribute identifies this file as being CF-compliant. Let's plot it, this time adding a title:

Running at verbosity level 4, we see the range of values:

If needed, run convert to reformat:

The next example extracts a specific time string from a precipitation analysis dataset:

Note the following from the header:

Let's request a specific time string, specify a title, and use the same plotting range as above, rather than defaulting to the range of input data values.

If running Plot-Data-Plane for multiple output times or data sources, consider using the -plot_range option to make their color scales comparable.

Lastly, let's plot output a NetCDF mask file that was created by an earlier version of MET's Gen-Vx-Mask tool.

We will plot the NCEP_Region_ID variable with only 2 dimensions:

Let's specify a different color table rather than using the default one.

Python Embedding

While the MET tools can read data from a few input gridded data file types, its ability to read data in memory from python greatly enhances its utility. Support for python embedding is optional, and must be enabled at compilation time as described in Appendix F of the MET User's Guide. MET supports three types of python embedding:

1. Reading a field of gridded data values.
2. Passing a list of point observations.
3. Passing a list of matched forecast/observation pair values.

On this page, we'll demonstrate only the first type, reading a field of gridded data values. When getting started with a new dataset via python, validating the logic by running Plot-Data-Plane is crucial. When MET reads data from flat files, important information, like the location and orientation of the data, is defined in the metadata. That is not the case for python embedding, and the responsibility for confirming those details falls to the user. So while python embedding provide extra flexibility, it also requires additional diligence.

Let's run the simplest of examples using sample data included with the MET release. The first step is to confirm that python script runs by itself, outside of MET.

This sample read_ascii_numpy.py script reads data from the input fcst.txt ASCII file and gives it a name, Forecast. Always run new python scripts on the command line first to confirm there aren't any syntax errors in the script itself. The required conventions for the python script are details in the Python Embedding for 2D Data section of the MET User's Guide.

Next, let's run Plot-Data-Plane using this python script to define the input data. As described in The Field String, this is done with the name configuration string and the level string does not apply.

Since there is no input_filename to be specified as the first required argument for Plot-Data-Plane, we provide the constant string PYTHON_NUMPY in that spot. This triggers Plot-Data-Plane to interpret the field string as a python embedding script to be run. Specifying PYTHON_XARRAY also works but requires slightly different conventions in the python embedding script.

When MET is compiled, it links to python libraries that it uses to instantiate a python interpreter at runtime. That compile time instance does have a few required packages, but will likely not include all packages that every user may want to load. You may find that your python script runs fine on the command line, but Plot-Data-Plane's call to python can't load a requested module. In that case, set the ${MET_PYTHON_EXE} environment variable to tell MET which instance of python you'd like to run.

The following command just uses that version of python that is already present in your path:

Rerunning the command from above should produce the same result, but if you look closely at the log messages, you'll see that your custom python version writes a temporary file, and MET's compile time python version reads data from it.

You can find several python embedding examples on the Sample Analysis Scripts page of the MET website. Each example includes both a python script and sample input data file. Please also see METplus Python Embedding use case examples. 

Start by making an output directory for Gen-Vx-Mask and changing directories:

Since Gen-Vx-Mask performs a simple masking step, no configuration file is needed.

We'll run the Gen-Vx-Mask tool to apply a polyline for the CONUS (Contiguous United States) to our model domain. Run Gen-Vx-Mask on the command line using the following command:

Re-run using verbosity level 3 and look closely at the log messages. How many grid points were included in this mask?

Gen-Vx-Mask should run very quickly since the grid is coarse (185x129 points) and there are 243 lat/lon points in the CONUS polyline. The more you increase the grid resolution and number of polyline points, the longer it will take to run. View the NetCDF bitmap file generated by executing the following command:

Notice that the bitmap has a value of 1 inside the CONUS polyline and 0 everywhere else. We'll use the CONUS mask we just defined in the next step.

You could try running plot_data_plane to create a PostScript image of this masking region. Can you remember how?

Notice that there are several ways that gen_vx_mask can be run to define regions of interest, some of which will be demonstrated over the next few pages.

Using a pre-defined NCEP grid from the NCEP ON388 Grid Identification Table, we'll create a latitude band, using the "lat" masking type, for the tropics region. Run Gen-Vx-Mask on the command line using the following command:

Most of the grids defined in ON388 Table B can be referenced in MET as "GNNN" where NNN is the 3-digit grid number. Run "ncdump -h" on the output file and notice that the mask variable is named "lat_mask":

Use the "-name" command line option, as shown below, to override this default.

To compute the intersection or union of two masks, use the output from the first run as input to the second. Run Gen-Vx-Mask on the command line using the following command, which uses the "lon" masking type:

Compare this intersection output to the union of the two masks, computed below:

Now we'll use the "grid" masking type to select a subgrid from a larger grid. Run Gen-Vx-Mask on the command line using the following command:

On the next page, we'll demonstrate using the "data" and "solar_alt" masking types.

On this page, we provide examples for land/sea mask and also a solar altitude to show where it is daytime on a global grid.

Run Gen-Vx-Mask on the command line using the following command:

A corresponding water mask could be created using two different methods. One way is to simply rerun the land mask command above using the -complement option:

Another way is to specify a different threshold (eq1 instead of eq0) rather than taking the complement:

Now we'll run Gen-Vx-Mask to show where it’s daytime on a global grid:

Next, combine the LAND output from the previous run with the solar altitude mask for daytime:

This creates a mask for grid points on land experiencing daylight at 18:30 UTC on 20170601.

On the next page, we'll demonstrate using the "track" and "circle" masking types.

On this page, we provide examples for using the "track" masking type using BEST track hurricane data and "circle" masking type.

Start by extracting the lat/lon locations for Hurricane Dorian:

Inspect the dorian.poly file that begins with "DORIAN" and has the lat/lon storm locations on subsequent lines.

Compute a mask for all grid points within 200 km of the Hurricane Dorian track. Run Gen-Vx-Mask using the following command:

Extract the first track point for Dorian and use the circle masking option to select all grid points within 500 km:

Run Gen-Vx-Mask on the command line using the following command:

The "circle" mask type draws a circle around each of the lat/lon points in the input poly file. The "circle" option may be useful when verifying within a certain radius of known radar locations. If you rerun without the "-thresh" command line option, Gen-Vx-Mask still runs but prints a warning message:

Can you figure out what this output file now contains? Run ncview or plot_data_plane to visualize it.

Gen-Vx-Mask also supports the "box", "solar_azi", and "shape" masking types, not covered in these exercises. Interested users can download Natural Earth shapefiles and run Gen-Vx-Mask using the "-type shape" option.

Next, we'll take a look at using the "shape" masking type with Gen-Vx-Mask.

We will demonstrate the Gen-Vx-Mask "shape" masking type using freely available shapefiles from Natural Earth.  While multiple resolutions are provided, we'll use the coarsest version for this example since it's the smallest in size.

Download the Natural Earth administrative shapefiles for countries boundaries.

Run gen_vx_mask to define the mask for the USA (index 4).

Re-run to compute the union with Canada (index 3).

Re-run to add Mexico (index 27).

The result is good but not perfect. There are a few missing grid points along the boundary. But this demonstrates how the tool works. Consider re-running all three commands again, but this time use the "-value" command line option to define the mask value to be written. Just make "-value" match the "-shapeno" option (.e.g. -value 4 for USA, -value 3 for Canada, and -value 27, for Mexico). What impact does that have on the result?

Next, we'll take a look at the functionality that Grid-Stat offers.

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

The behavior of Grid-Stat is controlled by the contents of the configuration file passed to it on the command line. The default Grid-Stat configuration file may be found in the data/config/GridStatConfig_default file. Prior to modifying the configuration file, users are advised to make a copy of the default:

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

The configurable items for Grid-Stat are used to specify how the verification is to be performed. The configurable items include specifications for the following:

• The forecast fields to be verified at the specified vertical level or accumulation interval
• The threshold values to be applied
• The areas over which to aggregate statistics - as predefined grids, configurable lat/lon polylines, or gridded data fields
• The confidence interval methods to be used
• The smoothing methods to be applied (as opposed to interpolation methods)
• The types of verification methods to be used

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

For this tutorial, we'll configure Grid-Stat to verify the 12-hour accumulated precipitation output of PCP-Combine. We'll be using Grid-Stat to verify a single field using NetCDF input for both the forecast and observation files. However, Grid-Stat may in general be used to verify an arbitrary number of fields. Edit the GridStatConfig_tutorial file as follows:

• Set:
  
  fcst = {
     field = [
        {
          name       = "APCP_12";
          level      = [ "(*,*)" ];
          cat_thresh = [ >0.0, >=5.0, >=10.0 ];
        }
     ];
  }
  obs = fcst;

  To verify the field of precipitation accumulated over 12 hours using the 3 thresholds specified.

• Set:
  
  mask = {
     grid = [];
     poly = [ "../gen_vx_mask/CONUS_mask.nc",
              "MET_BASE/poly/NWC.poly",
              "MET_BASE/poly/SWC.poly",
              "MET_BASE/poly/GRB.poly",
              "MET_BASE/poly/SWD.poly",
              "MET_BASE/poly/NMT.poly",
              "MET_BASE/poly/SMT.poly",
              "MET_BASE/poly/NPL.poly",
              "MET_BASE/poly/SPL.poly",
              "MET_BASE/poly/MDW.poly",
              "MET_BASE/poly/LMV.poly",
              "MET_BASE/poly/GMC.poly",
              "MET_BASE/poly/APL.poly",
              "MET_BASE/poly/NEC.poly",
              "MET_BASE/poly/SEC.poly" ];
  }

  To accumulate statistics over the Continental United States (CONUS) and the 14 NCEP verification regions in the United States defined by the polylines specified. To see a plot of these regions, execute the following command:

  gv ${MET_BUILD_BASE}/share/met/poly/ncep_vx_regions.pdf &

• In the boot dictionary, set:
  
  n_rep = 500;

  To turn on the computation of bootstrap confidence intervals using 500 replicates.

• In the nbrhd dictionary, set:
  
  width          = [ 3, 5 ];
  cov_thresh = [ >=0.5, >=0.75 ];

  To define two neighborhood sizes and two fractional coverage field thresholds.

• Set:
  
  output_flag = {
     fho     = NONE;
     ctc      = BOTH;
     cts      = BOTH;
     mctc   = NONE;
     mcts   = NONE;
     cnt     = BOTH;
     sl1l2   = BOTH;
     sal1l2 = NONE;
     vl1l2   = NONE;
     val1l2 = NONE;
     vcnt = NONE;
     pct     = NONE;
     pstd   = NONE;
     pjc     = NONE;
     prc     = NONE;
     eclv    = NONE;
     nbrctc = BOTH;
     nbrcts = BOTH;
     nbrcnt = BOTH;
     grad    = BOTH;
     dmap  = NONE;
     seeps  = NONE;
  }

  To compute contingency table counts (CTC), contingency table statistics (CTS), continuous statistics (CNT), scalar partial sums (SL1L2), neighborhood contingency table counts (NBRCTC), neighborhood contingency table statistics (NBRCTS), and neighborhood continuous statistics (NBRCNT).

Next, run Grid-Stat on the command line using the following command:

Grid-Stat is now performing the verification tasks we requested in the configuration file. It should take a minute or two to run. The status messages written to the screen indicate progress.

In this example, Grid-Stat performs several verification tasks in evaluating the 12-hour accumulated precipiation field:

• For continuous statistics and partial sums (CNT and SL1L2), 15 output lines each:
  (1 field * 15 masking regions)
• For contingency table counts and statistics (CTC and CTS), 45 output lines each:
  (1 field * 3 raw thresholds * 15 masking regions)
• For neighborhood methods (NBRCNT, NBRCTC, and NBRCTS), 90 output lines each:
  (1 field * 3 raw thresholds * 2 neighborhood sizes * 15 masking regions)

To greatly increase the runtime performance of Grid-Stat, you could disable the computation of bootstrap confidence intervals in the configuration file. Edit the GridStatConfig_tutorial file as follows:

• In the boot dictionary, set:
  
  n_rep = 0;

  To disable the computation of bootstrap confidence intervals.

Now, try rerunning the Grid-Stat command listed above and notice how much faster it runs. While bootstrap confidence intervals are nice to have, they take a long time to compute, especially for gridded data.

The output of Grid-Stat is one or more ASCII files containing statistics summarizing the verification performed and a NetCDF file containing difference fields. In this example, the output is written to the current directory, as we requested on the command line. It should now contain 10 Grid-Stat output files beginning with the grid_stat_ prefix, one each for the CTC, CTS, CNT, SL1L2, GRAD, NBRCTC, NBRCTS, and NBRCNT ASCII files, a STAT file, and a NetCDF matched pairs file.

The format of the CTC, CTS, CNT, and SL1L2 ASCII files will be covered for the Point-Stat tool. The neighborhood method and gradient output are unique to the Grid-Stat tool.

• Rather than comparing forecast/observation values at individual grid points, the neighborhood method compares areas of forecast values to areas of observation values. At each grid box, a fractional coverage value is computed for each field as the number of grid points within the neighborhood (centered on the current grid point) that exceed the specified raw threshold value. The forecast/observation fractional coverage values are then compared rather than the raw values themselves.
• Gradient statistics are computed on the forecast and observation gradients in the X and Y directions.

Since the lines of data in these ASCII files 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.

Execute the following command to view the NetCDF output of Grid-Stat:

Click through the 2d vars variable names in the ncview window to see plots of the forecast, observation, and difference fields for each masking region. If you see a warning message about the min/max values being zero, just click OK.

Now dump the NetCDF header:

View the NetCDF header to see how the variable names are defined.

Notice how *MANY* variables there are, separate output for each of the masking regions defined. Try editing the config file again by setting apply_mask = FALSE; and gradient = TRUE; in the nc_pairs_flag dictionary. Re-run Grid-Stat and inspect the output NetCDF file. What affect did these changes have?

We have now successfully run the PCP-Combine and Grid-Stat tools to verify 12-hourly accumulated preciptation for a single output time. We did the following steps:

• Identified our forecast and observation datasets.
• Constructed PCP-Combine commands to put them into a common accumulation interval.
• Configured and ran Grid-Stat to compute our desired verification statistics.

Now that we've defined the logic for a single run, the next step would be writing a script to automate these steps for many model initializations and forecast lead times. Rather than every MET user rewriting the same type of scripts, use METplus to automate these steps in a use case!

Answers to Exercises from Session 1

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

PB2NC Tool: Configure

The behavior of PB2NC is controlled by the contents of the configuration file passed to it on the command line. The default PB2NC configuration may be found in the data/config/PB2NCConfig_default file.

The configurable items for PB2NC are used to filter out the PrepBufr observations that should be retained or derived. You may find a complete description of the configurable items in the pb2nc configuration file section of the MET User's Guide or in the Configuration File Overview.

For this tutorial, edit the PB2NCConfig_tutorial_run1 file as follows:

• Set:
  
  message_type = [ "ADPUPA", "ADPSFC" ];

  to retain only those 2 message types. Message types are described in:
  http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_1.htm

• Set:
  
  obs_window = {
     beg = -1800;
     end =  1800;
  }

  so that only observations within 1800 second (30 minutes) of the file time will be retained.

• Set:
  
  mask = {
     grid = "G212";
     poly = "";
  }

  to retain only those observations residing within NCEP Grid 212, on which the forecast data resides.

• Set:
  
  obs_bufr_var = [ "QOB", "TOB", "UOB", "VOB", "D_WIND", "D_RH" ];

  to retain observations for specific humidity, temperature, the u-component of wind, and the v-component of wind and to derive observation values for wind speed and relative humidity.

  While we are request these observation variable names from the input file, the following corresponding strings will be written to the output file: SPFH, TMP, UGRD, VGRD, WIND, RH. This mapping of input PrepBufr variable names to output variable names is specified by the obs_prepbufr_map config file entry. This enables the new features in the current version of MET to be backward compatible with earlier versions.

PB2NC Tool: Run

PB2NC is now filtering the observations from the PrepBufr file using the configuration settings we specified and writing the output to the NetCDF file name we chose. This should take a few minutes to run. As it runs, you should see several status messages printed to the screen to indicate progress. You may use the -v command line option to turn off (-v 0) or change the amount of log information printed to the screen.

If you'd like to filter down the observations further, you may want to narrow the time window or modify other filtering criteria. We will do that after inspecting the resultant NetCDF file.

PB2NC Tool: Output

When PB2NC is finished, you may view the output NetCDF file it wrote using the ncdump utility.

In the NetCDF header, you'll see that the file contains nine dimensions and nine variables. The obs_arr variable contains the actual observation values. The obs_qty variable contains the corresponding quality flags. The four header variables (hdr_typ, hdr_sid, hdr_vld, hdr_arr) contain information about the observing locations.

The obs_var, obs_unit, and obs_desc variables describe the observation variables contained in the output. The second entry of the obs_arr variable (i.e. var_id) lists the index into these array for each observation. For example, for observations of temperature, you'd see TMP in obs_var, KELVIN in obs_unit, and TEMPERATURE OBSERVATION in obs_desc. For observations of temperature in obs_arr, the second entry (var_id) would list the index of that temperature information.

Plot-Point-Obs

The plot_point_obs tool plots the location of these NetCDF point observations. Just like plot_data_plane is useful to visualize gridded data, run plot_point_obs to make sure you have point observations where you expect.

Each red dot in the plot represents the location of at least one observation value. The plot_point_obs tool has additional command line options for filtering which observations get plotted and the area to be plotted.

By default, the points are plotted on the full globe.

MET extracts the grid information from the first record of that GRIB file and plots the points on that domain.

The plot_data_plane tool can be run on the NetCDF output of any of the MET point observation pre-processing tools (pb2nc, ascii2nc, madis2nc, and lidar2nc).

PB2NC Tool: Reconfigure and Rerun

Now we'll rerun PB2NC, but this time we'll tighten the observation acceptance criteria.

• Set:
  
  message_type = [];

  to retain all message types.

• Set:
  
  obs_window = {
     beg = -25*30;
     end =  25*30;
  }

  so that only observations 25 minutes before and 25 minutes after the top of the hour are retained.

• Set:
  
  quality_mark_thresh = 1;

  to retain only the observations marked "Good" by the NCEP quality control system.

The majority of the observations were rejected because their valid time no longer fell inside the tighter obs_window setting.

When configuring PB2NC for your own projects, you should err on the side of keeping more data rather than less. As you'll see, the grid-to-point verification tools (Point-Stat and Ensemble-Stat) allow you to further refine which point observations are actually used in the verification. However, keeping a lot of point observations that you'll never actually use will make the data files larger and slightly slow down the verification. For example, if you're using a Global Data Assimilation (GDAS) PREPBUFR file to verify a model over Europe, it would make sense to only keep those point observations that fall within your model domain.

ASCII2NC Tool: Run

Since ASCII2NC performs a simple reformatting step, typically no configuration file is needed. However, when processing high-frequency (1 or 3-minute) SURFRAD data, a configuration file may be used to define a time window and summary metric for each station. For example, you might compute the average observation value +/- 15 minutes at the top of each hour for each station. In this example, we will not use a configuration file.

The sample ASCII observations in the MET tarball are still identified by GRIB code rather than the newer variable name option.

For example, 52 corresponds to RH:

ASCII2NC should perform this reformatting step very quickly since the sample file only contains data for 5 stations.

ASCII2NC Tool: Output

When ASCII2NC is finished, you may view the output NetCDF file it wrote using the ncdump utility.

The NetCDF header should look nearly identical to the output of the NetCDF output of PB2NC. You can see the list of stations for which we have data by inspecting the hdr_sid_table variable:

This ASCII data only contains observations at a few locations.

Next, we'll use the NetCDF output of PB2NC and ASCII2NC to perform Grid-to-Point verification using the Point-Stat tool.

Point-Stat Tool: Configure

The behavior of Point-Stat is controlled by the contents of the configuration file passed to it on the command line. The default Point-Stat configuration file may be found in the data/config/PointStatConfig_default file.

The configurable items for Point-Stat are used to specify how the verification is to be performed. The configurable items include specifications for the following:

• The forecast fields to be verified at the specified vertical levels.
• The type of point observations to be matched to the forecasts.
• The threshold values to be applied.
• The areas over which to aggregate statistics - as predefined grids, lat/lon polylines, or individual stations.
• The confidence interval methods to be used.
• The interpolation methods to be used.
• The types of verification methods to be used.

• Set:
  
  fcst = {
    message_type = [ "ADPUPA" ];
    field = [
      {
        name     = "TMP";
        level      = [ "P850-1050", "P500-850" ];
        cat_thresh = [ <=273, >273 ];
      }
    ];
  }
  obs = fcst;

  to verify temperature over two different pressure ranges against ADPUPA observations using the thresholds specified.

• Set:
  
  ci_alpha = [ 0.05, 0.10 ];

  to compute confidence intervals using both a 5% and a 10% level of certainty.

• Set:
  
  output_flag = {
     fho    = BOTH;
     ctc    = BOTH;
     cts    = STAT;
     mctc   = NONE;
     mcts   = NONE;
     cnt    = BOTH;
     sl1l2  = STAT;
     sal1l2 = NONE;
     vl1l2  = NONE;
     val1l2 = NONE;
     pct    = NONE;
     pstd   = NONE;
     pjc    = NONE;
     prc    = NONE;
     ecnt   = NONE;
     eclv   = BOTH;
     mpr    = BOTH;
  }

  to indicate that the forecast-hit-observation (FHO) counts, contingency table counts (CTC), contingency table statistics (CTS), continuous statistics (CNT), partial sums (SL1L2), economic cost/loss value (ECLV), and the matched pair data (MPR) line types should be output. Setting SL1L2 and CTS to STAT causes those lines to only be written to the output .stat file, while setting others to BOTH causes them to be written to both the .stat file and the optional LINE_TYPE.txt file.

• Set:
  
  output_prefix = "run1";

  to customize the output file names for this run.

Note that in the mask dictionary, the grid entry is set to FULL. This instructs Point-Stat to compute statistics over the entire input model domain. Setting grid to FULL has this special meaning.

Point-Stat Tool: Run

Next, run Point-Stat to compare a GRIB forecast to the NetCDF point observation output of the ASCII2NC tool.

Point-Stat is now performing the verification tasks we requested in the configuration file. It should take less than a minute to run. You should see several status messages printed to the screen to indicate progress.

If you receive a syntax error such as the one listed below, review PointStatConfig_tutorial_run1 for an extra comma after the "}" on line number 59

Notice the more detailed information about which observations were used for each verification task. If you run Point-Stat and get fewer matched pairs than you expected, try using the -v 3 option to see why the observations were rejected.

Point-Stat Tool: Output

The output of Point-Stat is one or more ASCII files containing statistics summarizing the verification performed. Since we wrote output to the current directory, it should now contain 6 ASCII files that begin with the point_stat_ prefix, one each for the FHO, CTC, CNT, ECLV, and MPR types, and a sixth for the STAT file. The STAT file contains all of the output statistics while the other ASCII files contain the exact same data organized by line type.

• This is a simple ASCII file consisting of several rows of data.
• Each row contains data for a single verification task.
• The FCST_LEAD, FCST_VALID_BEG, and FCST_VALID_END columns indicate the timing information of the forecast field.
• The OBS_LEAD, OBS_VALID_BEG, and OBS_VALID_END columns indicate the timing information of the observation field.
• The FCST_VAR, FCST_UNITS, FCST_LEV, OBS_VAR, OBS_UNITS, and OBS_LEV columns indicate the two parts of the forecast and observation fields set in the configure file.
• The OBTYPE column indicates the PrepBufr message type used for this verification task.
• The VX_MASK column indicates the masking region over which the statistics were accumulated.
• The INTERP_MTHD and INTERP_PNTS columns indicate the method used to interpolate the forecast data to the observation location.
• The FCST_THRESH and OBS_THRESH columns indicate the thresholds applied to FCST_VAR and OBS_VAR.
• The COV_THRESH column is not applicable here and will always have NA when using point_stat.
• The ALPHA column indicates the alpha used for confidence intervals.
• The LINE_TYPE column indicates that these are CTC contingency table count lines.
• The TOTAL column indicates the total number of matched pairs.
• The remaining columns contain the counts for the contingency table computed by applying the threshold to the forecast/observation matched pairs. The FY_OY (forecast: yes, observation: yes), FY_ON (forecast: yes, observation: no), FN_OY (forecast: no, observation: yes), and FN_ON (forecast: no, observation: no) columns indicate those counts.

1. What do you notice about the structure of the contingency table counts with respect to the two thresholds used? Does this make sense?
2. Does the model appear to resolve relatively cold surface temperatures?
3. Based on these observations, are temperatures >273 relatively rare or common in the P850-500 range? How can this affect the ability to get a good score using contingency table statistics? What about temperatures <=273 at the surface?

• The columns prior to LINE_TYPE contain the same data as the previous file we viewed.
• The LINE_TYPE column indicates that these are CNT continuous lines.
• The remaining columns contain continuous statistics derived from the raw forecast/observation pairs. See the CNT OUTPUT FORMAT section in the Point-Stat section of the MET User's Guide for a thorough description of the output.
• Again, confidence intervals are given for each of these statistics as described above.

1. What conclusions can you draw about the model's performance at each level using continuous statistics? Justify your answer. Did you use a single metric in your evaluation? Why or why not?
2. Comparing the first line with an alpha value of 0.05 to the second line with an alpha value of 0.10, how does the level of confidence change the upper and lower bounds of the confidence intervals (CIs)?
3. Similarly, comparing the first line with few numbers of matched pairs in the TOTAL column to the third line with more, how does the sample size affect how you interpret your results?

• The columns prior to LINE_TYPE contain the same data as the previous file we viewed.
• The LINE_TYPE column indicates that these are FHO forecast-hit-observation rate lines.
• The remaining columns are similar to the contingency table output and contain the total number of matched pairs, the forecast rate, the hit rate, and observation rate.
• The forecast, hit, and observation rates should back up your answer to the third question about the contingency table output.

• The columns prior to LINE_TYPE contain the same data as the previous file we viewed.
• The LINE_TYPE column indicates that these are MPR matched pair lines.
• The remaining columns are similar to the contingency table output and contain the total number of matched pairs, the matched pair index, the latitude, longitude, and elevation of the observation, the forecasted value, the observed value, and the climatological value (if applicable).
• There is a lot of data here and it is recommended that the MPR line_type is used only to verify the tool is working properly.

Point-Stat Tool: Reconfigure

Now we'll reconfigure and rerun Point-Stat.

This time, we'll use two dictionary entries to specify the forecast field in order to set different thresholds for each vertical level. Point-Stat may be configured to verify as many or as few model variables and vertical levels as you desire.

• Set:
  
  fcst = {
     field = [
      {
        name       = "TMP";
        level      = [ "Z2" ];
        cat_thresh = [ >273, >278, >283, >288 ];
      },
      {
        name       = "TMP";
        level      = [ "P750-850" ];
        cat_thresh = [ >278 ];
      }
    ];
  }
  obs = fcst;

  to verify 2-meter temperature and temperature fields between 750hPa and 850hPa, using the thresholds specified.

• Set:
  
  message_type = ["ADPUPA","ADPSFC"];
  sid_inc = [];
  sid_exc = [];
  obs_quality = [];
  duplicate_flag = NONE;
  obs_summary = NONE;
  obs_perc_value = 50;

  to include the Upper Air (UPA) and Surface (SFC) observations in the evaluation

• Set:
  
  mask = {
     grid  = [ "G212" ];
     poly  = [ "MET_BASE/poly/EAST.poly",
               "MET_BASE/poly/WEST.poly" ];
     sid   = [];
     llpnt = [];
  }

  to compute statistics over the NCEP Grid 212 region and over the Eastern and Western United States, as defined by the polylines specified.

• Set:
  
  interp = {
     vld_thresh = 1.0;
     shape       = SQUARE;
     type = [
        {
           method = NEAREST;
           width  = 1;
        },
        {
           method = DW_MEAN;
           width  = 5;
        }
     ];
  }

  to indicate that the forecast values should be interpolated to the observation locations using the nearest neighbor method and by computing a distance-weighted average of the forecast values over the 5 by 5 box surrounding the observation location.

• Set:
  
  output_flag = {
     fho    = BOTH;
     ctc    = BOTH;
     cts    = BOTH;
     mctc   = NONE;
     mcts   = NONE;
     cnt    = BOTH;
     sl1l2  = BOTH;
     sal1l2 = NONE;
     vl1l2  = NONE;
     val1l2 = NONE;
     pct    = NONE;
     pstd   = NONE;
     pjc    = NONE;
     prc    = NONE;
     ecnt   = NONE;
     eclv   = BOTH;
     mpr    = BOTH;
  }

  to switch the SL1L2 and CTS output to BOTH and generate the optional ASCII output files for them.

• Set:
  
  output_prefix = "run2";

  to customize the output file names for this run.

• 2 fields: TMP/Z2 and TMP/P750-850
• 2 observing message types: ADPUPA and ADPSFC
• 3 masking regions: G212, EAST.poly, and WEST.poly
• 2 interpolations: UW_MEAN width 1 (nearest-neighbor) and DW_MEAN width 5

Multiplying 2 * 2 * 3 * 2 = 24. So in this example, Point-Stat will accumulate matched forecast/observation pairs into 24 groups. However, some of these groups will result in 0 matched pairs being found. To each non-zero group, the specified threshold(s) will be applied to compute contingency tables.

Point-Stat Tool: Rerun

Next, run Point-Stat to compare a GRIB forecast to the NetCDF point observation output of the PB2NC tool, as opposed to the much smaller ASCII2NC output we used in the first run.
 

Point-Stat is now performing the verification tasks we requested in the configuration file. It should take a minute or two to run. You should see several status messages printed to the screen to indicate progress. Note the number of matched pairs found for each verification task, some of which are 0.

Plot-Data-Plane Tool

In this step, we have verified 2-meter temperature. The Plot-Data-Plane tool within MET provides a way to visualize the gridded data fields that MET can read.

Plot-Data-Plane requires an input gridded data file, an output postscript image file name, and a configuration string defining which 2-D field is to be plotted.

1. Set the title to 2-m Temperature.
2. Set the plotting range as 250 to 305.
3. Use the color table named ${MET_BUILD_BASE}/share/met/colortables/NCL_colortables/wgne15.ctable

Next, we'll take a look at the Point-Stat output we just generated.

See the usage statement for all MET tools using the --help command line option or with no options at all.

Point-Stat Tool: Output

The format for the CTC, CTS, and CNT line types are the same. However, the numbers will be different as we used a different set of observations for the verification.

• The columns prior to LINE_TYPE contain header information.
• The LINE_TYPE column indicates that these are CTS lines.
• The remaining columns contain statistics derived from the threshold contingency table counts. See the point_stat output section of the MET User's Guide for a thorough description of the output.
• Confidence intervals are given for each of these statistics, computed using either one or two methods. The columns ending in _NCL(normal confidence lower) and _NCU (normal confidence upper) give lower and upper confidence limits computed using assumptions of normality. The columns ending in _BCL (bootstrap confidence lower) and _BCU (bootstrap confidence upper) give lower and upper confidence limits computed using bootstrapping.

• The columns prior to LINE_TYPE contain header information.
• The LINE_TYPE column indicates these are SL1L2 partial sums lines.

Lastly, the point_stat_run2_360000L_20070331_120000V.stat file contains all of the same data we just viewed but in a single file. The Stat-Analysis tool, which we'll use later in this tutorial, searches for the .stat output files by default but can also read the .txt output files.

Stat-Analysis Tool: Configure

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.

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.

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

The output for these two jobs are printed to the screen.

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

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

Stat-Analysis Tool: Output

On the previous page, we generated the output file aggr_ctc_lines.out by using the -out command line argument.

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.

This job should aggregate 2 CTC lines for 2-meter temperature across the EAST and WEST regions.

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.

    

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

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

Series-Analysis Tool: Configure

The behavior of Series-Analysis is controlled by the contents of the configuration file passed to it on the command line. The default Series-Analysis configuration file may be found in the data/config/SeriesAnalysisConfig_default file.

The configurable items for Series-Analysis are used to specify how the verification is to be performed. The configurable items include specifications for the following:

• The forecast fields to be verified at the specified vertical level or accumulation interval.
• The threshold values to be applied.
• The area over which to limit the computation of statistics - as predefined grids or configurable lat/lon polylines.
• The confidence interval methods to be used.
• The smoothing methods to be applied.
• The types of statistics to be computed.

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

• Set the fcst dictionary to
  fcst = {
     field = [
        {
          name  = "APCP";
          level = [ "A3" ];
        }
     ];
  }

  To request the GRIB abbreviation for precipitation (APCP) accumulated over 3 hours (A3).

• Delete obs = fcst; and insert
  obs = {
     field = [
        {
          name  = "APCP_03";
          level = [ "(*,*)" ];
        }
     ];
  }

  To request the NetCDF variable named APCP_03 where its two dimensions are the gridded dimensions (*,*).

• Look up a few lines above the fcst dictionary and set
  cat_thresh = [ >0.0, >=5.0 ];

  To define the categorical thresholds of interest. By defining this at the top level of config file context, these thresholds will be applied to both the fcst and obs settings.

• In the mask dictionary, set
  grid = "G212";

  To limit the computation of statistics to only those grid points falling inside the NCEP Grid 212 domain.

• Set
  block_size = 10000;

  To process 10,000 grid points in each pass through the data. Setting block_size larger should make the tool run faster but use more memory.

• In the output_stats dictionary, set
     fho    = [ "F_RATE", "O_RATE" ];
     ctc    = [ "FY_OY", "FN_ON" ];
     cts    = [ "CSI", "GSS" ];
     mctc   = [];
     mcts   = [];
     cnt    = [ "TOTAL", "RMSE" ];
     sl1l2  = [];
     pct    = [];
     pstd   = [];
     pjc    = [];
     prc    = [];

  For each line type, you can select statistics to be computed at each grid point over the series. These are the column names from those line types. Here, we select the forecast rate (FHO: F_RATE), observation rate (FHO: O_RATE), number of forecast yes and observation yes (CTC: FY_OY), number of forecast no and observation no (CTC: FN_ON), critical success index (CTS: CSI), and the Gilbert Skill Score (CTS: GSS) for each threshold, along with the root mean squared error (CNT: RMSE).

Series-Analysis Tool: Run

Note that the previous set of PCP-Combine commands could easily be run by looping through times in METplus Wrappers! The MET tools are often run using METplus Wrappers rather than typing individual commands by hand. You'll learn more about automation using the METplus Wrappers throughout the tutorial.

The statistics we requested in the configuration file will be computed separately for each grid location and accumulated over a time series of eight three-hour accumulations over a 24-hour period. Each grid point will have up to 8 matched pair values.

Note how long this command line is. Imagine how long it would be for a series of 100 files! Instead of listing all of the input files on the command line, you can list them in an ASCII file and pass that to Series-Analysis using the -fcst and -obs options.

Series-Analysis Tool: Output

The output of Series-Analysis is one NetCDF file containing the requested output statistics for each grid location on the same grid as the input files.

You may view the output NetCDF file that Series-Analysis wrote using the ncdump utility. Run the following command to view the header of the NetCDF output file:

In the NetCDF header, we see that the file contains many arrays of data. For each threshold (>0.0 and >=5.0), there are values for the requested statistics: F_RATE, O_RATE, FY_OY, FN_ON, CSI, and GSS. The file also contains the requested RMSE and TOTAL number of matched pairs for each grid location over the 24-hour period.

Next, run the ncview utility to display the contents of the NetCDF output file:

Click through the different variables to see how the performance varies over the domain. Looking at the series_cnt_RMSEvariable, are the errors larger in the south eastern or north western regions of the United States?

Why does the extent of missing data increase for CSI for the higher threshold? Compare series_cts_CSI_gt0.0 to series_cts_CSI_ge5.0. (Hint: Find the definition of Critical Success index (CSI) in the MET User's Guide and look closely at the denominator.)

Try running Plot-Data-Plane to visualize the observation rate variable for non-zero precipitation (i.e. series_fho_O_RATE_gt0.0). Since the valid range of values for this data is 0 to 1, use that to set the -plot_range option.

Setting block_size to 10000 still required 3 passes through our 185x129 grid (= 23865 grid points). What happens when you increase block_size to 24000 and re-run? Does it run slower or faster?

Similar to other MET tools, the behavior of Gen-Ens-Prod is controlled by the contents of the configuration file passed to it on the command line. The default Gen-Ens-Prod configuration file may be found in the data/config/GenEnsProdConfig_default file. The configurations used by the test script may be found in the scripts/config/GenEnsProdConfig* files.

The configurable items for Gen-Ens-Prod are a more limited version of those found in tools such as Grid-Stat and Point-Stat. The reason for this is tied to the purpose of the tool: Gen-Ens-Prod is for creating simple ensemble products and not for statistical verification with observation datasets. This is made clear in the configuration file's use of the ensemble ens dictionary, rather than a forecast and observation dictionary.

The configuration file also does not have the ability to request statistical line type output, instead allowing users to create ensemble products with the ensemble_flag dictionary.

For this tutorial, we'll configure Ensemble-Stat to 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.

Setting ens_thresh at 0.5 allows half of the ensemble member files being passed at runtime to contain invalid data and before GenEnsProd quits with an error. Similarly, vld_thresh at 0.5 allows each grid point to have invalid data for half of the ensemble member files passed at runtime and still compute the requested products for that grid point.

These products will be output to the netCDF designated at runtime with the -out flag. They will show a small sample of the options from Gen-Ens-Prod, as well as highlighting how the output changes when one of the files is entered as the control ensemble member, which will be done later this session.

Gen-Ens-Prod creates the products we requested in the configuration file, using the categorical thresholds specified. Note that we've passed the input ensemble data directly on the command line by specifying the ensemble member names using wildcards.

When Gen-Ens-Prod has completed running, there will be 1 netCDF output file, GenEnsProd_APCP24.nc

As mentioned previously, Gen-Ens-Prod only produces 1 netCDF output file. This file contains all of the requested products that were made in the configuration file.

Note that the file contains 9 variables: the latitude and longitudes, the ensemble mean and standard deviation, the ensemble valid data count, and the four categorical thresholds that were set, all with the prefix APCP_24_A24_ENS_FREQ_. These thresholds act as uncalibrated probability forecasts and can be verified against observational datasets with other MET tools.

Click through the variable names in the ncview window to see plots of the content we saw in the ncdump command.

Now that we've seen a successful run of the Gen-Ens-Prod tool, let's change the run command slightly to show how the -ctrl setting works.

Now that we've seen the output for Gen-Ens-Prod without any control ensemble members, let's change the run slightly by selecting one of the previous ensemble members as the control. Because the required changes will be performed in the run command, no edits will be made to the configuration file.

When MET has successfully finished running, 1 new netCDF file will be in the current directory containing the same ensemble products that were available in the first run. However, the output is now changed by using one of the original ensemble members as the control member.

You'll see as you click through the variables that the only difference between the two runs of Gen-Ens-Prod is APCP_24_A24_ENS_STDEV. The mean (APCP_24_A24_ENS_MEAN) and all of the categorical thresholds still use all 6 ensemble members. This is quickly confirmed by viewing APCP_24_A24_ENS_VLD for both, which shows 6 ensembles were used for both runs across the domain.

To better illustrate the difference in the APCP_24_A24_ENS_STDEV products, let's calculate the difference of the two fields using PCP-Combine, and plot the results using Plot-Data-Plane.

If the differences weren't apparent before, this plot makes it very clear how using a control ensemble member can drastically change a product. As long as it's desired, the -ctrl option is a quick way to run Gen-Ens-Prod without control members for some products, while including it in others.

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.

The configurable items for Ensemble-Stat are similar to those found in Grid-Stat and Point-Stat: forecast and observation dictionaries read in the the ensemble and observation fields, respectively, while the remaining dictionaries and options control the computed statistics and the area over which those statistics are calculated.

For this tutorial, we'll configure Ensemble-Stat to 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.

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

  To verify the 24-hour accumulated precipitation fields.

• In the observation filtering options, set:
     message_type = [ "ADPSFC" ];

  To verify against surface observations.

• In the field entry section, set:
     prob_cat_thresh= [ >=0, >=5.0, >=10.0 ];

  To specify thresholds to use for computation of the Ranked Probability Score (RPS).

• 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).

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 8 output files in the current directory: 7 ASCII statistics files (.stat, _ecnt.txt, _rhist.txt, _phist.txt, _orank.txt, _ssvar.txt , and _relp.txt ), and a NetCDF matched pairs file (_orank.nc).

The output from Ensemble-Stat is 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 12 output line types: ECNT, RPS, RHIST, PHIST, RELP, SSVAR, PCT, PSTD, PJC, PRC, ECLV, 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 RPS line type contains the Ranked Probability Score, as well as the number of ensembles that were used to calculate the score, the Ranked Probability Skill Score, and the Ranked Probability Score decomposed into its terms of reliability, resolution, and uncertainty.
3. 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, yielding 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.
4. 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. Similarly, when ensemble members drop out, RHIST lines cannot be aggregated together but PHIST lines can.
5. 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.
6. The PCT line type contains the contingency table counts for probabilistic forecasts.
7. The PSTD line type is the probabilistic statistics for dichotomous outcomes for derived ensemble relative frequencies.
8. The PJC line type contains joint and conditional factorization for derived ensemble relative frequencies
9. The PRC line type has the receiver operating characteristic for derived ensemble relative frequencies
10. The ECLV line type  is the economic cost/loss relative value for derived ensemble relative frequencies
11. 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, its rank, and the corresponding ensemble values for that point. When verifying against a griddedanalysis, the ranks can be written to the NetCDF output file.
12. 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 determined 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.txt files 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.

• 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 observation 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.

• 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).

• This file contains 1866 lines, 1 line for each observation value falling inside each verification region (VX_MASK).
• Each line contains 44 columns, including header information, the observation location and value, its 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.

This file is only created when you've verified using gridded observations and have requested its output using the output_flag parameter 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.

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 observation are tied (i.e. zero precip).  If you want to save it to a different file, make sure you set output_prefix to something meaningful, such as "run2", or "skip-constant".
• 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?

Answers to Exercises from Session 4

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

MODE Tool: Configure

The behavior of MODE is controlled by the contents of the configuration file passed to it on the command line. The default MODE configuration file may be found in the data/config/MODEConfig_default file.

We'll be using these three configuration files during this session.

The configuration items for MODE are used to specify how the object-based verification approach is to be performed. In MODE, as in the other MET statistics tools, you can compare any two fields. When necessary, the items in the configuration file are specified separately for the forecast and observation fields. In most cases though, users will be comparing the same forecast and observation fields. The configurable items include parameters for the following:

• The forecast and observation fields and vertical levels or accumulation intervals to be compared
• Options to mask out a portion of or threshold the raw fields
• The forecast and observation object definition parameters
• Options to filter out objects that don't meet a size or intensity criteria
• Flags to control the logic for matching/merging
• Weights to be applied for the fuzzy engine matching/merging algorithm
• Interest functions to be used for the fuzzy engine matching/merging algorithm
• Total interest threshold for matching/merging
• Various plotting options

We'll start here using by running the configuration files we copied over, as-is.

MODE Tool: Run