Feature #3304 power_spectra

data/config/ConfigConstants

@@ -214,3 +214,7 @@ CLIMO_STD_ANOM = 3;
 FCST_ANOM      = 4;
 FCST_STD_ANOM  = 5;
 
+// Power spectrum missing flag types
+NONE  = 1;
+MEAN  = 2;
+VALUE = 3;

data/config/GridDiagConfig_default

@@ -32,6 +32,15 @@ regrid = {
 censor_thresh = [];
 censor_val    = [];
 
+//
+// Specify the handling of missing data for power spectra
+//
+power_spectrum = {
+   missing_flag  = NONE;
+   missing_value = 0.0;
+   vld_thresh    = 0.5;
+}
+
 //
 // Data fields
 //
@@ -68,9 +77,10 @@ mask = {
 // Diagnostic output types
 //
 output_flag = {
-   histogram_1d = TRUE;
-   histogram_2d = TRUE;
-   info_theory  = FALSE;
+   histogram_1d   = TRUE;
+   histogram_2d   = TRUE;
+   info_theory    = FALSE;
+   power_spectrum = FALSE;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

docs/Users_Guide/grid-diag.rst

@@ -7,9 +7,13 @@ Grid-Diag Tool
 Introduction
 ============
 
-The Grid-Diag tool creates histograms (probability distributions when normalized) for an arbitrary collection of data fields and levels. Joint histograms are created for all possible pairs of variables. If no masking region is specified to subset the data fields spatially, then all points in the input domain are used. However, an arbitrary number of masking regions can be specified and output is created for each one. The histograms are accumulated over all of the input data files. Typically this tool is run with a time series of input data files, similar to Series-Analysis.
+Unlike the Grid-Stat and Point-Stat tools, which compute verification statistics for matched pairs of forecast and observation values, the Grid-Diag tool provides diagnostic information based on the distribution of values found in gridded input files.
 
-The Grid-Diag tool also uses the histograms to derive information theory statistics. Entropy is derived from each 1-dimensional histogram, and joint entropy and mutual information are derived from each 2-dimensional joint histogram. These statistics are defined using log base 2, rather than the natural logarithm which is also commonly used.
+Grid-Diag creates histograms (probability distributions when normalized) for an arbitrary collection of data fields and levels. Joint histograms are created for all possible pairs of variables. If no masking region is specified to subset the data fields spatially, then all points in the input domain are used. However, an arbitrary number of masking regions can be specified and output is created for each one. The histograms are accumulated over all of the input data files. Typically this tool is run with a time series of input data files, similar to Series-Analysis.
+
+Grid-Diag also uses the histograms to derive information theory statistics. Entropy is derived from each 1-dimensional histogram, and joint entropy and mutual information are derived from each 2-dimensional joint histogram. These statistics are defined using log base 2, rather than the natural logarithm which is also commonly used.
+
+Finally, Grid-Diag also derives power spectra for each input field using discrete cosine transforms, as described in :ref:`Denis et al. (2002) <Denis-2002>` and :ref:`Durran et al. (2017) <Durran-2017>`. However, power spectra are computed over the full model domain rather than being subset by spatial masking regions. Note that special logic is applied when processing kinetic energy. Grid-Diag reads the corresponding U and V-wind vector components and applies a discrete cosine transform to them separately prior to computing the power spectrum. If multiple inputs are provided with the **-data** command line option, the error power spectrum of the difference fields is also computed.
 
 Practical Information
 =====================
@@ -77,6 +81,20 @@ The configuration options listed above are common to many MET tools and are desc
 
 _____________________
 
+.. code-block:: none
+
+  power_spectrum = {
+     missing_flag  = NONE;
+     missing_value = 0.0;
+     vld_thresh    = 0.5;
+  }
+
+The **power_spectrum** dictionary defines options for computing power spectra and can be specified separately for each **data.field** entry below.
+
+The **missing_flag** and **missing_value** entries define how bad data values should be handled. For all other output types, bad data values are ignored but they are problematic for power spectra. Set **missing_flag** to **NONE** (default) to skip power spectrum when bad data is present, to **MEAN** to replace bad data with the mean of each input field, or to **VALUE** to replace bad data with the constant numeric value specified by **missing_value**. Set **vld_thresh** to a number beween 0 and 1 to define the required ratio of valid data to be present to compute power spectra output for that field.
+
+_____________________
+
 .. code-block:: none
 
   data = {
@@ -105,9 +123,10 @@ _____________________
 .. code-block:: none
 
    output_flag = {
-      histogram_1d = TRUE;
-      histogram_2d = TRUE;
-      info_theory  = FALSE;
+      histogram_1d   = TRUE;
+      histogram_2d   = TRUE;
+      info_theory    = FALSE;
+      power_spectrum = FALSE;
    }
 
 The **output_flag** dictionary controls the type of output that the Grid-Diag tool generates. Each flag should be set to **TRUE** or **FALSE** to enable the computation and writing of one or more variables to the output NetCDF file, as described below:
@@ -118,6 +137,8 @@ The **output_flag** dictionary controls the type of output that the Grid-Diag to
 
 3. **info_theory** for information theory metrics, including entropy for each **data.field** entry and mutual information and joint entropy for each pair of entries.
 
+4. **power_spectrum** for power spectrum output.
+
 grid_diag Output File
 ---------------------
 
@@ -133,3 +154,4 @@ If 2-dimensional joint historgrams are requested, a corresponding **hist_** vari
 
 If information theory output is requested, **entropy_**, **joint_entropy_**, and **mutual_information_** variables are written. Shannon entropy is derived from each 1-dimensional histogram, while joint entropy and mutual information are derived from each 2-dimensional joint histogram. These variables have one dimension for the number of masking regions and are computed using log base 2 rather than the natural logarithm. As such, their units are specified in the output as "bits" rather than "nats".
 
+If power spectrum output is requested, **wavenumber**, **wavelength**, and **power_spectrum_** variables are written. These variables have a single **wavenumber** dimension which is the minimum of the Nx and Ny dimensions of the verification grid. If multiple inputs are provided with the **-data** command line option, **error_power_spectrum_** variables are also written for the difference fields.

docs/Users_Guide/refs.rst

@@ -132,6 +132,22 @@ References
 |   *the Royal Statistical Society* A147, 278-292.
 |
 
+.. _Denis-2002:
+
+| Denis, B., J. Cote, R. Laprise, 2002:
+|    Spectral Decomposition of Two-Dimensional Atmospheric Fields on Limited-Area
+|    Domains Using the Discrete Cosine Transform (DCT). *Monthly Weather Review*, 130, 1812-1829.
+|    https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
+|
+
+.. _Durran-2017:
+
+| Durran, D., J. Weyn, M. Menchaca, 2017:
+|    Practical Considerations for Computing Dimensional Spectra from
+|    Gridded Data. *Monthly Weather Review*, 145, 3901-3910.
+|    https://doi.org/10.1175/MWR-D-17-0056.1
+|
+
 .. _Ebert-2008:
 
 | Ebert, E.E., 2008: Fuzzy verification of high-resolution gridded forecasts:

docs/Users_Guide/release-notes.rst

@@ -173,6 +173,9 @@ MET Version 13.0.0 Upgrade Instructions
      * The "mask.grid" and "mask.poly" entries have changed from strings to arrays of strings
        to support the processing of multiple masking regions.
 
+     * The "power_spectrum" dictionary is added to configure how missing data values should
+       handled when computing power spectra.
+
      * The new "output_flag" entry is a dictionary specifying the desired output types.
 
    * Gen-Ens-Prod configuration file
@@ -190,9 +193,7 @@ MET Version 13.0.0 Upgrade Instructions
        for a description of wind rotation and derivation and :numref:`config_wind_field_names` for the
        corresponding configuration options.
 
-.. dropdown:: Output format changes - NONE
-
-.. dropdown:: Output format changes - NONE
+.. dropdown:: Output format changes
 
    MET version 13.0.0 adds or modifies the following output file formats:
 
@@ -203,14 +204,19 @@ MET Version 13.0.0 Upgrade Instructions
 
      * Existing histogram variables are modified to include the "mask" dimension.
 
-     * New information theory variables are added for "entropy", "joint_entropy", and "mutual_information".
+     * New information theory variables are added for "entropy", "joint_entropy", and
+       "mutual_information".
+
+     * A new power spectrum "wavenumber" dimension is added, along with "wavenumber" and
+       "wavelength" variables. "power_spectrum" variables are written for each input,
+       and "error_power_spectrum" variables contain the differences for multiple inputs.
 
    * Gen-Ens-Prod output format
 
      * Adds new output variables with names include "EAS" and "EAS_WIDTH" for the Ensemble
        Agreement Scale algorithm.
 
-.. dropdown:: Output data changes - NONE
+.. dropdown:: Output data changes
 
    MET version 13.0.0 modifies existing output data values in the following ways:
 

internal/test_unit/config/GridDiagConfig

@@ -32,6 +32,15 @@ regrid = {
 censor_thresh = [];
 censor_val    = [];
 
+//
+// Specify the handling of missing data for power spectra
+//
+power_spectrum = {
+   missing_flag  = NONE;
+   missing_value = 0.0;
+   vld_thresh    = 0.5;
+}
+
 //
 // Data fields
 //

internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS → internal/test_unit/config/GridDiagConfig_FCST_OBS

@@ -32,15 +32,19 @@ regrid = {
 censor_thresh = [];
 censor_val    = [];
 
+//
+// Specify the handling of missing data for power spectra
+//
+power_spectrum = {
+   missing_flag  = VALUE;
+   missing_value = 0.0;
+   vld_thresh    = 0.5;
+}
+
 //
 // Data fields
 //
 data = {
-   name   = "APCP";
-   level  = ["A06"];
-   n_bins = 50;
-   range  = [0,25];
-
    field = [
       { ${FIELD1} },
       { ${FIELD2} }
@@ -62,10 +66,6 @@ mask = {
 //
 // Diagnostic output types
 //
-output_flag = {
-   histogram_1d = TRUE;
-   histogram_2d = TRUE;
-   info_theory  = FALSE;
-}
+output_flag = TRUE;
 
 ////////////////////////////////////////////////////////////////////////////////

internal/test_unit/config/GridDiagConfig_TMP

@@ -32,6 +32,15 @@ regrid = {
 censor_thresh = [];
 censor_val    = [];
 
+//
+// Specify the handling of missing data for power spectra
+//
+power_spectrum = {
+   missing_flag  = NONE;
+   missing_value = 0.0;
+   vld_thresh    = 0.5;
+}
+
 //
 // Data fields
 //
@@ -74,9 +83,10 @@ mask = {
 // Diagnostic output types
 //
 output_flag = {
-   histogram_1d = TRUE;
-   histogram_2d = TRUE;
-   info_theory  = FALSE;
+   histogram_1d   = TRUE;
+   histogram_2d   = TRUE;
+   info_theory    = FALSE;
+   power_spectrum = FALSE;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

internal/test_unit/xml/unit_grid_diag.xml

@@ -73,13 +73,13 @@
                 > &OUTPUT_DIR;/grid_diag/obs_file_list; \
          &MET_BIN;/grid_diag</exec>
     <env>
-      <pair><name>FIELD1</name> <value>set_attr_name = "FCST_APCP";</value></pair>
-      <pair><name>FIELD2</name> <value>set_attr_name = "OBS_APCP";</value></pair>
+      <pair><name>FIELD1</name> <value>set_attr_name = "FCST_APCP"; name = "APCP"; level = "A06"; n_bins = 50; range = [0,25];</value></pair>
+      <pair><name>FIELD2</name> <value>set_attr_name = "OBS_APCP"; name = "APCP"; level = "A06"; n_bins = 50; range = [0,25];</value></pair>
     </env>
     <param> \
       -data &OUTPUT_DIR;/grid_diag/fcst_file_list \
       -data &OUTPUT_DIR;/grid_diag/obs_file_list \
-      -config &CONFIG_DIR;/GridDiagConfig_APCP_06_FCST_OBS \
+      -config &CONFIG_DIR;/GridDiagConfig_FCST_OBS \
       -out &OUTPUT_DIR;/grid_diag/grid_diag_APCP_06_FCST_OBS.nc \
       -v 3
     </param>
@@ -109,13 +109,13 @@
                 > &OUTPUT_DIR;/grid_diag/obs_file_list; \
          &MET_BIN;/grid_diag</exec>
     <env>
-      <pair><name>FIELD1</name> <value>name = "APCP";</value></pair>
-      <pair><name>FIELD2</name> <value>name = "APCP";</value></pair>
+      <pair><name>FIELD1</name> <value>name = "APCP"; level = "A06"; n_bins = 50; range = [0,25];</value></pair>
+      <pair><name>FIELD2</name> <value>name = "APCP"; level = "A06"; n_bins = 50; range = [0,25];</value></pair>
     </env>
     <param> \
       -data &OUTPUT_DIR;/grid_diag/fcst_file_list \
       -data &OUTPUT_DIR;/grid_diag/obs_file_list \
-      -config &CONFIG_DIR;/GridDiagConfig_APCP_06_FCST_OBS \
+      -config &CONFIG_DIR;/GridDiagConfig_FCST_OBS \
       -out &OUTPUT_DIR;/grid_diag/grid_diag_APCP_06_VARN.nc \
       -v 3
     </param>
@@ -124,4 +124,60 @@
     </output>
   </test>
 
+  <test name="grid_diag_WIND_FCST_OBS">
+    <exec>echo "&DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F006.grib2 \
+                &DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F012.grib2 \
+                &DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F018.grib2 \
+                &DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F024.grib2" \
+                > &OUTPUT_DIR;/grid_diag/uv_wind_fcst_file_list; \
+          echo "&DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_0600_000.grb2 \
+                &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_1200_000.grb2 \
+                &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_1800_000.grb2 \
+                &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120410_0000_000.grb2" \
+                > &OUTPUT_DIR;/grid_diag/uv_wind_anly_file_list; \
+         &MET_BIN;/grid_diag</exec>
+    <env>
+      <pair><name>FIELD1</name> <value>name = "WIND"; level = "Z10"; n_bins = 200; range = [0,20];</value></pair>
+      <pair><name>FIELD2</name> <value>name = "WIND"; level = "Z10"; n_bins = 200; range = [0,20];</value></pair>
+    </env>
+    <param> \
+      -data &OUTPUT_DIR;/grid_diag/uv_wind_fcst_file_list \
+      -data &OUTPUT_DIR;/grid_diag/uv_wind_anly_file_list \
+      -config &CONFIG_DIR;/GridDiagConfig_FCST_OBS \
+      -out &OUTPUT_DIR;/grid_diag/grid_diag_WIND_Z10.nc \
+      -v 3
+    </param>
+    <output>
+      <grid_nc>&OUTPUT_DIR;/grid_diag/grid_diag_WIND_Z10.nc</grid_nc>
+    </output>
+  </test>
+
+  <test name="grid_diag_KENG_FCST_OBS">
+    <exec>echo "&DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F006.grib2 \
+                &DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F012.grib2 \
+                &DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F018.grib2 \
+                &DATA_DIR_MODEL;/grib2/gfs/gfs_2012040900_F024.grib2" \
+                > &OUTPUT_DIR;/grid_diag/keng_fcst_file_list; \
+          echo "&DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_0600_000.grb2 \
+                &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_1200_000.grb2 \
+                &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120409_1800_000.grb2 \
+                &DATA_DIR_MODEL;/grib2/gfsanl/gfsanl_4_20120410_0000_000.grb2" \
+                > &OUTPUT_DIR;/grid_diag/keng_anly_file_list; \
+         &MET_BIN;/grid_diag</exec>
+    <env>
+      <pair><name>FIELD1</name> <value>name = "KENG"; level = "Z10"; n_bins = 140; range = [0,140];</value></pair>
+      <pair><name>FIELD2</name> <value>name = "KENG"; level = "Z10"; n_bins = 140; range = [0,140];</value></pair>
+    </env>
+    <param> \
+      -data &OUTPUT_DIR;/grid_diag/keng_fcst_file_list \
+      -data &OUTPUT_DIR;/grid_diag/keng_anly_file_list \
+      -config &CONFIG_DIR;/GridDiagConfig_FCST_OBS \
+      -out &OUTPUT_DIR;/grid_diag/grid_diag_KENG_Z10.nc \
+      -v 3
+    </param>
+    <output>
+      <grid_nc>&OUTPUT_DIR;/grid_diag/grid_diag_KENG_Z10.nc</grid_nc>
+    </output>
+  </test>
+
 </met_test>