diff --git a/data/config/ConfigConstants b/data/config/ConfigConstants index f6edb9906a..c8d5c432c7 100644 --- a/data/config/ConfigConstants +++ b/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; diff --git a/data/config/GridDiagConfig_default b/data/config/GridDiagConfig_default index a732974c27..5cb4c79c22 100644 --- a/data/config/GridDiagConfig_default +++ b/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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/docs/Users_Guide/grid-diag.rst b/docs/Users_Guide/grid-diag.rst index cf11fcb44e..81a6baba45 100644 --- a/docs/Users_Guide/grid-diag.rst +++ b/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) ` and :ref:`Durran et al. (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. diff --git a/docs/Users_Guide/refs.rst b/docs/Users_Guide/refs.rst index 23f501303f..a1411e6d90 100644 --- a/docs/Users_Guide/refs.rst +++ b/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: diff --git a/docs/Users_Guide/release-notes.rst b/docs/Users_Guide/release-notes.rst index 067a5c93bd..c7b8849683 100644 --- a/docs/Users_Guide/release-notes.rst +++ b/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: diff --git a/internal/test_unit/config/GridDiagConfig b/internal/test_unit/config/GridDiagConfig index aa6553e03d..d0f83ce170 100644 --- a/internal/test_unit/config/GridDiagConfig +++ b/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 // diff --git a/internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS b/internal/test_unit/config/GridDiagConfig_FCST_OBS similarity index 87% rename from internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS rename to internal/test_unit/config/GridDiagConfig_FCST_OBS index f8d68d04ba..c65f376124 100644 --- a/internal/test_unit/config/GridDiagConfig_APCP_06_FCST_OBS +++ b/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; //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/GridDiagConfig_TMP b/internal/test_unit/config/GridDiagConfig_TMP index da7c02302e..f253d6c25a 100644 --- a/internal/test_unit/config/GridDiagConfig_TMP +++ b/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; } //////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/xml/unit_grid_diag.xml b/internal/test_unit/xml/unit_grid_diag.xml index 81c720a2b6..df6b62471e 100644 --- a/internal/test_unit/xml/unit_grid_diag.xml +++ b/internal/test_unit/xml/unit_grid_diag.xml @@ -73,13 +73,13 @@ > &OUTPUT_DIR;/grid_diag/obs_file_list; \ &MET_BIN;/grid_diag - FIELD1 set_attr_name = "FCST_APCP"; - FIELD2 set_attr_name = "OBS_APCP"; + FIELD1 set_attr_name = "FCST_APCP"; name = "APCP"; level = "A06"; n_bins = 50; range = [0,25]; + FIELD2 set_attr_name = "OBS_APCP"; name = "APCP"; level = "A06"; n_bins = 50; range = [0,25]; \ -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 @@ -109,13 +109,13 @@ > &OUTPUT_DIR;/grid_diag/obs_file_list; \ &MET_BIN;/grid_diag - FIELD1 name = "APCP"; - FIELD2 name = "APCP"; + FIELD1 name = "APCP"; level = "A06"; n_bins = 50; range = [0,25]; + FIELD2 name = "APCP"; level = "A06"; n_bins = 50; range = [0,25]; \ -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 @@ -124,4 +124,60 @@ + + 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 + + FIELD1 name = "WIND"; level = "Z10"; n_bins = 200; range = [0,20]; + FIELD2 name = "WIND"; level = "Z10"; n_bins = 200; range = [0,20]; + + \ + -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 + + + &OUTPUT_DIR;/grid_diag/grid_diag_WIND_Z10.nc + + + + + 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 + + FIELD1 name = "KENG"; level = "Z10"; n_bins = 140; range = [0,140]; + FIELD2 name = "KENG"; level = "Z10"; n_bins = 140; range = [0,140]; + + \ + -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 + + + &OUTPUT_DIR;/grid_diag/grid_diag_KENG_Z10.nc + + + diff --git a/src/basic/vx_config/config_constants.h b/src/basic/vx_config/config_constants.h index 4e3ee5b740..e62fc46d1a 100644 --- a/src/basic/vx_config/config_constants.h +++ b/src/basic/vx_config/config_constants.h @@ -607,6 +607,35 @@ enum class MatchType { NoMerge // Match with no additional merging }; +//////////////////////////////////////////////////////////////////////// + +// +// Enumeration for Grid-Diag Power Spectrum missing flag options +// + +enum class MissingDataType { + None, // No missing data type + Mean, // Replace missing data with the mean of the field + Value // Replace missing data with a constant value +}; + +// +// Struct to store power spectrum information +// + +struct PowerSpectrumInfo { + MissingDataType missing_flag; + double missing_value; + double vld_thresh; + bool skip; + + PowerSpectrumInfo() { clear(); } + ~PowerSpectrumInfo() { clear(); } + PowerSpectrumInfo(PowerSpectrumInfo const &i) { *this = i; } + PowerSpectrumInfo &operator=(const PowerSpectrumInfo &a) noexcept; + void clear(); +}; + //////////////////////////////////////////////////////////////////////// // // Constants used in configuartion files @@ -827,6 +856,7 @@ static const char conf_key_is_v_wind[] = "is_v_wind"; static const char conf_key_is_grid_relative[] = "is_grid_relative"; static const char conf_key_is_wind_speed[] = "is_wind_speed"; static const char conf_key_is_wind_direction[] = "is_wind_direction"; +static const char conf_key_is_kinetic_energy[] = "is_kinetic_energy"; static const char conf_key_is_prob[] = "is_prob"; // @@ -980,9 +1010,13 @@ static const char conf_val_beta[] = "BETA"; // Grid-Diag specific parameter key names // -static const char conf_key_hist1d_flag[] = "histogram_1d"; -static const char conf_key_hist2d_flag[] = "histogram_2d"; -static const char conf_key_info_theory_flag[] = "info_theory"; +static const char conf_key_hist1d_flag[] = "histogram_1d"; +static const char conf_key_hist2d_flag[] = "histogram_2d"; +static const char conf_key_info_theory_flag[] = "info_theory"; +static const char conf_key_power_spectrum_flag[] = "power_spectrum"; +static const char conf_key_power_spectrum[] = "power_spectrum"; +static const char conf_key_missing_flag[] = "missing_flag"; +static const char conf_key_missing_value[] = "missing_value"; // // STAT-Analysis and Pair-Stat specific parameter key names @@ -1500,6 +1534,14 @@ static const char conf_val_merge_both[] = "MERGE_BOTH"; static const char conf_val_merge_fcst[] = "MERGE_FCST"; static const char conf_val_no_merge[] = "NO_MERGE"; +// +// Grid-Diag specific parameter value names +// + +// Power spectrum missing flag values +static const char conf_val_mean[] = "MEAN"; +static const char conf_val_value[] = "VALUE"; + //////////////////////////////////////////////////////////////////////// // Default GRIB version 1 parameters diff --git a/src/basic/vx_config/config_util.cc b/src/basic/vx_config/config_util.cc index f4f17915a4..e1d3083603 100644 --- a/src/basic/vx_config/config_util.cc +++ b/src/basic/vx_config/config_util.cc @@ -432,7 +432,6 @@ map parse_conf_output_flag(Dictionary *dict, /////////////////////////////////////////////////////////////////////////////// map parse_conf_output_stats(Dictionary *dict) { - Dictionary *out_dict = (Dictionary *) nullptr; map output_map; STATLineType line_type; StringArray sa; @@ -444,7 +443,7 @@ map parse_conf_output_stats(Dictionary *dict) { } // Get the output flag dictionary - out_dict = dict->lookup_dictionary(conf_key_output_stats); + auto out_dict = dict->lookup_dictionary(conf_key_output_stats); // Loop over the output flag dictionary entries for(int i=0; in_entries(); i++) { @@ -1006,7 +1005,6 @@ NumArray parse_conf_eclv_points(Dictionary *dict) { /////////////////////////////////////////////////////////////////////////////// TimeSummaryInfo parse_conf_time_summary(Dictionary *dict) { - Dictionary *ts_dict = (Dictionary *) nullptr; TimeSummaryInfo info; bool is_correct_type = false; @@ -1017,7 +1015,7 @@ TimeSummaryInfo parse_conf_time_summary(Dictionary *dict) { } // Conf: time_summary - ts_dict = dict->lookup_dictionary(conf_key_time_summary); + auto ts_dict = dict->lookup_dictionary(conf_key_time_summary); // Conf: flag info.flag = ts_dict->lookup_bool(conf_key_flag); @@ -1160,7 +1158,6 @@ void parse_add_conf_key_values_map( map parse_conf_key_value_map( Dictionary *dict, const char *conf_key_map_name, const char *caller) { - Dictionary *map_dict = (Dictionary *) nullptr; map m; ConcatString key, val; const char *method_name = (nullptr != caller) ? caller : "parse_conf_key_value_map() -> "; @@ -1171,7 +1168,7 @@ map parse_conf_key_value_map( } // Conf: map_name: message_type_map, obs_var_map, etc - map_dict = dict->lookup_array(conf_key_map_name); + auto map_dict = dict->lookup_array(conf_key_map_name); // Loop through the array entries for(int i=0; in_entries(); i++) { @@ -1254,7 +1251,6 @@ map parse_conf_obs_to_qc_map(Dictionary *dict) { map parse_conf_key_convert_map( Dictionary *dict, const char *conf_key_map_name, const char *caller) { - Dictionary *map_dict = (Dictionary *) nullptr; StringArray sa; ConcatString key; UserFunc_1Arg fx; @@ -1267,7 +1263,7 @@ map parse_conf_key_convert_map( } // Conf: diag_convert_map - map_dict = dict->lookup_array(conf_key_map_name); + auto map_dict = dict->lookup_array(conf_key_map_name); // Loop through the array entries for(int i=0; in_entries(); i++) { @@ -1790,8 +1786,7 @@ InterpInfo & InterpInfo::operator=(const InterpInfo &a) noexcept { /////////////////////////////////////////////////////////////////////////////// InterpInfo parse_conf_interp(Dictionary *dict, const char *conf_key) { - Dictionary *interp_dict = (Dictionary *) nullptr; - Dictionary *type_dict = (Dictionary *) nullptr; + Dictionary *type_dict = nullptr; InterpInfo info; NumArray mthd_na; NumArray wdth_na; @@ -1809,7 +1804,7 @@ InterpInfo parse_conf_interp(Dictionary *dict, const char *conf_key) { } // Conf: interp - interp_dict = dict->lookup_dictionary(conf_key); + auto interp_dict = dict->lookup_dictionary(conf_key); // Conf: field - may be missing v = interp_dict->lookup_int(conf_key_field, false); @@ -2017,7 +2012,6 @@ ClimoCDFInfo &ClimoCDFInfo::operator=(const ClimoCDFInfo &a) noexcept { /////////////////////////////////////////////////////////////////////////////// ClimoCDFInfo parse_conf_climo_cdf(Dictionary *dict) { - Dictionary *cdf_dict = (Dictionary *) nullptr; ClimoCDFInfo info; NumArray bins; bool center; @@ -2030,7 +2024,7 @@ ClimoCDFInfo parse_conf_climo_cdf(Dictionary *dict) { } // Conf: climo_cdf - cdf_dict = dict->lookup_dictionary(conf_key_climo_cdf); + auto cdf_dict = dict->lookup_dictionary(conf_key_climo_cdf); // Conf: cdf_bins bins = cdf_dict->lookup_num_array(conf_key_cdf_bins); @@ -2120,7 +2114,6 @@ NbrhdInfo &NbrhdInfo::operator=(const NbrhdInfo &a) noexcept { /////////////////////////////////////////////////////////////////////////////// NbrhdInfo parse_conf_nbrhd(Dictionary *dict, const char *conf_key) { - Dictionary *nbrhd_dict = (Dictionary *) nullptr; NbrhdInfo info; int i; int v; @@ -2132,7 +2125,7 @@ NbrhdInfo parse_conf_nbrhd(Dictionary *dict, const char *conf_key) { } // Conf: nbrhd - nbrhd_dict = dict->lookup_dictionary(conf_key); + auto nbrhd_dict = dict->lookup_dictionary(conf_key); // Conf: field - may be missing @@ -2353,7 +2346,6 @@ HiRAInfo &HiRAInfo::operator=(const HiRAInfo &a) noexcept { /////////////////////////////////////////////////////////////////////////////// HiRAInfo parse_conf_hira(Dictionary *dict) { - Dictionary *hira_dict = (Dictionary *) nullptr; HiRAInfo info; int v; @@ -2364,7 +2356,7 @@ HiRAInfo parse_conf_hira(Dictionary *dict) { } // Conf: hira - hira_dict = dict->lookup_dictionary(conf_key_hira); + auto hira_dict = dict->lookup_dictionary(conf_key_hira); // Conf: flag info.flag = hira_dict->lookup_bool(conf_key_flag); @@ -3144,7 +3136,7 @@ void check_mctc_thresh(const ThreshArray &ta) { /////////////////////////////////////////////////////////////////////////////// const char * statlinetype_to_string(const STATLineType t) { - const char *s = (const char *) nullptr; + const char *s = nullptr; switch(t) { case STATLineType::sl1l2: s = stat_sl1l2_str; break; @@ -3265,7 +3257,7 @@ STATLineType string_to_statlinetype(const char *s) { /////////////////////////////////////////////////////////////////////////////// const char * bootintervaltype_to_string(const BootIntervalType t) { - auto s = (const char *) nullptr; + const char *s = nullptr; switch(t) { case BootIntervalType::BCA: s = conf_val_bca; break; @@ -3909,6 +3901,95 @@ NormalizeType parse_conf_normalize(Dictionary *dict) { return t; } +/////////////////////////////////////////////////////////////////////////////// + +MissingDataType int_to_missingdatatype(int v) { + MissingDataType t = MissingDataType::None; + + // Convert integer to enumerated MissingDataType + if(v == conf_const.lookup_int(conf_val_none)) t = MissingDataType::None; + else if(v == conf_const.lookup_int(conf_val_mean)) t = MissingDataType::Mean; + else if(v == conf_const.lookup_int(conf_val_value)) t = MissingDataType::Value; + else { + mlog << Error << "\nint_to_missingdatatype() -> " + << "Unexpected value of " << v << ".\n\n"; + exit(1); + } + + return t; +} + +/////////////////////////////////////////////////////////////////////////////// + +ConcatString missingdatatype_to_string(MissingDataType type) { + ConcatString s; + + // Convert enumerated MissingDataType to string + switch(type) { + case MissingDataType::None: s = conf_val_none; break; + case MissingDataType::Mean: s = conf_val_mean; break; + case MissingDataType::Value: s = conf_val_value; break; + default: + mlog << Error << "\nmissingdatatype_to_string() -> " + << "Unexpected MissingDataType value of " + << enum_class_as_int(type) << ".\n\n"; + exit(1); + } + + return s; +} + +/////////////////////////////////////////////////////////////////////////////// +// +// Code for PowerSpectrumInfo struct +// +/////////////////////////////////////////////////////////////////////////////// + +void PowerSpectrumInfo::clear() { + missing_flag = MissingDataType::None; + missing_value = bad_data_double; + vld_thresh = bad_data_double; + skip = false; +} + +/////////////////////////////////////////////////////////////////////////////// + +PowerSpectrumInfo &PowerSpectrumInfo::operator=(const PowerSpectrumInfo &a) noexcept { + if(this != &a) { + missing_flag = a.missing_flag; + missing_value = a.missing_value; + vld_thresh = a.vld_thresh; + skip = a.skip; + } + return *this; +} + +/////////////////////////////////////////////////////////////////////////////// + +PowerSpectrumInfo parse_conf_power_spectrum(Dictionary *dict) { + PowerSpectrumInfo info; + + if(!dict) { + mlog << Error << "\nparse_conf_power_spectrum() -> " + << "empty dictionary!\n\n"; + exit(1); + } + + // Conf: power_spectrum + auto ps_dict = dict->lookup_dictionary(conf_key_power_spectrum); + + // Conf: missing_flag + info.missing_flag = int_to_missingdatatype(ps_dict->lookup_int(conf_key_missing_flag)); + + // Conf: missing_value + info.missing_value = ps_dict->lookup_double(conf_key_missing_value); + + // Conf: vld_thresh + info.vld_thresh = ps_dict->lookup_double(conf_key_vld_thresh); + + return info; +} + /////////////////////////////////////////////////////////////////////////////// // // Print consistent error message and exit diff --git a/src/basic/vx_config/config_util.h b/src/basic/vx_config/config_util.h index 39ff66c2c5..a1ee1b4c60 100644 --- a/src/basic/vx_config/config_util.h +++ b/src/basic/vx_config/config_util.h @@ -202,7 +202,11 @@ extern ConcatString dist_to_string(DistType, const NumArray &); extern ConcatString griddecomptype_to_string(GridDecompType); extern ConcatString wavelettype_to_string(WaveletType); +extern MissingDataType int_to_missingdatatype(int); +extern ConcatString missingdatatype_to_string(MissingDataType); + extern int parse_conf_percentile(Dictionary *dict); +extern PowerSpectrumInfo parse_conf_power_spectrum(Dictionary *dict); extern void python_compile_error(const char *caller=nullptr); extern void ugrid_compile_error(const char *caller=nullptr); diff --git a/src/basic/vx_util/data_plane.cc b/src/basic/vx_util/data_plane.cc index 67912013fc..f828c41ead 100644 --- a/src/basic/vx_util/data_plane.cc +++ b/src/basic/vx_util/data_plane.cc @@ -458,6 +458,16 @@ void DataPlane::set_all(float *data, int nx, int ny) { /////////////////////////////////////////////////////////////////////////////// +bool DataPlane::has_bad_data() const { + + // Search for bad data + auto it = find(Data.begin(), Data.end(), bad_data_double); + + return it != Data.end(); +} + +/////////////////////////////////////////////////////////////////////////////// + bool DataPlane::is_all_bad_data() const { bool status = true; @@ -494,6 +504,27 @@ int DataPlane::n_good_data() const { /////////////////////////////////////////////////////////////////////////////// +double DataPlane::mean() const { + double sum = 0.0; + int n = 0; + +#pragma omp parallel default(none) \ + shared(sum, n) + { + + // Accumulate sum +#pragma omp for schedule(static) \ + reduction(+: sum, n) + for(int j=0; j 0 ? sum / n : bad_data_double); +} + +/////////////////////////////////////////////////////////////////////////////// + double DataPlane::get(int x, int y) const { int n; diff --git a/src/basic/vx_util/data_plane.h b/src/basic/vx_util/data_plane.h index 6b548a5b5d..883347bae0 100644 --- a/src/basic/vx_util/data_plane.h +++ b/src/basic/vx_util/data_plane.h @@ -98,8 +98,10 @@ class DataPlane { int nxy() const; bool is_empty() const; + bool has_bad_data() const; bool is_all_bad_data() const; int n_good_data() const; + double mean() const; unixtime init() const; unixtime valid() const; diff --git a/src/libcode/vx_data2d/data_class.cc b/src/libcode/vx_data2d/data_class.cc index 2a2462a419..d58032f802 100644 --- a/src/libcode/vx_data2d/data_class.cc +++ b/src/libcode/vx_data2d/data_class.cc @@ -263,24 +263,17 @@ mlog << Debug(3) << "Resetting grid definition from \"" << "\".\n"; // - // Make sure the grid dimensions do not change + // Check the grid dimensions // - if ( raw_nx() <= 0 && raw_ny() <= 0 ) { - - mlog << Warning << "\nMet2dDataFile::set_grid() -> " - << "When resetting the grid definition to \"" - << grid.serialize() << "\", the grid dimensions " - << "are changed (" << grid.nx() << ", " << grid.ny() - << ") != (" << raw_nx() << ", " << raw_ny() << ").\n\n"; - } - else if ( raw_nx() != grid.nx() || raw_ny() != grid.ny() ) { + if ( raw_nx() > 0 && raw_ny() > 0 && + ( raw_nx() != grid.nx() || raw_ny() != grid.ny() ) ) { mlog << Error << "\nMet2dDataFile::set_grid() -> " << "When resetting the grid definition to \"" << grid.serialize() << "\", the grid dimensions " - << "cannot change to (" << grid.nx() << ", " << grid.ny() - << ") from (" << raw_nx() << ", " << raw_ny() << ").\n\n"; + << "cannot change from (" << raw_nx() << ", " << raw_ny() + << ") to (" << grid.nx() << ", " << grid.ny() << ").\n\n"; exit ( 1 ); diff --git a/src/libcode/vx_data2d/var_info.cc b/src/libcode/vx_data2d/var_info.cc index 309f61ef4c..31d9e5a20d 100644 --- a/src/libcode/vx_data2d/var_info.cc +++ b/src/libcode/vx_data2d/var_info.cc @@ -142,6 +142,7 @@ void VarInfo::assign(const VarInfo &v) { SetAttrIsGridRelative = v.SetAttrIsGridRelative; SetAttrIsWindSpeed = v.SetAttrIsWindSpeed; SetAttrIsWindDirection = v.SetAttrIsWindDirection; + SetAttrIsKineticEnergy = v.SetAttrIsKineticEnergy; SetAttrIsProb = v.SetAttrIsProb; return; @@ -207,6 +208,7 @@ void VarInfo::clear() { SetAttrIsGridRelative = bad_data_int; SetAttrIsWindSpeed = bad_data_int; SetAttrIsWindDirection = bad_data_int; + SetAttrIsKineticEnergy = bad_data_int; SetAttrIsProb = bad_data_int; return; @@ -625,12 +627,14 @@ bool VarInfo::set_dict(Dictionary &dict, bool do_exit) { parse_set_attr_flag(dict, conf_key_is_u_wind); SetAttrIsVWind = parse_set_attr_flag(dict, conf_key_is_v_wind); + SetAttrIsGridRelative = + parse_set_attr_flag(dict, conf_key_is_grid_relative); SetAttrIsWindSpeed = parse_set_attr_flag(dict, conf_key_is_wind_speed); SetAttrIsWindDirection = parse_set_attr_flag(dict, conf_key_is_wind_direction); - SetAttrIsGridRelative = - parse_set_attr_flag(dict, conf_key_is_grid_relative); + SetAttrIsKineticEnergy = + parse_set_attr_flag(dict, conf_key_is_kinetic_energy); SetAttrIsProb = parse_set_attr_flag(dict, conf_key_is_prob); @@ -972,13 +976,17 @@ bool VarInfo::validate_wind_attributes(bool do_exit, const char *caller_name) co if(SetAttrIsVWind == 1) n++; if(SetAttrIsWindSpeed == 1) n++; if(SetAttrIsWindDirection == 1) n++; + if(SetAttrIsKineticEnergy == 1) n++; if(n > 1) { ConcatString msg; msg << "\n" << (caller_name == nullptr ? "" : caller_name) << "At most one wind attribute flag (" - << conf_key_is_u_wind << ", " << conf_key_is_v_wind << ", " - << conf_key_is_wind_speed << ", " << conf_key_is_wind_direction + << conf_key_is_u_wind << ", " + << conf_key_is_v_wind << ", " + << conf_key_is_wind_speed << ", " + << conf_key_is_wind_direction << ", " + << conf_key_is_kinetic_energy << ") can be set to true for each field.\n\n"; handle_config_error(msg, do_exit); return false; diff --git a/src/libcode/vx_data2d/var_info.h b/src/libcode/vx_data2d/var_info.h index d68e306058..1599405113 100644 --- a/src/libcode/vx_data2d/var_info.h +++ b/src/libcode/vx_data2d/var_info.h @@ -83,6 +83,7 @@ class VarInfo int SetAttrIsGridRelative; int SetAttrIsWindSpeed; int SetAttrIsWindDirection; + int SetAttrIsKineticEnergy; int SetAttrIsProb; void init_from_scratch(); diff --git a/src/libcode/vx_grid/latlon_grid.h b/src/libcode/vx_grid/latlon_grid.h index 81ccf0afe2..52df6ba48d 100644 --- a/src/libcode/vx_grid/latlon_grid.h +++ b/src/libcode/vx_grid/latlon_grid.h @@ -89,7 +89,7 @@ class LatLonGrid : public GridRep { //////////////////////////////////////////////////////////////////////// -inline double LatLonGrid::scale_km() const { return -1.0; } +inline double LatLonGrid::scale_km() const { return delta_lat * M_PI * grib_earth_radius_km / 180.0; } inline bool LatLonGrid::wrap_lon() const { return wrapLon; } diff --git a/src/libcode/vx_gsl_prob/Makefile.am b/src/libcode/vx_gsl_prob/Makefile.am index 29277933d0..1cf9ef08c7 100644 --- a/src/libcode/vx_gsl_prob/Makefile.am +++ b/src/libcode/vx_gsl_prob/Makefile.am @@ -14,6 +14,7 @@ noinst_LIBRARIES = libvx_gsl_prob.a libvx_gsl_prob_a_SOURCES = \ gsl_bvn.cc gsl_bvn.h \ gsl_cdf.cc gsl_cdf.h \ + gsl_fft2d.cc gsl_fft2d.h \ gsl_randist.cc gsl_randist.h \ gsl_statistics.cc gsl_statistics.h \ gsl_wavelet2d.cc gsl_wavelet2d.h \ diff --git a/src/libcode/vx_gsl_prob/Makefile.in b/src/libcode/vx_gsl_prob/Makefile.in index 3dbe1f7966..95c7d768c6 100644 --- a/src/libcode/vx_gsl_prob/Makefile.in +++ b/src/libcode/vx_gsl_prob/Makefile.in @@ -109,6 +109,7 @@ libvx_gsl_prob_a_AR = $(AR) $(ARFLAGS) libvx_gsl_prob_a_LIBADD = am_libvx_gsl_prob_a_OBJECTS = libvx_gsl_prob_a-gsl_bvn.$(OBJEXT) \ libvx_gsl_prob_a-gsl_cdf.$(OBJEXT) \ + libvx_gsl_prob_a-gsl_fft2d.$(OBJEXT) \ libvx_gsl_prob_a-gsl_randist.$(OBJEXT) \ libvx_gsl_prob_a-gsl_statistics.$(OBJEXT) \ libvx_gsl_prob_a-gsl_wavelet2d.$(OBJEXT) @@ -130,6 +131,7 @@ depcomp = $(SHELL) $(top_srcdir)/depcomp am__maybe_remake_depfiles = depfiles am__depfiles_remade = ./$(DEPDIR)/libvx_gsl_prob_a-gsl_bvn.Po \ ./$(DEPDIR)/libvx_gsl_prob_a-gsl_cdf.Po \ + ./$(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Po \ ./$(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Po \ ./$(DEPDIR)/libvx_gsl_prob_a-gsl_statistics.Po \ ./$(DEPDIR)/libvx_gsl_prob_a-gsl_wavelet2d.Po @@ -349,6 +351,7 @@ noinst_LIBRARIES = libvx_gsl_prob.a libvx_gsl_prob_a_SOURCES = \ gsl_bvn.cc gsl_bvn.h \ gsl_cdf.cc gsl_cdf.h \ + gsl_fft2d.cc gsl_fft2d.h \ gsl_randist.cc gsl_randist.h \ gsl_statistics.cc gsl_statistics.h \ gsl_wavelet2d.cc gsl_wavelet2d.h \ @@ -405,6 +408,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_gsl_prob_a-gsl_bvn.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_gsl_prob_a-gsl_cdf.Po@am__quote@ # am--include-marker +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_gsl_prob_a-gsl_statistics.Po@am__quote@ # am--include-marker @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libvx_gsl_prob_a-gsl_wavelet2d.Po@am__quote@ # am--include-marker @@ -457,6 +461,20 @@ libvx_gsl_prob_a-gsl_cdf.obj: gsl_cdf.cc @AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_gsl_prob_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_gsl_prob_a-gsl_cdf.obj `if test -f 'gsl_cdf.cc'; then $(CYGPATH_W) 'gsl_cdf.cc'; else $(CYGPATH_W) '$(srcdir)/gsl_cdf.cc'; fi` +libvx_gsl_prob_a-gsl_fft2d.o: gsl_fft2d.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_gsl_prob_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_gsl_prob_a-gsl_fft2d.o -MD -MP -MF $(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Tpo -c -o libvx_gsl_prob_a-gsl_fft2d.o `test -f 'gsl_fft2d.cc' || echo '$(srcdir)/'`gsl_fft2d.cc +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Tpo $(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='gsl_fft2d.cc' object='libvx_gsl_prob_a-gsl_fft2d.o' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_gsl_prob_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_gsl_prob_a-gsl_fft2d.o `test -f 'gsl_fft2d.cc' || echo '$(srcdir)/'`gsl_fft2d.cc + +libvx_gsl_prob_a-gsl_fft2d.obj: gsl_fft2d.cc +@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_gsl_prob_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_gsl_prob_a-gsl_fft2d.obj -MD -MP -MF $(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Tpo -c -o libvx_gsl_prob_a-gsl_fft2d.obj `if test -f 'gsl_fft2d.cc'; then $(CYGPATH_W) 'gsl_fft2d.cc'; else $(CYGPATH_W) '$(srcdir)/gsl_fft2d.cc'; fi` +@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Tpo $(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Po +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='gsl_fft2d.cc' object='libvx_gsl_prob_a-gsl_fft2d.obj' libtool=no @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_gsl_prob_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o libvx_gsl_prob_a-gsl_fft2d.obj `if test -f 'gsl_fft2d.cc'; then $(CYGPATH_W) 'gsl_fft2d.cc'; else $(CYGPATH_W) '$(srcdir)/gsl_fft2d.cc'; fi` + libvx_gsl_prob_a-gsl_randist.o: gsl_randist.cc @am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libvx_gsl_prob_a_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT libvx_gsl_prob_a-gsl_randist.o -MD -MP -MF $(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Tpo -c -o libvx_gsl_prob_a-gsl_randist.o `test -f 'gsl_randist.cc' || echo '$(srcdir)/'`gsl_randist.cc @am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Tpo $(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Po @@ -625,6 +643,7 @@ clean-am: clean-generic clean-noinstLIBRARIES mostlyclean-am distclean: distclean-am -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_bvn.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_cdf.Po + -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_statistics.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_wavelet2d.Po @@ -675,6 +694,7 @@ installcheck-am: maintainer-clean: maintainer-clean-am -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_bvn.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_cdf.Po + -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_fft2d.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_randist.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_statistics.Po -rm -f ./$(DEPDIR)/libvx_gsl_prob_a-gsl_wavelet2d.Po diff --git a/src/libcode/vx_gsl_prob/gsl_fft2d.cc b/src/libcode/vx_gsl_prob/gsl_fft2d.cc new file mode 100644 index 0000000000..76a338bfa3 --- /dev/null +++ b/src/libcode/vx_gsl_prob/gsl_fft2d.cc @@ -0,0 +1,166 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2026 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include "gsl/gsl_matrix.h" + +#include "vx_log.h" +#include "gsl_fft2d.h" + +using namespace std; + +//////////////////////////////////////////////////////////////////////// +// +// Apply a 1D orthonormal DCT Type-II transform in place, computed via +// a real FFT using the even/odd reordering and twiddle-factor rotation +// described by Makhoul (1980). This produces coefficients consistent +// with scipy.fftpack.dct(type=2, norm='ortho'): +// +// X[k] = 2 * sum_n x[n] * cos(pi*k*(2n+1)/(2N)), 0 <= k < N +// scaled by sqrt(1/(4N)) for k=0, sqrt(1/(2N)) for k>0. +// +// GSL has no native DCT, so this reconstructs true cosine +// coefficients from a real FFT rather than (incorrectly) applying +// DCT-II scale factors directly to half-complex FFT output. +// +// Equations to recover the unnormalized DCT-II coefficients: +// X[k] = 2 * Re( V[k] * exp(-i*pi*k/(2n)) ) +// = 2 * (Re(V[k])*cos(theta) + Im(V[k])*sin(theta)) +// +//////////////////////////////////////////////////////////////////////// + +static void dct_typeII_1d(double *x, size_t stride, int n, + const gsl_fft_real_wavetable *wt, + gsl_fft_real_workspace *ws) { + + // Reorder the input into a contiguous buffer: + // v[k] = x[2k] for the first half + // v[n-1-k] = x[2k+1] for the second half + // This is the standard permutation used to compute a DCT-II via a + // same-length FFT, valid for both even and odd n. + vector v(n); + int half = n / 2; + for(int k=0; k X(n); + for(int k=0; k + +#include "vx_util.h" +#include "gsl/gsl_fft_real.h" + +//////////////////////////////////////////////////////////////////////// + +extern void dct_typeII_2d(double *data, int ncol, int nrow); + +//////////////////////////////////////////////////////////////////////// + +#endif /* __MET_GSL_FFT2D_H__ */ + +//////////////////////////////////////////////////////////////////////// + diff --git a/src/libcode/vx_gsl_prob/vx_gsl_prob.h b/src/libcode/vx_gsl_prob/vx_gsl_prob.h index 2ee6f3b774..4f10c74ee7 100644 --- a/src/libcode/vx_gsl_prob/vx_gsl_prob.h +++ b/src/libcode/vx_gsl_prob/vx_gsl_prob.h @@ -16,6 +16,7 @@ #include "gsl_cdf.h" #include "gsl_bvn.h" +#include "gsl_fft2d.h" #include "gsl_randist.h" #include "gsl_statistics.h" #include "gsl_wavelet2d.h" diff --git a/src/libcode/vx_nav/nav.cc b/src/libcode/vx_nav/nav.cc index a7a1f5f1ea..db7e6147d4 100644 --- a/src/libcode/vx_nav/nav.cc +++ b/src/libcode/vx_nav/nav.cc @@ -23,12 +23,6 @@ using namespace std; ////////////////////////////////////////////////////////////////// -//static const double cf = 57.2957795130823208768; - - -////////////////////////////////////////////////////////////////// - - struct Vector3D { double x, y, z; }; @@ -58,31 +52,17 @@ double gc_dist(double lat1, double lon1, double lat2, double lon2) { -double x, dp, dl; -double lat1_radians, lat2_radians; -double lon1_radians, lon2_radians; +double lat1_radians = lat1 * rad_per_deg; +double lat2_radians = lat2 * rad_per_deg; +double lon1_radians = lon1 * rad_per_deg; +double lon2_radians = lon2 * rad_per_deg; -lat1_radians = lat1 * rad_per_deg; -lat2_radians = lat2 * rad_per_deg; -lon1_radians = lon1 * rad_per_deg; -lon2_radians = lon2 * rad_per_deg; -//lat1_radians = deg_to_rad(lat1); -//lat2_radians = deg_to_rad(lat2); -//lon1_radians = deg_to_rad(lon1); -//lon2_radians = deg_to_rad(lon2); +double dp = (lat1_radians - lat2_radians); +double dl = (lon1_radians - lon2_radians); -dp = (lat1_radians - lat2_radians); -dl = (lon1_radians - lon2_radians); +double x = haversine(dp) + cos(lat1_radians)*cos(lat2_radians)*haversine(dl); -//dp = (lat1 - lat2)/cf; -//dl = (lon1 - lon2)/cf; - -x = haversine(dp) + cos(lat1_radians)*cos(lat2_radians)*haversine(dl); -//x = haversine(dp) + cos(lat1/cf)*cos(lat2/cf)*haversine(dl); - -x = earth_radius_km*ahaversine(x); - -return x; +return earth_radius_km*ahaversine(x); } @@ -157,21 +137,19 @@ double rl_dist(double lat1, double lon1, double lat2, double lon2) { -double d, beta; +double d; if ( fabs(lat1 - lat2) < 0.0001 ) { d = earth_radius_km*cosd(lat1)*(lon1 - lon2) * rad_per_deg; -// d = earth_radius_km*cos(lat1/cf)*((lon1 - lon2)/cf); return fabs(d); } -beta = rl_bearing(lat1, lon1, lat2, lon2); +double beta = rl_bearing(lat1, lon1, lat2, lon2); d = earth_radius_km*( rad_per_deg * (lat1 - lat2)/cosd(beta) ); -//d = earth_radius_km*( ((lat1 - lat2)/cf)/cos(beta/cf) ); return fabs(d); @@ -199,13 +177,10 @@ double rl_bearing(double lat1, double lon1, double lat2, double lon2) { -double mp1, mp2, beta; - -mp1 = meridional_parts(lat1); -mp2 = meridional_parts(lat2); +double mp1 = meridional_parts(lat1); +double mp2 = meridional_parts(lat2); -beta = atan2d( rad_per_deg * (lon1 - lon2), mp2 - mp1); -//beta = cf*atan2((lon1 - lon2)/cf, mp2 - mp1); +double beta = atan2d( rad_per_deg * (lon1 - lon2), mp2 - mp1); return beta; @@ -231,7 +206,6 @@ double meridional_parts(double a) { a = a * rad_per_deg; -//a /= cf; return log(tan(piover4 + (0.5*a))); @@ -259,24 +233,13 @@ double gc_bearing(double lat1, double lon1, double lat2, double lon2) { -double x, y, dl, beta; +double dl = lon1 - lon2; -//lat1 /= cf; lon1 /= cf; +double x = cosd(lat1)*sind(lat2) - sind(lat1)*cosd(lat2)*cosd(dl); -//lat2 /= cf; lon2 /= cf; +double y = cosd(lat2)*sind(dl); -dl = lon1 - lon2; - -x = cosd(lat1)*sind(lat2) - sind(lat1)*cosd(lat2)*cosd(dl); -//x = cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dl); - -y = cosd(lat2)*sind(dl); -//y = cos(lat2)*sin(dl); - -beta = atan2d(y, x); -//beta = cf*atan2(y, x); - -return beta; +return atan2d(y, x); } @@ -306,51 +269,34 @@ void gc_point_v1(double lat1, double lon1, double lat2, double lon2, { -double sp1, sp2, sl1, sl2; -double cp1, cp2, cl1, cl2; -double t, st, sth, stmt, theta; -double x, y, z; - -theta = gc_dist(lat1, lon1, lat2, lon2)/earth_radius_km; - -t = dist/earth_radius_km; - -sth = sin(theta); +double theta = gc_dist(lat1, lon1, lat2, lon2)/earth_radius_km; -st = sin(t); +double t = dist/earth_radius_km; -stmt = sin(theta - t); +double sth = sin(theta); +double st = sin(t); +double stmt = sin(theta - t); -//lat1 /= cf; lon1 /= cf; +double sp1 = sind(lat1); +double sp2 = sind(lat2); +double sl1 = sind(lon1); +double sl2 = sind(lon2); -//lat2 /= cf; lon2 /= cf; +double cp1 = cosd(lat1); +double cp2 = cosd(lat2); +double cl1 = cosd(lon1); +double cl2 = cosd(lon2); -sp1 = sind(lat1); sp2 = sind(lat2); -//sp1 = sin(lat1); sp2 = sin(lat2); - -sl1 = sind(lon1); sl2 = sind(lon2); -//sl1 = sin(lon1); sl2 = sin(lon2); - -cp1 = cosd(lat1); cp2 = cosd(lat2); -//cp1 = cos(lat1); cp2 = cos(lat2); - -cl1 = cosd(lon1); cl2 = cosd(lon2); -//cl1 = cos(lon1); cl2 = cos(lon2); - -x = cp1*cl1*stmt + cp2*cl2*st; - -y = cp1*sl1*stmt + cp2*sl2*st; - -z = (sp1*stmt + sp2*st)/sth; +double x = cp1*cl1*stmt + cp2*cl2*st; +double y = cp1*sl1*stmt + cp2*sl2*st; +double z = (sp1*stmt + sp2*st)/sth; lat = asind(z); -//lat = cf*asin(z); if ( (fabs(x) + fabs(y)) < 1.0e-6 ) lon = 0.0; else lon = atan2d(y, x); -// lon = cf*atan2(y, x); return; @@ -384,41 +330,28 @@ void gc_point_v2(double lat1, double lon1, double bear, double dist, { -double x, y, z, t; -double sp, cp, cl, sl, sb, cb, st, ct; - -//lat1 /= cf; -//lon1 /= cf; - -//bear /= cf; - -t = dist/earth_radius_km; +double t = dist/earth_radius_km; -sp = sind(lat1); cp = cosd(lat1); -//sp = sin(lat1); cp = cos(lat1); +double sp = sind(lat1); +double cp = cosd(lat1); +double sl = sind(lon1); +double cl = cosd(lon1); -sl = sind(lon1); cl = cosd(lon1); -//sl = sin(lon1); cl = cos(lon1); +double sb = sind(bear); +double cb = cosd(bear); +double st = sin(t); +double ct = cos(t); -sb = sind(bear); cb = cosd(bear); -//sb = sin(bear); cb = cos(bear); - -st = sin(t); ct = cos(t); - -x = cp*sl*ct - sp*sl*cb*st - cl*sb*st; - -y = cp*cl*ct - sp*cl*cb*st + sl*sb*st; - -z = sp*ct + cp*cb*st; +double x = cp*sl*ct - sp*sl*cb*st - cl*sb*st; +double y = cp*cl*ct - sp*cl*cb*st + sl*sb*st; +double z = sp*ct + cp*cb*st; lat = asind(z); -//lat = cf*asin(z); if ( (fabs(x) + fabs(y)) < 1.0e-6 ) lon = 0.0; else lon = atan2d(x, y); -// lon = cf*atan2(x, y); return; @@ -485,39 +418,28 @@ void rl_point_v2(double lat1, double lon1, double bear, double dist, { -double t, sb, cb, tb, mp, mp1; +double sb = sind(bear); +double cb = cosd(bear); -//bear /= cf; - -sb = sind(bear); cb = cosd(bear); -//sb = sin(bear); cb = cos(bear); - -t = dist/earth_radius_km; +double t = dist/earth_radius_km; lat = lat1 + deg_per_rad * (t*cb); -//lat = lat1 + cf*t*cb; if ( fabs(cb) < 1.0e-5 ) lon = lon1 * rad_per_deg - t*( sb/cosd(lat1) ); -// lon = (lon1/cf) - t*( sb/cos(lat1/cf) ); else { - tb = tand(bear); -// tb = tan(bear); - - mp = meridional_parts(lat); - - mp1 = meridional_parts(lat1); + double tb = tand(bear); + double mp = meridional_parts(lat); + double mp1 = meridional_parts(lat1); lon = lon1 * rad_per_deg - tb*( mp - mp1 ); -// lon = (lon1/cf) - tb*( mp - mp1 ); } lon += twopi*floor( 0.5 - (lon/twopi) ); lon = lon * deg_per_rad; -//lon *= cf; return; @@ -546,31 +468,32 @@ return; double gc_dist_to_line(double lat1, double lon1, double lat2, double lon2, double lat3, double lon3) { - Vector3D a, b, c, g, f, t; - double lat4, lon4, dist12, dist14, dist24, dist; // Convert to cartesian coordinates - a = latlon_to_xyz(lat1, lon1); - b = latlon_to_xyz(lat2, lon2); - c = latlon_to_xyz(lat3, lon3); + Vector3D a = latlon_to_xyz(lat1, lon1); + Vector3D b = latlon_to_xyz(lat2, lon2); + Vector3D c = latlon_to_xyz(lat3, lon3); // Compute intersection of arc AB with perpendicular from C - g = cross_product(a, b); - f = cross_product(c, g); - t = cross_product(g, f); + Vector3D g = cross_product(a, b); + Vector3D f = cross_product(c, g); + Vector3D t = cross_product(g, f); // Convert intersection point from cartesian back to lat/lon + double lat4; + double lon4; xyz_to_latlon(t, lat4, lon4); // Length of the arc segment - dist12 = gc_dist(lat1, lon1, lat2, lon2); + double dist12 = gc_dist(lat1, lon1, lat2, lon2); // Distance from intersection point to end points - dist14 = gc_dist(lat1, lon1, lat4, lon4); - dist24 = gc_dist(lat2, lon2, lat4, lon4); + double dist14 = gc_dist(lat1, lon1, lat4, lon4); + double dist24 = gc_dist(lat2, lon2, lat4, lon4); // If intersection is between the end points, // use the perpendicular distance + double dist; if(dist14 <= dist12 && dist24 <= dist12) { dist = gc_dist(lat3, lon3, lat4, lon4); } diff --git a/src/tools/core/wavelet_stat/wavelet_stat.cc b/src/tools/core/wavelet_stat/wavelet_stat.cc index 50d89d0911..695e40fca0 100644 --- a/src/tools/core/wavelet_stat/wavelet_stat.cc +++ b/src/tools/core/wavelet_stat/wavelet_stat.cc @@ -91,7 +91,6 @@ static void setup_ps_file (unixtime, int); static void build_outfile_name(unixtime, int, const char *, ConcatString &); static double get_fill_value(const DataPlane &, int); -static void fill_bad_data(DataPlane &, double); static void pad_field(DataPlane &, double); static void get_tile(const DataPlane &, const DataPlane &, int, int, NumArray &, NumArray &); @@ -405,10 +404,10 @@ static void process_scores() { // Replace any bad data in the fields with a fill value mlog << Debug(2) << "Forecast field: "; - fill_bad_data(fcst_dp_fill, fcst_fill); + fcst_dp_fill.replace_bad_data(fcst_fill); mlog << Debug(2) << "Observation field: "; - fill_bad_data(obs_dp_fill, obs_fill); + obs_dp_fill.replace_bad_data(obs_fill); // Pad the fields out to the nearest power of two if requested if(conf_info.grid_decomp_flag == GridDecompType::Pad) { @@ -769,27 +768,8 @@ static double get_fill_value(const DataPlane &dp, int i_vx) { // Otherwise, fill them with the mean of the valid data. // else { - int count = 0; - double sum = 0.0; - int nxy = dp.nxy(); - -#pragma omp parallel default(none) \ - shared(nxy, dp, sum, count) - { - -#pragma omp for schedule(static) \ - reduction(+: sum, count) - for(int i=0; i 0) fill_val = sum/count; - else fill_val = 0.0; + fill_val = dp.mean(); + if(is_bad_data(fill_val)) fill_val = 0.0; } return fill_val; @@ -797,41 +777,6 @@ static double get_fill_value(const DataPlane &dp, int i_vx) { //////////////////////////////////////////////////////////////////////// -static void fill_bad_data(DataPlane &dp, double fill_val) { - - // - // Replace any bad data values with the fill value - // - int count = 0; - int nxy = dp.nxy(); - -#pragma omp parallel default(none) \ - shared(nxy, dp, fill_val, count) - { - -#pragma omp for schedule(static) \ - reduction(+: count) - for(int i=0; i 0) { - mlog << Debug(2) << "Replaced " << count - << " bad data values out of " << nxy - << " points with fill value of " - << fill_val << ".\n"; - } - - return; -} - -//////////////////////////////////////////////////////////////////////// - static void pad_field(DataPlane &dp, double pad_val) { DataPlane dp_pad; diff --git a/src/tools/other/grid_diag/Makefile.am b/src/tools/other/grid_diag/Makefile.am index 59e735bf05..c279ecc888 100644 --- a/src/tools/other/grid_diag/Makefile.am +++ b/src/tools/other/grid_diag/Makefile.am @@ -28,6 +28,7 @@ grid_diag_LDADD = -lvx_stat_out \ -lvx_data2d_nc_cf \ $(UGRID_MET_LIBS) \ $(PYTHON_MET_LIBS) \ + -lvx_nav \ -lvx_statistics \ -lvx_data2d \ -lvx_nc_util \ diff --git a/src/tools/other/grid_diag/Makefile.in b/src/tools/other/grid_diag/Makefile.in index 2e0f91a9c9..666078f5dc 100644 --- a/src/tools/other/grid_diag/Makefile.in +++ b/src/tools/other/grid_diag/Makefile.in @@ -343,6 +343,7 @@ grid_diag_LDADD = -lvx_stat_out \ -lvx_data2d_nc_cf \ $(UGRID_MET_LIBS) \ $(PYTHON_MET_LIBS) \ + -lvx_nav \ -lvx_statistics \ -lvx_data2d \ -lvx_nc_util \ diff --git a/src/tools/other/grid_diag/grid_diag.cc b/src/tools/other/grid_diag/grid_diag.cc index 624946d791..7ae0a45612 100644 --- a/src/tools/other/grid_diag/grid_diag.cc +++ b/src/tools/other/grid_diag/grid_diag.cc @@ -22,6 +22,7 @@ // 005 10/03/22 Prestopnik MET #2227 Remove using namespace std and netCDF from header files // 006 10/26/22 Linden MET #2232 Refine the Grid-Diag output variable names when specifying two input data sources // 007 01/07/26 Halley Gotway MET #3171 Multiple masks and information theory +// 008 02/12/26 Halley Gotway MET #3304 Power spectrum // //////////////////////////////////////////////////////////////////////// @@ -32,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -43,6 +45,7 @@ #include "series_data.h" #include "series_pdf.h" +#include "nav.h" #include "vx_statistics.h" #include "vx_nc_util.h" #include "vx_regrid.h" @@ -60,19 +63,42 @@ using namespace netCDF; static void process_command_line(int, char **); static void setup_diag_info(void); static void process_series(void); -static void process_hist1d(const vector &); -static void process_hist2d(const vector &); +static void process_hist1d(const vector &); +static void process_hist2d(const vector &); +static void prepare_power_spectrum_data(vector &); +static void process_power_spectrum(const vector &); +static void process_error_power_spectrum(const vector &); +static std::vector radial_spectral_variance(const DataPlane &); +static void sum_spectral_variance(vector &, const vector &); +static DataPlane dct_typeII(const DataPlane &); static void process_info_theory(void); static void setup_nc_file(void); +static void read_series_data(int, VarInfo *, const StringArray &, + const GrdFileType, InputDataInfo &); +static bool read_wind_series(const StringArray &, + int, const VarInfo *, const StringArray &, + const GrdFileType, DataPlane &); +static void regrid_data(const VarInfo *, const Grid &, DataPlane &); static ConcatString get_nc_var_str(const VarInfo *, int); +static ConcatString get_nc_att_str(const ConcatString &, + const ConcatString &); static void write_nc_var_int(const char *, const char *, int); -static void add_var_att_local(NcVar *, const char *, const ConcatString &); +static void add_var_data_atts(NcVar *, const ConcatString &, + const ConcatString &, + const ConcatString &); +static void add_var_att_local(NcVar *, const char *, + const ConcatString &); static void write_hist_bins(void); static void write_hist1d(void); static void write_hist2d(void); +static void write_wavelengths(void); +static void write_power_spectrum(void); +static void write_error_power_spectrum(void); static void write_info_theory(void); static void clean_up(void); +static double get_grid_res_km(const Grid &); + static Met2dDataFile *get_mtddf(const StringArray &, const int); static void usage(int exit_code=1); @@ -109,10 +135,17 @@ int met_main(int argc, char *argv[]) { // Write information theory output if(conf_info.nc_info.do_info_theory) write_info_theory(); + // Write power spectrum output for the full domain + if(conf_info.nc_info.do_power_spectrum) { + write_wavelengths(); + write_power_spectrum(); + if(multiple_data_sources) write_error_power_spectrum(); + } + // Write benchmarking metrics - #ifdef WITH_PROFILER +#ifdef WITH_PROFILER ctrack::result_print(); - #endif +#endif // Close files and deallocate memory clean_up(); @@ -173,8 +206,8 @@ static void process_command_line(int argc, char **argv) { // List the config files mlog << Debug(1) - << "Default Config File: " << default_config_file << "\n" - << "User Config File: " << config_file << "\n"; + << "Default Config File: " << default_config_file << "\n" + << "User Config File: " << config_file << "\n"; // Read the config files conf_info.read_config(default_config_file.c_str(), @@ -230,7 +263,7 @@ static void process_command_line(int argc, char **argv) { } // end for i // Process the configuration - conf_info.process_config(file_types); + conf_info.process_config(file_types, data_grid); // Determine the verification grid grid = parse_vx_grid(conf_info.data_info[0]->regrid(), @@ -259,21 +292,22 @@ string get_tool_name() { //////////////////////////////////////////////////////////////////////// static void setup_diag_info(void) { - #ifdef WITH_PROFILER +#ifdef WITH_PROFILER CTRACK; - #endif +#endif // Resize based on the number of variables and masks diag_info.resize(conf_info.get_n_data()); for(auto &info : diag_info) info.resize(conf_info.get_n_mask()); + power_info.resize(conf_info.get_n_data()); // Loop over variables for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { // Find bin ranges - const VarInfo *i_data = conf_info.data_info[i_var]; - NumArray range(i_data->range()); - int i_n_bins = i_data->n_bins(); + const VarInfo *i_vinfo = conf_info.data_info[i_var]; + NumArray range(i_vinfo->range()); + int i_n_bins = i_vinfo->n_bins(); double var_min = range[0]; double var_max = range[1]; double bin_delta = (var_max - var_min) / i_n_bins; @@ -290,27 +324,27 @@ static void setup_diag_info(void) { // 1D histogram mlog << Debug(2) - << "Initializing " << i_data->magic_str_attr() + << "Initializing " << i_vinfo->magic_str_attr() << " histogram with " << i_n_bins << " bins from " << var_min << " to " << var_max << ".\n"; vector hist1d; init_pdf(i_n_bins, hist1d); // Keep track of unique output variable names - if(nc_var_sa.has(i_data->magic_str_attr())) unique_variable_names = false; - nc_var_sa.add(i_data->magic_str_attr()); + if(nc_var_sa.has(i_vinfo->magic_str_attr())) unique_variable_names = false; + nc_var_sa.add(i_vinfo->magic_str_attr()); // 2D histograms map > hist2d; for(int j_var=i_var+1; j_var < conf_info.get_n_data(); j_var++) { - const VarInfo *j_data = conf_info.data_info[j_var]; + const VarInfo *j_vinfo = conf_info.data_info[j_var]; - int j_n_bins = j_data->n_bins(); + int j_n_bins = j_vinfo->n_bins(); mlog << Debug(2) - << "Initializing " << i_data->magic_str_attr() << "_" - << j_data->magic_str_attr() << " joint histogram with " + << "Initializing " << i_vinfo->magic_str_attr() << "_" + << j_vinfo->magic_str_attr() << " joint histogram with " << i_n_bins << " x " << j_n_bins << " bins.\n"; hist2d[j_var] = vector(); @@ -336,15 +370,14 @@ static void setup_diag_info(void) { //////////////////////////////////////////////////////////////////////// static void process_series(void) { - vector data_dp(conf_info.get_n_data()); + vector in_data(conf_info.get_n_data()); const StringArray *cur_files; const GrdFileType *cur_ftype; - Grid cur_grid; // List the lengths of the series options mlog << Debug(1) - << "Processing " << conf_info.get_n_data() << " data fields" - << " from " << n_series << " input file(s).\n"; + << "Processing " << conf_info.get_n_data() << " data fields" + << " from " << n_series << " input file(s).\n"; // Loop over the input files for(int i_series=0; i_series < n_series; i_series++) { @@ -357,7 +390,7 @@ static void process_series(void) { // Read the input data for this series entry for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - VarInfo *i_data = conf_info.data_info[i_var]; + VarInfo *i_vinfo = conf_info.data_info[i_var]; // Check for separate data files for each field if(data_files.size() > 1) { @@ -370,58 +403,34 @@ static void process_series(void) { } mlog << Debug(2) - << "Reading field " << i_data->magic_str_attr() + << "Reading field " << i_vinfo->magic_str_attr() << " data from file: " << (*cur_files)[i_series] << "\n"; - get_series_entry(i_series, i_data, *cur_files, *cur_ftype, - data_dp[i_var], cur_grid); - - // Regrid, if necessary - if(!(cur_grid == grid)) { - mlog << Debug(2) - << "Regridding field " << i_data->magic_str_attr() - << " to the verification grid using " - << i_data->regrid().get_str() << ".\n"; - data_dp[i_var] = met_regrid(data_dp[i_var], - cur_grid, grid, - i_data->regrid()); - } + // Read the series data + read_series_data(i_series, i_vinfo, *cur_files, *cur_ftype, + in_data[i_var]); - // Initialize time ranges - if(i_series == 0 && i_var == 0) { - init_beg = init_end = data_dp[i_var].init(); - valid_beg = valid_end = data_dp[i_var].valid(); - lead_beg = lead_end = data_dp[i_var].lead(); - } - // Update time ranges - else { - if(data_dp[i_var].init() < init_beg) { - init_beg = data_dp[i_var].init(); - } - if(data_dp[i_var].init() > init_end) { - init_end = data_dp[i_var].init(); - } - if(data_dp[i_var].valid() < valid_beg) { - valid_beg = data_dp[i_var].valid(); - } - if(data_dp[i_var].valid() > valid_end) { - valid_end = data_dp[i_var].valid(); - } - if(data_dp[i_var].lead() < lead_beg) { - lead_beg = data_dp[i_var].lead(); - } - if(data_dp[i_var].lead() > lead_end) { - lead_end = data_dp[i_var].lead(); - } - } } // end for i_var // Process the 1D histograms - process_hist1d(data_dp); + process_hist1d(in_data); // Process the 2D histograms - process_hist2d(data_dp); + process_hist2d(in_data); + + // Process the power spectrum + if(conf_info.nc_info.do_power_spectrum) { + + // Prepare power spectrum data + prepare_power_spectrum_data(in_data); + + // Process the power spectrum + process_power_spectrum(in_data); + + // Process the error power spectrum + if(multiple_data_sources) process_error_power_spectrum(in_data); + } } // end for i_series @@ -429,15 +438,15 @@ static void process_series(void) { if(conf_info.nc_info.do_info_theory) process_info_theory(); } - + //////////////////////////////////////////////////////////////////////// -static void process_hist1d(const vector &data_dp) { +static void process_hist1d(const vector &in_data) { // Update the 1D histogram counts for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; // Loop over the masks for(int i_mask=0; i_mask < conf_info.get_n_mask(); i_mask++) { @@ -445,7 +454,7 @@ static void process_hist1d(const vector &data_dp) { DiagInfo *i_diag = &diag_info[i_var][i_mask]; // Apply the mask before updating the data ranges - DataPlane dp(data_dp[i_var]); + DataPlane dp(in_data[i_var].dp); apply_mask(dp, conf_info.mask_mp[i_mask]); double min; double max; @@ -461,27 +470,27 @@ static void process_hist1d(const vector &data_dp) { update_pdf(i_diag->bin_min[0], i_diag->bin_delta, i_diag->hist1d, - data_dp[i_var], + in_data[i_var].dp, conf_info.mask_mp[i_mask]); mlog << Debug(2) - << "Processed " << i_data->magic_str_attr() + << "Processed " << i_vinfo->magic_str_attr() << " data over region " << conf_info.mask_name[i_mask] << " with range (" << i_diag->var_min << ", " << i_diag->var_max << ") into bins with range (" - << i_data->range()[0] << ", " - << i_data->range()[1] << ").\n"; + << i_vinfo->range()[0] << ", " + << i_vinfo->range()[1] << ").\n"; // Compare input data and bin ranges - if(i_diag->var_min < i_data->range()[0] || - i_diag->var_max > i_data->range()[1]) { + if(i_diag->var_min < i_vinfo->range()[0] || + i_diag->var_max > i_vinfo->range()[1]) { mlog << Warning << "\nprocess_hist1d() -> " - << "the range of the " << i_data->magic_str_attr() + << "the range of the " << i_vinfo->magic_str_attr() << " data over region " << conf_info.mask_name[i_mask] << " (" << i_diag->var_min << ", " << i_diag->var_max << ") falls outside the configuration file range (" - << i_data->range()[0] << ", " - << i_data->range()[1] << ")!\n\n"; + << i_vinfo->range()[0] << ", " + << i_vinfo->range()[1] << ")!\n\n"; } } // end for i_mask } // end for i_var @@ -489,16 +498,16 @@ static void process_hist1d(const vector &data_dp) { //////////////////////////////////////////////////////////////////////// -static void process_hist2d(const vector &data_dp) { +static void process_hist2d(const vector &in_data) { // Process the 2D joint histograms for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; for(int j_var=i_var+1; j_var < conf_info.get_n_data(); j_var++) { - const VarInfo *j_data = conf_info.data_info[j_var]; + const VarInfo *j_vinfo = conf_info.data_info[j_var]; for(int i_mask=0; i_mask < conf_info.get_n_mask(); i_mask++) { @@ -506,14 +515,14 @@ static void process_hist2d(const vector &data_dp) { DiagInfo *j_diag = &diag_info[j_var][i_mask]; // Update 2D histogram counts - update_joint_pdf(i_data->n_bins(), - j_data->n_bins(), + update_joint_pdf(i_vinfo->n_bins(), + j_vinfo->n_bins(), i_diag->bin_min[0], j_diag->bin_min[0], i_diag->bin_delta, j_diag->bin_delta, i_diag->hist2d[j_var], - data_dp[i_var], data_dp[j_var], + in_data[i_var].dp, in_data[j_var].dp, conf_info.mask_mp[i_mask]); } // end for i_mask } // end for j_var @@ -522,6 +531,262 @@ static void process_hist2d(const vector &data_dp) { //////////////////////////////////////////////////////////////////////// +static void prepare_power_spectrum_data(vector &in_data) { + + for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + + int nxy = in_data[i_var].dp.nxy(); + int n_vld = in_data[i_var].dp.n_good_data(); + + // Check for bad data + if(n_vld < nxy) { + + PowerSpectrumInfo *ps_ptr = &conf_info.ps_info[i_var]; + + // Skip power spectrum for missing data + if(ps_ptr->missing_flag == MissingDataType::None || + (double) n_vld/nxy < ps_ptr->vld_thresh) { + + mlog << Debug(3) << "Skipping " + << conf_info.data_info[i_var]->magic_str_attr() + << " power spectrum due to missing data.\n"; + ps_ptr->skip = true; + } + else { + + // Initialize to the config file value + double fill; + ConcatString desc; + + // Replace with the mean of the field + if(ps_ptr->missing_flag == MissingDataType::Mean) { + fill = in_data[i_var].dp.mean(); + desc = "mean"; + } + // Use config file value + else { + fill = ps_ptr->missing_value; + desc = "constant"; + } + + mlog << Debug(3) << "Replacing " + << conf_info.data_info[i_var]->magic_str_attr() + << " power spectrum missing data with a " << desc + << " value of " << fill << ".\n"; + in_data[i_var].dp.replace_bad_data(fill); + } + } + } // end for i_var +} + +//////////////////////////////////////////////////////////////////////// + +static void process_power_spectrum(const vector &in_data) { + + // Process the power spectrum for the full input domain + for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + + // Check skip + if(conf_info.ps_info[i_var].skip) continue; + + // Process U/V for kinetic energy + if(in_data[i_var].uv_flag) { + + // Apply the discrete cosine transforms + DataPlane u_dct_dp(dct_typeII(in_data[i_var].u_dp)); + DataPlane v_dct_dp(dct_typeII(in_data[i_var].v_dp)); + + // Compute the radial energy + vector u_re = radial_spectral_variance(u_dct_dp); + vector v_re = radial_spectral_variance(v_dct_dp); + + // Combine the components + vector re(u_re.size()); + for(size_t i=0; i re = radial_spectral_variance(dct_dp); + + // Sum the spectral variance + sum_spectral_variance(power_info[i_var].power, re); + } + } // end for i_var +} + +//////////////////////////////////////////////////////////////////////// + +static void process_error_power_spectrum(const vector &in_data) { + + // Process the power spectrum for the full input domain + for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + + // Check skip + if(conf_info.ps_info[i_var].skip) continue; + + for(int j_var=i_var+1; j_var < conf_info.get_n_data(); j_var++) { + + // Check skip + if(conf_info.ps_info[j_var].skip) continue; + + // Process U/V for kinetic energy + if(in_data[i_var].uv_flag && in_data[j_var].uv_flag) { + + // Compute difference fields + DataPlane u_diff_dp(subtract(in_data[i_var].u_dp, + in_data[j_var].u_dp)); + DataPlane v_diff_dp(subtract(in_data[i_var].v_dp, + in_data[j_var].v_dp)); + + // Apply the discrete cosine transforms + DataPlane u_dct_dp(dct_typeII(u_diff_dp)); + DataPlane v_dct_dp(dct_typeII(v_diff_dp)); + + // Compute the radial energy + vector u_re = radial_spectral_variance(u_dct_dp); + vector v_re = radial_spectral_variance(v_dct_dp); + + // Combine the components + vector re(u_re.size()); + for(size_t i=0; i re = radial_spectral_variance(dct_dp); + + // Sum the spectral variance + sum_spectral_variance(power_info[i_var].error_power[j_var], re); + } + } // end for j_var + } // end for i_var +} + +//////////////////////////////////////////////////////////////////////// +// +// Compute the radial spectral variance from a 2D array of DCT-II +// coefficients, following Denis et al. (2002). +// - Each axis is normalized independently by its own dimension +// (alph = sqrt((x/nx)^2 + (y/ny)^2)), giving an elliptical +// (rather than circular) truncation appropriate for rectangular +// domains, per section 2b of Denis et al. (2002). +// - Points with alph >= 1 fall outside the valid isotropic +// wavenumber range (the "corners" of the coefficient array) and +// are excluded. +// - The (0,0) DC coefficient and the entire near-DC shell +// (alph < 1/N) are excluded, per the k=1..N-1 range above. +// +//////////////////////////////////////////////////////////////////////// + +static vector radial_spectral_variance(const DataPlane &dp) { + + // Use the smaller dimension + int nx = dp.nx(); + int ny = dp.ny(); + int N = min(nx, ny); + + // Output spans k = 1 .. N-1 inclusive, per Denis et al. (2002). + // re[k-1] holds the energy for wavenumber k. + int n_bins = N - 1; + vector re(n_bins, 0.0); + + double inv_N = 1.0 / N; + + // Accumulate spectral variance for each DCT coefficient + for(int x=0; x= 1) -- neither + // corresponds to a valid k in {1, ..., N-1}. + if(alph < inv_N || alph >= 1.0) continue; + + // Map alph to wavenumber k in {1, ..., N-1}, then to the + // 0-based array index (k - 1). + auto k = (int) (alph * N); + if(k >= N) k = N - 1; + + re[k - 1] += dp(x,y) * dp(x,y) / dp.nxy(); + + } // end for y + } // end for x + + return re; +} + +//////////////////////////////////////////////////////////////////////// + +static void sum_spectral_variance(vector &sum, + const vector &cur) { + + // Initialize the sum if needed + if(sum.empty()) { + sum = cur; + } + // Otherwise, accumulate values + else { + + // Must be the same size + if(sum.size() != cur.size()) { + mlog << Error << "\nsum_spectral_variance() -> " + << "vector lengths do not match (" << sum.size() + << " != " << cur.size() << ")!\n\n"; + exit(1); + } + + // Increment each element + for(int i=0; i hist2d_i_sum(i_data->n_bins(), 0); - vector hist2d_j_sum(j_data->n_bins(), 0); + vector hist2d_i_sum(i_vinfo->n_bins(), 0); + vector hist2d_j_sum(j_vinfo->n_bins(), 0); - for(int i=0; in_bins(); i++) { - for(int j=0; jn_bins(); j++) { + for(int i=0; in_bins(); i++) { + for(int j=0; jn_bins(); j++) { - int n = i * j_data->n_bins() + j; + int n = i * j_vinfo->n_bins() + j; // Increment sums hist2d_ij_sum += i_diag->hist2d[j_var][n]; @@ -584,15 +849,15 @@ static void process_info_theory() { } // end for i // Compute probabilities and acccumulate mutual information - for(int i=0; in_bins(); i++) { + for(int i=0; in_bins(); i++) { auto p_i = (double) hist2d_i_sum[i] / (double) hist2d_ij_sum; - for(int j=0; jn_bins(); j++) { + for(int j=0; jn_bins(); j++) { auto p_j = (double) hist2d_j_sum[j] / (double) hist2d_ij_sum; - int n = i * j_data->n_bins() + j; + int n = i * j_vinfo->n_bins() + j; auto p_ij = (double) i_diag->hist2d[j_var][n] / (double) hist2d_ij_sum; @@ -613,6 +878,105 @@ static void process_info_theory() { //////////////////////////////////////////////////////////////////////// +static void read_series_data(int i_series, VarInfo *i_vinfo, + const StringArray &in_files, + const GrdFileType in_ftype, + InputDataInfo &in_data) { + Grid cur_grid; + + // Read the requested field + get_series_entry(i_series, i_vinfo, in_files, in_ftype, + in_data.dp, cur_grid); + + // Regrid, if needed + regrid_data(i_vinfo, cur_grid, in_data.dp); + + // Check for kinetic energy + in_data.uv_flag = i_vinfo->is_kinetic_energy(); + + // Read U/V components for kinetic energy + if(in_data.uv_flag && + (!read_wind_series(i_vinfo->wind_info().u_wind, + i_series, i_vinfo, in_files, in_ftype, in_data.u_dp) || + !read_wind_series(i_vinfo->wind_info().v_wind, + i_series, i_vinfo, in_files, in_ftype, in_data.v_dp))) { + mlog << Error << "\nread_series_data() -> " + << "trouble reading U/V wind components " + << "for kinetic energy field \"" + << i_vinfo->magic_str() << "\"!\n\n"; + exit(1); + } + + // Initialize time ranges + auto dp = &in_data.dp; + if(valid_beg == 0 && valid_end == 0) { + init_beg = init_end = dp->init(); + valid_beg = valid_end = dp->valid(); + lead_beg = lead_end = dp->lead(); + } + // Update time ranges + else { + if(dp->init() < init_beg) init_beg = dp->init(); + if(dp->init() > init_end) init_end = dp->init(); + if(dp->valid() < valid_beg) valid_beg = dp->valid(); + if(dp->valid() > valid_end) valid_end = dp->valid(); + if(dp->lead() < lead_beg) lead_beg = dp->lead(); + if(dp->lead() > lead_end) lead_end = dp->lead(); + } +} + +//////////////////////////////////////////////////////////////////////// + +static bool read_wind_series(const StringArray &wind_names, + int i_series, const VarInfo *i_vinfo, + const StringArray &in_files, + const GrdFileType in_ftype, + DataPlane &wind_dp) { + bool found = false; + auto wind_vinfo = i_vinfo->clone(); + Grid cur_grid; + + // Loop over possible wind names + for(int i=0; ireset_dict_with_name(wind_names[i].c_str()) && + get_series_entry(i_series, wind_vinfo.get(), in_files, in_ftype, + wind_dp, cur_grid)) { + found = true; + + mlog << Debug(3) << "Found matching wind field \"" + << wind_vinfo->magic_str() << "\".\n"; + + // Regrid, if needed + regrid_data(i_vinfo, cur_grid, wind_dp); + + break; + } + } + + return found; +} + +//////////////////////////////////////////////////////////////////////// + +static void regrid_data(const VarInfo *vinfo, const Grid &cur_grid, + DataPlane &dp) { + + // Check for grid match + if(cur_grid == grid) return; + + mlog << Debug(2) + << "Regridding field " << vinfo->magic_str_attr() + << " to the verification grid using " + << vinfo->regrid().get_str() << ".\n"; + dp = met_regrid(dp, cur_grid, grid, vinfo->regrid()); + + return; +} + +//////////////////////////////////////////////////////////////////////// + static ConcatString get_nc_var_str(const VarInfo *info, int index) { ConcatString cs; @@ -630,6 +994,21 @@ static ConcatString get_nc_var_str(const VarInfo *info, int index) { //////////////////////////////////////////////////////////////////////// +static ConcatString get_nc_att_str(const ConcatString &cs1, + const ConcatString &cs2) { + + // Return one if equal + if(cs1 == cs2) return cs1; + + // Otherwise, include both strings + ConcatString cs(cs1); + cs << " and " << cs2; + + return cs; +} + +//////////////////////////////////////////////////////////////////////// + static void setup_nc_file(void) { // Create NetCDF file @@ -685,6 +1064,12 @@ static void setup_nc_file(void) { int mask_size = conf_info.mask_mp[i_mask].count(); mask_size_var.putVar(offsets, counts, &mask_size); } + + // Add the power spectra dimension + if(conf_info.nc_info.do_power_spectrum) { + wavenumber_dim = add_dim(nc_out, "wavenumber", + (long) min(grid.nx(), grid.ny())-1); + } } //////////////////////////////////////////////////////////////////////// @@ -705,6 +1090,21 @@ static void write_nc_var_int(const char *var_name, //////////////////////////////////////////////////////////////////////// +static void add_var_data_atts(NcVar *var, + const ConcatString &long_name_cs, + const ConcatString &level_cs, + const ConcatString &units_cs) { + + // Add variable attributes for long_name, level, and units + if(var) { + add_var_att_local(var, "long_name", long_name_cs); + add_var_att_local(var, "level", level_cs); + add_var_att_local(var, "units", units_cs); + } +} + +//////////////////////////////////////////////////////////////////////// + static void add_var_att_local(NcVar *var, const char *att_name, const ConcatString &att_value) { if(att_value.nonempty()) add_att(var, att_name, att_value.c_str()); @@ -717,15 +1117,15 @@ static void write_hist_bins(void) { for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; DiagInfo *i_diag = &diag_info[i_var][0]; // Define NetCDF variable name - ConcatString var_str(get_nc_var_str(i_data, i_var+1)); + ConcatString var_str(get_nc_var_str(i_vinfo, i_var+1)); // Define NetCDF dimensions NcDim var_dim = add_dim(nc_out, var_str, - (long) i_data->n_bins()); + (long) i_vinfo->n_bins()); data_var_dims.emplace_back(var_dim); // Create NetCDF variable @@ -745,16 +1145,11 @@ static void write_hist_bins(void) { // Add variable attributes ConcatString cs; cs << cs_erase << "Minimum value of " << var_str << " bin"; - add_var_att_local(&var_min, "long_name", cs); - add_var_att_local(&var_min, "units", i_data->units_attr()); - + add_var_data_atts(&var_min, cs, i_vinfo->level_attr(), i_vinfo->units_attr()); cs << cs_erase << "Maximum value of " << var_str << " bin"; - add_var_att_local(&var_max, "long_name", cs); - add_var_att_local(&var_max, "units", i_data->units_attr()); - + add_var_data_atts(&var_max, cs, i_vinfo->level_attr(), i_vinfo->units_attr()); cs << cs_erase << "Midpoint value of " << var_str << " bin"; - add_var_att_local(&var_mid, "long_name", cs); - add_var_att_local(&var_mid, "units", i_data->units_attr()); + add_var_data_atts(&var_mid, cs, i_vinfo->level_attr(), i_vinfo->units_attr()); // Write bin values for the current variable var_min.putVar(i_diag->bin_min.data()); @@ -773,10 +1168,10 @@ static void write_hist1d(void) { // Define and write 1D histograms for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; // Define NetCDF variable name - ConcatString var_str(get_nc_var_str(i_data, i_var+1)); + ConcatString var_str(get_nc_var_str(i_vinfo, i_var+1)); ConcatString var_name("hist_"); var_name << var_str; @@ -790,7 +1185,7 @@ static void write_hist1d(void) { // Add variable attributes ConcatString cs; cs << "Histogram of " << var_str << " values"; - add_var_att_local(&var, "long_name", cs); + add_var_data_atts(&var, cs, i_vinfo->level_attr(), i_vinfo->units_attr()); // Write 1D histogram for each mask for(int i_mask=0; i_mask < conf_info.get_n_mask(); i_mask++) { @@ -800,7 +1195,7 @@ static void write_hist1d(void) { offsets[0] = i_mask; offsets[1] = 0; counts[0] = 1; - counts[1] = i_data->n_bins(); + counts[1] = i_vinfo->n_bins(); var.putVar(offsets, counts, hist); @@ -817,17 +1212,17 @@ static void write_hist2d(void) { // Define and write 2D joint histograms for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; for(int j_var=i_var+1; j_var < conf_info.get_n_data(); j_var++) { - const VarInfo *j_data = conf_info.data_info[j_var]; + const VarInfo *j_vinfo = conf_info.data_info[j_var]; // Define NetCDF variable name ConcatString var_str; - var_str << get_nc_var_str(i_data, i_var+1) << "_" - << get_nc_var_str(j_data, j_var+1); - ConcatString var_name("hist_"); + var_str << get_nc_var_str(i_vinfo, i_var+1) << "_" + << get_nc_var_str(j_vinfo, j_var+1); + ConcatString var_name("hist_"); var_name << var_str; // Create NetCDF variable @@ -841,7 +1236,11 @@ static void write_hist2d(void) { // Add variable attributes ConcatString cs; cs << "Joint histogram of " << var_str << " values"; - add_var_att_local(&var, "long_name", cs); + add_var_data_atts(&var, cs, + get_nc_att_str(i_vinfo->level_attr(), + j_vinfo->level_attr()), + get_nc_att_str(i_vinfo->units_attr(), + j_vinfo->units_attr())); // Write 2D histogram for each mask for(int i_mask=0; i_mask < conf_info.get_n_mask(); i_mask++) { @@ -852,8 +1251,8 @@ static void write_hist2d(void) { offsets[1] = 0; offsets[2] = 0; counts[0] = 1; - counts[1] = i_data->n_bins(); - counts[2] = j_data->n_bins(); + counts[1] = i_vinfo->n_bins(); + counts[2] = j_vinfo->n_bins(); var.putVar(offsets, counts, hist); @@ -870,10 +1269,10 @@ static void write_info_theory(void) { // Write entropy for each 1D histogram for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; // Define NetCDF variable name - ConcatString var_str(get_nc_var_str(i_data, i_var+1)); + ConcatString var_str(get_nc_var_str(i_vinfo, i_var+1)); ConcatString var_name("entropy_"); var_name << var_str; @@ -884,8 +1283,7 @@ static void write_info_theory(void) { // Add variable attributes ConcatString cs; cs << "Entropy value for " << var_str; - add_var_att_local(&var, "long_name", cs); - add_var_att_local(&var, "units", units_cs); + add_var_data_atts(&var, cs, i_vinfo->level_attr(), units_cs); // Store the data vector data(conf_info.get_n_mask()); @@ -901,18 +1299,18 @@ static void write_info_theory(void) { // Write joint entropy and mutual information for each 2D joint histogram for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { - const VarInfo *i_data = conf_info.data_info[i_var]; + const VarInfo *i_vinfo = conf_info.data_info[i_var]; for(int j_var=i_var+1; j_var < conf_info.get_n_data(); j_var++) { - const VarInfo *j_data = conf_info.data_info[j_var]; + const VarInfo *j_vinfo = conf_info.data_info[j_var]; ConcatString var_str; - var_str << get_nc_var_str(i_data, i_var+1) << "_" - << get_nc_var_str(j_data, j_var+1); + var_str << get_nc_var_str(i_vinfo, i_var+1) << "_" + << get_nc_var_str(j_vinfo, j_var+1); // Define NetCDF variable names - ConcatString je_var_name("joint_entropy_"); + ConcatString je_var_name("joint_entropy_"); je_var_name << var_str; ConcatString mi_var_name("mutual_information_"); mi_var_name << var_str; @@ -923,26 +1321,165 @@ static void write_info_theory(void) { NcVar mi_var = add_var(nc_out, mi_var_name, ncFloat, mask_dim, deflate_level); + // Level attribute + ConcatString level_cs(get_nc_att_str(i_vinfo->level_attr(), + j_vinfo->level_attr())); + // Add variable attributes ConcatString cs; cs << "Joint entropy value for " << var_str; - add_var_att_local(&je_var, "long_name", cs); - add_var_att_local(&je_var, "units", units_cs); + add_var_data_atts(&je_var, cs, level_cs, units_cs); cs << cs_erase << "Mutual information value for " << var_str; - add_var_att_local(&mi_var, "long_name", cs); - add_var_att_local(&mi_var, "units", units_cs); + add_var_data_atts(&mi_var, cs, level_cs, units_cs); // Store the data vector je_data(conf_info.get_n_mask()); - vector mi_data(conf_info.get_n_mask()); + vector mi_vinfo(conf_info.get_n_mask()); for(int i_mask=0; i_mask < conf_info.get_n_mask(); i_mask++) { je_data[i_mask] = diag_info[i_var][i_mask].joint_entropy[j_var]; - mi_data[i_mask] = diag_info[i_var][i_mask].mutual_information[j_var]; + mi_vinfo[i_mask] = diag_info[i_var][i_mask].mutual_information[j_var]; } // Write the data je_var.putVar(je_data.data()); - mi_var.putVar(mi_data.data()); + mi_var.putVar(mi_vinfo.data()); + + } // end for j_var + } // end for i_var +} + +//////////////////////////////////////////////////////////////////////// + +static void write_wavelengths(void) { + + // Define wavenumbers and wavelengths based on the number of bins + int N = min(grid.nx(), grid.ny()); + int n_bins = N - 1; + double grid_res_km = get_grid_res_km(grid); + + vector wavenumber(n_bins); + vector wavelength(n_bins); + + for(int i=0; iunits_attr() << ")^2"; + add_var_data_atts(&var, cs, i_vinfo->level_attr(), units_cs); + add_var_att_local(&var, "mask", full_domain_str); + + // Write power spectrum + PowerInfo *p_diag = &power_info[i_var]; + + // Divide sums by the series length + for(auto &x : p_diag->power) x /= n_series; + + // Write the mean power data + var.putVar(p_diag->power.data()); + + } // end for i_var +} + +//////////////////////////////////////////////////////////////////////// + +static void write_error_power_spectrum(void) { + + // Define and write the error power spectrum + for(int i_var=0; i_var < conf_info.get_n_data(); i_var++) { + + // Check skip + if(conf_info.ps_info[i_var].skip) continue; + + const VarInfo *i_vinfo = conf_info.data_info[i_var]; + + for(int j_var=i_var+1; j_var < conf_info.get_n_data(); j_var++) { + + // Check skip + if(conf_info.ps_info[j_var].skip) continue; + + const VarInfo *j_vinfo = conf_info.data_info[j_var]; + + // Define NetCDF variable name + ConcatString i_var_str(get_nc_var_str(i_vinfo, i_var+1)); + ConcatString j_var_str(get_nc_var_str(j_vinfo, j_var+1)); + ConcatString var_name("error_power_spectrum_"); + var_name << i_var_str << "_" << j_var_str; + + // Create NetCDF variable + NcVar var = add_var(nc_out, var_name, ncFloat, wavenumber_dim, + deflate_level); + + // Level attribute + ConcatString level_cs(get_nc_att_str(i_vinfo->level_attr(), + j_vinfo->level_attr())); + + // Units attribute + ConcatString units_cs("("); + units_cs << get_nc_att_str(i_vinfo->units_attr(), + j_vinfo->units_attr()) << ")^2"; + + // Add variable attributes + ConcatString cs("Power spectrum of errors for "); + cs << i_var_str << " minus " << j_var_str; + add_var_data_atts(&var, cs, level_cs, units_cs); + add_var_att_local(&var, "mask", full_domain_str); + + // Write power spectrum + PowerInfo *p_diag = &power_info[i_var]; + + // Divide sums by the series length + for(auto &x : p_diag->error_power[j_var]) x /= n_series; + + // Write the mean error power + var.putVar(p_diag->error_power[j_var].data()); } // end for j_var } // end for i_var @@ -952,20 +1489,16 @@ static void write_info_theory(void) { static Met2dDataFile *get_mtddf(const StringArray &file_list, const int i_field) { - Met2dDataFile *mtddf = nullptr; - Dictionary *dict = nullptr; - Dictionary i_dict; - GrdFileType file_type; int i; // Conf: data.field - dict = conf_info.conf.lookup_array(conf_key_data_field); + Dictionary *dict = conf_info.conf.lookup_array(conf_key_data_field); // Get the i-th data.field entry - i_dict = parse_conf_i_vx_dict(dict, i_field); + Dictionary i_dict = parse_conf_i_vx_dict(dict, i_field); // Look for file_type in the i-th data.field entry - file_type = parse_conf_file_type(&i_dict); + GrdFileType file_type = parse_conf_file_type(&i_dict); // Find the first file that actually exists for(i=0; i < file_list.n(); i++) { @@ -980,6 +1513,7 @@ static Met2dDataFile *get_mtddf(const StringArray &file_list, } // Read first valid file + Met2dDataFile *mtddf = nullptr; if(!(mtddf = Met2dDataFileFactory::new_met_2d_data_file( file_list[i].c_str(), file_type))) { mlog << Error << "\nget_mtddf() -> " @@ -1008,6 +1542,57 @@ static void clean_up(void) { return; } +//////////////////////////////////////////////////////////////////////// +// +// Estimate the grid spacking in km +// +//////////////////////////////////////////////////////////////////////// + +static double get_grid_res_km(const Grid &g) { + double res_km; + + // Use the grid scale, if well-defined + if(g.scale_km() > 0) { + res_km = g.scale_km(); + } + // Otherwise, determine the scale from the grid points + else { + + int nx = grid.nx(); + int ny = grid.ny(); + + // X-spacing at the center + double lat1; + double lon1; + g.xy_to_latlon(0, (ny - 1)/2.0, lat1, lon1); + double lat2; + double lon2; + g.xy_to_latlon(nx - 1, (ny - 1)/2.0, lat2, lon2); + double dx_km = gc_dist(lat1, lon1, lat2, lon2) / (nx - 1); + + // Y-spacing at the center + g.xy_to_latlon((nx - 1)/2.0, 0, lat1, lon1); + g.xy_to_latlon((nx - 1)/2.0, ny - 1, lat2, lon2); + double dy_km = gc_dist(lat1, lon1, lat2, lon2) / (ny - 1); + + // Log message when grid spacing differs + if(!is_eq(dx_km, dy_km, 0.1)) { + mlog << Debug(3) << "Grid spacing in the X (" << dx_km + << " km) and Y (" << dy_km << " km) dimensions differ.\n"; + } + res_km = min(dx_km, dy_km); + } + + // Attempt to round to the nearest integer + if(is_eq(res_km, (double) nint(res_km), 0.1)) { + res_km = (double) nint(res_km); + } + + mlog << Debug(3) << "Using grid spacing of " << res_km << " km.\n"; + + return res_km; +} + //////////////////////////////////////////////////////////////////////// __attribute__((noreturn)) static void usage(int exit_code) { @@ -1053,8 +1638,8 @@ __attribute__((noreturn)) static void usage(int exit_code) { //////////////////////////////////////////////////////////////////////// static void set_data_files(const StringArray & a) { - data_files.emplace_back(a); if(!data_files.empty()) multiple_data_sources = true; + data_files.emplace_back(a); } //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/other/grid_diag/grid_diag.h b/src/tools/other/grid_diag/grid_diag.h index 810b9aee50..6b0701b2c6 100644 --- a/src/tools/other/grid_diag/grid_diag.h +++ b/src/tools/other/grid_diag/grid_diag.h @@ -87,6 +87,7 @@ static GridDiagConfInfo conf_info; // Output NetCDF file static netCDF::NcFile *nc_out = nullptr; netCDF::NcDim mask_dim; +netCDF::NcDim wavenumber_dim; std::vector data_var_dims; int deflate_level; @@ -108,6 +109,14 @@ static Grid grid; // Input files static Met2dDataFile *data_mtddf = nullptr; +// Struct to store scalar and vector input data +struct InputDataInfo { + DataPlane dp; + DataPlane u_dp; + DataPlane v_dp; + bool uv_flag = false; +}; + // Struct to store diagnostic info for each field and masking region struct DiagInfo { @@ -117,7 +126,7 @@ struct DiagInfo { std::vector bin_mid; double bin_delta; - // Input data info + // Range of input data values double var_min; double var_max; @@ -131,9 +140,20 @@ struct DiagInfo { std::map mutual_information; }; -// DiagInfo objects [n_data][n_mask] +// DiagInfo objects, sized as n_data by n_mask std::vector > diag_info; +// Struct to store power spectrum information +struct PowerInfo { + + // Power spectrum + std::vector power; + std::map > error_power; +}; + +// PowerInfo objects, sized as n_data +std::vector power_info; + // Series length static int n_series = bad_data_int; diff --git a/src/tools/other/grid_diag/grid_diag_conf_info.cc b/src/tools/other/grid_diag/grid_diag_conf_info.cc index 068e9f2d3d..a41e1b8568 100644 --- a/src/tools/other/grid_diag/grid_diag_conf_info.cc +++ b/src/tools/other/grid_diag/grid_diag_conf_info.cc @@ -46,7 +46,8 @@ void GridDiagNcOutInfo::clear() { bool GridDiagNcOutInfo::all_false() const { - bool status = do_hist1d || do_hist2d || do_info_theory; + bool status = do_hist1d || do_hist2d || + do_info_theory || do_power_spectrum; return !status; } @@ -55,9 +56,10 @@ bool GridDiagNcOutInfo::all_false() const { void GridDiagNcOutInfo::set_all_false() { - do_hist1d = false; - do_hist2d = false; - do_info_theory = false; + do_hist1d = false; + do_hist2d = false; + do_info_theory = false; + do_power_spectrum = false; return; } @@ -66,9 +68,10 @@ void GridDiagNcOutInfo::set_all_false() { void GridDiagNcOutInfo::set_all_true() { - do_hist1d = true; - do_hist2d = true; - do_info_theory = true; + do_hist1d = true; + do_hist2d = true; + do_info_theory = true; + do_power_spectrum = true; return; } @@ -113,6 +116,7 @@ void GridDiagConfInfo::clear() { if(info) { delete info; info = nullptr; } } data_info.clear(); + ps_info.clear(); return; } @@ -153,11 +157,13 @@ void GridDiagConfInfo::set_n_data() { // Allocate space based on the number of verification tasks data_info.resize(n_data, nullptr); + ps_info.resize(n_data); } //////////////////////////////////////////////////////////////////////// -void GridDiagConfInfo::process_config(vector file_types) { +void GridDiagConfInfo::process_config(vector file_types, + Grid &data_grid) { ConcatString s; StringArray sa; Dictionary i_dict; @@ -176,7 +182,7 @@ void GridDiagConfInfo::process_config(vector file_types) { Dictionary *dict = conf.lookup_array(conf_key_data_field); // Parse the data field information - for(int i=0; i<(int) data_info.size(); i++) { + for(int i=0; i 1 ? @@ -191,6 +197,11 @@ void GridDiagConfInfo::process_config(vector file_types) { // Set the current dictionaries data_info[i]->set_dict(i_dict); + // Update the grid definition, if needed + if(!data_grid.is_set() && data_info[i]->grid_attr().is_set()) { + data_grid = data_info[i]->grid_attr(); + } + // Dump the contents of the current VarInfo if(mlog.verbosity_level() >= 5) { mlog << Debug(5) @@ -209,6 +220,9 @@ void GridDiagConfInfo::process_config(vector file_types) { exit(1); } + // Conf: power_spectrum + ps_info[i] = parse_conf_power_spectrum(&i_dict); + } // end for i // Conf: output_flag @@ -249,9 +263,10 @@ void GridDiagConfInfo::parse_output_flag() { // Parse the various entries auto d = e->dict_value(); - nc_info.do_hist1d = d->lookup_bool(conf_key_hist1d_flag); - nc_info.do_hist2d = d->lookup_bool(conf_key_hist2d_flag); - nc_info.do_info_theory = d->lookup_bool(conf_key_info_theory_flag); + nc_info.do_hist1d = d->lookup_bool(conf_key_hist1d_flag); + nc_info.do_hist2d = d->lookup_bool(conf_key_hist2d_flag); + nc_info.do_info_theory = d->lookup_bool(conf_key_info_theory_flag); + nc_info.do_power_spectrum = d->lookup_bool(conf_key_power_spectrum_flag); return; } diff --git a/src/tools/other/grid_diag/grid_diag_conf_info.h b/src/tools/other/grid_diag/grid_diag_conf_info.h index d0012d92a5..2fb1393743 100644 --- a/src/tools/other/grid_diag/grid_diag_conf_info.h +++ b/src/tools/other/grid_diag/grid_diag_conf_info.h @@ -32,6 +32,7 @@ struct GridDiagNcOutInfo { bool do_hist1d; bool do_hist2d; bool do_info_theory; + bool do_power_spectrum; GridDiagNcOutInfo(); @@ -61,6 +62,9 @@ class GridDiagConfInfo { std::vector data_info; // VarInfo pointer array [n_data] + // Power spectrum options + std::vector ps_info; // [n_data] + // Masking region names and MaskPlanes StringArray mask_name; std::vector mask_mp; @@ -74,7 +78,7 @@ class GridDiagConfInfo { void read_config(const char *, const char *); void set_n_data(); - void process_config(std::vector); + void process_config(std::vector, Grid &); void parse_output_flag(); void process_masks(const Grid &);