diff --git a/met/README b/met/README index 91965d9cd5..11a9bce6b7 100644 --- a/met/README +++ b/met/README @@ -191,6 +191,7 @@ sub-directory, or the MET Online Tutorial: command line) will contain 36 executables: - ascii2nc - ensemble_stat + - gen_ens_prod - gen_vx_mask - grid_stat - gis_dump_dbf diff --git a/met/configure.ac b/met/configure.ac index cbe955c398..2fc599e47a 100644 --- a/met/configure.ac +++ b/met/configure.ac @@ -485,6 +485,27 @@ else AC_MSG_NOTICE([ensemble_stat will not be compiled]) fi +# gen_ens_prod + +AC_ARG_ENABLE(gen_ens_prod, + [AS_HELP_STRING([--disable-gen_ens_prod], [Disable compilation of gen_ens_prod])], + [case "${enableval}" in + yes | no ) ENABLE_GEN_ENS_PROD="${enableval}" ;; + *) AC_MSG_ERROR(bad value ${enableval} for --disable-gen_ens_prod) ;; + esac], + [ENABLE_GEN_ENS_PROD="yes"] +) + +AM_CONDITIONAL([ENABLE_GEN_ENS_PROD], [test "x$ENABLE_GEN_ENS_PROD" = "xyes"]) + +if test "x$ENABLE_GEN_ENS_PROD" = "xyes"; then + AC_DEFINE([ENABLE_GEN_ENS_PROD], [], ["build gen_ens_prod"]) + AC_MSG_NOTICE([gen_ens_prod will be compiled]) +else + AC_MSG_NOTICE([gen_ens_prod will not be compiled]) +fi + + # gen_vx_mask AC_ARG_ENABLE(gen_vx_mask, @@ -1226,6 +1247,7 @@ AC_CONFIG_FILES([Makefile src/tools/other/Makefile src/tools/other/ascii2nc/Makefile src/tools/other/lidar2nc/Makefile + src/tools/other/gen_ens_prod/Makefile src/tools/other/gen_vx_mask/Makefile src/tools/other/gis_utils/Makefile src/tools/other/ioda2nc/Makefile diff --git a/met/data/config/GenEnsProdConfig_default b/met/data/config/GenEnsProdConfig_default new file mode 100644 index 0000000000..fbe81a00f4 --- /dev/null +++ b/met/data/config/GenEnsProdConfig_default @@ -0,0 +1,141 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Gen-Ens-Prod configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Output model name to be written +// +model = "WRF"; + +// +// Output description to be written +// May be set separately in each "obs.field" entry +// +desc = "NA"; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Verification grid +// May be set separately in each "field" entry +// +regrid = { + to_grid = NONE; + method = NEAREST; + width = 1; + vld_thresh = 0.5; + shape = SQUARE; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// May be set separately in each "field" entry +// +censor_thresh = []; +censor_val = []; +cat_thresh = []; +nc_var_str = ""; + +// +// Ensemble fields to be processed +// +ens = { + ens_thresh = 1.0; + vld_thresh = 1.0; + + field = [ + { + name = "APCP"; + level = "A03"; + cat_thresh = [ >0.0, >=5.0 ]; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Neighborhood ensemble probabilities +// +nbrhd_prob = { + width = [ 5 ]; + shape = CIRCLE; + vld_thresh = 0.0; +} + +// +// NMEP smoothing methods +// +nmep_smooth = { + vld_thresh = 0.0; + shape = CIRCLE; + gaussian_dx = 81.27; + gaussian_radius = 120; + type = [ + { + method = GAUSSIAN; + width = 1; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Climatology data +// +climo_mean = { + + file_name = []; + field = []; + + regrid = { + method = NEAREST; + width = 1; + vld_thresh = 0.5; + shape = SQUARE; + } + + time_interp_method = DW_MEAN; + day_interval = 31; + hour_interval = 6; +} + +climo_stdev = climo_mean; +climo_stdev = { + file_name = []; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Ensemble product output types +// May be set separately in each "ens.field" entry +// +ensemble_flag = { + latlon = TRUE; + mean = TRUE; + stdev = TRUE; + minus = TRUE; + plus = TRUE; + min = TRUE; + max = TRUE; + range = TRUE; + vld_count = TRUE; + frequency = TRUE; + nep = FALSE; + nmep = FALSE; + climo = FALSE; + climo_cdp = FALSE; +} + +//////////////////////////////////////////////////////////////////////////////// + +version = "V10.1.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/met/data/config/Makefile.am b/met/data/config/Makefile.am index c73bb4cadc..5ab35344a9 100644 --- a/met/data/config/Makefile.am +++ b/met/data/config/Makefile.am @@ -23,6 +23,7 @@ config_DATA = \ ConfigMapData \ Ascii2NcConfig_default \ EnsembleStatConfig_default \ + GenEnsProdConfig_default \ GridStatConfig_default \ GridDiagConfig_default \ IODA2NCConfig_default \ diff --git a/met/data/config/STATAnalysisConfig_CBS_Index b/met/data/config/STATAnalysisConfig_CBS_Index index 989b1de4ea..e65ff69483 100644 --- a/met/data/config/STATAnalysisConfig_CBS_Index +++ b/met/data/config/STATAnalysisConfig_CBS_Index @@ -106,6 +106,13 @@ weight = [ 6.4, 6.4, 6.4, 6.4, 6.4, //////////////////////////////////////////////////////////////////////////////// +// +// Array of STAT-Analysis jobs to be performed on the filtered data +// +jobs = []; + +//////////////////////////////////////////////////////////////////////////////// + // // Confidence interval settings // diff --git a/met/data/config/STATAnalysisConfig_GO_Index b/met/data/config/STATAnalysisConfig_GO_Index index a82daa1710..2277304894 100644 --- a/met/data/config/STATAnalysisConfig_GO_Index +++ b/met/data/config/STATAnalysisConfig_GO_Index @@ -114,6 +114,13 @@ weight = [ 4.0, 3.0, 2.0, 1.0, //////////////////////////////////////////////////////////////////////////////// +// +// Array of STAT-Analysis jobs to be performed on the filtered data +// +jobs = []; + +//////////////////////////////////////////////////////////////////////////////// + // // Confidence interval settings // diff --git a/met/data/config/STATAnalysisConfig_default b/met/data/config/STATAnalysisConfig_default index 193faf7d3b..9ac2467ccb 100644 --- a/met/data/config/STATAnalysisConfig_default +++ b/met/data/config/STATAnalysisConfig_default @@ -73,9 +73,7 @@ weight = []; // // Array of STAT-Analysis jobs to be performed on the filtered data // -jobs = [ - "-job filter -dump_row ./filter_job.stat" -]; +jobs = []; //////////////////////////////////////////////////////////////////////////////// diff --git a/met/docs/Flowchart/MET_flowchart.pptx b/met/docs/Flowchart/MET_flowchart.pptx index 575eaada14..5e4f364225 100644 Binary files a/met/docs/Flowchart/MET_flowchart.pptx and b/met/docs/Flowchart/MET_flowchart.pptx differ diff --git a/met/docs/Flowchart/MET_flowchart_v10.1.0.png b/met/docs/Flowchart/MET_flowchart_v10.1.0.png new file mode 100644 index 0000000000..e9f287d189 Binary files /dev/null and b/met/docs/Flowchart/MET_flowchart_v10.1.0.png differ diff --git a/met/docs/Users_Guide/gen-ens-prod.rst b/met/docs/Users_Guide/gen-ens-prod.rst new file mode 100644 index 0000000000..4b1d556d69 --- /dev/null +++ b/met/docs/Users_Guide/gen-ens-prod.rst @@ -0,0 +1,228 @@ +.. _gen-ens-prod: + +Gen-Ens-Prod Tool +================= + +Introduction +____________ + +The Gen-Ens-Prod tool generates simple ensemble products (mean, spread, probability, etc) from gridded ensemble member input files. While it processes model inputs, but it does not compare them to observations or compute statistics. However, the output products can be passed as input to the MET statistics tools for comparison against observations. Climatological mean and standard deviation data may also be provided to define thresholds based on the climatological distribution at each grid point. + +Note that this ensemble product generation step was provided by the Ensemble-Stat tool in earlier versions of MET. The Gen-Ens-Prod tool replaces and extends that functionality. Users are strongly encouraged to migrate ensemble product generation from Ensemble-Stat to Gen-Ens-Prod, as new features will only be added to Gen-Ens-Prod and the existing Ensemble-Stat functionality will be deprecated in a future version. + +Scientific and statistical aspects +__________________________________ + +Ensemble forecasts derived from a set of deterministic ensemble members +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Ensemble forecasts are often created as a set of deterministic forecasts. The ensemble members are rarely used separately. Instead, they can be combined in various ways to produce a forecast. MET can combine the ensemble members into some type of summary forecast according to user specifications. Ensemble means are the most common, and can be paired with the ensemble variance or spread. Maximum, minimum and other summary values are also available, with details in the practical information section. + +The -ctrl command line option specifies an input file for the ensemble control member. The fields specified in the configuration file are read from the control member file. Those fields are included in the computation of the ensemble mean and probabilities but excluded from the ensemble spread. + +The ensemble relative frequency is the simplest method for turning a set of deterministic forecasts into something resembling a probability forecast. MET will create the ensemble relative frequency as the proportion of ensemble members forecasting some event. For example, if 5 out of 10 ensemble members predict measurable precipitation at a grid location, then the ensemble relative frequency of precipitation will be :math:`5/10=0.5`. If the ensemble relative frequency is calibrated (unlikely) then this could be thought of as a probability of precipitation. + +The neighborhood ensemble probability (NEP) and neighborhood maximum ensemble probability (NMEP) methods are described in :ref:`Schwartz and Sobash (2017) `. They are an extension of the ensemble relative frequencies described above. The NEP value is computed by averaging the relative frequency of the event within the neighborhood over all ensemble members. The NMEP value is computed as the fraction of ensemble members for which the event is occurring somewhere within the surrounding neighborhood. The NMEP output is typically smoothed using a Gaussian kernel filter. The neighborhood sizes and smoothing options can be customized in the configuration file. + +The Gen-Ens-Prod tool writes the gridded relative frequencies, NEP, and NMEP fields to a NetCDF output file. Probabilistic verification methods can then be applied to those fields by evaluating them with the Grid-Stat and/or Point-Stat tools. + +Climatology data +~~~~~~~~~~~~~~~~ + +The ensemble relative frequencies derived by Gen-Ens-Prod are computed by applying threshold(s) to the input ensemble member data. Those thresholds can be simple and remain constant over the entire domain (e.g. >0) or can be defined relative to the climatological distribution at each grid point (e.g. >CDP90, for exceeding the 90-th percentile of climatology). When using climatological distribution percentile (CDP) thresholds, the climatological mean and standard deviation must be provided in the configuration file. + +Practical Information +_____________________ + +This section contains information about configuring and running the Gen-Ens-Prod tool. The Gen-Ens-Prod tool writes a NetCDF output file containing the requested ensemble product fields for each input field specified. If provided, the climatology data files must be gridded. All input gridded model and climatology datasets must be on the same grid. However, users may leverage the automated regridding feature in MET if the desired output grid is specified in the configuration file. + +gen_ens_prod usage +~~~~~~~~~~~~~~~~~~~ + +The usage statement for the Ensemble Stat tool is shown below: + +.. code-block:: none + + Usage: gen_ens_prod + -ens file_1 ... file_n | ens_file_list + -out file + -config file + [-ctrl file] + [-log file] + [-v level] + +gen_ens_prod has three required arguments and accepts several optional ones. + +Required arguments gen_ens_prod +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +1. The **-ens file_1 ... file_n** option specifies the ensemble member file names. This argument is not required when ensemble files are specified in the **ens_file_list**, detailed below. + +2. The **ens_file_list** option is an ASCII file containing a list of ensemble member file names. This is not required when a file list is included on the command line, as described above. + +3. The **-out file** option specifies the NetCDF output file name to be written. + +4. The **-config file** option is a **GenEnsProdConfig** file containing the desired configuration settings. + +Optional arguments for gen_ens_prod +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +4. The **-ctrl file** option specifies the input file for the ensemble control member. Data for this member is included in the computation of the ensemble mean, but excluded from the spread. + +5. The **-log** file outputs log messages to the specified file. + +6. The **-v level** option indicates the desired level of verbosity. The value of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity will increase the amount of logging. + +An example of the gen_ens_prod calling sequence is shown below: + +.. code-block:: none + + gen_ens_prod \ + -ens sample_fcst/2009123112/*gep*/d01_2009123112_02400.grib \ + -out out/gen_ens_prod/config/GenEnsProdConfig \ + -config config/GenEnsProdConfig -v 2 + +In this example, the Gen-Ens-Prod tool derives products from the input ensemble members listed on the command line. + +gen_ens_prod configuration file +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The default configuration file for the Gen-Ens-Prod tool named **GenEnsProdConfig_default** can be found in the installed *share/met/config* directory. Another version is located in *scripts/config*. We encourage users to make a copy of these files prior to modifying their contents. The contents of the configuration file are described in the subsections below. + +Note that environment variables may be used when editing configuration files, as described in :numref:`pb2nc configuration file` for the PB2NC tool. + +____________________ + +.. code-block:: none + + model = "WRF"; + desc = "NA"; + regrid = { ... } + censor_thresh = []; + censor_val = []; + nc_var_str = ""; + climo_mean = { ... } // Corresponding to ens.field entries + climo_stdev = { ... } // Corresponding to ens.field entries + rng = { ... } + version = "VN.N"; + +The configuration options listed above are common to many MET tools and are described in :numref:`config_options`. + +_____________________ + +.. code-block:: none + + ens = { + ens_thresh = 1.0; + vld_thresh = 1.0; + field = [ + { + name = "APCP"; + level = "A03"; + cat_thresh = [ >0.0, >=5.0 ]; + } + ]; + } + +The **ens** dictionary defines which ensemble fields should be processed. + +When summarizing the ensemble, compute a ratio of the number of valid ensemble fields to the total number of ensemble members. If this ratio is less than the **ens_thresh**, then quit with an error. This threshold must be between 0 and 1. Setting this threshold to 1 requires that all ensemble members input files exist and all requested data be present. + +When summarizing the ensemble, for each grid point compute a ratio of the number of valid data values to the number of ensemble members. If that ratio is less than **vld_thresh**, write out bad data for that grid point. This threshold must be between 0 and 1. Setting this threshold to 1 requires that each grid point contain valid data for all ensemble members in order to compute ensemble product values for that grid point. + +For each dictionary entry in the **field** array, give the name and vertical or accumulation level, plus one or more categorical thresholds in the **cat_thresh** entry. The formatting for threshold are described in :numref:`config_options`. It is the user's responsibility to know the units for each model variable and choose appropriate threshold values. The thresholds are used to define ensemble relative frequencies. For example, a threshold of >=5 is used to define the proportion of ensemble members predicting precipitation of at least 5mm at each grid point. + +_______________________ + +.. code-block:: none + + nbrhd_prob = { + width = [ 5 ]; + shape = CIRCLE; + vld_thresh = 0.0; + } + +The **nbrhd_prob** dictionary defines the neighborhoods used to compute NEP and NMEP output. + +The neighborhood **shape** is a **SQUARE** or **CIRCLE** centered on the current point, and the **width** array specifies the width of the square or diameter of the circle as an odd integer. The **vld_thresh** entry is a number between 0 and 1 specifying the required ratio of valid data in the neighborhood for an output value to be computed. + +If **ensemble_flag.nep** is set to TRUE, NEP output is created for each combination of the categorical threshold (**cat_thresh**) and neighborhood width specified. + +_____________________ + +.. code-block:: none + + nmep_smooth = { + vld_thresh = 0.0; + shape = CIRCLE; + gaussian_dx = 81.27; + gaussian_radius = 120; + type = [ + { + method = GAUSSIAN; + width = 1; + } + ]; + } + +Similar to the **interp** dictionary, the **nmep_smooth** dictionary includes a **type** array of dictionaries to define one or more methods for smoothing the NMEP data. Setting the interpolation method to nearest neighbor (**NEAREST**) effectively disables this smoothing step. + +If **ensemble_flag.nmep** is set to TRUE, NMEP output is created for each combination of the categorical threshold (**cat_thresh**), neighborhood width (**nbrhd_prob.width**), and smoothing method(**nmep_smooth.type**) specified. + +_____________________ + +.. code-block:: none + + ensemble_flag = { + latlon = TRUE; + mean = TRUE; + stdev = TRUE; + minus = TRUE; + plus = TRUE; + min = TRUE; + max = TRUE; + range = TRUE; + vld_count = TRUE; + frequency = TRUE; + nep = FALSE; + nmep = FALSE; + climo = FALSE; + climo_cdp = FALSE; + } + +The **ensemble_flag** specifies which derived ensemble fields should be calculated and output. Setting the flag to TRUE produces output of the specified field, while FALSE produces no output for that field type. The flags correspond to the following output line types: + +1. Grid Latitude and Longitude Fields + +2. Ensemble Mean Field + +3. Ensemble Standard Deviation Field + +4. Ensemble Mean - One Standard Deviation Field + +5. Ensemble Mean + One Standard Deviation Field + +6. Ensemble Minimum Field + +7. Ensemble Maximum Field + +8. Ensemble Range Field + +9. Ensemble Valid Data Count + +10. Ensemble Relative Frequency (i.e. uncalibrate probability forecast) for each categorical threshold (**cat_thresh**) specified + +11. Neighborhood Ensemble Probability for each categorical threshold (**cat_thresh**) and neighborhood width (**nbrhd_prob.width**) specified + +12. Neighborhood Maximum Ensemble Probability for each categorical threshold (**cat_thresh**), neighborhood width (**nbrhd_prob.width**), and smoothing method (**nmep_smooth.type**) specified + +13. Climatology mean (**climo_mean**) and standard deviation (**climo_stdev**) data regridded to the model domain + +14. Climatological Distribution Percentile field for each CDP threshold specified + +gen_ens_prod output +~~~~~~~~~~~~~~~~~~~~ + +The Gen-Ens-Prod tools writes a gridded NetCDF output file whose file name is specified using the -out command line option. The contents of that file depend on the contents of the **ens.field** array, the **ensemble_flag** options selected, and the presence of climatology data. The NetCDF variable names are self-describing and include the name/level of the field being processed, the type of ensemble product, and any relevant threshold information. If **nc_var_str** is defined for an **ens.field** array entry, that string is included in the corresponding NetCDF output variable names. + +The Gen-Ens-Prod NetCDF output can be passed as input to the MET statistics tools, like Point-Stat and Grid-Stat, for futher processing and comparison against observations. diff --git a/met/docs/Users_Guide/grid-stat.rst b/met/docs/Users_Guide/grid-stat.rst index ef1a465763..b96ed7a724 100644 --- a/met/docs/Users_Guide/grid-stat.rst +++ b/met/docs/Users_Guide/grid-stat.rst @@ -431,7 +431,7 @@ _____________________ } -The **nc_pairs_flag** entry may either be set to a boolean value or a dictionary specifying which fields should be written. Setting it to TRUE indicates the output NetCDF matched pairs file should be created with all available output fields, while setting all to FALSE disables its creation. This is done regardless of if **output_flag** dictionary indicates any statistics should be computed. The **latlon, raw**, and **diff** entries control the creation of output variables for the latitude and longitude, the raw forecast and observed fields, and the forecast minus observation difference fields. The **climo, weight**, and **nbrhd** entries control the creation of output variables for the climatological mean and standard deviation fields, the grid area weights applied, and the fractional coverage fields computed for neighborhood verification methods. Setting these entries to TRUE indicates that they should be written, while setting them to FALSE disables their creation. +The **nc_pairs_flag** entry may either be set to a boolean value or a dictionary specifying which fields should be written. Setting it to TRUE indicates the output NetCDF matched pairs file should be created with all available output fields, while setting all to FALSE disables its creation. This is done regardless of if **output_flag** dictionary indicates any statistics should be computed. The **latlon, raw**, and **diff** entries control the creation of output variables for the latitude and longitude, the forecast and observed fields after they have been modified by any user-defined regridding, censoring, and conversion, and the forecast minus observation difference fields, respectively. The **climo, weight**, and **nbrhd** entries control the creation of output variables for the climatological mean and standard deviation fields, the grid area weights applied, and the fractional coverage fields computed for neighborhood verification methods. Setting these entries to TRUE indicates that they should be written, while setting them to FALSE disables their creation. Setting the **climo_cdp** entry to TRUE enables the creation of an output variable for each climatological distribution percentile (CDP) threshold requested in the configuration file. Note that enabling **nbrhd** output may lead to very large output files. The **gradient** entry controls the creation of output variables for the FCST and OBS gradients in the grid-x and grid-y directions. The **distance_map** entry controls the creation of output variables for the FCST and OBS distance maps for each categorical threshold. The **apply_mask** entry controls whether to create the FCST, OBS, and DIFF output variables for all defined masking regions. Setting this to TRUE will create the FCST, OBS, and DIFF output variables for all defined masking regions. Setting this to FALSE will create the FCST, OBS, and DIFF output variables for only the FULL verification domain. diff --git a/met/docs/Users_Guide/index.rst b/met/docs/Users_Guide/index.rst index 213fcae61d..5ee52862e0 100644 --- a/met/docs/Users_Guide/index.rst +++ b/met/docs/Users_Guide/index.rst @@ -50,6 +50,7 @@ The National Center for Atmospheric Research (NCAR) is sponsored by NSF. The DTC config_options_tc reformat_point reformat_grid + gen-ens-prod masking point-stat grid-stat diff --git a/met/docs/Users_Guide/mode.rst b/met/docs/Users_Guide/mode.rst index defbb38603..ac42611764 100644 --- a/met/docs/Users_Guide/mode.rst +++ b/met/docs/Users_Guide/mode.rst @@ -479,6 +479,8 @@ _____________________ Each component of the pairs information in the NetCDF file can be turned on or off. The old syntax is still supported: **TRUE** means accept the defaults, **FALSE** means no NetCDF output is generated. NetCDF output can also be turned off by setting all the individual dictionary flags to false. +The nc_pairs_flag is described in :numref:`grid_stat-configuration-file` + _____________________ diff --git a/met/export.mk b/met/export.mk index 93bd734461..eeb741edab 100644 --- a/met/export.mk +++ b/met/export.mk @@ -41,6 +41,7 @@ export HDF5_BASE_DIR export ENABLE_ASCII2NC export ENABLE_ENSEMBLE_STAT +export ENABLE_GEN_ENS_PROD export ENABLE_GEN_VX_MASK export ENABLE_GRID_STAT export ENABLE_IODA2NC diff --git a/met/out/gen_ens_prod/.gitignore b/met/out/gen_ens_prod/.gitignore new file mode 100644 index 0000000000..d6b7ef32c8 --- /dev/null +++ b/met/out/gen_ens_prod/.gitignore @@ -0,0 +1,2 @@ +* +!.gitignore diff --git a/met/scripts/.gitignore b/met/scripts/.gitignore index 2c32bbda55..1263b3675e 100644 --- a/met/scripts/.gitignore +++ b/met/scripts/.gitignore @@ -1,6 +1,7 @@ ascii2nc ensemble_stat gen_vx_mask +gen_ens_prod grib2 grid_stat gis_utils diff --git a/met/scripts/Makefile b/met/scripts/Makefile index e9bd0c73ff..2081261e28 100755 --- a/met/scripts/Makefile +++ b/met/scripts/Makefile @@ -51,6 +51,7 @@ all: ${TESTS} include $(MK_DIR)/ascii2nc.mk include $(MK_DIR)/ensemble_stat.mk +include $(MK_DIR)/gen_ens_prod.mk include $(MK_DIR)/gen_vx_mask.mk include $(MK_DIR)/grib2.mk include $(MK_DIR)/grid_stat.mk @@ -94,6 +95,7 @@ clean: rm -f ${TESTS} rm -f ${TEST_OUT_DIR}/ascii2nc/* rm -f ${TEST_OUT_DIR}/ensemble_stat/* + rm -f ${TEST_OUT_DIR}/gen_ens_prod/* rm -f ${TEST_OUT_DIR}/gen_vx_mask/* rm -f ${TEST_OUT_DIR}/grid_stat/* rm -f ${TEST_OUT_DIR}/madis2nc/* diff --git a/met/scripts/config/GenEnsProdConfig b/met/scripts/config/GenEnsProdConfig new file mode 100644 index 0000000000..a0a52794c9 --- /dev/null +++ b/met/scripts/config/GenEnsProdConfig @@ -0,0 +1,163 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Gen-Ens-Prod configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Output model name to be written +// +model = "WRF"; + +// +// Output description to be written +// May be set separately in each "obs.field" entry +// +desc = "NA"; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Verification grid +// May be set separately in each "field" entry +// +regrid = { + to_grid = NONE; + method = NEAREST; + width = 1; + vld_thresh = 0.5; + shape = SQUARE; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// May be set separately in each "field" entry +// +censor_thresh = []; +censor_val = []; +cat_thresh = []; +nc_var_str = ""; + +// +// Ensemble fields to be processed +// +ens = { + ens_thresh = 1.0; + vld_thresh = 1.0; + + field = [ + { + name = "APCP"; + level = [ "A24" ]; + cat_thresh = [ >0.0, >=10.0 ]; + ensemble_flag = TRUE; + }, + { + name = "REFC"; + level = [ "L0" ]; + cat_thresh = [ >=35.0 ]; + GRIB1_ptv = 129; + }, + { + name = "UGRD"; + level = [ "Z10" ]; + cat_thresh = [ >=5.0 ]; + }, + { + name = "VGRD"; + level = [ "Z10" ]; + cat_thresh = [ >=5.0 ]; + }, + { + name = "WIND"; + level = [ "Z10" ]; + cat_thresh = [ >=5.0 ]; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Neighborhood ensemble probabilities +// +nbrhd_prob = { + width = [ 5 ]; + shape = CIRCLE; + vld_thresh = 0.0; +} + +// +// NMEP smoothing methods +// +nmep_smooth = { + vld_thresh = 0.0; + shape = CIRCLE; + gaussian_dx = 81.27; + gaussian_radius = 120; + type = [ + { + method = GAUSSIAN; + width = 1; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Climatology data +// +climo_mean = { + + file_name = []; + field = []; + + regrid = { + method = NEAREST; + width = 1; + vld_thresh = 0.5; + shape = SQUARE; + } + + time_interp_method = DW_MEAN; + day_interval = 31; + hour_interval = 6; +} + +climo_stdev = climo_mean; +climo_stdev = { + file_name = []; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Ensemble product output types +// May be set separately in each "ens.field" entry +// +ensemble_flag = { + latlon = TRUE; + mean = TRUE; + stdev = TRUE; + minus = FALSE; + plus = FALSE; + min = FALSE; + max = FALSE; + range = FALSE; + vld_count = FALSE; + frequency = TRUE; + nep = FALSE; + nmep = FALSE; + climo = FALSE; + climo_cdp = FALSE; +} + +//////////////////////////////////////////////////////////////////////////////// + +version = "V10.1.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/met/scripts/examples/test_gen_ens_prod.sh b/met/scripts/examples/test_gen_ens_prod.sh new file mode 100755 index 0000000000..c30d3cd273 --- /dev/null +++ b/met/scripts/examples/test_gen_ens_prod.sh @@ -0,0 +1,8 @@ +#/bin/sh + +echo +echo "*** Running Gen-Ens-Prod using GRIB forecasts ***" +gen_ens_prod \ + -ens ../data/sample_fcst/2009123112/*gep*/d01_2009123112_02400.grib \ + -config config/GenEnsProdConfig \ + -out ${TEST_OUT_DIR}/gen_ens_prod/gen_ens_prod_20100101_120000V_ens.nc -v 2 diff --git a/met/scripts/mk/gen_ens_prod.mk b/met/scripts/mk/gen_ens_prod.mk new file mode 100644 index 0000000000..7c9f0f448f --- /dev/null +++ b/met/scripts/mk/gen_ens_prod.mk @@ -0,0 +1,33 @@ + + + +######################################################################## + + +GEN_ENS_PROD_EXEC = ${OTHER_DIR}/gen_ens_prod/gen_ens_prod + + +######################################################################## + + + ## + ## gen_ens_prod + ## + ## prerequisites: + ## + + +gen_ens_prod: ${GEN_ENS_PROD_EXEC} + @ echo + @ echo "*** Running Gen-Ens-Prod using GRIB forecasts ***" + ${GEN_ENS_PROD_EXEC} \ + -ens ../data/sample_fcst/2009123112/*gep*/d01_2009123112_02400.grib \ + -config config/GenEnsProdConfig \ + -out ${TEST_OUT_DIR}/gen_ens_prod/gen_ens_prod_20100101_120000V_ens.nc -v 2 + @ + @ touch gen_ens_prod + + +######################################################################## + + diff --git a/met/src/basic/vx_math/nint.cc b/met/src/basic/vx_math/nint.cc index f027903500..d013916ca3 100644 --- a/met/src/basic/vx_math/nint.cc +++ b/met/src/basic/vx_math/nint.cc @@ -46,3 +46,13 @@ return ( a ); //////////////////////////////////////////////////////////////////////// +int positive_modulo(int i, int n) + +{ + +return ( i % n + n ) % n; + +} + + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_math/nint.h b/met/src/basic/vx_math/nint.h index 3c466858d6..712823ecba 100644 --- a/met/src/basic/vx_math/nint.h +++ b/met/src/basic/vx_math/nint.h @@ -11,8 +11,8 @@ //////////////////////////////////////////////////////////////////////// -#ifndef __PS_DITHER_NINT_H__ -#define __PS_DITHER_NINT_H__ +#ifndef __VX_MATH_NINT_H__ +#define __VX_MATH_NINT_H__ //////////////////////////////////////////////////////////////////////// @@ -24,7 +24,13 @@ extern int nint(double); //////////////////////////////////////////////////////////////////////// -#endif // __PS_DITHER_NINT_H__ +extern int positive_modulo(int, int); + + +//////////////////////////////////////////////////////////////////////// + + +#endif // __VX_MATH_NINT_H__ //////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/CircularTemplate.cc b/met/src/basic/vx_util/CircularTemplate.cc index a49f40c62f..2b2bac65d6 100644 --- a/met/src/basic/vx_util/CircularTemplate.cc +++ b/met/src/basic/vx_util/CircularTemplate.cc @@ -1,47 +1,25 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 -// ** University Corporation for Atmospheric Research (UCAR) -// ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -// RCS info: dixon $ -// $Locker: $ -// $Date: 2016/03/03 18:19:27 $ -// $Id: CircularTemplate.cc,v 1.6 2016/03/03 18:19:27 dixon Exp $ -// $Revision: 1.6 $ -// $State: Exp $ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ -/********************************************************************* - * CircularTemplate.cc: class implementing a circular template to be - * applied on gridded data. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - *********************************************************************/ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: CircularTemplate.cc +// +// Description: +// Class implementing a Circular template to be +// applied on gridded data. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// 001 09-07-21 Halley Gotway Add wrap_lon. +// +/////////////////////////////////////////////////////////////////////////////// #include @@ -56,126 +34,92 @@ using namespace std; +/////////////////////////////////////////////////////////////////////////////// -/********************************************************************** - * Constructor - */ +CircularTemplate::CircularTemplate(const int width, bool wrap_lon) : + GridTemplate(), _width(width) { -CircularTemplate::CircularTemplate(const int width) : - GridTemplate(), - _width(width) -{ + _wrapLon = wrap_lon; - //width of 2 is not supported - if (width == 2){ - mlog << Error << "\nCircularTemplate::CircularTemplate() -> " - << "unsupported width of " << width - << " for circles.\n\n"; + // width of 2 is not supported + if (width == 2) { + mlog << Error << "\nCircularTemplate::CircularTemplate() -> " + << "unsupported width of " << width << " for circles.\n\n"; exit(1); } - - bool evenWidth = ((width % 2) == 0); - // if the width is even, that means we are dealing with a point interpolation - // because grid interpolation has to be odd. - // for an ODD WIDTH the reference point is the same as the center point and is the nearest grid point - - // for an EVEN WIDTH, we move the "reference" point, to the lower left grid point, - // this means offsets are stored relative to the lower left corner of the true center. - // but we find distances based on the the true center location when determining if an - // offset is within the circle. - - double radius = (width-1)/2.0; + bool evenWidth = ((width % 2) == 0); + + // if the width is even, that means we are dealing with a point interpolation + // because grid interpolation has to be odd. + // for an ODD WIDTH the reference point is the same as the center point and is the nearest grid point + // for an EVEN WIDTH, we move the "reference" point, to the lower left grid point, + // this means offsets are stored relative to the lower left corner of the true center. + // but we find distances based on the the true center location when determining if an + // offset is within the circle. + + double radius = (width-1)/2.0; - // Create the offsets list. - - // If the radius is small, then just make the - // template cover the current grid square. - // radius < 1 no longer supported. - /* - if (radius < 1.0) - { - _addOffset(0, 0); - return; - } - */ + // Create the offsets list. - // need to increase the area we look at if the width is even, because - // some valid offset points will actually be farther from the reference point - // than the radius, because the reference point is offset from the true - // center of the circle. - int maxOffset = static_cast(floor(radius)); - if (evenWidth){ - maxOffset++; - } + // Need to increase the area we look at if the width is even, because + // some valid offset points will actually be farther from the reference point + // than the radius, because the reference point is offset from the true + // center of the circle. + + int maxOffset = static_cast(floor(radius)); + if(evenWidth) maxOffset++; - int minOffset = static_cast(floor(-1 * radius)); + int minOffset = static_cast(floor(-1 * radius)); - for (int y = minOffset; y <= maxOffset; y++) - { - for (int x = minOffset; x <= maxOffset; x++) - { - double double_x = (double)x; - double double_y = (double)y; - - if (evenWidth){ - // if width is even, the reference point is actually shifted 1/2 a grid spacing down and to the left, - // from the true center of the circle. - // - // so when we calculate distance, we need to subtract .5 so that the distance reflects the distance from the center - // of the circle, instead of the distance from the reference. - // - // for example - a circle with width == 4. The reference point is the lower left corner of the center square. - // the point directly below that is at (0,-1), but it's actually (-.5, -1.5) from the center of the circle. - // - // another example - same circle. The point directly to the right of the reference point is (1,0), but it's - // actually (.5,-.5) from the center. + for(int y = minOffset; y <= maxOffset; y++) { + for(int x = minOffset; x <= maxOffset; x++) { + double double_x = (double)x; + double double_y = (double)y; + + if(evenWidth) { + // if width is even, the reference point is actually shifted 1/2 a grid spacing down and to the left, + // from the true center of the circle. + // + // so when we calculate distance, we need to subtract .5 so that the distance reflects the distance from the center + // of the circle, instead of the distance from the reference. + // + // for example - a circle with width == 4. The reference point is the lower left corner of the center square. + // the point directly below that is at (0,-1), but it's actually (-.5, -1.5) from the center of the circle. + // + // another example - same circle. The point directly to the right of the reference point is (1,0), but it's + // actually (.5,-.5) from the center. - double_x -= 0.5; - double_y -= 0.5; - } - double distance= sqrt((double_x * double_x) + (double_y * double_y)); - - if (distance <= radius) - { - _addOffset(x, y); - } - - } /* endfor - x */ - - } /* endfor - y */ - - _setEdgeOffsets(); - -} + double_x -= 0.5; + double_y -= 0.5; + } + double distance= sqrt((double_x * double_x) + (double_y * double_y)); + + if(distance <= radius) _addOffset(x, y); + } // end for x + } // end for y -/********************************************************************** - * Destructor - */ + _setEdgeOffsets(); -CircularTemplate::~CircularTemplate(void) -{ - // Do nothing } - -/********************************************************************** - * printOffsetList() - Print the offset list to the given stream. This - * is used for debugging. - */ +/////////////////////////////////////////////////////////////////////////////// -void CircularTemplate::printOffsetList(FILE *stream) -{ - fprintf(stream, "\n\n"); - fprintf(stream, "Circular template with width %d grid spaces:\n", - _width); - - GridTemplate::printOffsetList(stream); +CircularTemplate::~CircularTemplate(void) { + // Do nothing } +/////////////////////////////////////////////////////////////////////////////// + +void CircularTemplate::printOffsetList(FILE *stream) { + fprintf(stream, "\n\n"); + fprintf(stream, "Circular template:"); + fprintf(stream, " width = %d\n", _width); + fprintf(stream, " wrap_lon = %d\n", _wrapLon); + fprintf(stream, " grid points:\n"); + GridTemplate::printOffsetList(stream); +} -/********************************************************************** - * Private Member Functions * - **********************************************************************/ +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/CircularTemplate.h b/met/src/basic/vx_util/CircularTemplate.h index fc92c7db8b..835aa99ed5 100644 --- a/met/src/basic/vx_util/CircularTemplate.h +++ b/met/src/basic/vx_util/CircularTemplate.h @@ -1,53 +1,30 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 -// ** University Corporation for Atmospheric Research (UCAR) -// ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -/* RCS info - * $Author: dixon $ - * $Locker: $ - * $Date: 2016/03/03 19:21:31 $ - * $Id: CircularTemplate.hh,v 1.5 2016/03/03 19:21:31 dixon Exp $ - * $Revision: 1.5 $ - * $State: Exp $ - */ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ - -/************************************************************************ - * CircularTemplate.hh: class implementing a circular template to be - * applied on gridded data. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - ************************************************************************/ - -#ifndef CircularTemplate_HH -#define CircularTemplate_HH +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: CircularTemplate.h +// +// Description: +// Class implementing a Circular template to be +// applied on gridded data. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// 001 09-07-21 Halley Gotway Add wrap_lon. +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __CIRCULAR_TEMPLATE_H__ +#define __CIRCULAR_TEMPLATE_H__ + +/////////////////////////////////////////////////////////////////////////////// #include @@ -57,46 +34,39 @@ #include #include +/////////////////////////////////////////////////////////////////////////////// -class CircularTemplate : public GridTemplate -{ - public: +class CircularTemplate : public GridTemplate { - // Constructor + public: - CircularTemplate(int width = 0); - - // Destructor + CircularTemplate(int width, bool wrap_lon); + virtual ~CircularTemplate(void); - virtual ~CircularTemplate(void); + void printOffsetList(FILE *stream); - // Print the offset list to the given stream. This is used for debugging. + // Access methods - void printOffsetList(FILE *stream); - - // Access methods + int getWidth(void) const { + return _width; + } - int getWidth(void) const - { - return _width; - } - - - virtual const char* getClassName(void) const{ - return _className(); - } + virtual const char* getClassName(void) const { + return _className(); + } - private: + private: - int _width; - - // Return the class name for error messages. - static const char* _className(void) - { - return("CircleTemplate"); - } + int _width; + // Return the class name for error messages. + static const char* _className(void) { + return("CircleTemplate"); + } }; +/////////////////////////////////////////////////////////////////////////////// + +#endif // __CIRCULAR_TEMPLATE_H__ -#endif +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/GridOffset.cc b/met/src/basic/vx_util/GridOffset.cc index 4ae65fdf57..f21fd39636 100644 --- a/met/src/basic/vx_util/GridOffset.cc +++ b/met/src/basic/vx_util/GridOffset.cc @@ -1,48 +1,24 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -// RCS info -// $Author: dixon $ -// $Locker: $ -// $Date: 2016/03/03 18:19:27 $ -// $Id: GridOffset.cc,v 1.6 2016/03/03 18:19:27 dixon Exp $ -// $Revision: 1.6 $ -// $State: Exp $ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ -/********************************************************************* - * GridOffset.cc: class implementing grid points as described as offsets - * from an origin. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - *********************************************************************/ +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: GridOffset.cc +// +// Description: +// Class implementing grid index points as described as +// offsets from an origin. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// +/////////////////////////////////////////////////////////////////////////////// #include @@ -53,9 +29,7 @@ using namespace std; -/********************************************************************** - * Constructors - */ +/////////////////////////////////////////////////////////////////////////////// GridOffset::GridOffset(int cur_x_offset, int cur_y_offset) { @@ -63,6 +37,7 @@ GridOffset::GridOffset(int cur_x_offset, int cur_y_offset) this->y_offset = cur_y_offset; } +/////////////////////////////////////////////////////////////////////////////// GridOffset::GridOffset(const GridOffset& rhs) { @@ -70,6 +45,7 @@ GridOffset::GridOffset(const GridOffset& rhs) y_offset = rhs.y_offset; } +/////////////////////////////////////////////////////////////////////////////// GridOffset::GridOffset(const GridOffset* rhs) { @@ -77,16 +53,10 @@ GridOffset::GridOffset(const GridOffset* rhs) y_offset = rhs->y_offset; } - -/********************************************************************** - * Destructor - */ +/////////////////////////////////////////////////////////////////////////////// GridOffset::~GridOffset(void) { } - -/********************************************************************** - * Private Member Functions * - **********************************************************************/ +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/GridOffset.h b/met/src/basic/vx_util/GridOffset.h index 8e4826011b..dd04a36bcc 100644 --- a/met/src/basic/vx_util/GridOffset.h +++ b/met/src/basic/vx_util/GridOffset.h @@ -1,56 +1,34 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 -// ** University Corporation for Atmospheric Research (UCAR) -// ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -/* RCS info - * $Author: dixon $ - * $Locker: $ - * $Date: 2016/03/03 19:21:31 $ - * $Id: GridOffset.hh,v 1.5 2016/03/03 19:21:31 dixon Exp $ - * $Revision: 1.5 $ - * $State: Exp $ - */ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ - -/************************************************************************ - * GridOffset.hh: class implementing grid index points as described as - * offsets from an origin. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - ************************************************************************/ - -#ifndef GridOffset_HH -#define GridOffset_HH +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: GridOffset.h +// +// Description: +// Class implementing grid index points as described as +// offsets from an origin. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __GRID_OFFSET_H__ +#define __GRID_OFFSET_H__ + +/////////////////////////////////////////////////////////////////////////////// #include +/////////////////////////////////////////////////////////////////////////////// + using namespace std; class GridOffset @@ -74,5 +52,8 @@ class GridOffset }; +/////////////////////////////////////////////////////////////////////////////// + +#endif // __GRID_OFFSET_H__ -#endif +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/GridPoint.cc b/met/src/basic/vx_util/GridPoint.cc index 82581f2aeb..bbba9ee87a 100644 --- a/met/src/basic/vx_util/GridPoint.cc +++ b/met/src/basic/vx_util/GridPoint.cc @@ -1,47 +1,23 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -// RCS info -// $Author: dixon $ -// $Locker: $ -// $Date: 2016/03/03 18:19:27 $ -// $Id: GridPoint.cc,v 1.7 2016/03/03 18:19:27 dixon Exp $ -// $Revision: 1.7 $ -// $State: Exp $ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ -/********************************************************************* - * GridPoint.cc: class implementing grid index points. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - *********************************************************************/ +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: GridPoint.cc +// +// Description: +// Class implementing grid index points. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// +/////////////////////////////////////////////////////////////////////////////// #include @@ -52,42 +28,39 @@ #include using namespace std; -/********************************************************************** - * Constructors - */ +/////////////////////////////////////////////////////////////////////////////// GridPoint::GridPoint(int cur_x, int cur_y) { setPoint(cur_x, cur_y); } +/////////////////////////////////////////////////////////////////////////////// GridPoint::GridPoint(GridPoint *point) { setPoint(point); } +/////////////////////////////////////////////////////////////////////////////// GridPoint::GridPoint(GridPoint *point, GridOffset *offset) { setPoint(point, offset); } - -/********************************************************************** - * Destructor - */ +/////////////////////////////////////////////////////////////////////////////// GridPoint::~GridPoint(void) { } - -/********************************************************************** - * rotate() - Rotate the point about the origin by the given angle. - * The angle value must be given in degrees. - */ - +/////////////////////////////////////////////////////////////////////////////// +// +// Rotate the point about the origin by the given angle. +// The angle value must be given in degrees. +// +/////////////////////////////////////////////////////////////////////////////// void GridPoint::rotate(const double angle) { @@ -103,7 +76,4 @@ void GridPoint::rotate(const double angle) y = (int)(new_y + 0.5); } - -/********************************************************************** - * Private Member Functions * - **********************************************************************/ +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/GridPoint.h b/met/src/basic/vx_util/GridPoint.h index 84876e37fe..9e66aceb0e 100644 --- a/met/src/basic/vx_util/GridPoint.h +++ b/met/src/basic/vx_util/GridPoint.h @@ -1,57 +1,35 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -/* RCS info - * $Author: dixon $ - * $Locker: $ - * $Date: 2016/03/03 19:21:31 $ - * $Id: GridPoint.hh,v 1.10 2016/03/03 19:21:31 dixon Exp $ - * $Revision: 1.10 $ - * $State: Exp $ - */ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ - -/************************************************************************ - * GridPoint.hh: class implementing grid index points. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - ************************************************************************/ - -#ifndef GridPoint_HH -#define GridPoint_HH +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: GridPoint.h +// +// Description: +// Class implementing grid index points. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __GRID_POINT_H__ +#define __GRID_POINT_H__ + +/////////////////////////////////////////////////////////////////////////////// #include #include "GridOffset.h" +/////////////////////////////////////////////////////////////////////////////// + class GridPoint { public: @@ -161,5 +139,8 @@ class GridPoint }; +/////////////////////////////////////////////////////////////////////////////// + +#endif // __GRID_POINT_H__ -#endif +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/GridTemplate.cc b/met/src/basic/vx_util/GridTemplate.cc index f3b83b3252..486ff28b44 100644 --- a/met/src/basic/vx_util/GridTemplate.cc +++ b/met/src/basic/vx_util/GridTemplate.cc @@ -1,48 +1,25 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -// RCS info -// $Author: dixon $ -// $Locker: $ -// $Date: 2016/03/03 18:19:27 $ -// $Id: GridTemplate.cc,v 1.14 2016/03/03 18:19:27 dixon Exp $ -// $Revision: 1.14 $ -// $State: Exp $ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ -/********************************************************************* - * GridTemplate.cc: class implementing a template to be applied to - * gridded data. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - *********************************************************************/ +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: GridTemplate.cc +// +// Description: +// Class implementing a template to be applied to +// gridded data. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// 001 09-07-21 Halley Gotway Add wrap_lon. +// +/////////////////////////////////////////////////////////////////////////////// #include #include @@ -51,6 +28,7 @@ #include #include "vx_log.h" +#include "nint.h" #include "GridTemplate.h" #include "GridOffset.h" @@ -60,18 +38,19 @@ using namespace std; -/********************************************************************** - * Constructors - */ +/////////////////////////////////////////////////////////////////////////////// -GridTemplate::GridTemplate(void) -{ - // Do nothing +GridTemplate::GridTemplate(void) : + _wrapLon(false) { + // Do nothing } -GridTemplate::GridTemplate(const GridTemplate& rhs) -{ - vector< GridOffset* >::const_iterator offset_iter; +/////////////////////////////////////////////////////////////////////////////// + +GridTemplate::GridTemplate(const GridTemplate& rhs) { + _wrapLon = rhs._wrapLon; + + vector::const_iterator offset_iter; for (offset_iter = rhs._offsetList.begin(); offset_iter != rhs._offsetList.end(); ++offset_iter) @@ -83,536 +62,558 @@ GridTemplate::GridTemplate(const GridTemplate& rhs) _pointInGridReturn = rhs._pointInGridReturn; } -/********************************************************************** - * Destructor - */ +/////////////////////////////////////////////////////////////////////////////// -GridTemplate::~GridTemplate(void) -{ - // Reclaim the space for the offset list +GridTemplate::~GridTemplate(void) { - vector< GridOffset* >::iterator list_iter; - for (list_iter = _offsetList.begin(); list_iter != _offsetList.end(); + // Reclaim the space for the offset list + vector::iterator list_iter; + for(list_iter = _offsetList.begin(); list_iter != _offsetList.end(); ++list_iter) - delete *list_iter; - - _offsetList.erase(_offsetList.begin(), _offsetList.end()); + delete *list_iter; + _offsetList.erase(_offsetList.begin(), _offsetList.end()); } -/********************************************************************** - * getFirstInGrid() - Get the first template grid point within the given - * grid. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the first template grid point within the given grid. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// GridPoint *GridTemplate::getFirstInGrid( - const int &base_x, const int &base_y, - const int &nx, const int &ny) const -{ - // Set the grid information + const int &base_x, const int &base_y, + const int &nx, const int &ny) const { - setGrid(base_x, base_y, nx, ny); + // Set the grid information, applying the global wrap logic + setGrid((_wrapLon ? base_x % nx : base_x), base_y, nx, ny); - // Send back the first point + // Send back the first point - return getNextInGrid(); + return getNextInGrid(); } -/********************************************************************** - * getNextInGrid() - Get the next template grid point within the grid. - * Returns NULL when there are no more points in the - * grid. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ - -GridPoint *GridTemplate::getNextInGrid(void) const -{ - while (_pointInGridIterator != _offsetList.end()) - { - GridOffset *offset = *_pointInGridIterator; - - _pointInGridIterator++; - - _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; - _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; - - if (_pointInGridReturn.x >= 0 && - _pointInGridReturn.x < _pointInGridNumX && - _pointInGridReturn.y >= 0 && - _pointInGridReturn.y < _pointInGridNumY) - { - return &_pointInGridReturn; - } - - } - - return (GridPoint *)NULL; +/////////////////////////////////////////////////////////////////////////////// +// +// Get the next template grid point within the grid. +// Returns NULL when there are no more points in the grid. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// + +GridPoint *GridTemplate::getNextInGrid(void) const { + + while (_pointInGridIterator != _offsetList.end()) { + GridOffset *offset = *_pointInGridIterator; + + _pointInGridIterator++; + + _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; + _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; + + // Apply global wrap logic + _pointInGridReturn.x = (_wrapLon ? + positive_modulo(_pointInGridReturn.x, _pointInGridNumX) : + _pointInGridReturn.x); + + if(_pointInGridReturn.x >= 0 && + _pointInGridReturn.x < _pointInGridNumX && + _pointInGridReturn.y >= 0 && + _pointInGridReturn.y < _pointInGridNumY) { + return &_pointInGridReturn; + } + } // end while + + return (GridPoint *)NULL; } -/********************************************************************** - * getFirst() - Get the first template grid point without checking - * the grid bounds. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the first template grid point without checking the grid +// bounds. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// GridPoint *GridTemplate::getFirst(const int &base_x, const int &base_y, - const int &nx, const int &ny) const -{ - // Set the grid information - setGrid(base_x, base_y, nx, ny); + const int &nx, const int &ny) const { + + // Set the grid information, applying the global wrap logic + setGrid((_wrapLon ? positive_modulo(base_x, nx) : base_x), base_y, nx, ny); - // Send back the first point - return getNext(); + // Send back the first point + return getNext(); } -/********************************************************************** - * getNext() - Get the next template grid point without checking the - * grid bounds. Returns NULL when there are no more points. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the next template grid point without checking the grid bounds. +// Returns NULL when there are no more points. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// -GridPoint *GridTemplate::getNext(void) const -{ - GridPoint *next_point = (GridPoint *)NULL; - if (_pointInGridIterator != _offsetList.end()) - { - GridOffset *offset = *_pointInGridIterator; +GridPoint *GridTemplate::getNext(void) const { - _pointInGridIterator++; + GridPoint *next_point = (GridPoint *)NULL; + if(_pointInGridIterator != _offsetList.end()) { + GridOffset *offset = *_pointInGridIterator; - _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; - _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; + _pointInGridIterator++; - next_point = &_pointInGridReturn; + _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; + _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; - } + // Apply global wrap logic + _pointInGridReturn.x = (_wrapLon ? + positive_modulo(_pointInGridReturn.x, _pointInGridNumX) : + _pointInGridReturn.x); + + next_point = &_pointInGridReturn; + } - return next_point; + return next_point; } -/********************************************************************** - * getFirstInLftEdge() - Get the first template grid point in the left - * column. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the first template grid point in the left column. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// -GridPoint *GridTemplate::getFirstInLftEdge(void) const -{ - // Reset the iterator and send back the first point +GridPoint *GridTemplate::getFirstInLftEdge(void) const { - _pointInLftEdgeIterator = _offsetLftEdge.begin(); + // Reset the iterator and send back the first point + _pointInLftEdgeIterator = _offsetLftEdge.begin(); - return getNextInLftEdge(); + return getNextInLftEdge(); } -/********************************************************************** - * getNextInLftEdge() - Get the next template grid point in the first - * column. Returns NULL when there are no more - * points in the first column. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ - -GridPoint *GridTemplate::getNextInLftEdge(void) const -{ - while (_pointInLftEdgeIterator != _offsetLftEdge.end()) - { - GridOffset *offset = *_pointInLftEdgeIterator; - - _pointInLftEdgeIterator++; - - _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; - _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; - - if (_pointInGridReturn.x >= 0 && - _pointInGridReturn.x < _pointInGridNumX && - _pointInGridReturn.y >= 0 && - _pointInGridReturn.y < _pointInGridNumY) - { - return &_pointInGridReturn; - } - - } - - return (GridPoint *)NULL; +/////////////////////////////////////////////////////////////////////////////// +// +// Get the next template grid point in the first column. +// Returns NULL when there are no more points in the first column. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// + +GridPoint *GridTemplate::getNextInLftEdge(void) const { + + while(_pointInLftEdgeIterator != _offsetLftEdge.end()) { + + GridOffset *offset = *_pointInLftEdgeIterator; + + _pointInLftEdgeIterator++; + + _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; + _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; + + // Apply global wrap logic + _pointInGridReturn.x = (_wrapLon ? + positive_modulo(_pointInGridReturn.x, _pointInGridNumX) : + _pointInGridReturn.x); + + if(_pointInGridReturn.x >= 0 && + _pointInGridReturn.x < _pointInGridNumX && + _pointInGridReturn.y >= 0 && + _pointInGridReturn.y < _pointInGridNumY) { + return &_pointInGridReturn; + } + } // end while + + return (GridPoint *)NULL; } -/********************************************************************** - * getFirstInTopEdge() - Get the first template grid point in the - * top row. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the first template grid point in the top row. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// -GridPoint *GridTemplate::getFirstInTopEdge(void) const -{ - // Reset the iterator and send back the first point +GridPoint *GridTemplate::getFirstInTopEdge(void) const { - _pointInTopEdgeIterator = _offsetTopEdge.begin(); + // Reset the iterator and send back the first point + _pointInTopEdgeIterator = _offsetTopEdge.begin(); - return getNextInTopEdge(); + return getNextInTopEdge(); } -/********************************************************************** - * getNextInTopEdge() - Get the next template grid point in the top - * row. Returns NULL when there are no more - * points in the top row. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ - -GridPoint *GridTemplate::getNextInTopEdge(void) const -{ - while (_pointInTopEdgeIterator != _offsetTopEdge.end()) - { - GridOffset *offset = *_pointInTopEdgeIterator; - - _pointInTopEdgeIterator++; - - _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; - _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; - - if (_pointInGridReturn.x >= 0 && - _pointInGridReturn.x < _pointInGridNumX && - _pointInGridReturn.y >= 0 && - _pointInGridReturn.y < _pointInGridNumY) - { - return &_pointInGridReturn; - } - - } - - return (GridPoint *)NULL; +/////////////////////////////////////////////////////////////////////////////// +// +// Get the next template grid point in the top row. +// Returns NULL when there are no more points in the top row. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// + +GridPoint *GridTemplate::getNextInTopEdge(void) const { + + while(_pointInTopEdgeIterator != _offsetTopEdge.end()) { + GridOffset *offset = *_pointInTopEdgeIterator; + + _pointInTopEdgeIterator++; + + _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; + _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; + + // Apply global wrap logic + _pointInGridReturn.x = (_wrapLon ? + positive_modulo(_pointInGridReturn.x, _pointInGridNumX) : + _pointInGridReturn.x); + + if(_pointInGridReturn.x >= 0 && + _pointInGridReturn.x < _pointInGridNumX && + _pointInGridReturn.y >= 0 && + _pointInGridReturn.y < _pointInGridNumY) { + return &_pointInGridReturn; + } + } // end while + + return (GridPoint *)NULL; } -/********************************************************************** - * getFirstInRgtEdge() - Get the first template grid point in the - * right column. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the first template grid point in the right column. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// -GridPoint *GridTemplate::getFirstInRgtEdge(void) const -{ - // Reset the iterator and send back the first point +GridPoint *GridTemplate::getFirstInRgtEdge(void) const { - _pointInRgtEdgeIterator = _offsetRgtEdge.begin(); + // Reset the iterator and send back the first point + _pointInRgtEdgeIterator = _offsetRgtEdge.begin(); - return getNextInRgtEdge(); + return getNextInRgtEdge(); } -/********************************************************************** - * getNextInRgtEdge() - Get the next template grid point in the right - * column. Returns NULL when there are no more - * points in the right column. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ - -GridPoint *GridTemplate::getNextInRgtEdge(void) const -{ - while (_pointInRgtEdgeIterator != _offsetRgtEdge.end()) - { - GridOffset *offset = *_pointInRgtEdgeIterator; - - _pointInRgtEdgeIterator++; - - _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; - _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; - - if (_pointInGridReturn.x >= 0 && - _pointInGridReturn.x < _pointInGridNumX && - _pointInGridReturn.y >= 0 && - _pointInGridReturn.y < _pointInGridNumY) - { - return &_pointInGridReturn; - } - - } - - return (GridPoint *)NULL; +/////////////////////////////////////////////////////////////////////////////// +// +// Get the next template grid point in the right column. +// Returns NULL when there are no more points in the right column. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// + +GridPoint *GridTemplate::getNextInRgtEdge(void) const { + + while(_pointInRgtEdgeIterator != _offsetRgtEdge.end()) { + GridOffset *offset = *_pointInRgtEdgeIterator; + + _pointInRgtEdgeIterator++; + + _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; + _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; + + // Apply global wrap logic + _pointInGridReturn.x = (_wrapLon ? + positive_modulo(_pointInGridReturn.x, _pointInGridNumX) : + _pointInGridReturn.x); + + if(_pointInGridReturn.x >= 0 && + _pointInGridReturn.x < _pointInGridNumX && + _pointInGridReturn.y >= 0 && + _pointInGridReturn.y < _pointInGridNumY) { + return &_pointInGridReturn; + } + } // end while + + return (GridPoint *)NULL; } -/********************************************************************** - * getFirstInBotEdge() - Get the first template grid point in the - * bottom row. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Get the first template grid point in the bottom row. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// -GridPoint *GridTemplate::getFirstInBotEdge(void) const -{ - // Reset the iterator and send back the first point +GridPoint *GridTemplate::getFirstInBotEdge(void) const { - _pointInBotEdgeIterator = _offsetBotEdge.begin(); + // Reset the iterator and send back the first point + _pointInBotEdgeIterator = _offsetBotEdge.begin(); - return getNextInBotEdge(); + return getNextInBotEdge(); } -/********************************************************************** - * getNextInBotEdge() - Get the next template grid point in the bottom - * row. Returns NULL when there are no more points - * in the bottom row. - * - * Returns a pointer to a static object which must NOT be deleted by the - * calling routine. - */ - -GridPoint *GridTemplate::getNextInBotEdge(void) const -{ - while (_pointInBotEdgeIterator != _offsetBotEdge.end()) - { - GridOffset *offset = *_pointInBotEdgeIterator; - - _pointInBotEdgeIterator++; - - _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; - _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; - - if (_pointInGridReturn.x >= 0 && - _pointInGridReturn.x < _pointInGridNumX && - _pointInGridReturn.y >= 0 && - _pointInGridReturn.y < _pointInGridNumY) - { - return &_pointInGridReturn; - } - - } - - return (GridPoint *)NULL; +/////////////////////////////////////////////////////////////////////////////// +// +// Get the next template grid point in the bottom row. +// Returns NULL when there are no more points in the bottom row. +// Initialize the grid dimensions and base location. +// +// Returns a pointer to a static object which must NOT be deleted +// by the calling routine. +// +/////////////////////////////////////////////////////////////////////////////// + +GridPoint *GridTemplate::getNextInBotEdge(void) const { + + while(_pointInBotEdgeIterator != _offsetBotEdge.end()) { + GridOffset *offset = *_pointInBotEdgeIterator; + + _pointInBotEdgeIterator++; + + _pointInGridReturn.x = _pointInGridBase.x + offset->x_offset; + _pointInGridReturn.y = _pointInGridBase.y + offset->y_offset; + + // Apply global wrap logic + _pointInGridReturn.x = (_wrapLon ? + positive_modulo(_pointInGridReturn.x, _pointInGridNumX) : + _pointInGridReturn.x); + + if(_pointInGridReturn.x >= 0 && + _pointInGridReturn.x < _pointInGridNumX && + _pointInGridReturn.y >= 0 && + _pointInGridReturn.y < _pointInGridNumY) { + return &_pointInGridReturn; + } + } // end while + + return (GridPoint *)NULL; } -/********************************************************************** - * setGrid() - Initialize the grid dimensions and base location. - * - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Initialize the grid dimensions and base location. +// +/////////////////////////////////////////////////////////////////////////////// void GridTemplate::setGrid(const int &base_x, const int &base_y, - const int &nx, const int &ny) const -{ - // Set up the iterators and save the grid information + const int &nx, const int &ny) const { - _pointInGridIterator = _offsetList.begin(); + // Set up the iterators and save the grid information + _pointInGridIterator = _offsetList.begin(); - _pointInLftEdgeIterator = _offsetLftEdge.begin(); - _pointInRgtEdgeIterator = _offsetRgtEdge.begin(); - _pointInTopEdgeIterator = _offsetTopEdge.begin(); - _pointInBotEdgeIterator = _offsetBotEdge.begin(); + _pointInLftEdgeIterator = _offsetLftEdge.begin(); + _pointInRgtEdgeIterator = _offsetRgtEdge.begin(); + _pointInTopEdgeIterator = _offsetTopEdge.begin(); + _pointInBotEdgeIterator = _offsetBotEdge.begin(); - _pointInGridBase.x = base_x; - _pointInGridBase.y = base_y; + _pointInGridBase.x = base_x; + _pointInGridBase.y = base_y; - _pointInGridNumX = nx; - _pointInGridNumY = ny; + _pointInGridNumX = nx; + _pointInGridNumY = ny; - return; + return; } -/********************************************************************** - * incBaseX() - Increment the base_x location and reset the offset - * iterators. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Increment the base_x location and reset the offset iterators. +// +/////////////////////////////////////////////////////////////////////////////// + +void GridTemplate::incBaseX(const int &x_inc) const { -void GridTemplate::incBaseX(const int &x_inc) const -{ + _pointInGridBase.x += x_inc; - _pointInGridBase.x += x_inc; + // Apply global wrap logic + _pointInGridBase.x = (_wrapLon ? + positive_modulo(_pointInGridBase.x, _pointInGridNumX) : + _pointInGridBase.x); - if (_pointInGridBase.x < 0 || + if(_pointInGridBase.x < 0 || _pointInGridBase.x >= _pointInGridNumX ) { - mlog << Error << "\nGridTemplate::incBaseX() -> " - << "x (" << _pointInGridBase.x << ") out of range (0, " - << _pointInGridNumX << ").\n\n"; - exit(1); - } + mlog << Error << "\nGridTemplate::incBaseX() -> " + << "x (" << _pointInGridBase.x << ") out of range (0, " + << _pointInGridNumX << ").\n\n"; + exit(1); + } - return; + return; } -/********************************************************************** - * incBaseY() - Increment the base_y location and reset the offset - * iterators. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Increment the base_y location and reset the offset iterators. +// +/////////////////////////////////////////////////////////////////////////////// -void GridTemplate::incBaseY(const int &y_inc) const -{ +void GridTemplate::incBaseY(const int &y_inc) const { - _pointInGridBase.y += y_inc; + _pointInGridBase.y += y_inc; - if (_pointInGridBase.y < 0 || + if(_pointInGridBase.y < 0 || _pointInGridBase.y >= _pointInGridNumY ) { - mlog << Error << "\nGridTemplate::incBaseY() -> " - << "y (" << _pointInGridBase.y << ") out of range (0, " - << _pointInGridNumY << ").\n\n"; - exit(1); - } + mlog << Error << "\nGridTemplate::incBaseY() -> " + << "y (" << _pointInGridBase.y << ") out of range (0, " + << _pointInGridNumY << ").\n\n"; + exit(1); + } - return; + return; } -/********************************************************************** - * printOffsetList() - Print the offset list to the given stream. This - * is used for debugging. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Print the offset list to the given stream. +// This is used for debugging. +// +/////////////////////////////////////////////////////////////////////////////// -void GridTemplate::printOffsetList(FILE *stream) -{ - vector< GridOffset* >::iterator ol_iterator; +void GridTemplate::printOffsetList(FILE *stream) { + vector< GridOffset* >::iterator ol_iterator; - for (ol_iterator = _offsetList.begin(); + for(ol_iterator = _offsetList.begin(); ol_iterator != _offsetList.end(); - ol_iterator++) - { - GridOffset *offset = *ol_iterator; - - double x = (double)offset->x_offset; - double y = (double)offset->y_offset; - - double distance = sqrt((x * x) + (y * y)); + ol_iterator++) { + GridOffset *offset = *ol_iterator; - fprintf(stream, " %4d %4d %f\n", - offset->x_offset, offset->y_offset, distance); + double x = (double)offset->x_offset; + double y = (double)offset->y_offset; - } + double distance = sqrt((x * x) + (y * y)); + fprintf(stream, " %4d %4d %f\n", + offset->x_offset, offset->y_offset, distance); + } } -/********************************************************************** - * Private Member Functions * - **********************************************************************/ - -/********************************************************************** - * _addOffset() - Add the given offset to the offset list. - */ +/////////////////////////////////////////////////////////////////////////////// +// +// Add the given offset to the offset list. +// +/////////////////////////////////////////////////////////////////////////////// -void GridTemplate::_addOffset(int x_offset, int y_offset) -{ - GridOffset *offset = new GridOffset(x_offset, y_offset); +void GridTemplate::_addOffset(int x_offset, int y_offset) { + GridOffset *offset = new GridOffset(x_offset, y_offset); - _offsetList.push_back(offset); + _offsetList.push_back(offset); - return; + return; } -/********************************************************************** - * _setEdgeOffsets() - Process the current offset list and find the - * edges. - */ - -void GridTemplate::_setEdgeOffsets() -{ - vector< GridOffset* >::iterator v_iterator; - map< int, GridOffset* >::iterator m_iterator; - map< int, GridOffset* > min_x_by_y; - map< int, GridOffset* > max_x_by_y; - map< int, GridOffset* > min_y_by_x; - map< int, GridOffset* > max_y_by_x; - int x, y; - - // Loop over the offsets. - // For each row, find the min/max col. - // For each col, find the min/max row. - - for (v_iterator = _offsetList.begin(); +/////////////////////////////////////////////////////////////////////////////// +// +// Process the current offset list and find the edges. +// +/////////////////////////////////////////////////////////////////////////////// + +void GridTemplate::_setEdgeOffsets() { + vector::iterator v_iterator; + map::iterator m_iterator; + map min_x_by_y; + map max_x_by_y; + map min_y_by_x; + map max_y_by_x; + int x, y; + + // Loop over the offsets. + // For each row, find the min/max col. + // For each col, find the min/max row. + + for(v_iterator = _offsetList.begin(); v_iterator != _offsetList.end(); - v_iterator++) - { - GridOffset *offset = *v_iterator; - x = offset->x_offset; - y = offset->y_offset; - - // Min x for each y - if(min_x_by_y.count(y) == 0) { min_x_by_y[y] = offset; } - else if(x < min_x_by_y[y]->x_offset) { min_x_by_y[y] = offset; } - - // Max x for each y - if(max_x_by_y.count(y) == 0) { max_x_by_y[y] = offset; } - else if(x > max_x_by_y[y]->x_offset) { max_x_by_y[y] = offset; } - - // Min y for each x - if(min_y_by_x.count(x) == 0) { min_y_by_x[x] = offset; } - else if(y < min_y_by_x[x]->y_offset) { min_y_by_x[x] = offset; } - - // Max y for each x - if(max_y_by_x.count(x) == 0) { max_y_by_x[x] = offset; } - else if(y > max_y_by_x[x]->y_offset) { max_y_by_x[x] = offset; } - - } + v_iterator++) { + GridOffset *offset = *v_iterator; + x = offset->x_offset; + y = offset->y_offset; + + // Min x for each y + if(min_x_by_y.count(y) == 0) { min_x_by_y[y] = offset; } + else if(x < min_x_by_y[y]->x_offset) { min_x_by_y[y] = offset; } + + // Max x for each y + if(max_x_by_y.count(y) == 0) { max_x_by_y[y] = offset; } + else if(x > max_x_by_y[y]->x_offset) { max_x_by_y[y] = offset; } + + // Min y for each x + if(min_y_by_x.count(x) == 0) { min_y_by_x[x] = offset; } + else if(y < min_y_by_x[x]->y_offset) { min_y_by_x[x] = offset; } + + // Max y for each x + if(max_y_by_x.count(x) == 0) { max_y_by_x[x] = offset; } + else if(y > max_y_by_x[x]->y_offset) { max_y_by_x[x] = offset; } + } - // Store min_x_by_y map as _offsetLftEdge vector - for (m_iterator = min_x_by_y.begin(); + // Store min_x_by_y map as _offsetLftEdge vector + for(m_iterator = min_x_by_y.begin(); m_iterator != min_x_by_y.end(); - m_iterator++) - { - _offsetLftEdge.push_back(m_iterator->second); - } + m_iterator++) { + _offsetLftEdge.push_back(m_iterator->second); + } - // Store max_x_by_y map as _offsetRgtEdge vector - for (m_iterator = max_x_by_y.begin(); + // Store max_x_by_y map as _offsetRgtEdge vector + for(m_iterator = max_x_by_y.begin(); m_iterator != max_x_by_y.end(); - m_iterator++) - { - _offsetRgtEdge.push_back(m_iterator->second); - } + m_iterator++) { + _offsetRgtEdge.push_back(m_iterator->second); + } - // Store max_y_by_x map as _offsetTopEdge vector - for (m_iterator = max_y_by_x.begin(); + // Store max_y_by_x map as _offsetTopEdge vector + for(m_iterator = max_y_by_x.begin(); m_iterator != max_y_by_x.end(); - m_iterator++) - { - _offsetTopEdge.push_back(m_iterator->second); - } + m_iterator++) { + _offsetTopEdge.push_back(m_iterator->second); + } - // Store min_y_by_x map as _offsetBotEdge vector - for (m_iterator = min_y_by_x.begin(); + // Store min_y_by_x map as _offsetBotEdge vector + for(m_iterator = min_y_by_x.begin(); m_iterator != min_y_by_x.end(); - m_iterator++) - { - _offsetBotEdge.push_back(m_iterator->second); - } + m_iterator++) { + _offsetBotEdge.push_back(m_iterator->second); + } - return; + return; } -//////////////////////////////////////////////////////////////// -// GridTemplateFactory -//////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +// +// Code for class GridTemplateFactory. +// +/////////////////////////////////////////////////////////////////////////////// + GridTemplateFactory::GridTemplateFactory() { enum_to_string.resize(GridTemplate_NUM_TEMPLATES); enum_to_string[GridTemplate_None] = ""; enum_to_string[GridTemplate_Square] = "SQUARE"; enum_to_string[GridTemplate_Circle] = "CIRCLE"; - //enum_to_string[GridTemplate_Rectangle] = "RECTANGLE"; } +/////////////////////////////////////////////////////////////////////////////// + GridTemplateFactory::~GridTemplateFactory() { enum_to_string.clear(); } -///////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +// // log & exit on failure +// +/////////////////////////////////////////////////////////////////////////////// + GridTemplateFactory::GridTemplates GridTemplateFactory::string2Enum(string target) { - for (unsigned int ix = 0; ix < GridTemplate_NUM_TEMPLATES; ix++){ - if (enum_to_string[ix] == target){ + for(unsigned int ix = 0; ix < GridTemplate_NUM_TEMPLATES; ix++) { + if(enum_to_string[ix] == target) { return static_cast(ix); } } @@ -621,13 +622,15 @@ GridTemplateFactory::GridTemplates GridTemplateFactory::string2Enum(string targe exit(1); } - -///////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +// // log & exit on failure -string GridTemplateFactory::enum2String(GridTemplates target) { +// +/////////////////////////////////////////////////////////////////////////////// - if( static_cast(target) > enum_to_string.size() - 1){ +string GridTemplateFactory::enum2String(GridTemplates target) { + if(static_cast(target) > enum_to_string.size() - 1) { mlog << Error << "\nGridTemplateFactory::enum2String() -> " << "target out of range " << target << " > " << (static_cast(enum_to_string.size()) - 1) @@ -637,26 +640,36 @@ string GridTemplateFactory::enum2String(GridTemplates target) { return enum_to_string[static_cast(target)]; } -///////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +// // Caller assumes ownership of the returned pointer. -GridTemplate* GridTemplateFactory::buildGT(string gt, int width) { - return buildGT(string2Enum(gt), width); +// +/////////////////////////////////////////////////////////////////////////////// + +GridTemplate* GridTemplateFactory::buildGT(string gt, int width, bool wrap_lon) { + return buildGT(string2Enum(gt), width, wrap_lon); } -///////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// +// // Caller assumes ownership of the returned pointer. -GridTemplate* GridTemplateFactory::buildGT(GridTemplates gt, int width) { +// +/////////////////////////////////////////////////////////////////////////////// + +GridTemplate* GridTemplateFactory::buildGT(GridTemplates gt, int width, bool wrap_lon) { switch (gt) { case(GridTemplate_Square): - return new RectangularTemplate(width,width); + return new RectangularTemplate(width, width, wrap_lon); case(GridTemplate_Circle): - return new CircularTemplate(width); + return new CircularTemplate(width, wrap_lon); default: mlog << Error << "\nbuildGT() -> " - << "Unexpected gt value of " << gt << ".\n\n"; + << "Unexpected GridTemplates value (" << gt << ").\n\n"; exit(1); } } + +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/GridTemplate.h b/met/src/basic/vx_util/GridTemplate.h index 28cfffec62..fd356c5f36 100644 --- a/met/src/basic/vx_util/GridTemplate.h +++ b/met/src/basic/vx_util/GridTemplate.h @@ -1,53 +1,30 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -/* RCS info - * $Author: dixon $ - * $Locker: $ - * $Date: 2016/03/03 19:21:31 $ - * $Id: GridTemplate.hh,v 1.9 2016/03/03 19:21:31 dixon Exp $ - * $Revision: 1.9 $ - * $State: Exp $ - */ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ - -/************************************************************************ - * GridTemplate.hh: class implementing a template to be applied to - * gridded data. - * - * RAP, NCAR, Boulder CO - * - * January 1999 - * - * Nancy Rehak - * - ************************************************************************/ - -#ifndef GridTemplate_HH -#define GridTemplate_HH +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: GridTemplate.h +// +// Description: +// Class implementing a template to be applied to +// gridded data. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-99 Rehak Initial version. +// 001 09-07-21 Halley Gotway Add wrap_lon. +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __GRID_TEMPLATE_H__ +#define __GRID_TEMPLATE_H__ + +/////////////////////////////////////////////////////////////////////////////// #include #include @@ -55,177 +32,154 @@ #include "GridOffset.h" #include "GridPoint.h" -using namespace std; - -class GridTemplate -{ - public: - - // Constructors - - GridTemplate(void); - GridTemplate(const GridTemplate& rhs); - - // Destructor - - virtual ~GridTemplate(void); +/////////////////////////////////////////////////////////////////////////////// - // Methods for iterating through the template within the grid centered - // on the given point. To use these methods, first call getFirstInGrid() - // to get the first point. Then call getNextInGrid() to get all remaining - // points until a (GridPoint *)NULL is returned. Any time getFirstInGrid() - // is called, the iteration will be cleared and will start over again. - // - // base_x and base_y give the coordinates of the point around which the - // template is to be applied. - // - // These routines return a point to a static object which must NOT be - // deleted by the calling routine. +using namespace std; - GridPoint *getFirstInGrid(const int &base_x, const int &base_y, - const int &nx, const int &ny) const; - GridPoint *getNextInGrid(void) const; +class GridTemplate { - GridPoint *getFirst(const int &base_x, const int &base_y, - const int &nx, const int &ny) const; - GridPoint *getNext(void) const; + public: - GridPoint *getFirstInLftEdge(void) const; - GridPoint *getNextInLftEdge(void) const; + GridTemplate(void); + GridTemplate(const GridTemplate& rhs); + virtual ~GridTemplate(void); - GridPoint *getFirstInRgtEdge(void) const; - GridPoint *getNextInRgtEdge(void) const; + // Methods for iterating through the template within the grid centered + // on the given point. To use these methods, first call getFirstInGrid() + // to get the first point. Then call getNextInGrid() to get all remaining + // points until a (GridPoint *)NULL is returned. Any time getFirstInGrid() + // is called, the iteration will be cleared and will start over again. + // + // base_x and base_y give the coordinates of the point around which the + // template is to be applied. + // + // These routines return a point to a static object which must NOT be + // deleted by the calling routine. - GridPoint *getFirstInTopEdge(void) const; - GridPoint *getNextInTopEdge(void) const; + GridPoint *getFirstInGrid(const int &base_x, const int &base_y, + const int &nx, const int &ny) const; + GridPoint *getNextInGrid(void) const; - GridPoint *getFirstInBotEdge(void) const; - GridPoint *getNextInBotEdge(void) const; + GridPoint *getFirst(const int &base_x, const int &base_y, + const int &nx, const int &ny) const; + GridPoint *getNext(void) const; - void setGrid(const int &base_x, const int &base_y, - const int &nx, const int &ny) const; + GridPoint *getFirstInLftEdge(void) const; + GridPoint *getNextInLftEdge(void) const; - void incBaseX(const int &x_inc) const; - void incBaseY(const int &y_inc) const; + GridPoint *getFirstInRgtEdge(void) const; + GridPoint *getNextInRgtEdge(void) const; - // Printing methods + GridPoint *getFirstInTopEdge(void) const; + GridPoint *getNextInTopEdge(void) const; - void printOffsetList(FILE *stream); + GridPoint *getFirstInBotEdge(void) const; + GridPoint *getNextInBotEdge(void) const; - // Access methods + void setGrid(const int &base_x, const int &base_y, + const int &nx, const int &ny) const; - inline void addOffset(const GridOffset &offset) - { - _offsetList.push_back(new GridOffset(offset.x_offset, - offset.y_offset)); - } + void incBaseX(const int &x_inc) const; + void incBaseY(const int &y_inc) const; - inline void addOffset(const int x_offset, const int y_offset) - { - _offsetList.push_back(new GridOffset(x_offset, y_offset)); - } + // Printing methods - int size(void) const - { - return _offsetList.size(); - } + void printOffsetList(FILE *stream); - int getNumPts(void) const - { - return _offsetList.size(); - } + // Access methods - virtual const char* getClassName(void) const = 0; + inline void addOffset(const GridOffset &offset) { + _offsetList.push_back(new GridOffset(offset.x_offset, + offset.y_offset)); + } - virtual int getWidth() const = 0; + inline void addOffset(const int x_offset, const int y_offset) { + _offsetList.push_back(new GridOffset(x_offset, y_offset)); + } - protected: + int size(void) const { + return _offsetList.size(); + } - // The offsets that make up the circle + int getNumPts(void) const { + return _offsetList.size(); + } - vector< GridOffset* > _offsetList; + int getWrapLon(void) const { + return _wrapLon; + } - // The offsets that define the first and last rows and columns + virtual const char* getClassName(void) const = 0; - vector< GridOffset* > _offsetLftEdge; // not allocated - vector< GridOffset* > _offsetRgtEdge; // not allocated - vector< GridOffset* > _offsetTopEdge; // not allocated - vector< GridOffset* > _offsetBotEdge; // not allocated + virtual int getWidth() const = 0; - // Iterator for finding points within a grid + protected: - mutable vector< GridOffset* >::const_iterator _pointInGridIterator; - mutable vector< GridOffset* >::const_iterator _pointInLftEdgeIterator; - mutable vector< GridOffset* >::const_iterator _pointInRgtEdgeIterator; - mutable vector< GridOffset* >::const_iterator _pointInTopEdgeIterator; - mutable vector< GridOffset* >::const_iterator _pointInBotEdgeIterator; + bool _wrapLon; - mutable GridPoint _pointInGridBase; - mutable int _pointInGridNumX; - mutable int _pointInGridNumY; - mutable GridPoint _pointInGridReturn; + // The offsets that make up the circle + vector _offsetList; - // Add the given offset to the offset list + // The offsets that define the first and last rows and columns + vector _offsetLftEdge; // not allocated + vector _offsetRgtEdge; // not allocated + vector _offsetTopEdge; // not allocated + vector _offsetBotEdge; // not allocated - void _addOffset(int x_offset, int y_offset); + // Iterator for finding points within a grid + mutable vector::const_iterator _pointInGridIterator; + mutable vector::const_iterator _pointInLftEdgeIterator; + mutable vector::const_iterator _pointInRgtEdgeIterator; + mutable vector::const_iterator _pointInTopEdgeIterator; + mutable vector::const_iterator _pointInBotEdgeIterator; - // Determine the offsets for the First/Last Row/Column + mutable GridPoint _pointInGridBase; + mutable int _pointInGridNumX; + mutable int _pointInGridNumY; + mutable GridPoint _pointInGridReturn; - void _setEdgeOffsets(); + // Add the given offset to the offset list + void _addOffset(int x_offset, int y_offset); + // Determine the offsets for the First/Last Row/Column + void _setEdgeOffsets(); }; +/////////////////////////////////////////////////////////////////////////////// + // The factory knows about child classes and can create them for you class GridTemplateFactory { - public: - - // constructor - GridTemplateFactory(); + public: - // destructor - ~GridTemplateFactory(); - // do not assign specific values to these enumes. - // other code requires them to start at zero and increase by 1 - // make sure GridTemplate_NUM_TEMPLATES is always last. - // TODO: rename this Shape - enum GridTemplates { - GridTemplate_None, - GridTemplate_Square, - GridTemplate_Circle, - //GridTemplate_Rectangle, - GridTemplate_NUM_TEMPLATES - }; + GridTemplateFactory(); + ~GridTemplateFactory(); - // String corresponding to the enumerated values above - vector enum_to_string; + // do not assign specific values to these enumes. + // other code requires them to start at zero and increase by 1 + // make sure GridTemplate_NUM_TEMPLATES is always last. + enum GridTemplates { + GridTemplate_None, + GridTemplate_Square, + GridTemplate_Circle, + GridTemplate_NUM_TEMPLATES + }; -///////////////////////////////////////////////////////////////// -// exit on failure - GridTemplates string2Enum(string target); + // String corresponding to the enumerated values above + vector enum_to_string; -////////////////////////////////////////////////////////////// -// exit on failure - string enum2String(GridTemplates gt); + GridTemplates string2Enum(string target); + string enum2String(GridTemplates gt); -///////////////////////////////////////////////////////////////// -// Caller assumes ownership of the returned pointer. - GridTemplate* buildGT(string gt, int width); + GridTemplate* buildGT(string gt, int width, bool wrap_lon); + GridTemplate* buildGT(GridTemplates gt, int width, bool wrap_lon); -///////////////////////////////////////////////////////////////// -// Caller assumes ownership of the returned pointer. - GridTemplate* buildGT(GridTemplates gt, int width); - - //FUTURE DEV NOTE: - // If we ever decide to support rectangles or elipses, we are going to - // need to modify these methods to take more than just a single width & shape - // to define the grid template. I suggest building a structure/union that holds - // all the defining characteristics of a grid template (shape, height, width, radius, etc.) - // and passing that around everywhere that we currently pass around a shape & width. +}; +/////////////////////////////////////////////////////////////////////////////// -}; +#endif // __GRID_TEMPLATE_H__ -#endif +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/RectangularTemplate.cc b/met/src/basic/vx_util/RectangularTemplate.cc index 59c06901a0..7368e2125c 100644 --- a/met/src/basic/vx_util/RectangularTemplate.cc +++ b/met/src/basic/vx_util/RectangularTemplate.cc @@ -1,48 +1,25 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -// RCS info -// $Author: dixon $ -// $Locker: $ -// $Date: 2016/03/03 18:19:27 $ -// $Id: RectangularTemplate.cc,v 1.4 2016/03/03 18:19:27 dixon Exp $ -// $Revision: 1.4 $ -// $State: Exp $ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ -/********************************************************************* - * RectangularTemplate.cc: class implementing a Rectangular template - * to be applied on gridded data. - * - * RAP, NCAR, Boulder CO - * - * January 2007 - * - * Dan Megenhardt - * - *********************************************************************/ +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: RectangularTemplate.cc +// +// Description: +// Class implementing a Rectangular template to be +// applied on gridded data. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-07 Megenhardt Initial version. +// 001 09-07-21 Halley Gotway Add wrap_lon. +// +/////////////////////////////////////////////////////////////////////////////// #include @@ -52,82 +29,63 @@ #include #include #include + using namespace std; -/********************************************************************** - * Constructor - */ - -RectangularTemplate::RectangularTemplate(int height, int width) : - GridTemplate(), - _height(height), - _width(width) -{ - - /* - CONVENTION: - If width is even, we push the center to the lower left corner - and have the "extra" in the right & above. - */ - - bool evenWidth = ((width % 2) == 0); - bool evenHeight = ((height % 2) == 0); - - // equivalent of w/2 for even & (w-1)/2 for odd - int halfwidth = width / 2; - int halfheight = height / 2; - - int xmin = halfwidth; - int ymin = halfheight; - - if (evenWidth){ - xmin--; - } - if (evenHeight){ - ymin--; - } - - // Create the offsets list - for (int y = -ymin; y <= halfheight; y++) - { - for (int x = -xmin; x <= halfheight; x++) - { - _addOffset(x, y); - } /* endfor x = 0 */ - } /* endfor y = 0*/ - - _setEdgeOffsets(); -} +/////////////////////////////////////////////////////////////////////////////// +RectangularTemplate::RectangularTemplate(int height, int width, bool wrap_lon) : + GridTemplate(), _height(height), _width(width) { -/********************************************************************** - * Destructor - */ + _wrapLon = wrap_lon; -RectangularTemplate::~RectangularTemplate(void) -{ - // Do nothing -} + // + // CONVENTION: + // If width is even, we push the center to the lower left corner + // and have the "extra" in the right & above. + // + bool evenWidth = ((width % 2) == 0); + bool evenHeight = ((height % 2) == 0); -/********************************************************************** - * printOffsetList() - Print the offset list to the given stream. This - * is used for debugging. - */ + // equivalent of w/2 for even & (w-1)/2 for odd + int halfwidth = width / 2; + int halfheight = height / 2; -void RectangularTemplate::printOffsetList(FILE *stream) -{ - fprintf(stream, "\n\n"); - fprintf(stream, "Rectangular template:"); - fprintf(stream, " height = %d\n", _height); - fprintf(stream, " width = %d\n", _width); + int xmin = halfwidth; + int ymin = halfheight; - fprintf(stream, " grid points:\n"); + if(evenWidth) xmin--; + if(evenHeight) ymin--; + + // Create the offsets list + for(int y = -ymin; y <= halfheight; y++) { + for(int x = -xmin; x <= halfwidth; x++) { + _addOffset(x, y); + } // end for x + } // end for y + + _setEdgeOffsets(); - GridTemplate::printOffsetList(stream); } +/////////////////////////////////////////////////////////////////////////////// + +RectangularTemplate::~RectangularTemplate(void) { + // Do nothing +} + +/////////////////////////////////////////////////////////////////////////////// + +void RectangularTemplate::printOffsetList(FILE *stream) { + fprintf(stream, "\n\n"); + fprintf(stream, "Rectangular template:"); + fprintf(stream, " height = %d\n", _height); + fprintf(stream, " width = %d\n", _width); + fprintf(stream, " wrap_lon = %d\n", _wrapLon); + fprintf(stream, " grid points:\n"); + + GridTemplate::printOffsetList(stream); +} -/********************************************************************** - * Private Member Functions * - **********************************************************************/ +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/RectangularTemplate.h b/met/src/basic/vx_util/RectangularTemplate.h index b45dcb6a77..507286adbc 100644 --- a/met/src/basic/vx_util/RectangularTemplate.h +++ b/met/src/basic/vx_util/RectangularTemplate.h @@ -1,57 +1,30 @@ -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -// ** Copyright UCAR (c) 1990 - 2021 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) -// ** Boulder, Colorado, USA -// ** BSD licence applies - redistribution and use in source and binary -// ** forms, with or without modification, are permitted provided that -// ** the following conditions are met: -// ** 1) If the software is modified to produce derivative works, -// ** such modified software should be clearly marked, so as not -// ** to confuse it with the version available from UCAR. -// ** 2) Redistributions of source code must retain the above copyright -// ** notice, this list of conditions and the following disclaimer. -// ** 3) Redistributions in binary form must reproduce the above copyright -// ** notice, this list of conditions and the following disclaimer in the -// ** documentation and/or other materials provided with the distribution. -// ** 4) Neither the name of UCAR nor the names of its contributors, -// ** if any, may be used to endorse or promote products derived from -// ** this software without specific prior written permission. -// ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS -// ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED -// ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* -/*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/ - -/* RCS info - * $Author: dixon $ - * $Locker: $ - * $Date: 2016/03/03 19:21:31 $ - * $Id: RectangularTemplate.hh,v 1.3 2016/03/03 19:21:31 dixon Exp $ - * $Revision: 1.3 $ - * $State: Exp $ - */ - -/**-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**/ - -/************************************************************************ - * RectangularTemplate.hh: class implementing a Rectangular template to be - * applied on gridded data. - * - * RAP, NCAR, Boulder CO - * - * January 2007 - * - * Dan Megenhardt - * - ************************************************************************/ - -#ifndef RectangularTemplate_HH -#define RectangularTemplate_HH - -/* - **************************** includes ********************************** - */ +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +/////////////////////////////////////////////////////////////////////////////// +// +// Filename: RectangularTemplate.h +// +// Description: +// Class implementing a Rectangular template to be +// applied on gridded data. +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 01-01-07 Megenhardt Initial version. +// 001 09-07-21 Halley Gotway Add wrap_lon. +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef __RECTANGULAR_TEMPLATE_H__ +#define __RECTANGULAR_TEMPLATE_H__ + +/////////////////////////////////////////////////////////////////////////////// #include @@ -61,71 +34,44 @@ #include #include -/* - ******************************* defines ******************************** - */ +/////////////////////////////////////////////////////////////////////////////// -/* - ******************************* structures ***************************** - */ +class RectangularTemplate : public GridTemplate { -/* - ************************* global variables ***************************** - */ + public: -/* - ***************************** function prototypes ********************** - */ + RectangularTemplate(int height, int width, bool wrap_lon); + virtual ~RectangularTemplate(void); -/* - ************************* class definitions **************************** - */ + void printOffsetList(FILE *stream); -class RectangularTemplate : public GridTemplate -{ - public: + // Access methods - // Constructor - RectangularTemplate(int height, int width); + int getHeight(void) const { + return _height; + } - // Destructor + int getWidth(void) const { + return _width; + } - virtual ~RectangularTemplate(void); + const char* getClassName(void) const { + return RectangularTemplate::_className(); + } - // Print the offset list to the given stream. This is used for debugging. + private: - void printOffsetList(FILE *stream); - - // Access methods - - int getHeight(void) const - { - return _height; - } - - int getWidth(void) const - { - return _width; - } - - const char* getClassName(void) const{ - return RectangularTemplate::_className(); - } - - private: - - // The box dimensions - - int _height; - int _width; - - // Return the class name for error messages. - static const char* _className(void) - { - return("RectangularTemplate"); - } + int _height; + int _width; + // Return the class name for error messages. + static const char* _className(void) { + return("RectangularTemplate"); + } }; +/////////////////////////////////////////////////////////////////////////////// + +#endif // __RECTANGULAR_TEMPLATE_H__ -#endif +/////////////////////////////////////////////////////////////////////////////// diff --git a/met/src/basic/vx_util/data_plane_util.cc b/met/src/basic/vx_util/data_plane_util.cc index 0004abe0f8..3f5b5f7278 100644 --- a/met/src/basic/vx_util/data_plane_util.cc +++ b/met/src/basic/vx_util/data_plane_util.cc @@ -95,7 +95,7 @@ void rescale_probability(DataPlane &dp) { void smooth_field(const DataPlane &dp, DataPlane &smooth_dp, InterpMthd mthd, int width, const GridTemplateFactory::GridTemplates shape, - double t, const GaussianInfo &gaussian) { + bool wrap_lon, double t, const GaussianInfo &gaussian) { double v = 0.0; int x, y; @@ -107,10 +107,11 @@ void smooth_field(const DataPlane &dp, DataPlane &smooth_dp, // build the grid template GridTemplateFactory gtf; - GridTemplate* gt = gtf.buildGT(shape, width); + GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon); mlog << Debug(3) - << "Smoothing field using the " << interpmthd_to_string(mthd) + << "Smoothing " << (wrap_lon ? "global" : "non-global") + << " field using the " << interpmthd_to_string(mthd) << "(" << gt->size() << ") " << gt->getClassName() << " interpolation method.\n"; @@ -185,10 +186,10 @@ void smooth_field(const DataPlane &dp, DataPlane &smooth_dp, DataPlane smooth_field(const DataPlane &dp, InterpMthd mthd, int width, const GridTemplateFactory::GridTemplates shape, - double t, const GaussianInfo &gaussian) { + bool wrap_lon, double t, const GaussianInfo &gaussian) { DataPlane smooth_dp; - smooth_field(dp, smooth_dp, mthd, width, shape, t, gaussian); + smooth_field(dp, smooth_dp, mthd, width, shape, wrap_lon, t, gaussian); return(smooth_dp); } @@ -202,23 +203,53 @@ DataPlane smooth_field(const DataPlane &dp, void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp, int width, const GridTemplateFactory::GridTemplates shape, - SingleThresh t, double vld_t) { + bool wrap_lon, SingleThresh t, + const DataPlane *cmn, const DataPlane *csd, double vld_t) { GridPoint *gp = NULL; int x, y; int n_vld = 0; int n_thr = 0; double v; + double bad = bad_data_double; + bool use_climo = false; // Check that width is set to 1 or greater if(width < 1) { mlog << Error << "\nfractional_coverage() -> " - << "Grid must have at least one point in it. \n\n"; + << "grid must have at least one point in it. \n\n"; exit(1); } + // Check climatology data, if needed + if(cmn && csd) { + if(!cmn->is_empty() && !csd->is_empty()) use_climo = true; + } + + // Check climatology dimensions + if(use_climo) { + + // Check dimensions + if(cmn->nx() != dp.nx() || cmn->ny() != dp.ny()) { + mlog << Error << "\nfractional_coverage() -> " + << "climatology mean dimension (" + << cmn->nx() << ", " << cmn->ny() + << ") does not match the data dimenion (" + << dp.nx() << ", " << dp.ny() << ")!\n\n"; + exit(1); + } + if(csd->nx() != dp.nx() || csd->ny() != dp.ny()) { + mlog << Error << "\nfractional_coverage() -> " + << "climatology standard deviation dimension (" + << csd->nx() << ", " << csd->ny() + << ") does not match the data dimenion (" + << dp.nx() << ", " << dp.ny() << ")!\n\n"; + exit(1); + } + } + // Build the grid template GridTemplateFactory gtf; - GridTemplate* gt = gtf.buildGT(shape, width); + GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon); mlog << Debug(3) << "Computing fractional coverage field using the " @@ -246,7 +277,9 @@ void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp, gp = gt->getNextInGrid()) { if(is_bad_data(v = dp.get(gp->x, gp->y))) continue; n_vld++; - if(t.check(v)) n_thr++; + if(t.check(v, + (use_climo ? cmn->get(gp->x, gp->y) : bad), + (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++; } } // Subtract off the bottom edge, shift up, and add the top. @@ -258,7 +291,9 @@ void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp, gp = gt->getNextInBotEdge()) { if(is_bad_data(v = dp.get(gp->x, gp->y))) continue; n_vld--; - if(t.check(v)) n_thr--; + if(t.check(v, + (use_climo ? cmn->get(gp->x, gp->y) : bad), + (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr--; } // Increment Y @@ -270,7 +305,9 @@ void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp, gp = gt->getNextInTopEdge()) { if(is_bad_data(v = dp.get(gp->x, gp->y))) continue; n_vld++; - if(t.check(v)) n_thr++; + if(t.check(v, + (use_climo ? cmn->get(gp->x, gp->y) : bad), + (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++; } } @@ -291,161 +328,6 @@ void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp, return; } -//////////////////////////////////////////////////////////////////////// -// -// Convert the DataPlane field to the corresponding fractional coverage -// using the threshold critea specified. -// -//////////////////////////////////////////////////////////////////////// - -void fractional_coverage_square(const DataPlane &dp, DataPlane &frac_dp, - int width, SingleThresh t, double vld_t) { - int i, j, k, n, x, y, x_ll, y_ll, y_ur, xx, yy, half_width; - double v; - int count_vld = 0; - int count_thr = 0; - NumArray box_na; - - mlog << Debug(3) - << "Computing fractional coverage field using the " - << t.get_str() << " threshold and the " - << interpmthd_to_string(InterpMthd_Nbrhd) - << "(" << width*width << ") interpolation method.\n"; - - // Check that width is set to 1 or greater - if(width < 1) { - mlog << Error << "\nfractional_coverage_square() -> " - << "width must be set to a value of 1 or greater.\n\n"; - exit(1); - } - - // Initialize the fractional coverage field - frac_dp = dp; - frac_dp.set_constant(bad_data_double); - - // Compute the box half-width - half_width = (width - 1)/2; - - // Initialize the box - for(i=0; i= dp.nx() || - yy < 0 || yy >= dp.ny()) { - k = bad_data_int; - } - // Check v to see if it meets the threshold criteria - else { - v = dp.get(xx, yy); - if(is_bad_data(v)) k = bad_data_int; - else if(t.check(v)) k = 1; - else k = 0; - } - box_na.set(n, k); - - // Increment the counts - if(!is_bad_data(k)) { - count_vld += 1; - count_thr += k; - } - - } // end for j - } // end for i - } // end if - - // Otherwise, update one row of the box - else { - - // Compute the row of the neighborhood box to be updated - j = (y - 1) % width; - - for(i=0; i= dp.nx() || - yy < 0 || yy >= dp.ny()) { - k = bad_data_int; - } - // Check v to see if it meets the threshold criteria - else { - v = dp.get(xx, yy); - if(is_bad_data(v)) k = bad_data_int; - else if(t.check(v)) k = 1; - else k = 0; - } - box_na.set(n, k); - - // Increment the counts - if(!is_bad_data(k)) { - count_vld += 1; - count_thr += k; - } - - } // end for i - } // end else - - // Check whether enough valid grid points were found - if((double) count_vld/(width*width) < vld_t || - count_vld == 0) { - v = bad_data_double; - } - // Compute the fractional coverage - else { - v = (double) count_thr/count_vld; - } - - // Store the fractional coverage value - frac_dp.set(v, x, y); - - } // end for y - } // end for x - - return; - -} - //////////////////////////////////////////////////////////////////////// // // Select points inside the mask and write them to a NumArray. diff --git a/met/src/basic/vx_util/data_plane_util.h b/met/src/basic/vx_util/data_plane_util.h index 401a22c2f7..c5b8a4b65a 100644 --- a/met/src/basic/vx_util/data_plane_util.h +++ b/met/src/basic/vx_util/data_plane_util.h @@ -43,19 +43,17 @@ extern void rescale_probability(DataPlane &); extern void smooth_field(const DataPlane &dp, DataPlane &smooth_dp, InterpMthd mthd, int width, const GridTemplateFactory::GridTemplates shape, - double t, const GaussianInfo &gaussian); + bool wrap_lon, double t, const GaussianInfo &gaussian); extern DataPlane smooth_field(const DataPlane &dp, InterpMthd mthd, int width, const GridTemplateFactory::GridTemplates shape, - double t, const GaussianInfo &gaussian); + bool wrap_lon, double t, const GaussianInfo &gaussian); extern void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp, int width, const GridTemplateFactory::GridTemplates shape, - SingleThresh t, double vld_t); - -extern void fractional_coverage_square(const DataPlane &dp, DataPlane &frac_dp, - int width, SingleThresh t, double vld_t); + bool wrap_lon, SingleThresh t, + const DataPlane *cmn, const DataPlane *csd, double vld_t); extern void apply_mask(const DataPlane &, const MaskPlane &, NumArray &); extern void apply_mask(DataPlane &, const MaskPlane &); diff --git a/met/src/basic/vx_util/interp_util.cc b/met/src/basic/vx_util/interp_util.cc index 1f1de88bea..ac37311a9e 100644 --- a/met/src/basic/vx_util/interp_util.cc +++ b/met/src/basic/vx_util/interp_util.cc @@ -18,7 +18,6 @@ using namespace std; #include #include -//#include "interp_mthd.h" #include "interp_util.h" #include "GridTemplate.h" #include "RectangularTemplate.h" @@ -739,26 +738,30 @@ double interp_nbrhd(const DataPlane &dp, const GridTemplate >, int x, int y, // //////////////////////////////////////////////////////////////////////// -double interp_bilin(const DataPlane &dp, double obs_x, double obs_y, - const MaskPlane *mp) { - int x, y; +double interp_bilin(const DataPlane &dp, bool wrap_lon, + double obs_x, double obs_y, const MaskPlane *mp) { + int x, xp1, y; double bilin_v, dx, dy; double wtsw, wtse, wtnw, wtne; x = nint(floor(obs_x)); y = nint(floor(obs_y)); + // Apply global wrap logic to x and x+1 + x = (wrap_lon ? positive_modulo(x, dp.nx()) : x); + xp1 = (wrap_lon ? positive_modulo(x+1, dp.nx()) : x+1); + // Check the optional mask if(mp) { if(!(*mp)(x, y ) || - !(*mp)(x+1, y ) || + !(*mp)(xp1, y ) || !(*mp)(x, y+1) || - !(*mp)(x+1, y+1)) return(bad_data_double); + !(*mp)(xp1, y+1)) return(bad_data_double); } // Compute dx and dy - dx = obs_x - x; - dy = obs_y - y; + dx = obs_x - nint(floor(obs_x)); + dy = obs_y - nint(floor(obs_y)); // Compute weights for 4 corner points wtsw = (1.0-dx) * (1.0-dy); @@ -767,7 +770,7 @@ double interp_bilin(const DataPlane &dp, double obs_x, double obs_y, wtne = dx * dy; // On the grid boundary, check each corner point - if(x < 0 || x+1 >= dp.nx() || y < 0 || y+1 >= dp.ny()) { + if(x < 0 || xp1 >= dp.nx() || y < 0 || y+1 >= dp.ny()) { // Initialize weighted sum bilin_v = 0.0; @@ -784,12 +787,12 @@ double interp_bilin(const DataPlane &dp, double obs_x, double obs_y, // Southeast Corner if(!is_bad_data(bilin_v) && !is_eq(wtse, 0.0)) { - if(x+1 < 0 || x+1 >= dp.nx() || y < 0 || y >= dp.ny()) + if(xp1 < 0 || xp1 >= dp.nx() || y < 0 || y >= dp.ny()) bilin_v = bad_data_double; - else if(is_bad_data(dp.get(x+1, y))) + else if(is_bad_data(dp.get(xp1, y))) bilin_v = bad_data_double; else - bilin_v += wtse * dp.get(x+1, y); + bilin_v += wtse * dp.get(xp1, y); } // Northwest Corner @@ -804,27 +807,27 @@ double interp_bilin(const DataPlane &dp, double obs_x, double obs_y, // Northeast Corner if(!is_bad_data(bilin_v) && !is_eq(wtne, 0.0)) { - if(x+1 < 0 || x+1 >= dp.nx() || y+1 < 0 || y+1 >= dp.ny()) + if(xp1 < 0 || xp1 >= dp.nx() || y+1 < 0 || y+1 >= dp.ny()) bilin_v = bad_data_double; - else if(is_bad_data(dp.get(x+1, y+1))) + else if(is_bad_data(dp.get(xp1, y+1))) bilin_v = bad_data_double; else - bilin_v += wtne * dp.get(x+1, y+1); + bilin_v += wtne * dp.get(xp1, y+1); } } // On the grid interior, compute weighted average else { if(is_bad_data(dp.get(x, y )) || - is_bad_data(dp.get(x+1, y )) || + is_bad_data(dp.get(xp1, y )) || is_bad_data(dp.get(x, y+1)) || - is_bad_data(dp.get(x+1, y+1))) { + is_bad_data(dp.get(xp1, y+1))) { bilin_v = bad_data_double; } else { bilin_v = wtsw * dp.get(x, y) + - wtse * dp.get(x+1, y) + + wtse * dp.get(xp1, y) + wtnw * dp.get(x, y+1) + - wtne * dp.get(x+1, y+1); + wtne * dp.get(xp1, y+1); } } @@ -837,10 +840,13 @@ double interp_bilin(const DataPlane &dp, double obs_x, double obs_y, // //////////////////////////////////////////////////////////////////////// -double interp_xy(const DataPlane &dp, int x, int y, +double interp_xy(const DataPlane &dp, bool wrap_lon, int x, int y, const MaskPlane *mp) { double v; + // Apply global wrap logic + x = (wrap_lon ? positive_modulo(x, dp.nx()) : x); + // Check the optional mask if(mp) { if(!(*mp)(x, y)) return(bad_data_double); @@ -932,8 +938,8 @@ double compute_sfc_interp(const DataPlane &dp, double obs_elv, double obs_v, const InterpMthd mthd, const int width, const GridTemplateFactory::GridTemplates shape, - double interp_thresh, const SurfaceInfo &sfc_info, - bool is_land_obs) { + bool wrap_lon, double interp_thresh, + const SurfaceInfo &sfc_info, bool is_land_obs) { double v = bad_data_double; int x = nint(obs_x); int y = nint(obs_y); @@ -945,7 +951,7 @@ double compute_sfc_interp(const DataPlane &dp, } GridTemplateFactory gtf; - const GridTemplate* gt = gtf.buildGT(shape,width); + const GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon); MaskPlane sfc_mask = compute_sfc_mask(*gt, x, y, sfc_info, is_land_obs, obs_elv); @@ -979,11 +985,11 @@ double compute_sfc_interp(const DataPlane &dp, break; case(InterpMthd_Bilin): // Bilinear interpolation - v = interp_bilin(dp, obs_x, obs_y, &sfc_mask); + v = interp_bilin(dp, wrap_lon, obs_x, obs_y, &sfc_mask); break; case(InterpMthd_Nearest): // Nearest Neighbor - v = interp_xy(dp, x, y, &sfc_mask); + v = interp_xy(dp, wrap_lon, x, y, &sfc_mask); break; case(InterpMthd_Best): // Best Match @@ -991,19 +997,19 @@ double compute_sfc_interp(const DataPlane &dp, break; case(InterpMthd_Upper_Left): // Upper Left corner of the grid box - v = interp_xy(dp, floor(obs_x), ceil(obs_y), &sfc_mask); + v = interp_xy(dp, wrap_lon, floor(obs_x), ceil(obs_y), &sfc_mask); break; case(InterpMthd_Upper_Right): // Upper Right corner of the grid box - v = interp_xy(dp, ceil(obs_x), ceil(obs_y), &sfc_mask); + v = interp_xy(dp, wrap_lon, ceil(obs_x), ceil(obs_y), &sfc_mask); break; case(InterpMthd_Lower_Right): // Lower Right corner of the grid box - v = interp_xy(dp, ceil(obs_x), floor(obs_y), &sfc_mask); + v = interp_xy(dp, wrap_lon, ceil(obs_x), floor(obs_y), &sfc_mask); break; case(InterpMthd_Lower_Left): // Lower Left corner of the grid box - v = interp_xy(dp, floor(obs_x), floor(obs_y), &sfc_mask); + v = interp_xy(dp, wrap_lon, floor(obs_x), floor(obs_y), &sfc_mask); break; case(InterpMthd_Geog_Match): // Geography Match for surface point verification @@ -1088,9 +1094,11 @@ double compute_horz_interp(const DataPlane &dp, double obs_x, double obs_y, double obs_v, const InterpMthd mthd, const int width, const GridTemplateFactory::GridTemplates shape, - double interp_thresh, const SingleThresh *cat_thresh) { + bool wrap_lon, double interp_thresh, + const SingleThresh *cat_thresh) { return(compute_horz_interp(dp, obs_x, obs_y, obs_v, bad_data_double, - bad_data_double, mthd, width, shape, interp_thresh, cat_thresh)); + bad_data_double, mthd, width, shape, wrap_lon, + interp_thresh, cat_thresh)); } //////////////////////////////////////////////////////////////////////// @@ -1100,19 +1108,20 @@ double compute_horz_interp(const DataPlane &dp, double obs_v, double cmn, double csd, const InterpMthd mthd, const int width, const GridTemplateFactory::GridTemplates shape, - double interp_thresh, const SingleThresh *cat_thresh) { + bool wrap_lon, double interp_thresh, + const SingleThresh *cat_thresh) { double v = bad_data_double; int x = nint(obs_x); int y = nint(obs_y); // if width is even, push center to lower left point instead of nearest - if((width % 2 ) == 0) { + if((width % 2) == 0) { x = static_cast(floor(obs_x)); y = static_cast(floor(obs_y)); } GridTemplateFactory gtf; - const GridTemplate* gt = gtf.buildGT(shape,width); + const GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon); // Compute the interpolated value for the fields above and below switch(mthd) { @@ -1149,11 +1158,11 @@ double compute_horz_interp(const DataPlane &dp, break; case(InterpMthd_Bilin): // Bilinear interpolation - v = interp_bilin(dp, obs_x, obs_y); + v = interp_bilin(dp, wrap_lon, obs_x, obs_y); break; case(InterpMthd_Nearest): // Nearest Neighbor - v = interp_xy(dp, x, y); + v = interp_xy(dp, wrap_lon, x, y); break; case(InterpMthd_Best): // Best Match @@ -1161,19 +1170,19 @@ double compute_horz_interp(const DataPlane &dp, break; case(InterpMthd_Upper_Left): // Upper Left corner of the grid box - v = interp_xy(dp, floor(obs_x), ceil(obs_y)); + v = interp_xy(dp, wrap_lon, floor(obs_x), ceil(obs_y)); break; case(InterpMthd_Upper_Right): // Upper Right corner of the grid box - v = interp_xy(dp, ceil(obs_x), ceil(obs_y)); + v = interp_xy(dp, wrap_lon, ceil(obs_x), ceil(obs_y)); break; case(InterpMthd_Lower_Right): // Lower Right corner of the grid box - v = interp_xy(dp, ceil(obs_x), floor(obs_y)); + v = interp_xy(dp, wrap_lon, ceil(obs_x), floor(obs_y)); break; case(InterpMthd_Lower_Left): // Lower Left corner of the grid box - v = interp_xy(dp, floor(obs_x), floor(obs_y)); + v = interp_xy(dp, wrap_lon, floor(obs_x), floor(obs_y)); break; case(InterpMthd_Geog_Match): // Geography Match for surface point verification diff --git a/met/src/basic/vx_util/interp_util.h b/met/src/basic/vx_util/interp_util.h index 4369557be4..c8e6dc3f21 100644 --- a/met/src/basic/vx_util/interp_util.h +++ b/met/src/basic/vx_util/interp_util.h @@ -17,6 +17,7 @@ // Mod# Date Name Description // ---- ---- ---- ----------- // 000 11-17-08 Halley Gotway +// 001 09-07-21 Halley Gotway Add wrap_lon. // /////////////////////////////////////////////////////////////////////////////// @@ -84,8 +85,8 @@ extern double interp_geog_match(const DataPlane &, const GridTemplate >, dou extern double interp_nbrhd (const DataPlane &, const GridTemplate >, int x, int y, double t, const SingleThresh *, double cmn, double csd, const MaskPlane *mp = 0); -extern double interp_bilin (const DataPlane &, double obs_x, double obs_y, const MaskPlane *mp = 0); -extern double interp_xy (const DataPlane &, int x, int y, const MaskPlane *mp = 0); +extern double interp_bilin (const DataPlane &, bool wrap_lon, double obs_x, double obs_y, const MaskPlane *mp = 0); +extern double interp_xy (const DataPlane &, bool wrap_lon, int x, int y, const MaskPlane *mp = 0); extern double interp_best (const DataPlane &dp, const GridTemplate >, int x, int y, double obs_v, double t, const MaskPlane *mp = 0); @@ -102,8 +103,8 @@ extern double compute_sfc_interp(const DataPlane &dp, double obs_elv, double obs_v, const InterpMthd mthd, const int width, const GridTemplateFactory::GridTemplates shape, - double interp_thresh, const SurfaceInfo &sfc_info, - bool is_land_obs); + bool wrap_lon, double interp_thresh, + const SurfaceInfo &sfc_info, bool is_land_obs); extern MaskPlane compute_sfc_mask(const GridTemplate >, int x, int y, const SurfaceInfo &sfc_info, @@ -113,14 +114,16 @@ extern double compute_horz_interp(const DataPlane &dp, double obs_x, double obs_y, double obs_v, const InterpMthd mthd, const int width, const GridTemplateFactory::GridTemplates shape, - double interp_thresh, const SingleThresh *cat_thresh = 0); + bool wrap_lon, double interp_thresh, + const SingleThresh *cat_thresh = 0); extern double compute_horz_interp(const DataPlane &dp, double obs_x, double obs_y, double obs_v, double cmn, double csd, const InterpMthd mthd, const int width, const GridTemplateFactory::GridTemplates shape, - double interp_thresh, const SingleThresh *cat_thresh = 0); + bool wrap_lon, double interp_thresh, + const SingleThresh *cat_thresh = 0); extern double compute_vert_pinterp(double, double, double, double, double); extern double compute_vert_zinterp(double, double, double, double, double); diff --git a/met/src/basic/vx_util/num_array.cc b/met/src/basic/vx_util/num_array.cc index 3a6b10e5a1..a792109883 100644 --- a/met/src/basic/vx_util/num_array.cc +++ b/met/src/basic/vx_util/num_array.cc @@ -570,6 +570,28 @@ return; //////////////////////////////////////////////////////////////////////// +void NumArray::set_const(double v, int n) + +{ + +erase(); + +add_const(v, n); + + // + // a constant array is sorted + // + +Sorted = true; + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + void NumArray::sort_array() { diff --git a/met/src/basic/vx_util/num_array.h b/met/src/basic/vx_util/num_array.h index eaea0f7990..d1446bb114 100644 --- a/met/src/basic/vx_util/num_array.h +++ b/met/src/basic/vx_util/num_array.h @@ -83,6 +83,7 @@ class NumArray { void set(double); void set(int, int); void set(int, double); + void set_const(double, int); // Increment value void inc(int, int); diff --git a/met/src/libcode/vx_grid/gaussian_grid.cc b/met/src/libcode/vx_grid/gaussian_grid.cc index 6a175077a4..61edbbb7f8 100644 --- a/met/src/libcode/vx_grid/gaussian_grid.cc +++ b/met/src/libcode/vx_grid/gaussian_grid.cc @@ -431,7 +431,7 @@ return ( 0.0 ); //////////////////////////////////////////////////////////////////////// -bool GaussianGrid::is_global() const +bool GaussianGrid::wrap_lon() const { diff --git a/met/src/libcode/vx_grid/gaussian_grid.h b/met/src/libcode/vx_grid/gaussian_grid.h index 404b546097..1825f853e0 100644 --- a/met/src/libcode/vx_grid/gaussian_grid.h +++ b/met/src/libcode/vx_grid/gaussian_grid.h @@ -81,7 +81,7 @@ class GaussianGrid : public GridRep { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); diff --git a/met/src/libcode/vx_grid/goes_grid.cc b/met/src/libcode/vx_grid/goes_grid.cc index f7cc8e96fa..71e6ff61d5 100644 --- a/met/src/libcode/vx_grid/goes_grid.cc +++ b/met/src/libcode/vx_grid/goes_grid.cc @@ -376,7 +376,7 @@ return ( 0.0 ); //////////////////////////////////////////////////////////////////////// -bool GoesImagerGrid::is_global() const +bool GoesImagerGrid::wrap_lon() const { diff --git a/met/src/libcode/vx_grid/goes_grid.h b/met/src/libcode/vx_grid/goes_grid.h index a43869b1f7..291428fa33 100644 --- a/met/src/libcode/vx_grid/goes_grid.h +++ b/met/src/libcode/vx_grid/goes_grid.h @@ -72,7 +72,7 @@ class GoesImagerGrid : public GridRep { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); diff --git a/met/src/libcode/vx_grid/grid_base.cc b/met/src/libcode/vx_grid/grid_base.cc index 7b1fa99790..ece0332ce1 100644 --- a/met/src/libcode/vx_grid/grid_base.cc +++ b/met/src/libcode/vx_grid/grid_base.cc @@ -957,19 +957,19 @@ return ( rep->rot_grid_to_earth(x, y) ); //////////////////////////////////////////////////////////////////////// -bool Grid::is_global() const +bool Grid::wrap_lon() const { if ( !rep ) { - mlog << Error << "\nGrid::is_global() const -> empty grid!\n\n"; + mlog << Error << "\nGrid::wrap_lon() const -> empty grid!\n\n"; exit ( 1 ); } -return ( rep->is_global() ); +return ( rep->wrap_lon() ); } diff --git a/met/src/libcode/vx_grid/grid_base.h b/met/src/libcode/vx_grid/grid_base.h index c2894f0df1..5a225f6275 100644 --- a/met/src/libcode/vx_grid/grid_base.h +++ b/met/src/libcode/vx_grid/grid_base.h @@ -136,7 +136,7 @@ class GridInterface { // pure abstract class for grid public interface virtual double rot_grid_to_earth(int x, int y) const = 0; - virtual bool is_global() const = 0; + virtual bool wrap_lon() const = 0; virtual void shift_right(int) = 0; @@ -168,7 +168,7 @@ class GridRep : public GridInterface { virtual double rot_grid_to_earth(int x, int y) const = 0; - virtual bool is_global() const = 0; + virtual bool wrap_lon() const = 0; virtual void shift_right(int) = 0; @@ -246,7 +246,7 @@ class Grid : public GridInterface { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); diff --git a/met/src/libcode/vx_grid/latlon_grid.cc b/met/src/libcode/vx_grid/latlon_grid.cc index 6176424ed1..11e779cc66 100644 --- a/met/src/libcode/vx_grid/latlon_grid.cc +++ b/met/src/libcode/vx_grid/latlon_grid.cc @@ -1,5 +1,3 @@ - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) @@ -9,8 +7,6 @@ // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - - //////////////////////////////////////////////////////////////////////// @@ -77,6 +73,8 @@ lat_ll = lon_ll = 0.0; delta_lat = delta_lon = 0.0; +wrapLon = false; + memset(&Data, 0, sizeof(Data)); return; @@ -120,7 +118,29 @@ Name = data.name; Data = data; + // + // wrap longitudes when: + // - all 360 degrees are included + // - 360 is divisible by delta_lon + // +bool lon_rng_flag = fabs((Nx + 1)*delta_lon) >= 360.0; +double lon_div = 360.0 / delta_lon; +bool lon_div_flag = is_eq(lon_div, floor(lon_div)); + +wrapLon = (lon_rng_flag && lon_div_flag); + + // + // inform users about unexpected delta longitudes + // + +if ( lon_rng_flag && !lon_div_flag ) { + + mlog << Debug(4) << "Cannot wrap longitudes since 360 is not " + << "evenly divisible by delta_lon = " + << delta_lon << ".\n"; + +} return; @@ -133,6 +153,7 @@ return; void LatLonGrid::latlon_to_xy(double lat, double lon, double &x, double &y) const { + double n; y = (lat - lat_ll)/delta_lat; @@ -249,6 +270,8 @@ out << prefix << "lon_ll = " << lon_ll << "\n"; out << prefix << "delta_lat_ll = " << delta_lat << "\n"; out << prefix << "delta_lon_ll = " << delta_lon << "\n"; +out << prefix << "wrapLon = " << bool_to_string(wrapLon) << "\n"; + // // done // @@ -282,6 +305,8 @@ snprintf(junk, sizeof(junk), " lon_ll: %.3f", lon_ll); a << junk; snprintf(junk, sizeof(junk), " delta_lat: %.3f", delta_lat); a << junk; snprintf(junk, sizeof(junk), " delta_lon: %.3f", delta_lon); a << junk; +snprintf(junk, sizeof(junk), " wrapLon: %s", bool_to_string(wrapLon)); a << junk; + // // done // @@ -328,33 +353,13 @@ return ( 0.0 ); //////////////////////////////////////////////////////////////////////// -bool LatLonGrid::is_global() const - -{ - -const double lon_range = fabs((Nx + 1)*delta_lon); -const double lat_range = fabs((Ny + 1)*delta_lat); - -const bool full_range_lat = (lat_range >= 180.0); -const bool full_range_lon = (lon_range >= 360.0); - -const bool answer = full_range_lat && full_range_lon; - -return ( answer ); - -} - - -//////////////////////////////////////////////////////////////////////// - - void LatLonGrid::shift_right(int N) { if ( N == 0 ) return; -if ( ! is_global() ) { +if ( ! wrapLon ) { mlog << Error << "\n\n LatLonGrid::shift_right(int) -> " @@ -450,5 +455,3 @@ return; //////////////////////////////////////////////////////////////////////// - - diff --git a/met/src/libcode/vx_grid/latlon_grid.h b/met/src/libcode/vx_grid/latlon_grid.h index 4f0f01cb5b..21c86d6e2c 100644 --- a/met/src/libcode/vx_grid/latlon_grid.h +++ b/met/src/libcode/vx_grid/latlon_grid.h @@ -1,5 +1,3 @@ - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) @@ -9,8 +7,6 @@ // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - - //////////////////////////////////////////////////////////////////////// @@ -46,6 +42,8 @@ class LatLonGrid : public GridRep { int Nx; int Ny; + bool wrapLon; + ConcatString Name; LatLonData Data; @@ -79,7 +77,7 @@ class LatLonGrid : public GridRep { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); @@ -91,7 +89,8 @@ class LatLonGrid : public GridRep { //////////////////////////////////////////////////////////////////////// -inline double LatLonGrid::scale_km() const { return ( -1.0 ); } +inline double LatLonGrid::scale_km() const { return ( -1.0 ); } +inline bool LatLonGrid::wrap_lon() const { return ( wrapLon ); } //////////////////////////////////////////////////////////////////////// @@ -101,6 +100,3 @@ inline double LatLonGrid::scale_km() const { return ( -1.0 ); } //////////////////////////////////////////////////////////////////////// - - - diff --git a/met/src/libcode/vx_grid/lc_grid.cc b/met/src/libcode/vx_grid/lc_grid.cc index aab40b015f..41a2e375e9 100644 --- a/met/src/libcode/vx_grid/lc_grid.cc +++ b/met/src/libcode/vx_grid/lc_grid.cc @@ -638,7 +638,7 @@ return ( angle ); //////////////////////////////////////////////////////////////////////// -bool LambertGrid::is_global() const +bool LambertGrid::wrap_lon() const { diff --git a/met/src/libcode/vx_grid/lc_grid.h b/met/src/libcode/vx_grid/lc_grid.h index 6b67eabb0b..0b3fefdaaa 100644 --- a/met/src/libcode/vx_grid/lc_grid.h +++ b/met/src/libcode/vx_grid/lc_grid.h @@ -118,7 +118,7 @@ class LambertGrid : public GridRep { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); diff --git a/met/src/libcode/vx_grid/merc_grid.cc b/met/src/libcode/vx_grid/merc_grid.cc index c45dee4d6d..d89a7898e4 100644 --- a/met/src/libcode/vx_grid/merc_grid.cc +++ b/met/src/libcode/vx_grid/merc_grid.cc @@ -572,7 +572,7 @@ return ( 0.0 ); //////////////////////////////////////////////////////////////////////// -bool MercatorGrid::is_global() const +bool MercatorGrid::wrap_lon() const { diff --git a/met/src/libcode/vx_grid/merc_grid.h b/met/src/libcode/vx_grid/merc_grid.h index 94dd842bda..f46a9c3bf8 100644 --- a/met/src/libcode/vx_grid/merc_grid.h +++ b/met/src/libcode/vx_grid/merc_grid.h @@ -101,7 +101,7 @@ class MercatorGrid : public GridRep { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); diff --git a/met/src/libcode/vx_grid/rot_latlon_grid.cc b/met/src/libcode/vx_grid/rot_latlon_grid.cc index 2e97cba498..12b6286135 100644 --- a/met/src/libcode/vx_grid/rot_latlon_grid.cc +++ b/met/src/libcode/vx_grid/rot_latlon_grid.cc @@ -1,5 +1,3 @@ - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) @@ -321,6 +319,8 @@ snprintf(junk, sizeof(junk), " rot_lon_ll: %.3f", RData.rot_lon_ll); a << junk snprintf(junk, sizeof(junk), " delta_rot_lat: %.3f", RData.delta_rot_lat); a << junk; snprintf(junk, sizeof(junk), " delta_rot_lon: %.3f", RData.delta_rot_lon); a << junk; +snprintf(junk, sizeof(junk), " wrapLon: %s", bool_to_string(wrapLon)); a << junk; + snprintf(junk, sizeof(junk), " true_lat_south_pole: %.3f", RData.true_lat_south_pole); a << junk; snprintf(junk, sizeof(junk), " true_lon_south_pole: %.3f", RData.true_lon_south_pole); a << junk; @@ -372,26 +372,6 @@ return ( 0.0 ); //////////////////////////////////////////////////////////////////////// -bool RotatedLatLonGrid::is_global() const - -{ - -const double lon_range = fabs((Nx + 1)*delta_lon); -const double lat_range = fabs((Ny + 1)*delta_lat); - -const bool full_range_lat = (lat_range >= 180.0); -const bool full_range_lon = (lon_range >= 360.0); - -const bool answer = full_range_lat && full_range_lon; - -return ( answer ); - -} - - -//////////////////////////////////////////////////////////////////////// - - void RotatedLatLonGrid::shift_right(int N) { @@ -472,5 +452,3 @@ return; //////////////////////////////////////////////////////////////////////// - - diff --git a/met/src/libcode/vx_grid/rot_latlon_grid.h b/met/src/libcode/vx_grid/rot_latlon_grid.h index eea580d285..79a75e1869 100644 --- a/met/src/libcode/vx_grid/rot_latlon_grid.h +++ b/met/src/libcode/vx_grid/rot_latlon_grid.h @@ -1,5 +1,3 @@ - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) @@ -72,7 +70,7 @@ class RotatedLatLonGrid : public LatLonGrid { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); @@ -84,7 +82,8 @@ class RotatedLatLonGrid : public LatLonGrid { //////////////////////////////////////////////////////////////////////// -inline double RotatedLatLonGrid::scale_km() const { return ( -1.0 ); } +inline double RotatedLatLonGrid::scale_km() const { return ( -1.0 ); } +inline bool RotatedLatLonGrid::wrap_lon() const { return ( wrapLon ); } //////////////////////////////////////////////////////////////////////// @@ -94,6 +93,3 @@ inline double RotatedLatLonGrid::scale_km() const { return ( -1.0 ); } //////////////////////////////////////////////////////////////////////// - - - diff --git a/met/src/libcode/vx_grid/st_grid.cc b/met/src/libcode/vx_grid/st_grid.cc index 65246ae5fc..3b64e4a14d 100644 --- a/met/src/libcode/vx_grid/st_grid.cc +++ b/met/src/libcode/vx_grid/st_grid.cc @@ -553,7 +553,7 @@ return ( angle ); //////////////////////////////////////////////////////////////////////// -bool StereographicGrid::is_global() const +bool StereographicGrid::wrap_lon() const { diff --git a/met/src/libcode/vx_grid/st_grid.h b/met/src/libcode/vx_grid/st_grid.h index e80f3aadff..9fecdba445 100644 --- a/met/src/libcode/vx_grid/st_grid.h +++ b/met/src/libcode/vx_grid/st_grid.h @@ -99,7 +99,7 @@ class StereographicGrid : public GridRep { double rot_grid_to_earth(int x, int y) const; - bool is_global() const; + bool wrap_lon() const; void shift_right(int); diff --git a/met/src/libcode/vx_regrid/vx_regrid.cc b/met/src/libcode/vx_regrid/vx_regrid.cc index 0a7da68e3e..6be4189fa7 100644 --- a/met/src/libcode/vx_regrid/vx_regrid.cc +++ b/met/src/libcode/vx_regrid/vx_regrid.cc @@ -133,6 +133,7 @@ to_data.set_accum (from_data.accum()); // // copy data // + for (xt=0; xt<(to_grid.nx()); ++xt) { for (yt=0; yt<(to_grid.ny()); ++yt) { @@ -144,14 +145,14 @@ for (xt=0; xt<(to_grid.nx()); ++xt) { xf = nint(x_from); yf = nint(y_from); - if ( (xf < 0) || (xf >= from_grid.nx()) || (yf < 0) || (yf >= from_grid.ny()) ) { - + if ( ( (xf < 0 || xf >= from_grid.nx()) && !from_grid.wrap_lon() ) || + yf < 0 || yf >= from_grid.ny() ) { value = bad_data_float; - - } else { - value = compute_horz_interp(from_data, x_from, y_from, bad_data_double, - info.method, info.width, info.shape, info.vld_thresh); - + } + else { + value = compute_horz_interp(from_data, x_from, y_from, + bad_data_double, info.method, info.width, + info.shape, from_grid.wrap_lon(), info.vld_thresh); } to_data.put(value, xt, yt); @@ -327,6 +328,7 @@ to_data.set_accum (from_data.accum()); // // copy data // + for (xt=0; xt<(to_grid.nx()); ++xt) { for (yt=0; yt<(to_grid.ny()); ++yt) { @@ -338,15 +340,14 @@ for (xt=0; xt<(to_grid.nx()); ++xt) { xf = nint(x_from); yf = nint(y_from); - if ( (xf < 0) || (xf >= from_grid.nx()) || - (yf < 0) || (yf >= from_grid.ny()) ) { - + if ( ( (xf < 0 || xf >= from_grid.nx()) && !from_grid.wrap_lon() ) || + yf < 0 || yf >= from_grid.ny() ) { value = bad_data_float; } else { value = compute_horz_interp(from_data, x_from, y_from, bad_data_double, InterpMthd_Max, info.width, - info.shape, info.vld_thresh); + info.shape, from_grid.wrap_lon(), info.vld_thresh); } to_data.put(value, xt, yt); diff --git a/met/src/libcode/vx_regrid/vx_regrid_budget.cc b/met/src/libcode/vx_regrid/vx_regrid_budget.cc index 881e56ef1c..42157cb729 100644 --- a/met/src/libcode/vx_regrid/vx_regrid_budget.cc +++ b/met/src/libcode/vx_regrid/vx_regrid_budget.cc @@ -77,7 +77,7 @@ for (ixt=0; ixt<(to_grid.nx()); ++ixt) { from_grid.latlon_to_xy(lat, lon, dxf, dyf); - value = interp_bilin(from_data, dxf, dyf); + value = interp_bilin(from_data, from_grid.wrap_lon(), dxf, dyf); if ( value != bad_data_double ) { sum += value; ++count; } diff --git a/met/src/libcode/vx_statistics/apply_mask.cc b/met/src/libcode/vx_statistics/apply_mask.cc index 1aea2fe836..01d88f601e 100644 --- a/met/src/libcode/vx_statistics/apply_mask.cc +++ b/met/src/libcode/vx_statistics/apply_mask.cc @@ -185,7 +185,7 @@ void parse_grid_mask(const ConcatString &mask_grid_str, const Grid &grid, // // Check to make sure that we're not using the full domain // - if( full_domain_str != mask_grid_str) { + if(full_domain_str != mask_grid_str) { // // Search for the grid name in the predefined grids diff --git a/met/src/libcode/vx_statistics/pair_base.cc b/met/src/libcode/vx_statistics/pair_base.cc index d90e3167a6..aaf73755b9 100644 --- a/met/src/libcode/vx_statistics/pair_base.cc +++ b/met/src/libcode/vx_statistics/pair_base.cc @@ -876,6 +876,7 @@ double compute_interp(const DataPlaneArray &dpa, const double obs_v, const double cmn, const double csd, const InterpMthd method, const int width, const GridTemplateFactory::GridTemplates shape, + const bool wrap_lon, const double thresh, const bool spfh_flag, const LevelType lvl_typ, const double to_lvl, const int i_blw, const int i_abv, @@ -886,14 +887,16 @@ double compute_interp(const DataPlaneArray &dpa, if(dpa.n_planes() == 0) return(bad_data_double); v_blw = compute_horz_interp(dpa[i_blw], obs_x, obs_y, obs_v, cmn, csd, - method, width, shape, thresh, cat_thresh); + method, width, shape, wrap_lon, + thresh, cat_thresh); if(i_blw == i_abv) { v = v_blw; } else { v_abv = compute_horz_interp(dpa[i_abv], obs_x, obs_y, obs_v, cmn, csd, - method, width, shape, thresh, cat_thresh); + method, width, shape, wrap_lon, + thresh, cat_thresh); // Check for bad data prior to vertical interpolation if(is_bad_data(v_blw) || is_bad_data(v_abv)) { @@ -932,6 +935,7 @@ void get_interp_points(const DataPlaneArray &dpa, const double obs_x, const double obs_y, const InterpMthd method, const int width, const GridTemplateFactory::GridTemplates shape, + const bool wrap_lon, const double thresh, const bool spfh_flag, const LevelType lvl_typ, const double to_lvl, const int i_blw, const int i_abv, @@ -947,7 +951,7 @@ void get_interp_points(const DataPlaneArray &dpa, int i, n_vld; NumArray pts_blw, pts_abv; GridTemplateFactory gtf; - const GridTemplate* gt = gtf.buildGT(shape, width); + const GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon); // Get interpolation points below the observation pts_blw = interp_points(dpa[i_blw], *gt, obs_x, obs_y); diff --git a/met/src/libcode/vx_statistics/pair_base.h b/met/src/libcode/vx_statistics/pair_base.h index ebe0f845fc..465432f16f 100644 --- a/met/src/libcode/vx_statistics/pair_base.h +++ b/met/src/libcode/vx_statistics/pair_base.h @@ -197,6 +197,7 @@ extern double compute_interp(const DataPlaneArray &dpa, const double obs_v, const double cmn, const double csd, const InterpMthd method, const int width, const GridTemplateFactory::GridTemplates shape, + const bool wrap_lon, const double thresh, const bool spfh_flag, const LevelType lvl_typ, const double to_lvl, const int i_blw, const int i_abv, @@ -206,6 +207,7 @@ extern void get_interp_points(const DataPlaneArray &dpa, const double obs_x, const double obs_y, const InterpMthd method, const int width, const GridTemplateFactory::GridTemplates shape, + const bool wrap_lon, const double thresh, const bool spfh_flag, const LevelType lvl_typ, const double to_lvl, const int i_blw, const int i_abv, diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.cc b/met/src/libcode/vx_statistics/pair_data_ensemble.cc index 0fcf84fc99..683aa499b9 100644 --- a/met/src/libcode/vx_statistics/pair_data_ensemble.cc +++ b/met/src/libcode/vx_statistics/pair_data_ensemble.cc @@ -1408,8 +1408,8 @@ void VxPairDataEnsemble::add_point_obs(float *hdr_arr, int *hdr_typ_arr, y = nint(obs_y); // Check if the observation's lat/lon is on the grid - if(x < 0 || x >= gr.nx() || - y < 0 || y >= gr.ny()) return; + if(((x < 0 || x >= gr.nx()) && !gr.wrap_lon()) || + y < 0 || y >= gr.ny()) return; // For pressure levels, check if the observation pressure level // falls in the requsted range. @@ -1556,7 +1556,7 @@ void VxPairDataEnsemble::add_point_obs(float *hdr_arr, int *hdr_typ_arr, cmn_v = compute_interp(climo_mn_dpa, obs_x, obs_y, obs_v, bad_data_double, bad_data_double, pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, + pd[0][0][k].interp_shape, gr.wrap_lon(), interp_thresh, spfh_flag, fcst_info->level().type(), to_lvl, cmn_lvl_blw, cmn_lvl_abv); @@ -1583,7 +1583,7 @@ void VxPairDataEnsemble::add_point_obs(float *hdr_arr, int *hdr_typ_arr, csd_v = compute_interp(climo_sd_dpa, obs_x, obs_y, obs_v, bad_data_double, bad_data_double, pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, + pd[0][0][k].interp_shape, gr.wrap_lon(), interp_thresh, spfh_flag, fcst_info->level().type(), to_lvl, csd_lvl_blw, csd_lvl_abv); @@ -1612,7 +1612,7 @@ void VxPairDataEnsemble::add_point_obs(float *hdr_arr, int *hdr_typ_arr, //////////////////////////////////////////////////////////////////////// -void VxPairDataEnsemble::add_ens(int member, bool mn) { +void VxPairDataEnsemble::add_ens(int member, bool mn, Grid &gr) { int i, j, k, l; int f_lvl_blw, f_lvl_abv; double to_lvl, fcst_v; @@ -1654,6 +1654,7 @@ void VxPairDataEnsemble::add_ens(int member, bool mn) { pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, pd[0][0][k].interp_shape, + gr.wrap_lon(), interp_thresh, spfh_flag, fcst_info->level().type(), to_lvl, f_lvl_blw, f_lvl_abv); diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.h b/met/src/libcode/vx_statistics/pair_data_ensemble.h index 84aa45796b..7d59896593 100644 --- a/met/src/libcode/vx_statistics/pair_data_ensemble.h +++ b/met/src/libcode/vx_statistics/pair_data_ensemble.h @@ -279,7 +279,7 @@ class VxPairDataEnsemble { unixtime, const char *, float *, Grid &, const char * = 0, const DataPlane * = 0); - void add_ens(int, bool mn); + void add_ens(int, bool mn, Grid &); int get_n_pair() const; diff --git a/met/src/libcode/vx_statistics/pair_data_point.cc b/met/src/libcode/vx_statistics/pair_data_point.cc index a77749abba..7bfda5fa42 100644 --- a/met/src/libcode/vx_statistics/pair_data_point.cc +++ b/met/src/libcode/vx_statistics/pair_data_point.cc @@ -967,8 +967,9 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, y = nint(obs_y); // Check if the observation's lat/lon is on the grid - if(x < 0 || x >= gr.nx() || - y < 0 || y >= gr.ny()) { + if(((x < 0 || x >= gr.nx()) && !gr.wrap_lon()) || + y < 0 || y >= gr.ny()) { + mlog << Debug(4) << "For " << fcst_info->magic_str() << " versus " << obs_info->magic_str() @@ -989,7 +990,8 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, double topo = compute_horz_interp( *sfc_info.topo_ptr, obs_x, obs_y, hdr_elv, InterpMthd_Bilin, 2, - GridTemplateFactory::GridTemplate_Square, 1.0); + GridTemplateFactory::GridTemplate_Square, + gr.wrap_lon(), 1.0); // Skip bad topography values if(is_bad_data(hdr_elv) || is_bad_data(topo)) { @@ -1165,7 +1167,7 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, cmn_v = compute_interp(climo_mn_dpa, obs_x, obs_y, obs_v, bad_data_double, bad_data_double, pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, + pd[0][0][k].interp_shape, gr.wrap_lon(), interp_thresh, spfh_flag, fcst_info->level().type(), to_lvl, cmn_lvl_blw, cmn_lvl_abv); @@ -1192,8 +1194,8 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, // Compute the interpolated climatology standard deviation csd_v = compute_interp(climo_sd_dpa, obs_x, obs_y, obs_v, bad_data_double, bad_data_double, - pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, + pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, + pd[0][0][k].interp_shape, gr.wrap_lon(), interp_thresh, spfh_flag, fcst_info->level().type(), to_lvl, csd_lvl_blw, csd_lvl_abv); @@ -1222,14 +1224,14 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, fcst_v = compute_sfc_interp(fcst_dpa[0], obs_x, obs_y, hdr_elv, obs_v, pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, interp_thresh, sfc_info, - is_land); + pd[0][0][k].interp_shape, gr.wrap_lon(), + interp_thresh, sfc_info, is_land); } // Otherwise, compute interpolated value else { fcst_v = compute_interp(fcst_dpa, obs_x, obs_y, obs_v, cmn_v, csd_v, pd[0][0][k].interp_mthd, pd[0][0][k].interp_wdth, - pd[0][0][k].interp_shape, + pd[0][0][k].interp_shape, gr.wrap_lon(), interp_thresh, spfh_flag, fcst_info->level().type(), to_lvl, f_lvl_blw, f_lvl_abv); @@ -1239,8 +1241,10 @@ void VxPairDataPoint::add_point_obs(float *hdr_arr, const char *hdr_typ_str, mlog << Debug(4) << "For " << fcst_info->magic_str() << " versus " << obs_info->magic_str() - << ", skipping observation due to bad data in the interpolated " - << "forecast value:\n" + << ", skipping observation due to bad data in the " + << interpmthd_to_string(pd[0][0][k].interp_mthd) << "(" + << pd[0][0][k].interp_wdth * pd[0][0][k].interp_wdth + << ") interpolated forecast value:\n" << point_obs_to_string(hdr_arr, hdr_typ_str, hdr_sid_str, hdr_ut, obs_qty, obs_arr, var_name) << "\n"; diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc index d87c377bdf..6dfc98dee3 100644 --- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc +++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc @@ -1123,7 +1123,7 @@ int process_point_ens(int i_ens, int &n_miss) { conf_info.vx_opt[i].vx_pd.set_fcst_dpa(fcst_dpa); // Compute forecast values for this ensemble member - conf_info.vx_opt[i].vx_pd.add_ens(i_ens-n_miss, is_ens_mean); + conf_info.vx_opt[i].vx_pd.add_ens(i_ens-n_miss, is_ens_mean, grid); } // end for i @@ -1588,13 +1588,13 @@ void process_grid_vx() { // Smooth the ensemble mean field, if requested if(ens_mean_flag && (field == FieldType_Fcst || field == FieldType_Both)) { - emn_dp = smooth_field(emn_dp, mthd, wdth, shape, + emn_dp = smooth_field(emn_dp, mthd, wdth, shape, grid.wrap_lon(), vld_thresh, gaussian); } // Smooth the observation field, if requested if(field == FieldType_Obs || field == FieldType_Both) { - obs_dp = smooth_field(obs_dp, mthd, wdth, shape, + obs_dp = smooth_field(obs_dp, mthd, wdth, shape, grid.wrap_lon(), vld_thresh, gaussian); } @@ -1617,7 +1617,7 @@ void process_grid_vx() { // Smooth the forecast field, if requested if(field == FieldType_Fcst || field == FieldType_Both) { - fcst_dp[k] = smooth_field(fcst_dp[k], mthd, wdth, shape, + fcst_dp[k] = smooth_field(fcst_dp[k], mthd, wdth, shape, grid.wrap_lon(), vld_thresh, gaussian); } @@ -2060,7 +2060,9 @@ void track_counts(int i_vx, const DataPlane &dp) { fractional_coverage(dp, frac_dp, conf_info.nbrhd_prob.width[j], conf_info.nbrhd_prob.shape, - ThreshBuf[i], conf_info.nbrhd_prob.vld_thresh); + grid.wrap_lon(), ThreshBuf[i], + (const DataPlane *) 0, (const DataPlane *) 0, + conf_info.nbrhd_prob.vld_thresh); // Increment counts const double *Frac = frac_dp.data(); @@ -2496,6 +2498,7 @@ void write_ens_nc(int i_ens, DataPlane &dp) { smooth_dp = smooth_field(prob_dp, InterpMthd_UW_Mean, conf_info.nbrhd_prob.width[j], conf_info.nbrhd_prob.shape, + grid.wrap_lon(), conf_info.nbrhd_prob.vld_thresh, info); for(k=0; kshape, interp->width[j]); + GridTemplate* gt = gtf.buildGT(interp->shape, interp->width[j], grid.wrap_lon()); shc.set_interp_mthd(interp_mthd, interp->shape); int interp_pnts = gt->size(); @@ -769,8 +769,8 @@ void process_scores() { if(interp->field == FieldType_Fcst || interp->field == FieldType_Both) { smooth_field(fcst_dp, fcst_dp_smooth, interp_mthd, - interp->width[j], interp->shape, interp->vld_thresh, - interp->gaussian); + interp->width[j], interp->shape, grid.wrap_lon(), + interp->vld_thresh, interp->gaussian); } // Do not smooth the forecast field else { @@ -781,8 +781,8 @@ void process_scores() { if(interp->field == FieldType_Obs || interp->field == FieldType_Both) { smooth_field(obs_dp, obs_dp_smooth, interp_mthd, - interp->width[j], interp->shape, interp->vld_thresh, - interp->gaussian); + interp->width[j], interp->shape, grid.wrap_lon(), + interp->vld_thresh, interp->gaussian); } // Do not smooth the observation field else { @@ -970,8 +970,8 @@ void process_scores() { if(interp->field == FieldType_Fcst || interp->field == FieldType_Both) { smooth_field(fu_dp, fu_dp_smooth, interp_mthd, - interp->width[j], interp->shape, interp->vld_thresh, - interp->gaussian); + interp->width[j], interp->shape, grid.wrap_lon(), + interp->vld_thresh, interp->gaussian); } // Do not smooth the forecast field else { @@ -984,8 +984,8 @@ void process_scores() { interp->field == FieldType_Both) { smooth_field(ou_dp, ou_dp_smooth, interp_mthd, interp->width[j], - interp->shape, interp->vld_thresh, - interp->gaussian); + interp->shape, grid.wrap_lon(), + interp->vld_thresh, interp->gaussian); } // Do not smooth the observation field else { @@ -1392,7 +1392,9 @@ void process_scores() { // Compute fractional coverage fractional_coverage(fcst_dp, fcst_dp_smooth, nbrhd->width[j], nbrhd->shape, + grid.wrap_lon(), conf_info.vx_opt[i].fcat_ta[k], + &cmn_dp, &csd_dp, nbrhd->vld_thresh); // Compute the binary threshold field @@ -1430,7 +1432,9 @@ void process_scores() { // Compute fractional coverage fractional_coverage(obs_dp, obs_dp_smooth, nbrhd->width[j], nbrhd->shape, + grid.wrap_lon(), conf_info.vx_opt[i].ocat_ta[k], + &cmn_dp, &csd_dp, nbrhd->vld_thresh); // Compute the binary threshold field diff --git a/met/src/tools/core/point_stat/point_stat.cc b/met/src/tools/core/point_stat/point_stat.cc index 5f609feeff..549e6fce9f 100644 --- a/met/src/tools/core/point_stat/point_stat.cc +++ b/met/src/tools/core/point_stat/point_stat.cc @@ -1676,7 +1676,8 @@ void do_hira_ens(int i_vx, const PairDataPoint *pd_ptr) { // Determine the number of points in the area GridTemplateFactory gtf; GridTemplate* gt = gtf.buildGT(conf_info.vx_opt[i_vx].hira_info.shape, - conf_info.vx_opt[i_vx].hira_info.width[i]); + conf_info.vx_opt[i_vx].hira_info.width[i], + grid.wrap_lon()); // Initialize hira_pd.clear(); @@ -1696,7 +1697,7 @@ void do_hira_ens(int i_vx, const PairDataPoint *pd_ptr) { get_interp_points(conf_info.vx_opt[i_vx].vx_pd.fcst_dpa, pd_ptr->x_na[j], pd_ptr->y_na[j], InterpMthd_Nbrhd, conf_info.vx_opt[i_vx].hira_info.width[i], - conf_info.vx_opt[i_vx].hira_info.shape, + conf_info.vx_opt[i_vx].hira_info.shape, grid.wrap_lon(), conf_info.vx_opt[i_vx].hira_info.vld_thresh, spfh_flag, conf_info.vx_opt[i_vx].vx_pd.fcst_info->level().type(), pd_ptr->lvl_na[j], lvl_blw, lvl_abv, f_ens); @@ -1873,7 +1874,7 @@ void do_hira_prob(int i_vx, const PairDataPoint *pd_ptr) { pd_ptr->x_na[k], pd_ptr->y_na[k], pd_ptr->o_na[k], pd_ptr->cmn_na[k], pd_ptr->csd_na[k], InterpMthd_Nbrhd, conf_info.vx_opt[i_vx].hira_info.width[j], - conf_info.vx_opt[i_vx].hira_info.shape, + conf_info.vx_opt[i_vx].hira_info.shape, grid.wrap_lon(), conf_info.vx_opt[i_vx].hira_info.vld_thresh, spfh_flag, conf_info.vx_opt[i_vx].vx_pd.fcst_info->level().type(), pd_ptr->lvl_na[k], lvl_blw, lvl_abv, &cat_thresh); @@ -1892,7 +1893,7 @@ void do_hira_prob(int i_vx, const PairDataPoint *pd_ptr) { pd_ptr->x_na[k], pd_ptr->y_na[k], pd_ptr->o_na[k], pd_ptr->cmn_na[k], pd_ptr->csd_na[k], InterpMthd_Nbrhd, conf_info.vx_opt[i_vx].hira_info.width[j], - conf_info.vx_opt[i_vx].hira_info.shape, + conf_info.vx_opt[i_vx].hira_info.shape, grid.wrap_lon(), conf_info.vx_opt[i_vx].hira_info.vld_thresh, spfh_flag, conf_info.vx_opt[i_vx].vx_pd.fcst_info->level().type(), pd_ptr->lvl_na[k], lvl_blw, lvl_abv, &cat_thresh); diff --git a/met/src/tools/core/point_stat/point_stat_conf_info.cc b/met/src/tools/core/point_stat/point_stat_conf_info.cc index 87183cb5d0..6117a0f7fc 100644 --- a/met/src/tools/core/point_stat/point_stat_conf_info.cc +++ b/met/src/tools/core/point_stat/point_stat_conf_info.cc @@ -1051,6 +1051,7 @@ void PointStatVxOpt::set_vx_pd(PointStatConfInfo *conf_info) { for(i=0; i 0) { - + if(default_job.job_type != no_stat_job_type) { + process_job(command_line_job_options.c_str(), 1); + } + // + // Or process config file jobs. + // + else if(config_file.length() > 0) { jobs_sa = conf.lookup_string_array(conf_key_jobs); + // + // At least one job in the config file. + // + if(jobs_sa.n() == 0) { + mlog << Error << "\nmain() -> " + << "no jobs defined in \"" << config_file << "\"!\n\n"; + throw(1); + } + for(i=0; i " + << "at least one job must be specified on the command line " + << "with \"-job\" or in a configuration file with \"-config\"!\n\n"; + throw(1); } - } - catch(int j) { // Catch an error + catch(int j) { // Catch errors mlog << Error << "\nmain() -> " << "encountered an error value of " << j diff --git a/met/src/tools/other/Makefile.am b/met/src/tools/other/Makefile.am index 3f56c89a18..47c8da59ce 100644 --- a/met/src/tools/other/Makefile.am +++ b/met/src/tools/other/Makefile.am @@ -22,6 +22,10 @@ if ENABLE_GIS_UTILS SUBDIRS += gis_utils endif +if ENABLE_GEN_ENS_PROD + SUBDIRS += gen_ens_prod +endif + if ENABLE_GEN_VX_MASK SUBDIRS += gen_vx_mask endif diff --git a/met/src/tools/other/gen_ens_prod/.gitignore b/met/src/tools/other/gen_ens_prod/.gitignore new file mode 100644 index 0000000000..3ad703829a --- /dev/null +++ b/met/src/tools/other/gen_ens_prod/.gitignore @@ -0,0 +1,7 @@ +gen_ens_prod +*.o +*.a +.deps +Makefile +Makefile.in +*.dSYM diff --git a/met/src/tools/other/gen_ens_prod/Makefile.am b/met/src/tools/other/gen_ens_prod/Makefile.am new file mode 100644 index 0000000000..5c1191f123 --- /dev/null +++ b/met/src/tools/other/gen_ens_prod/Makefile.am @@ -0,0 +1,44 @@ +## @start 1 +## Makefile.am -- Process this file with automake to produce Makefile.in +## @end 1 + +MAINTAINERCLEANFILES = Makefile.in + +# Include the project definitions + +include ${top_srcdir}/Make-include + +# The program + +bin_PROGRAMS = gen_ens_prod +gen_ens_prod_SOURCES = gen_ens_prod.cc \ + gen_ens_prod_conf_info.cc +gen_ens_prod_CPPFLAGS = ${MET_CPPFLAGS} +gen_ens_prod_LDFLAGS = ${MET_LDFLAGS} +gen_ens_prod_LDADD = -lvx_stat_out \ + -lvx_statistics \ + -lvx_shapedata \ + -lvx_gsl_prob \ + -lvx_analysis_util \ + -lvx_data2d_factory \ + -lvx_data2d_nc_met \ + -lvx_data2d_grib $(GRIB2_LIBS) \ + -lvx_data2d_nc_pinterp \ + -lvx_data2d_nccf \ + $(PYTHON_LIBS) \ + -lvx_data2d \ + -lvx_nc_obs \ + -lvx_nc_util \ + -lvx_regrid \ + -lvx_grid \ + -lvx_config \ + -lvx_color \ + -lvx_util \ + -lvx_math \ + -lvx_cal \ + -lvx_log \ + $(PYTHON_LIBS) \ + -lm -lnetcdf_c++4 -lnetcdf -lgsl -lgslcblas + +EXTRA_DIST = gen_ens_prod.h \ + gen_ens_prod_conf_info.h diff --git a/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc b/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc new file mode 100644 index 0000000000..ce956467af --- /dev/null +++ b/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc @@ -0,0 +1,1053 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// +// +// Filename: gen_ens_prod.cc +// +// Description: +// +// Mod# Date Name Description +// ---- ---- ---- ----------- +// 000 09/10/21 Halley Gotway Initial version (MET #1904). +// +//////////////////////////////////////////////////////////////////////// + +using namespace std; + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "gen_ens_prod.h" + +#include "vx_nc_util.h" +#include "vx_data2d_nc_met.h" +#include "vx_regrid.h" +#include "vx_log.h" + +#include "nc_obs_util.h" +#include "nc_point_obs_in.h" + +//////////////////////////////////////////////////////////////////////// + +static void process_command_line(int, char **); +static void process_grid(const Grid &); +static void process_ensemble(); +static bool get_data_plane(const char *, GrdFileType, VarInfo *, DataPlane &); + +static void clear_counts(); +static void track_counts(int, const DataPlane &, bool, + const DataPlane &, const DataPlane &); + +static void setup_nc_file(); +static void write_ens_nc(int, int, const DataPlane &, + const DataPlane &, const DataPlane &); +static void write_ens_var_float(int, float *, const DataPlane &, + const char *, const char *); +static void write_ens_var_int(int, int *, const DataPlane &, + const char *, const char *); +static void write_ens_data_plane(int, const DataPlane &, const DataPlane &, + const char *, const char *); + +static void add_var_att_local(VarInfo *, NcVar *, bool is_int, + const DataPlane &, const char *, const char *); + +static void clean_up(); +static void usage(); + +static void set_ens_files (const StringArray &); +static void set_out_file (const StringArray &); +static void set_config_file(const StringArray &); +static void set_ctrl_file (const StringArray &); + +//////////////////////////////////////////////////////////////////////// + +int main(int argc, char *argv[]) { + + // Set handler to be called for memory allocation error + set_new_handler(oom); + + // Process the command line arguments + process_command_line(argc, argv); + + // Process the ensemble fields + process_ensemble(); + + // Close output files and deallocate memory + clean_up(); + + return(0); +} + +//////////////////////////////////////////////////////////////////////// + +void process_command_line(int argc, char **argv) { + int i; + CommandLine cline; + ConcatString default_config_file; + Met2dDataFile *ens_mtddf = (Met2dDataFile *) 0; + + // + // Check for zero arguments + // + if(argc == 1) usage(); + + // + // Parse the command line into tokens + // + cline.set(argc, argv); + + // + // Set the usage function + // + cline.set_usage(usage); + + // + // Add the options function calls + // + cline.add(set_ens_files, "-ens", -1); + cline.add(set_out_file, "-out", 1); + cline.add(set_config_file, "-config", 1); + cline.add(set_ctrl_file, "-ctrl", 1); + + // + // Parse the command line + // + cline.parse(); + + // Check for error, there should be zero arguments left + if(cline.n() != 0) usage(); + + // Check that the required arguments have been set + n_ens = ens_files.n(); + if(ens_files.n() == 0) { + mlog << Error << "\nprocess_command_line() -> " + << "the ensemble file list must be set using the " + << "\"-ens\" option.\n\n"; + exit(1); + } + if(config_file.length() == 0) { + mlog << Error << "\nprocess_command_line() -> " + << "the output file must be set using the " + << "\"-out\" option.\n\n"; + exit(1); + } + if(config_file.length() == 0) { + mlog << Error << "\nprocess_command_line() -> " + << "the configuration file must be set using the " + << "\"-config\" option.\n\n"; + exit(1); + } + + // Create the default config file name + default_config_file = replace_path(default_config_filename); + + // List the config files + mlog << Debug(1) + << "Default Config File: " << default_config_file << "\n" + << "User Config File: " << config_file << "\n"; + + // Read the config files + conf_info.read_config(default_config_file, config_file); + + // Get the ensemble file type from config, if present + etype = parse_conf_file_type(conf_info.conf.lookup_dictionary(conf_key_ens)); + + // Read the first input ensemble file + if(!(ens_mtddf = mtddf_factory.new_met_2d_data_file(ens_files[0].c_str(), etype))) { + mlog << Error << "\nprocess_command_line() -> " + << "trouble reading ensemble file \"" << ens_files[0] << "\"\n\n"; + exit(1); + } + + // Store the input ensemble file type + etype = ens_mtddf->file_type(); + + // Process the configuration + conf_info.process_config(etype); + + // Allocate arrays to store threshold counts + thresh_cnt_na = new NumArray [conf_info.get_max_n_cat()]; + thresh_nbrhd_cnt_na = new NumArray * [conf_info.get_max_n_cat()]; + + for(i=0; i " + << "cannot open input ensemble file: " + << ens_files[i] << "\n\n"; + ens_file_vld.add(0); + } + else { + ens_file_vld.add(1); + } + } + + // Deallocate memory for data files + if(ens_mtddf) { delete ens_mtddf; ens_mtddf = (Met2dDataFile *) 0; } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void process_grid(const Grid &fcst_grid) { + Grid obs_grid; + + // Parse regridding logic + RegridInfo ri; + if(conf_info.get_n_var() > 0) { + ri = conf_info.ens_info[0]->regrid(); + } + else { + mlog << Error << "\nprocess_grid() -> " + << "at least one ensemble field must be specified!\n\n"; + exit(1); + } + + // Determine the verification grid + grid = parse_vx_grid(ri, &fcst_grid, &fcst_grid); + nxy = grid.nx() * grid.ny(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +bool get_data_plane(const char *infile, GrdFileType ftype, + VarInfo *info, DataPlane &dp) { + bool found; + Met2dDataFile *mtddf = (Met2dDataFile *) 0; + + // Read the current ensemble file + if(!(mtddf = mtddf_factory.new_met_2d_data_file(infile, ftype))) { + mlog << Error << "\nget_data_plane() -> " + << "trouble reading file \"" << infile << "\"\n\n"; + exit(1); + } + + // Read the gridded data field + if((found = mtddf->data_plane(*info, dp))) { + + // Setup the verification grid, if necessary + if(nxy == 0) process_grid(mtddf->grid()); + + // Create the output file, if necessary + if(nc_out == (NcFile *) 0) setup_nc_file(); + + // Regrid, if requested and necessary + if(!(mtddf->grid() == grid)) { + mlog << Debug(1) + << "Regridding field \"" << info->magic_str() + << "\" to the verification grid.\n"; + dp = met_regrid(dp, mtddf->grid(), grid, info->regrid()); + } + + // Store the valid time, if not already set + if(ens_valid_ut == (unixtime) 0) { + ens_valid_ut = dp.valid(); + } + // Check to make sure that the valid time doesn't change + else if(ens_valid_ut != dp.valid()) { + mlog << Warning << "\nget_data_plane() -> " + << "The valid time has changed, " + << unix_to_yyyymmdd_hhmmss(ens_valid_ut) + << " != " << unix_to_yyyymmdd_hhmmss(dp.valid()) + << " in \"" << infile << "\"\n\n"; + } + + } // end if found + + // Deallocate the data file pointer, if necessary + if(mtddf) { delete mtddf; mtddf = (Met2dDataFile *) 0; } + + return(found); +} + +//////////////////////////////////////////////////////////////////////// + +void process_ensemble() { + int i_var, i_file, n_ens_vld; + bool need_reset; + DataPlane ens_dp, ctrl_dp, cmn_dp, csd_dp; + unixtime max_init_ut = bad_data_ll; + + // Loop through each of the ensemble fields to be processed + for(i_var=0; i_varmagic_str() << "\n"; + + // Need to reinitialize counts and sums for each ensemble field + need_reset = true; + + // Loop through the ensemble member files + for(i_file=0, n_ens_vld=0; i_file " + << "ensemble field \"" << conf_info.ens_info[i_var]->magic_str() + << "\" not found in file \"" << ens_files[i_file] << "\"\n\n"; + continue; + } + else { + n_ens_vld++; + } + + // Reinitialize for the current variable + if(need_reset) { + + need_reset = false; + + // Reset the running sums and counts + clear_counts(); + + // Read climatology data for this field + cmn_dp = read_climo_data_plane( + conf_info.conf.lookup_array(conf_key_climo_mean_field, false), + i_var, ens_valid_ut, grid); + + csd_dp = read_climo_data_plane( + conf_info.conf.lookup_array(conf_key_climo_stdev_field, false), + i_var, ens_valid_ut, grid); + + // Read ensemble control member data, if provided + if(ctrl_file.nonempty()) { + + // Error out if missing + if(!get_data_plane(ctrl_file.c_str(), etype, + conf_info.ens_info[i_var], ctrl_dp)) { + mlog << Error << "\nprocess_ensemble() -> " + << "control member ensemble field \"" + << conf_info.ens_info[i_var]->magic_str() + << "\" not found in file \"" << ctrl_file << "\"\n\n"; + exit(1); + } + + // Apply current data to the running sums and counts + track_counts(i_var, ctrl_dp, true, cmn_dp, csd_dp); + + } // end if ctrl_file + + mlog << Debug(3) + << "Found " << (ctrl_dp.is_empty() ? 0 : 1) + << " control member, " << (cmn_dp.is_empty() ? 0 : 1) + << " climatology mean, and " << (csd_dp.is_empty() ? 0 : 1) + << " climatology standard deviation field(s) for \"" + << conf_info.ens_info[i_var]->magic_str() << "\".\n"; + + } // end if need_reset + + // Apply current data to the running sums and counts + track_counts(i_var, ens_dp, false, cmn_dp, csd_dp); + + // Keep track of the maximum initialization time + if(is_bad_data(max_init_ut) || ens_dp.init() > max_init_ut) { + max_init_ut = ens_dp.init(); + } + + } // end for i_file + + // Check for too much missing data + if(((double) n_ens_vld/n_ens) < conf_info.vld_ens_thresh) { + mlog << Error << "\nprocess_ensemble() -> " + << n_ens - n_ens_vld << " of " << n_ens + << " missing fields for \"" << conf_info.ens_info[i_var]->magic_str() + << "\" exceeds the maximum allowable specified by \"" + << conf_key_ens_ens_thresh << "\" (" << conf_info.vld_ens_thresh + << ") in the configuration file.\n\n"; + exit(1); + } + + // Write out the ensemble information to a NetCDF file + ens_dp.set_init(max_init_ut); + write_ens_nc(i_var, n_ens_vld, ens_dp, cmn_dp, csd_dp); + + } // end for i_var + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void clear_counts() { + int i, j; + + cnt_na.set_const(0.0, nxy); + min_na.set_const(bad_data_double, nxy); + max_na.set_const(bad_data_double, nxy); + sum_na.set_const(0.0, nxy); + + stdev_cnt_na.set_const(0.0, nxy); + stdev_sum_na.set_const(0.0, nxy); + stdev_ssq_na.set_const(0.0, nxy); + + for(i=0; i= max_na.buf()[i] || is_bad_data(max_na.buf()[i])) max_na.buf()[i] = ens; + + // Standard deviation sum, sum of squares, and count, excluding control member + if(!is_ctrl) { + stdev_sum_na.buf()[i] += ens; + stdev_ssq_na.buf()[i] += ens*ens; + stdev_cnt_na.buf()[i] += 1; + } + + // Event frequency + for(j=0; j 0 + if(conf_info.nc_info[i_var].do_nmep) { + DataPlane frac_dp; + + // Loop over thresholds + for(i=0; i 0) thresh_nbrhd_cnt_na[i][j].inc(k, 1); + } // end for k + + } // end for j + } // end for i + } // end if do_nmep + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void setup_nc_file() { + + // Create a new NetCDF file and open it + nc_out = open_ncfile(out_file.c_str(), true); + + if(IS_INVALID_NC_P(nc_out)) { + mlog << Error << "\nsetup_nc_file() -> " + << "trouble opening output NetCDF file " + << out_file << "\n\n"; + exit(1); + } + + // Add global attributes + write_netcdf_global(nc_out, out_file.text(), program_name, + conf_info.model.c_str()); + + // Add the projection information + write_netcdf_proj(nc_out, grid); + + // Define Dimensions + lat_dim = add_dim(nc_out, "lat", (long) grid.ny()); + lon_dim = add_dim(nc_out, "lon", (long) grid.nx()); + + // Add the lat/lon variables + if(conf_info.nc_info[0].do_latlon) { + write_netcdf_latlon(nc_out, &lat_dim, &lon_dim, grid); + } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void write_ens_nc(int i_var, int n_ens_vld, + const DataPlane &ens_dp, + const DataPlane &cmn_dp, + const DataPlane &csd_dp) { + int i, j, k; + double t; + char type_str[max_str_len]; + DataPlane prob_dp, nbrhd_dp; + + // Allocate memory for storing ensemble data + float *ens_mean = new float [nxy]; + float *ens_stdev = new float [nxy]; + float *ens_minus = new float [nxy]; + float *ens_plus = new float [nxy]; + float *ens_min = new float [nxy]; + float *ens_max = new float [nxy]; + float *ens_range = new float [nxy]; + int *ens_vld = new int [nxy]; + + // Store the threshold for the ratio of valid data points + t = conf_info.vld_data_thresh; + + // Store the data + for(i=0; i simp; + conf_info.cat_ta[i_var].get_simple_nodes(simp); + + // Process all CDP thresholds except 0 and 100 + for(vector::iterator it = simp.begin(); + it != simp.end(); it++) { + if(it->ptype() == perc_thresh_climo_dist && + !is_eq(it->pvalue(), 0.0) && + !is_eq(it->pvalue(), 100.0)) { + snprintf(type_str, sizeof(type_str), "CLIMO_CDP%i", + nint(it->pvalue())); + cdp_dp = normal_cdf_inv(it->pvalue()/100.0, cmn_dp, csd_dp); + write_ens_data_plane(i_var, cdp_dp, ens_dp, + type_str, + "Climatology distribution percentile"); + } + } // end for it + } + + // Deallocate and clean up + if(ens_mean) { delete [] ens_mean; ens_mean = (float *) 0; } + if(ens_stdev) { delete [] ens_stdev; ens_stdev = (float *) 0; } + if(ens_minus) { delete [] ens_minus; ens_minus = (float *) 0; } + if(ens_plus) { delete [] ens_plus; ens_plus = (float *) 0; } + if(ens_min) { delete [] ens_min; ens_min = (float *) 0; } + if(ens_max) { delete [] ens_max; ens_max = (float *) 0; } + if(ens_range) { delete [] ens_range; ens_range = (float *) 0; } + if(ens_vld) { delete [] ens_vld; ens_vld = (int *) 0; } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void write_ens_var_float(int i_var, float *ens_data, const DataPlane &dp, + const char *type_str, + const char *long_name_str) { + NcVar ens_var; + ConcatString ens_var_name, var_str, name_str, cs; + + // Append nc_var_str config file entry + cs = conf_info.nc_var_str[i_var]; + if(cs.length() > 0) var_str << "_" << cs; + + // Construct the variable name + ens_var_name << cs_erase + << conf_info.ens_info[i_var]->name_attr() << "_" + << conf_info.ens_info[i_var]->level_attr() + << var_str << "_" << type_str; + + // Skip variable names that have already been written + if(nc_ens_var_sa.has(ens_var_name)) return; + + // Otherwise, add to the list of previously defined variables + nc_ens_var_sa.add(ens_var_name); + + ens_var = add_var(nc_out, (string)ens_var_name, ncFloat, lat_dim, lon_dim); + + // + // Construct the variable name attribute + // For the ensemble mean, just use the variable name. + // For all other fields, append the field type. + // + if(strcmp(type_str, "ENS_MEAN") == 0) { + name_str << cs_erase + << conf_info.ens_info[i_var]->name_attr(); + } + else { + name_str << cs_erase + << conf_info.ens_info[i_var]->name_attr() << "_" + << type_str; + } + + // Add the variable attributes + add_var_att_local(conf_info.ens_info[i_var], &ens_var, false, dp, + name_str.c_str(), long_name_str); + + // Write the data + if(!put_nc_data_with_dims(&ens_var, &ens_data[0], grid.ny(), grid.nx())) { + mlog << Error << "\nwrite_ens_var_float() -> " + << "error in ens_var->put for the " << ens_var_name + << " field.\n\n"; + exit(1); + } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void write_ens_var_int(int i_var, int *ens_data, const DataPlane &dp, + const char *type_str, + const char *long_name_str) { + NcVar ens_var; + ConcatString ens_var_name, var_str, name_str, cs; + + // Append nc_var_str config file entry + cs = conf_info.nc_var_str[i_var]; + if(cs.length() > 0) var_str << "_" << cs; + + // Construct the variable name + ens_var_name << cs_erase + << conf_info.ens_info[i_var]->name_attr() << "_" + << conf_info.ens_info[i_var]->level_attr() + << var_str << "_" << type_str; + + // Skip variable names that have already been written + if(nc_ens_var_sa.has(ens_var_name)) return; + + // Otherwise, add to the list of previously defined variables + nc_ens_var_sa.add(ens_var_name); + + int deflate_level = conf_info.get_compression_level(); + ens_var = add_var(nc_out, (string)ens_var_name, ncInt, lat_dim, lon_dim, deflate_level); + + // Construct the variable name attribute + name_str << cs_erase + << conf_info.ens_info[i_var]->name_attr() << "_" + << type_str; + + // Add the variable attributes + add_var_att_local(conf_info.ens_info[i_var], &ens_var, true, dp, + name_str.c_str(), long_name_str); + + // Write the data + if(!put_nc_data_with_dims(&ens_var, &ens_data[0], grid.ny(), grid.nx())) { + mlog << Error << "\nwrite_ens_var_int() -> " + << "error in ens_var->put for the " << ens_var_name + << " field.\n\n"; + exit(1); + } + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void write_ens_data_plane(int i_var, const DataPlane &ens_dp, const DataPlane &dp, + const char *type_str, + const char *long_name_str) { + + // Allocate memory for this data + float *ens_data = new float [nxy]; + + // Store the data in an array of floats + for(int i=0; iname_attr() << " at " + << info->level_attr() << " " + << long_name_str; + + // Add variable attributes + add_att(nc_var, "name", name_str); + add_att(nc_var, "long_name", (string)att_str); + add_att(nc_var, "level", (string)info->level_attr()); + add_att(nc_var, "units", (string)info->units_attr()); + + if(is_int) add_att(nc_var, "_FillValue", bad_data_int); + else add_att(nc_var, "_FillValue", bad_data_float); + + // Write out times + write_netcdf_var_times(nc_var, dp); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void clean_up() { + int i, j; + + mlog << Debug(2) << "\n" << sep_str << "\n\n"; + + // List the output NetCDF file + mlog << Debug(1) << "Output file: " << out_file << "\n"; + + // Close the output NetCDF file + if(nc_out) { delete nc_out; nc_out = (NcFile *) 0; } + + // Deallocate threshold count arrays + if(thresh_cnt_na) { + for(i=0; i +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +using namespace netCDF; + +#include "gen_ens_prod_conf_info.h" + +#include "vx_data2d_factory.h" +#include "vx_grid.h" +#include "vx_util.h" +#include "vx_stat_out.h" + +//////////////////////////////////////////////////////////////////////// +// +// Constants +// +//////////////////////////////////////////////////////////////////////// + +static const char * program_name = "gen_ens_prod"; + +// Default configuration file name +static const char * default_config_filename = + "MET_BASE/config/GenEnsProdConfig_default"; + +//////////////////////////////////////////////////////////////////////// +// +// Variables for Command Line Arguments +// +//////////////////////////////////////////////////////////////////////// + +static StringArray ens_files; +static IntArray ens_file_vld; +static GrdFileType etype = FileType_None; +static int n_ens; +static GenEnsProdConfInfo conf_info; +static ConcatString config_file; +static ConcatString out_file; +static ConcatString ctrl_file; +static unixtime ens_valid_ut = (unixtime) 0; + +//////////////////////////////////////////////////////////////////////// +// +// Variables for Output Files +// +//////////////////////////////////////////////////////////////////////// + +// Output NetCDF file +static NcFile *nc_out = (NcFile *) 0; +static NcDim lat_dim; +static NcDim lon_dim; + +// List of output NetCDF variable names +static StringArray nc_ens_var_sa; + +//////////////////////////////////////////////////////////////////////// +// +// Miscellaneous Variables +// +//////////////////////////////////////////////////////////////////////// + +// Grid variables +static Grid grid; +static int nxy = 0; + +// Data file factory and input files +static Met2dDataFileFactory mtddf_factory; + +// Arrays to store running sums and counts +static NumArray cnt_na, min_na, max_na, sum_na; +static NumArray stdev_cnt_na, stdev_sum_na, stdev_ssq_na; +static NumArray *thresh_cnt_na = (NumArray *) 0; // [n_thresh] +static NumArray **thresh_nbrhd_cnt_na = (NumArray **) 0; // [n_thresh][n_nbrhd] + +//////////////////////////////////////////////////////////////////////// + +#endif // __GEN_ENS_PROD_H__ + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/other/gen_ens_prod/gen_ens_prod_conf_info.cc b/met/src/tools/other/gen_ens_prod/gen_ens_prod_conf_info.cc new file mode 100644 index 0000000000..ef7a274802 --- /dev/null +++ b/met/src/tools/other/gen_ens_prod/gen_ens_prod_conf_info.cc @@ -0,0 +1,362 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + +using namespace std; + +#include +#include +#include +#include +#include +#include +#include + +#include "gen_ens_prod_conf_info.h" +#include "configobjecttype_to_string.h" + +#include "vx_data2d_factory.h" +#include "vx_data2d.h" +#include "vx_log.h" + +#include "GridTemplate.h" + +//////////////////////////////////////////////////////////////////////// +// +// Code for class GenEnsProdConfInfo +// +//////////////////////////////////////////////////////////////////////// + +GenEnsProdConfInfo::GenEnsProdConfInfo() { + init_from_scratch(); +} + +//////////////////////////////////////////////////////////////////////// + +GenEnsProdConfInfo::~GenEnsProdConfInfo() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdConfInfo::init_from_scratch() { + + clear(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdConfInfo::clear() { + vector::const_iterator it; + + // Clear, erase, and initialize members + model.clear(); + desc.clear(); + for(it = ens_info.begin(); it != ens_info.end(); it++) { + if(*it) { delete *it; } + } + ens_info.clear(); + cdf_info.clear(); + cat_ta.clear(); + nc_var_str.clear(); + nbrhd_prob.clear(); + nmep_smooth.clear(); + vld_ens_thresh = bad_data_double; + vld_data_thresh = bad_data_double; + nc_info.clear(); + version.clear(); + + // Reset counts + n_var = 0; + max_n_cat = 0; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdConfInfo::read_config(const ConcatString default_file_name, + const ConcatString user_file_name) { + + // Read the config file constants + conf.read(replace_path(config_const_filename).c_str()); + + // Read the default config file + conf.read(default_file_name.c_str()); + + // Read the user-specified config file + conf.read(user_file_name.c_str()); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdConfInfo::process_config(GrdFileType etype) { + int i; + VarInfoFactory info_factory; + Dictionary *edict = (Dictionary *) 0; + Dictionary i_edict; + InterpMthd mthd; + + // Dump the contents of the config file + if(mlog.verbosity_level() >= 5) conf.dump(cout); + + // Initialize + clear(); + + // Conf: version + version = parse_conf_version(&conf); + + // Conf: model + model = parse_conf_string(&conf, conf_key_model); + + // Conf: desc + desc = parse_conf_string(&conf, conf_key_desc); + + // Conf: ens.field + edict = conf.lookup_array(conf_key_ens_field); + + // Determine the number of ensemble fields to be processed + if((n_var = parse_conf_n_vx(edict)) == 0) { + mlog << Error << "\nGenEnsProdConfInfo::process_config() -> " + << "At least one field must be specified in the \"" + << conf_key_ens_field << "\" dictionary!\n\n"; + exit(1); + } + + // Parse the ensemble field information + for(i=0,max_n_cat=0; iset_dict(i_edict); + + // Dump the contents of the current VarInfo + if(mlog.verbosity_level() >= 5) { + mlog << Debug(5) + << "Parsed ensemble field number " << i+1 << ":\n"; + ens_info[i]->dump(cout); + } + + // Conf: nc_var_str + nc_var_str.add(parse_conf_string(&i_edict, conf_key_nc_var_str, false)); + + // Conf: cat_thresh + cat_ta.push_back(i_edict.lookup_thresh_array(conf_key_cat_thresh)); + + // Dump the contents of the current thresholds + if(mlog.verbosity_level() >= 5) { + mlog << Debug(5) + << "Parsed thresholds for ensemble field number " << i+1 << ":\n"; + cat_ta[i].dump(cout); + } + + // Keep track of the maximum number of thresholds + if(cat_ta[i].n() > max_n_cat) max_n_cat = cat_ta[i].n(); + + // Conf: ensemble_flag + nc_info.push_back(parse_nc_info(&i_edict)); + } + + // Conf: ens.ens_thresh + vld_ens_thresh = conf.lookup_double(conf_key_ens_ens_thresh); + + // Check that the valid ensemble threshold is between 0 and 1. + if(vld_ens_thresh < 0.0 || vld_ens_thresh > 1.0) { + mlog << Error << "\nGenEnsProdConfInfo::process_config() -> " + << "The \"" << conf_key_ens_ens_thresh << "\" parameter (" + << vld_ens_thresh << ") must be set between 0 and 1.\n\n"; + exit(1); + } + + // Conf: ens.vld_thresh + vld_data_thresh = conf.lookup_double(conf_key_ens_vld_thresh); + + // Check that the valid data threshold is between 0 and 1. + if(vld_data_thresh < 0.0 || vld_data_thresh > 1.0) { + mlog << Error << "\nGenEnsProdConfInfo::process_config() -> " + << "The \"" << conf_key_ens_vld_thresh << "\" parameter (" + << vld_data_thresh << ") must be set between 0 and 1.\n\n"; + exit(1); + } + + // Conf: nbrhd_prob + nbrhd_prob = parse_conf_nbrhd(edict, conf_key_nbrhd_prob); + + // Conf: nmep_smooth + nmep_smooth = parse_conf_interp(edict, conf_key_nmep_smooth); + + // Loop through the neighborhood probability smoothing options + for(i=0; i " + << "Neighborhood probability smoothing methods DW_MEAN, " + << "LS_FIT, and BILIN are not supported for \"" + << conf_key_nmep_smooth << "\".\n\n"; + exit(1); + } + + // Check for valid neighborhood probability interpolation widths + if(nmep_smooth.width[i] < 1 || nmep_smooth.width[i]%2 == 0) { + mlog << Error << "\nGenEnsProdConfInfo::process_config() -> " + << "Neighborhood probability smoothing widths must be set " + << "to odd values greater than or equal to 1 (" + << nmep_smooth.width[i] << ") for \"" + << conf_key_nmep_smooth << "\".\n\n"; + exit(1); + } + } // end for i + + return; +} + +//////////////////////////////////////////////////////////////////////// + +GenEnsProdNcOutInfo GenEnsProdConfInfo::parse_nc_info(Dictionary *dict) { + GenEnsProdNcOutInfo cur; + + // Parse the ensemble flag + const DictionaryEntry *e = dict->lookup(conf_key_ensemble_flag); + + if(!e) { + mlog << Error + << "\nGenEnsProdConfInfo::parse_nc_info() -> " + << "lookup failed for key \"" << conf_key_ensemble_flag + << "\"\n\n"; + exit(1); + } + + // Boolean type + if(e->type() == BooleanType) { + if(e->b_value()) cur.set_all_true(); + else cur.set_all_false(); + } + // Dictionary type + else { + + // Check the type + if(e->type() != DictionaryType) { + mlog << Error << "\nGenEnsProdConfInfo::parse_nc_info() -> " + << "bad type (" << configobjecttype_to_string(e->type()) + << ") for key \"" << conf_key_ensemble_flag << "\"\n\n"; + exit(1); + } + + // Parse the various entries + Dictionary * d = e->dict_value(); + + cur.do_latlon = d->lookup_bool(conf_key_latlon_flag); + cur.do_mean = d->lookup_bool(conf_key_mean_flag); + cur.do_stdev = d->lookup_bool(conf_key_stdev_flag); + cur.do_minus = d->lookup_bool(conf_key_minus_flag); + cur.do_plus = d->lookup_bool(conf_key_plus_flag); + cur.do_min = d->lookup_bool(conf_key_min_flag); + cur.do_max = d->lookup_bool(conf_key_max_flag); + cur.do_range = d->lookup_bool(conf_key_range_flag); + cur.do_vld = d->lookup_bool(conf_key_vld_count_flag); + cur.do_freq = d->lookup_bool(conf_key_frequency_flag); + cur.do_nep = d->lookup_bool(conf_key_nep_flag); + cur.do_nmep = d->lookup_bool(conf_key_nmep_flag); + cur.do_climo = d->lookup_bool(conf_key_climo_flag); + cur.do_climo_cdp = d->lookup_bool(conf_key_climo_cdp_flag); + } + + return(cur); +} + +//////////////////////////////////////////////////////////////////////// +// +// Code for struct GenEnsProdNcOutInfo +// +//////////////////////////////////////////////////////////////////////// + +GenEnsProdNcOutInfo::GenEnsProdNcOutInfo() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdNcOutInfo::clear() { + + set_all_true(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +bool GenEnsProdNcOutInfo::all_false() const { + + bool status = do_latlon || do_mean || do_stdev || do_minus || + do_plus || do_min || do_max || do_range || + do_vld || do_freq || do_nep || do_nmep || + do_climo || do_climo_cdp; + + return(!status); +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdNcOutInfo::set_all_false() { + + do_latlon = false; + do_mean = false; + do_stdev = false; + do_minus = false; + do_plus = false; + do_min = false; + do_max = false; + do_range = false; + do_vld = false; + do_freq = false; + do_nep = false; + do_nmep = false; + do_climo = false; + do_climo_cdp = false; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenEnsProdNcOutInfo::set_all_true() { + + do_latlon = true; + do_mean = true; + do_stdev = true; + do_minus = true; + do_plus = true; + do_min = true; + do_max = true; + do_range = true; + do_vld = true; + do_freq = true; + do_nep = true; + do_nmep = true; + do_climo = true; + do_climo_cdp = true; + + return; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/other/gen_ens_prod/gen_ens_prod_conf_info.h b/met/src/tools/other/gen_ens_prod/gen_ens_prod_conf_info.h new file mode 100644 index 0000000000..2ec6014824 --- /dev/null +++ b/met/src/tools/other/gen_ens_prod/gen_ens_prod_conf_info.h @@ -0,0 +1,121 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** 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 +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + +#ifndef __GEN_ENS_PROD_CONF_INFO_H__ +#define __GEN_ENS_PROD_CONF_INFO_H__ + +//////////////////////////////////////////////////////////////////////// + +#include + +#include "vx_config.h" +#include "vx_data2d.h" +#include "vx_grid.h" +#include "vx_util.h" +#include "vx_cal.h" +#include "vx_math.h" + +//////////////////////////////////////////////////////////////////////// + +struct GenEnsProdNcOutInfo { + + bool do_latlon; + bool do_mean; + bool do_stdev; + bool do_minus; + bool do_plus; + bool do_min; + bool do_max; + bool do_range; + bool do_vld; + bool do_freq; + bool do_nep; + bool do_nmep; + bool do_climo; + bool do_climo_cdp; + + GenEnsProdNcOutInfo(); + + void clear(); + + bool all_false() const; + + void set_all_false(); + void set_all_true(); +}; + +//////////////////////////////////////////////////////////////////////// + +class GenEnsProdConfInfo { + + private: + + void init_from_scratch(); + + // Ensemble processing + int n_var; // Number of ensemble fields to be processed + int max_n_cat; // Maximum number of ensemble thresholds + + public: + + GenEnsProdConfInfo(); + ~GenEnsProdConfInfo(); + + ////////////////////////////////////////////////////////////////// + + // Gen-Ens-Prod configuration object + MetConfig conf; + + // Data parsed from the Gen-Ens-Prod configuration object + ConcatString model; // Model name + ConcatString desc; // Description + + vector ens_info; // Array of VarInfo pointers (allocated) + vector cdf_info; // Array of climo CDF info objects + vector cat_ta; // Array for ensemble categorical thresholds + StringArray nc_var_str; // Array of ensemble variable name strings + vector nc_info; // Array of ensemble product outputs + + NbrhdInfo nbrhd_prob; // Neighborhood probability definition + InterpInfo nmep_smooth; // Neighborhood maximum smoothing information + + double vld_ens_thresh; // Required ratio of valid input files + double vld_data_thresh; // Required ratio of valid data for each point + + ConcatString version; // Config file version + + ////////////////////////////////////////////////////////////////// + + void clear(); + + void read_config (const ConcatString, const ConcatString); + void process_config(GrdFileType); + + GenEnsProdNcOutInfo parse_nc_info(Dictionary *); + + // Accessor functions + int get_n_var() const; + int get_max_n_cat() const; + int get_n_nbrhd() const; + int get_compression_level(); +}; + +//////////////////////////////////////////////////////////////////////// + +inline int GenEnsProdConfInfo::get_n_var() const { return(n_var); } +inline int GenEnsProdConfInfo::get_max_n_cat() const { return(max_n_cat); } +inline int GenEnsProdConfInfo::get_n_nbrhd() const { return(nbrhd_prob.width.n()); } +inline int GenEnsProdConfInfo::get_compression_level() { return(conf.nc_compression()); } + +//////////////////////////////////////////////////////////////////////// + +#endif /* __GEN_ENS_PROD_CONF_INFO_H__ */ + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/other/shift_data_plane/shift_data_plane.cc b/met/src/tools/other/shift_data_plane/shift_data_plane.cc index 979babbfe1..5534ec3a48 100644 --- a/met/src/tools/other/shift_data_plane/shift_data_plane.cc +++ b/met/src/tools/other/shift_data_plane/shift_data_plane.cc @@ -258,7 +258,7 @@ void process_data_file() { for(x=0; x=35.0 ]; + GRIB1_ptv = 129; + ensemble_flag = TRUE; + }, + { + name = "UGRD"; + level = [ "Z10" ]; + cat_thresh = [ CDP75 ]; + }, + { + name = "WIND"; + level = [ "Z10" ]; + cat_thresh = [ >=CDP25&&<=CDP75 ]; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Neighborhood ensemble probabilities +// +nbrhd_prob = { + width = [ 5 ]; + shape = CIRCLE; + vld_thresh = 0.0; +} + +// +// NMEP smoothing methods +// +nmep_smooth = { + vld_thresh = 0.0; + shape = CIRCLE; + gaussian_dx = 81.27; + gaussian_radius = 120; + type = [ + { + method = GAUSSIAN; + width = 1; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Climatology data +// +climo_mean = ens; +climo_mean = { + + file_name = [ "${CLIMO_MEAN_FILE}" ]; + + regrid = { + method = NEAREST; + width = 1; + vld_thresh = 0.5; + shape = SQUARE; + } + + time_interp_method = DW_MEAN; + day_interval = 1; + hour_interval = 6; +} + +climo_stdev = climo_mean; +climo_stdev = { + file_name = [ "${CLIMO_STDEV_FILE}" ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Ensemble product output types +// May be set separately in each "ens.field" entry +// +ensemble_flag = { + latlon = TRUE; + mean = TRUE; + stdev = TRUE; + minus = FALSE; + plus = FALSE; + min = FALSE; + max = FALSE; + range = FALSE; + vld_count = FALSE; + frequency = TRUE; + nep = TRUE; + nmep = TRUE; + climo = TRUE; + climo_cdp = TRUE; +} + +//////////////////////////////////////////////////////////////////////////////// + +version = "V10.1.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/test/config/STATAnalysisConfig_SFC_SS_Index b/test/config/STATAnalysisConfig_SFC_SS_Index index ca08d2565c..994f091052 100644 --- a/test/config/STATAnalysisConfig_SFC_SS_Index +++ b/test/config/STATAnalysisConfig_SFC_SS_Index @@ -80,6 +80,13 @@ weight = [ 4.0, 3.0, 2.0, 1.0, //////////////////////////////////////////////////////////////////////////////// +// +// Array of STAT-Analysis jobs to be performed on the filtered data +// +jobs = []; + +//////////////////////////////////////////////////////////////////////////////// + // // Confidence interval settings // diff --git a/test/xml/unit_gen_ens_prod.xml b/test/xml/unit_gen_ens_prod.xml new file mode 100644 index 0000000000..9d3f2b0c75 --- /dev/null +++ b/test/xml/unit_gen_ens_prod.xml @@ -0,0 +1,78 @@ + + + + + + + + + +]> + + + + &TEST_DIR; + true + + + + + + + + + + echo "&DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep2/arw-sch-gep2_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep3/arw-tom-gep3_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/nmm-fer-gep4/nmm-fer-gep4_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-fer-gep5/arw-fer-gep5_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep6/arw-sch-gep6_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep7/arw-tom-gep7_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/nmm-fer-gep8/nmm-fer-gep8_2012040912_F024.grib" \ + > &OUTPUT_DIR;/gen_ens_prod/input_file_list; \ + &MET_BIN;/gen_ens_prod + + CLIMO_MEAN_FILE &DATA_DIR_CLIMO;/NCEP_1.0deg/cmean_1d.19790410 + CLIMO_STDEV_FILE &DATA_DIR_CLIMO;/NCEP_1.0deg/cstdv_1d.19790410 + + \ + -ens &OUTPUT_DIR;/gen_ens_prod/input_file_list \ + -config &CONFIG_DIR;/GenEnsProdConfig \ + -out &OUTPUT_DIR;/gen_ens_prod/gen_ens_prod_NO_CTRL_20120410_120000V.nc \ + -v 2 + + + &OUTPUT_DIR;/gen_ens_prod/gen_ens_prod_NO_CTRL_20120410_120000V.nc + + + + + echo "&DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep2/arw-sch-gep2_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep3/arw-tom-gep3_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/nmm-fer-gep4/nmm-fer-gep4_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-fer-gep5/arw-fer-gep5_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep6/arw-sch-gep6_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep7/arw-tom-gep7_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/nmm-fer-gep8/nmm-fer-gep8_2012040912_F024.grib" \ + > &OUTPUT_DIR;/gen_ens_prod/input_file_list; \ + &MET_BIN;/gen_ens_prod + + CLIMO_MEAN_FILE &DATA_DIR_CLIMO;/NCEP_1.0deg/cmean_1d.19790410 + CLIMO_STDEV_FILE &DATA_DIR_CLIMO;/NCEP_1.0deg/cstdv_1d.19790410 + + \ + -ens &OUTPUT_DIR;/gen_ens_prod/input_file_list \ + -ctrl &DATA_DIR_MODEL;/grib1/arw-tom-gep0/arw-tom-gep0_2012040912_F024.grib \ + -config &CONFIG_DIR;/GenEnsProdConfig \ + -out &OUTPUT_DIR;/gen_ens_prod/gen_ens_prod_WITH_CTRL_20120410_120000V.nc \ + -v 2 + + + &OUTPUT_DIR;/gen_ens_prod/gen_ens_prod_WITH_CTRL_20120410_120000V.nc + + + + diff --git a/test/xml/unit_met_test_scripts.xml b/test/xml/unit_met_test_scripts.xml index 2812eeeaaa..1d73214f86 100644 --- a/test/xml/unit_met_test_scripts.xml +++ b/test/xml/unit_met_test_scripts.xml @@ -23,7 +23,7 @@ true - + @@ -52,6 +52,23 @@ + + + + + + &MET_BIN;/gen_ens_prod + \ + -ens &MET_DATA;/sample_fcst/2009123112/*gep*/d01_2009123112_02400.grib \ + -config &MET_SCRIPTS;/config/GenEnsProdConfig \ + -out &OUTPUT_DIR;/gen_ens_prod/gen_ens_prod_20100101_120000V_ens.nc \ + -v 2 + + + &OUTPUT_DIR;/gen_ens_prod/gen_ens_prod_20100101_120000V_ens.nc + + + diff --git a/test/xml/unit_regrid.xml b/test/xml/unit_regrid.xml index 64f31916ed..8fdc3cb7a6 100644 --- a/test/xml/unit_regrid.xml +++ b/test/xml/unit_regrid.xml @@ -331,7 +331,7 @@ &DATA_DIR_MODEL;/grib1/gfs/gfs_2012040900_F036.grib \ &DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040912_F024.grib \ &OUTPUT_DIR;/regrid/regrid_data_plane_GFS_TO_HMT_MAX_5_SQUARE.nc \ - -field 'name="TMP"; level="Z2";' \ + -field 'name="TMP"; level="Z2";' \ -method MAX -width 5 -shape SQUARE \ -v 1 @@ -346,7 +346,7 @@ &DATA_DIR_MODEL;/grib1/gfs/gfs_2012040900_F036.grib \ G212 \ &OUTPUT_DIR;/regrid/regrid_data_plane_GFS_TO_G212_CONVERT_CENSOR.nc \ - -field 'name="TMP"; level="Z2";' \ + -field 'name="TMP"; level="Z2";' \ -convert 'convert(x) = x - 273.15;' \ -censor lt0 -9999 @@ -355,4 +355,20 @@ + + + + &MET_BIN;/regrid_data_plane + \ + &DATA_DIR_MODEL;/nccf/GloTEC_TEC_2015_03_17.nc \ + "latlon 360 91 -45.0 0 1.0 1.0" \ + &OUTPUT_DIR;/regrid/regrid_data_plane_WRAP_LON.nc \ + -field 'name="TEC"; level="(0,*,*)"; file_type=NETCDF_NCCF;' \ + -method MAX -width 5 -shape CIRCLE + + + &OUTPUT_DIR;/regrid/regrid_data_plane_WRAP_LON.nc + + +