diff --git a/data/config/GridStatConfig_default b/data/config/GridStatConfig_default index 2a2050f689..9ac3423238 100644 --- a/data/config/GridStatConfig_default +++ b/data/config/GridStatConfig_default @@ -234,6 +234,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -246,6 +247,7 @@ nc_pairs_flag = { diff = TRUE; climo = TRUE; climo_cdp = FALSE; + seeps = FALSE; weight = FALSE; nbrhd = FALSE; fourier = FALSE; @@ -254,6 +256,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/data/config/PointStatConfig_default b/data/config/PointStatConfig_default index f7e60a8e32..89676eeeac 100644 --- a/data/config/PointStatConfig_default +++ b/data/config/PointStatConfig_default @@ -290,6 +290,10 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; //////////////////////////////////////////////////////////////////////////////// tmp_dir = "/tmp"; diff --git a/data/config/TCPairsConfig_default b/data/config/TCPairsConfig_default index a768ac23e6..85a71cc97f 100644 --- a/data/config/TCPairsConfig_default +++ b/data/config/TCPairsConfig_default @@ -138,44 +138,65 @@ watch_warn = { } // -// Diagnostics to be extracted -// -diag_name = []; +// Metadata for diagnostic sources +// +diag_info_map = [ + { + diag_source = "CIRA_DIAG_RT"; + track_source = "GFS"; + field_source = "GFS_0p50"; + match_to_track = [ "GFS" ]; + diag_name = []; + }, + { + diag_source = "SHIPS_DIAG_RT"; + track_source = "SHIPS_TRK"; + field_source = "GFS_0p50"; + match_to_track = [ "OFCL" ]; + diag_name = []; + } +]; // -// Unit conversions to be applied based on diagnostic names and units +// Unit conversions for diagnostic sources, names, and units // diag_convert_map = [ - { source = "TCDIAG"; - key = [ "(10C)", "(10KT)", "(10M/S)" ]; - convert(x) = x / 10; }, - - { source = "LSDIAG_RT"; - key = [ "LAT", "LON", "CSST", "RSST", "DSST", "DSTA", "XDST", "XNST", "NSST", "NSTA", - "NTMX", "NTFR", "U200", "U20C", "V20C", "E000", "EPOS", "ENEG", "EPSS", "ENSS", - "T000", "TLAT", "TLON", "TWAC", "TWXC", "G150", "G200", "G250", "V000", "V850", - "V500", "V300", "SHDC", "SHGC", "T150", "T200", "T250", "SHRD", "SHRS", "SHRG", - "HE07", "HE05", "PW01", "PW02", "PW03", "PW04", "PW05", "PW06", "PW07", "PW08", - "PW09", "PW10", "PW11", "PW12", "PW13", "PW14", "PW15", "PW16", "PW17", "PW18", - "PW20", "PW21" ]; - convert(x) = x / 10; }, - - { source = "LSDIAG_RT"; - key = [ "VVAV", "VMFX", "VVAC" ]; - convert(x) = x / 100; }, - - { source = "LSDIAG_RT"; - key = [ "TADV" ]; - convert(x) = x / 1000000; }, - - { source = "LSDIAG_RT"; - key = [ "Z850", "D200", "TGRD", "DIVC" ]; - convert(x) = x / 10000000; }, - - { source = "LSDIAG_RT"; - key = [ "PENC", "PENV" ]; - convert(x) = x / 10 + 1000; } - + { + diag_source = "CIRA_DIAG"; + key = [ "(10C)", "(10KT)", "(10M/S)" ]; + convert(x) = x / 10; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "LAT", "LON", "CSST", "RSST", "DSST", "DSTA", "XDST", "XNST", "NSST", "NSTA", + "NTMX", "NTFR", "U200", "U20C", "V20C", "E000", "EPOS", "ENEG", "EPSS", "ENSS", + "T000", "TLAT", "TLON", "TWAC", "TWXC", "G150", "G200", "G250", "V000", "V850", + "V500", "V300", "SHDC", "SHGC", "T150", "T200", "T250", "SHRD", "SHRS", "SHRG", + "HE07", "HE05", "PW01", "PW02", "PW03", "PW04", "PW05", "PW06", "PW07", "PW08", + "PW09", "PW10", "PW11", "PW12", "PW13", "PW14", "PW15", "PW16", "PW17", "PW18", + "PW20", "PW21" ]; + convert(x) = x / 10; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "VVAV", "VMFX", "VVAC" ]; + convert(x) = x / 100; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "TADV" ]; + convert(x) = x / 1000000; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "Z850", "D200", "TGRD", "DIVC" ]; + convert(x) = x / 10000000; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "PENC", "PENV" ]; + convert(x) = x / 10 + 1000; + } ]; // diff --git a/data/table_files/met_header_columns_V11.0.txt b/data/table_files/met_header_columns_V11.0.txt index c991a2847f..c393863057 100644 --- a/data/table_files/met_header_columns_V11.0.txt +++ b/data/table_files/met_header_columns_V11.0.txt @@ -36,6 +36,6 @@ V11.0 : STAT : SSIDX : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID V11.0 : MODE : OBJ : VERSION MODEL N_VALID GRID_RES DESC FCST_LEAD FCST_VALID FCST_ACCUM OBS_LEAD OBS_VALID OBS_ACCUM FCST_RAD FCST_THR OBS_RAD OBS_THR FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE OBJECT_ID OBJECT_CAT CENTROID_X CENTROID_Y CENTROID_LAT CENTROID_LON AXIS_ANG LENGTH WIDTH AREA AREA_THRESH CURVATURE CURVATURE_X CURVATURE_Y COMPLEXITY INTENSITY_10 INTENSITY_25 INTENSITY_50 INTENSITY_75 INTENSITY_90 INTENSITY_USER INTENSITY_SUM CENTROID_DIST BOUNDARY_DIST CONVEX_HULL_DIST ANGLE_DIFF ASPECT_DIFF AREA_RATIO INTERSECTION_AREA UNION_AREA SYMMETRIC_DIFF INTERSECTION_OVER_AREA CURVATURE_RATIO COMPLEXITY_RATIO PERCENTILE_INTENSITY_RATIO INTEREST V11.0 : MODE : CTS : VERSION MODEL N_VALID GRID_RES DESC FCST_LEAD FCST_VALID FCST_ACCUM OBS_LEAD OBS_VALID OBS_ACCUM FCST_RAD FCST_THR OBS_RAD OBS_THR FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE FIELD TOTAL FY_OY FY_ON FN_OY FN_ON BASER FMEAN ACC FBIAS PODY PODN POFD FAR CSI GSS HK HSS ODDS -V11.0 : TCST : TCMPR : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE TOTAL INDEX LEVEL WATCH_WARN INITIALS ALAT ALON BLAT BLON TK_ERR X_ERR Y_ERR ALTK_ERR CRTK_ERR ADLAND BDLAND AMSLP BMSLP AMAX_WIND BMAX_WIND AAL_WIND_34 BAL_WIND_34 ANE_WIND_34 BNE_WIND_34 ASE_WIND_34 BSE_WIND_34 ASW_WIND_34 BSW_WIND_34 ANW_WIND_34 BNW_WIND_34 AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH NUM_MEMBERS TRACK_SPREAD DIST_MEAN MSLP_SPREAD MAX_WIND_SPREAD -V11.0 : TCST : TCDIAG : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE TOTAL INDEX DIAG_SOURCE (N_DIAG) DIAG_[0-9]* VALUE_[0-9]* +V11.0 : TCST : TCMPR : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE TOTAL INDEX LEVEL WATCH_WARN INITIALS ALAT ALON BLAT BLON TK_ERR X_ERR Y_ERR ALTK_ERR CRTK_ERR ADLAND BDLAND AMSLP BMSLP AMAX_WIND BMAX_WIND AAL_WIND_34 BAL_WIND_34 ANE_WIND_34 BNE_WIND_34 ASE_WIND_34 BSE_WIND_34 ASW_WIND_34 BSW_WIND_34 ANW_WIND_34 BNW_WIND_34 AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH NUM_MEMBERS TRACK_SPREAD TRACK_STDEV MSLP_STDEV MAX_WIND_STDEV +V11.0 : TCST : TCDIAG : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE TOTAL INDEX DIAG_SOURCE TRACK_SOURCE FIELD_SOURCE (N_DIAG) DIAG_[0-9]* VALUE_[0-9]* V11.0 : TCST : PROBRIRW : VERSION AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT LEAD VALID INIT_MASK VALID_MASK LINE_TYPE ALAT ALON BLAT BLON INITIALS TK_ERR X_ERR Y_ERR ADLAND BDLAND RIRW_BEG RIRW_END RIRW_WINDOW AWIND_END BWIND_BEG BWIND_END BDELTA BDELTA_MAX BLEVEL_BEG BLEVEL_END (N_THRESH) THRESH_[0-9]* PROB_[0-9]* diff --git a/docs/Users_Guide/config_options.rst b/docs/Users_Guide/config_options.rst index b01f079c64..f9cf4b903a 100644 --- a/docs/Users_Guide/config_options.rst +++ b/docs/Users_Guide/config_options.rst @@ -1456,6 +1456,14 @@ In this way, the number of bins impacts the resolution of the climatological probabilities. These derived probability values are used to compute the climatological Brier Score and Brier Skill Score. + +The "seeps_p1_thresh" option controls the threshold of p1 (probability of being dry) values. + +.. code-block:: none + + seeps_p1_thresh = >=0.1&&<=0.85; + + mask_missing_flag ^^^^^^^^^^^^^^^^^ diff --git a/docs/Users_Guide/grid-stat.rst b/docs/Users_Guide/grid-stat.rst index 7b550ae0b7..2887cd77a6 100644 --- a/docs/Users_Guide/grid-stat.rst +++ b/docs/Users_Guide/grid-stat.rst @@ -362,6 +362,7 @@ _____________________ nbrcnt = BOTH; grad = BOTH; dmap = BOTH; + seeps = NONE; } @@ -410,9 +411,14 @@ The **output_flag** array controls the type of output that the Grid-Stat tool ge 21. **DMAP** for Distance Map Statistics +22. **SEEPS** for SEEPS (Stable Equitable Error in Probability Space) score. It's described in :numref:`table_PS_format_info_SEEPS`. The SEEPS score of matched pair data is saved into the NetCDF. + Note that the first two line types are easily derived from one another. The user is free to choose which measure is most desired. The output line types are described in more detail in :numref:`grid_stat-output`. +The SEEPS climo file is not distributed with MET tools because of the file size. It should be configured by using the environment variable, MET_SEEPS_GRID_CLIMO_NAME. + + _____________________ .. code-block:: none diff --git a/docs/Users_Guide/point-stat.rst b/docs/Users_Guide/point-stat.rst index f52db68945..c5bb22b5e9 100644 --- a/docs/Users_Guide/point-stat.rst +++ b/docs/Users_Guide/point-stat.rst @@ -487,6 +487,8 @@ Note that writing out matched pair data (MPR lines) for a large number of cases If all line types corresponding to a particular verification method are set to NONE, the computation of those statistics will be skipped in the code and thus make the Point-Stat tool run more efficiently. For example, if FHO, CTC, and CTS are all set to NONE, the Point-Stat tool will skip the categorical verification step. +The default SEEPS climo file exists at MET_BASE/share/met/climo/seeps/PPT24_seepsweights.nc. It can be overridden by using the environment variable, MET_SEEPS_POINT_CLIMO_NAME. + .. _point_stat-output: point_stat output diff --git a/docs/Users_Guide/stat-analysis.rst b/docs/Users_Guide/stat-analysis.rst index 30ad22187f..1c1f1db4c0 100644 --- a/docs/Users_Guide/stat-analysis.rst +++ b/docs/Users_Guide/stat-analysis.rst @@ -40,12 +40,12 @@ The **-derive** job command option can be used to perform the derivation of stat Aggregated values from multiple STAT lines ------------------------------------------ -The Stat-Analysis "aggregate" job aggregates values from multiple STAT lines of the same type. The user may specify the specific line type of interest and any other relevant search criteria. The Stat-Analysis tool then creates sums of each of the values in all lines matching the search criteria. The aggregated data are output as the same line type as the user specified. The STAT line types which may be aggregated in this way are the contingency table (FHO, CTC, PCT, MCTC, NBRCTC), partial sums (SL1L2, SAL1L2, VL1L2, and VAL1L2), and other (ISC, ECNT, RPS, RHIST, PHIST, RELP, NBRCNT, SSVAR, and GRAD) line types. +The Stat-Analysis "aggregate" job aggregates values from multiple STAT lines of the same type. The user may specify the specific line type of interest and any other relevant search criteria. The Stat-Analysis tool then creates sums of each of the values in all lines matching the search criteria. The aggregated data are output as the same line type as the user specified. The STAT line types which may be aggregated in this way are the contingency table (FHO, CTC, PCT, MCTC, NBRCTC), partial sums (SL1L2, SAL1L2, VL1L2, and VAL1L2), and other (ISC, ECNT, RPS, RHIST, PHIST, RELP, NBRCNT, SSVAR, GRAD, and SEEPS) line types. Aggregate STAT lines and produce aggregated statistics ------------------------------------------------------ -The Stat-Analysis "aggregate-stat" job aggregates multiple STAT lines of the same type together and produces relevant statistics from the aggregated line. This may be done in the same manner listed above in :numref:`StA_Aggregated-values-from`. However, rather than writing out the aggregated STAT line itself, the relevant statistics generated from that aggregated line are provided in the output. Specifically, if a contingency table line type (FHO, CTC, PCT, MCTC, or NBRCTC) has been aggregated, contingency table statistics (CTS, ECLV, PSTD, MCTS, or NBRCTS) line types can be computed. If a partial sums line type (SL1L2 or SAL1L2) has been aggregated, the continuous statistics (CNT) line type can be computed. If a vector partial sums line type (VL1L2 or VAL1L2) has been aggregated, the vector continuous statistics (VCNT) line type can be computed. For ensembles, the ORANK line type can be accumulated into ECNT, RPS, RHIST, PHIST, RELP, or SSVAR output. If the matched pair line type (MPR) has been aggregated, may output line types (FHO, CTC, CTS, CNT, MCTC, MCTS, SL1L2, SAL1L2, VL1L2, VCNT, WDIR, PCT, PSTD, PJC, PRC, or ECLV) can be computed. Multiple output line types may be specified for each "aggregate-stat" job, as long as each output is derivable from the input. +The Stat-Analysis "aggregate-stat" job aggregates multiple STAT lines of the same type together and produces relevant statistics from the aggregated line. This may be done in the same manner listed above in :numref:`StA_Aggregated-values-from`. However, rather than writing out the aggregated STAT line itself, the relevant statistics generated from that aggregated line are provided in the output. Specifically, if a contingency table line type (FHO, CTC, PCT, MCTC, or NBRCTC) has been aggregated, contingency table statistics (CTS, ECLV, PSTD, MCTS, or NBRCTS) line types can be computed. If a partial sums line type (SL1L2 or SAL1L2) has been aggregated, the continuous statistics (CNT) line type can be computed. If a vector partial sums line type (VL1L2 or VAL1L2) has been aggregated, the vector continuous statistics (VCNT) line type can be computed. For ensembles, the ORANK line type can be accumulated into ECNT, RPS, RHIST, PHIST, RELP, or SSVAR output. If a SEEPS matched pair line type (SEEPS_MPR) has been aggregated, the aggregated SEEPS line type (SEEPS) can be computed. If the matched pair line type (MPR) has been aggregated, may output line types (FHO, CTC, CTS, CNT, MCTC, MCTS, SL1L2, SAL1L2, VL1L2, VCNT, WDIR, PCT, PSTD, PJC, PRC, or ECLV) can be computed. Multiple output line types may be specified for each "aggregate-stat" job, as long as each output is derivable from the input. When aggregating the matched pair line type (MPR), additional required job command options are determined by the requested output line type(s). For example, the "-out_thresh" (or "-out_fcst_thresh" and "-out_obs_thresh" options) are required to compute contingnecy table counts (FHO, CTC) or statistics (CTS). Those same job command options can also specify filtering thresholds when computing continuous partial sums (SL1L2, SAL1L2) or statistics (CNT). Output is written for each threshold specified. diff --git a/docs/Users_Guide/tc-pairs.rst b/docs/Users_Guide/tc-pairs.rst index 02e1c2cb2d..40162dee4d 100644 --- a/docs/Users_Guide/tc-pairs.rst +++ b/docs/Users_Guide/tc-pairs.rst @@ -14,10 +14,39 @@ Scientific and statistical aspects .. _TC-Pairs_Diagnostics: -Storm Diagnostics +TC Diagnostics ----------------- -TODO: Add a paragraph about storm diagnostics, describing what they are, why they are important, and how they can be generated. +TC diagnostics provide information about a TC's structure or its environment. Each TC diagnostic is a single-valued measure that corresponds to some aspect of the storm itself or the surrounding large-scale environment. TC diagnostics can be derived from observational analyses, model fields, or even satellite observations. Examples include: + + * Inner core diagnostics provide information about the structure of the storm near the storm center. Examples include the intensity of the storm and the radius of maximum winds. + + * Large scale diagnostics provide information about quantities that characterize its environment. Examples include environmental vertical wind shear, total precipitable water, the average relative humidity, measures of convective instability, and the upper bound of intensity that a storm may be expected to achieve in its current environment. These diagnostics are typically derived from model fields as an average of the quantity of interest over either a circular area or an annulus centered on the storm center. Often, the storm center is taken to be the underlying model's storm center. In other cases, the diagnostics may be computed along some other specified track. + + * Ocean-based diagnostics provide information about the sea's thermal characteristics in the vicinity of the storm center. Examples include the sea surface temperature, ocean heat content, and the depth of warm water of a given temperature. + + * Satellite-based diagnostics provide information about the storm structure as observed by geostationary satellite infrared imagery. Examples include information about the shape and extent of the cold-cirrus canopy of the TC and whether patterns are present that may portend intensification. + +Diagnostics are critically important for training and running statistical-dynamical models that predict a TC's intensity or size. One of the most well-known diagnostics sets is that of the Statistical Hurricane Intensity Prediction Scheme (SHIPS), which supports predictions of TC intensity. A large 30-year development dataset of TC diagnostics has been retrospectively derived to support the training of the SHIPS intensity model as well as other related models such as the Logistic Growth Equation Model (LGEM), SHIPS Rapid Intensification Index (SHIPS-RII), and others. These diagnostics, called *lsdiag* for "large scale" environment, are computed using a *perfect prog* approach in which the diagnostics are computed on the reference model's verifying analyses to generate a set of time-dependent diagnostics from t=0 out to the desired maximum forecast lead time. This is repeated for each initialization, building up a full history of diagnostics for each storm. By using the subsequent verifying analysis for later lead times, the model is taken to be "perfect", removing the impact of model forecast errors. The resulting developmental dataset is ideal for training statistical-dynamical models such as SHIPS. To generate forecasts in real-time, the diagnostics are computed along a forecast track (often taken to be the National Hurricane Center's official forecast) using the fields of the underlying NWP model (e.g, the Global Forecast System, or GFS model). The resulting diagnostics are then used as *predictors* in models like SHIPS and LGEM to predict a TC's future intensity or probability of undergoing rapid intensification. + +Beside their use in TC prediction, TC diagnostics can be very useful to forecasters to understand the forecast scenario. They are also useful to model developers for evaluation of model errors and understanding model performance under different environmental conditions. For instance, a modeler may wish to understand their model's track biases under conditions of high vertical wind shear. TC diagnostics can also be used to understand the sensitivity of the model's intensity predictions to oceanic conditions such as upwelling. The TC-Pairs tool allows filtering and subsetting based on the values of one or several TC diagnostics. + +As of MET v11.0.0, two types of TC diagnostics are supported in TC-Pairs: + + .. SHIPS_DIAG_DEV: Includes a plethora of inner core, environmental, oceanic, and satellite-based diagnostics. These diagnostics are computed using the *perfect prog* approach. + + * SHIPS_DIAG_RT: Real-time SHIPS diagnostics computed from a NWP model such as the Global Forecast System (GFS) model along a GFS forecast track defined by a SHIPS-specific tracker. Note that these SHIPS-derived forecast tracks do not appear in the NHC adeck data. + + * CIRA_DIAG_RT: Real-time model-based diagnostics computed along the model's predicted track. + +Diagnostics from the SHIPS Development Datasets (SHIPS_DIAG_DEV) will be supported in a future release of MET. + +A future version of MET will also allow the CIRA model diagnostics to be computed directly from model forecast fields. Until then, users may obtain the SHIPS diagnostics at the following locations: + + * SHIPS_DIAG_DEV: https://rammb2.cira.colostate.edu/research/tropical-cyclones/ships/#DevelopmentalData + + * SHIPS_DIAG_RT: https://ftp.nhc.noaa.gov/atcf/lsdiag/ + .. _TC-Pairs_Practical-information: @@ -58,7 +87,7 @@ Required arguments for tc_pairs Optional arguments for tc_pairs ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -5. The **-diag source path** argument indicates the TC-Pairs acceptable format data containing the tropical cyclone diagnostics dataset corresponding to the adeck tracks. The **source** can be set to TCDIAG, LSDIAG_RT, or LSDIAG_DEV to indicate the input diagnostics data source. The **path** argument specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. +5. The **-diag source path** argument indicates the TC-Pairs acceptable format data containing the tropical cyclone diagnostics dataset corresponding to the adeck tracks. The **source** can be set to CIRA_DIAG_RT or SHIPS_DIAG_RT to indicate the input diagnostics data source. The **path** argument specifies the name of a TC-Pairs acceptable format file or top-level directory containing TC-Pairs acceptable format files ending in ".dat" to be processed. Support for additional diagnostic sources will be added in future releases. 6. The -**out base** argument indicates the path of the output file base. This argument overrides the default output file base (**./out_tcmpr**). @@ -70,8 +99,6 @@ This tool currently only supports the rapid intensification (**RI**) edeck proba At least one **-adeck** or **-edeck** option must be specified. The **-adeck, -edeck**, and **-bdeck** options may optionally be followed with **suffix=string** to append that string to all model names found within that data source. This option may be useful when processing track data from two different sources which reuse the same model names. -The **-diag** option may optionally be followed with **model=string** to override the model name of the tracks to which those diagnostics correspond. The **string** specifies a comma-separated list of one or more ATCF ID's to which these diagnostics should be paired (e.g. **model=OFCL,SHIP**). - An example of the tc_pairs calling sequence is shown below: .. code-block:: none @@ -175,7 +202,7 @@ ____________________ } ]; -The **consensus** field allows the user to generate a user-defined consensus forecasts from any number of models. All models used in the consensus forecast need to be included in the **model** field (first entry in **TCPairsConfig_default**). The name field is the desired consensus model name. The **members** field is a comma-separated list of model IDs that make up the members of the consensus. The **required** field is a comma-separated list of true/false values associated with each consensus member. If a member is designated as true, the member is required to be present in order for the consensus to be generated. If a member is false, the consensus will be generated regardless of whether the member is present. The length of the required array must be the same length as the members array. The **min_req** field is the number of members required in order for the consensus to be computed. The required and min_req field options are applied at each forecast lead time. If any member of the consensus has a non-valid position or intensity value, the consensus for that valid time will not be generated. If a consensus model is indicated in the configuration file there will be non-missing output for the consensus track variables in the output file (NUM_MEMBERS, TRACK_SPREAD, DIST_MEAN, MSLP_SPREAD, MAX_WIND_SPREAD). See the TCMPR line type definitions below. +The **consensus** field allows the user to generate a user-defined consensus forecasts from any number of models. All models used in the consensus forecast need to be included in the **model** field (first entry in **TCPairsConfig_default**). The name field is the desired consensus model name. The **members** field is a comma-separated list of model IDs that make up the members of the consensus. The **required** field is a comma-separated list of true/false values associated with each consensus member. If a member is designated as true, the member is required to be present in order for the consensus to be generated. If a member is false, the consensus will be generated regardless of whether the member is present. The length of the required array must be the same length as the members array. The **min_req** field is the number of members required in order for the consensus to be computed. The required and min_req field options are applied at each forecast lead time. If any member of the consensus has a non-valid position or intensity value, the consensus for that valid time will not be generated. If a consensus model is indicated in the configuration file there will be non-missing output for the consensus track variables in the output file (NUM_MEMBERS, TRACK_SPREAD, TRACK_STDEV, MSLP_STDEV, MAX_WIND_STDEV). See the TCMPR line type definitions below. ____________________ @@ -258,52 +285,75 @@ ____________________ .. code-block:: none - diag_name = []; - -The **diag_name** entry specifies a comma-separated list of strings for the tropical cyclone diagnostics of interest. This applies when the **-tcdiag** and/or **-lsdiag** command line options have been used to provide storm diagnostics data. If a non-zero list of diagnostic names is specified, only those diagnostics appearing in the list are written to the TCDIAG output line type. If defined as an empty list (default), all diagnostics found in the input are written to the TCDIAG output lines. - -A TCMPR line is written to the output for each track point. If diagnostics data is also defined for that track point, a TCDIAG line is written immediately after the corresponding TCMPR line. The contents of that TCDIAG line is determined by diagnostic names requested in the **diag_name** entry. + diag_info_map = [ + { + diag_source = "CIRA_DIAG_RT"; + track_source = "GFS"; + field_source = "GFS_0p50"; + match_to_track = [ "GFS" ]; + diag_name = []; + }, + { + diag_source = "SHIPS_DIAG_RT"; + track_source = "SHIPS_TRK"; + field_source = "GFS_0p50"; + match_to_track = [ "OFCL" ]; + diag_name = []; + } + ]; + +A TCMPR line is written to the output for each track point. If diagnostics data is also defined for that track point, a TCDIAG line is written immediately after the corresponding TCMPR line. The contents of that TCDIAG line is determined by the **diag_info_map** entry. + +The **diag_info_map** entries define how the diagnostics read with the **-diag** command line option should be used. Each array element is a dictionary consisting of entries for **diag_source**, **track_source**, **field_source**, **match_to_track**, and **diag_name**. + +The **diag_source** entry is one of the supported diagnostics data sources. The **track_source** entry is a string defining the ATCF ID of the track data used to define the locations at which diagnostics are computed. This string is written to the **TRACK_SOURCE** column of the TCDIAG output line. The **field_source** entry is a string describing the gridded model data from which the diagnostics are computed. This string is written to the **FIELD_SOURCE** column of the TCDIAG output line type. The **match_to_track** entry specifies a comma-separated list of strings defining the ATCF ID(s) of the tracks to which these diagnostic values should be matched. The **diag_name** entry specifies a comma-separated list of strings for the tropical cyclone diagnostics of interest. If a non-zero list of diagnostic names is specified, only those diagnostics appearing in the list are written to the TCDIAG output line type. If defined as an empty list (default), all diagnostics found in the input are written to the TCDIAG output lines. ____________________ .. code-block:: none - diag_convert_map = [ - { source = "TCDIAG"; - key = [ "(10C)", "(10KT)", "(10M/S)" ]; - convert(x) = x / 10; }, - - { source = "LSDIAG_RT"; - key = [ "LAT", "LON", "CSST", "RSST", "DSST", "DSTA", "XDST", "XNST", "NSST", "NSTA", - "NTMX", "NTFR", "U200", "U20C", "V20C", "E000", "EPOS", "ENEG", "EPSS", "ENSS", - "T000", "TLAT", "TLON", "TWAC", "TWXC", "G150", "G200", "G250", "V000", "V850", - "V500", "V300", "SHDC", "SHGC", "T150", "T200", "T250", "SHRD", "SHRS", "SHRG", - "HE07", "HE05", "PW01", "PW02", "PW03", "PW04", "PW05", "PW06", "PW07", "PW08", - "PW09", "PW10", "PW11", "PW12", "PW13", "PW14", "PW15", "PW16", "PW17", "PW18", - "PW20", "PW21" ]; - convert(x) = x / 10; }, - - { source = "LSDIAG_RT"; - key = [ "VVAV", "VMFX", "VVAC" ]; - convert(x) = x / 100; }, - - { source = "LSDIAG_RT"; - key = [ "TADV" ]; - convert(x) = x / 1000000; }, - - { source = "LSDIAG_RT"; - key = [ "Z850", "D200", "TGRD", "DIVC" ]; - convert(x) = x / 10000000; }, - - { source = "LSDIAG_RT"; - key = [ "PENC", "PENV" ]; - convert(x) = x / 10 + 1000; } - - ]; - -The **diag_convert_map** entries define conversion functions to be applied to diagnostics data read with the **-diag** command line option. Each array element is a dictionary consisting of a **source**, **key**, and **convert(x)** entry. - -The **source** is one of the supported diagnostics data sources. The **key** is an array of strings. The strings can specify diagnostic names or units, although units are only checked for **TCDIAG** sources. If both the name and units are specified, the conversion function for the name takes precedence. **convert(x)** is a function of one variable which defines how the diagnostic data should be converted. The defined function is applied to any diagnostic value whose name or units appears in the **key**. + diag_convert_map = [ + { + diag_source = "CIRA_DIAG"; + key = [ "(10C)", "(10KT)", "(10M/S)" ]; + convert(x) = x / 10; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "LAT", "LON", "CSST", "RSST", "DSST", "DSTA", "XDST", "XNST", "NSST", "NSTA", + "NTMX", "NTFR", "U200", "U20C", "V20C", "E000", "EPOS", "ENEG", "EPSS", "ENSS", + "T000", "TLAT", "TLON", "TWAC", "TWXC", "G150", "G200", "G250", "V000", "V850", + "V500", "V300", "SHDC", "SHGC", "T150", "T200", "T250", "SHRD", "SHRS", "SHRG", + "HE07", "HE05", "PW01", "PW02", "PW03", "PW04", "PW05", "PW06", "PW07", "PW08", + "PW09", "PW10", "PW11", "PW12", "PW13", "PW14", "PW15", "PW16", "PW17", "PW18", + "PW20", "PW21" ]; + convert(x) = x / 10; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "VVAV", "VMFX", "VVAC" ]; + convert(x) = x / 100; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "TADV" ]; + convert(x) = x / 1000000; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "Z850", "D200", "TGRD", "DIVC" ]; + convert(x) = x / 10000000; + }, + { + diag_source = "SHIPS_DIAG"; + key = [ "PENC", "PENV" ]; + convert(x) = x / 10 + 1000; + } + ]; + +The **diag_convert_map** entries define conversion functions to be applied to diagnostics data read with the **-diag** command line option. Each array element is a dictionary consisting of a **diag_source**, **key**, and **convert(x)** entry. + +The **diag_source** entry is one of the supported diagnostics data sources. Partial string matching logic is applied, so **SHIPS_DIAG** entries are matched to both **SHIPS_DIAG_RT** and **SHIPS_DIAG_DEV** diagnostic sources. The **key** entry is an array of strings. The strings can specify diagnostic names or units, although units are only checked for **CIRA_DIAG** sources. If both the name and units are specified, the conversion function for the name takes precedence. The **convert(x)** entry is a function of one variable which defines how the diagnostic data should be converted. The defined function is applied to any diagnostic value whose name or units appears in the **key**. ____________________ @@ -549,15 +599,15 @@ TC-Pairs produces output in TCST format. The default output file name can be ove - consensus variable: number of models (or ensemble members) that were used to build the consensus track * - 81 - TRACK_SPREAD - - consensus variable: the standard deviation of the distances from the member locations to the consensus track location (nm) - * - 82 - - DIST_MEAN - consensus variable: the mean of the distances from the member location to the consensus track location (nm) + * - 82 + - TRACK_STDEV + - consensus variable: the standard deviation of the distances from the member locations to the consensus track location (nm) * - 83 - - MSLP_SPREAD + - MSLP_STDEV - consensus variable: the standard deviation of the member's mean sea level pressure values * - 84 - - MAX_WIND_SPREAD + - MAX_WIND_STDEV - consensus variable: the standard deviation of the member's maximum wind speed values .. _TCDIAG Line Type: @@ -583,14 +633,20 @@ TC-Pairs produces output in TCST format. The default output file name can be ove - Index of the current track pair * - 16 - DIAG_SOURCE - - Diagnostics data source + - Diagnostics data source indicated by the `-diag` command line option * - 17 + - TRACK_SOURCE + - ATCF ID of the track data used to define the diagnostics + * - 18 + - FIELD_SOURCE + - Description of gridded field data source used to define the diagnostics + * - 19 - N_DIAG - Number of storm diagnostic name and value columns to follow - * - 18 + * - 20 - DIAG_i - Name of the of the ith storm diagnostic (repeated) - * - 19 + * - 21 - VALUE_i - Value of the ith storm diagnostic (repeated) diff --git a/internal/test_unit/config/GridStatConfig_APCP_regrid b/internal/test_unit/config/GridStatConfig_APCP_regrid index c9e10b1d6f..c0a4ed6ca5 100644 --- a/internal/test_unit/config/GridStatConfig_APCP_regrid +++ b/internal/test_unit/config/GridStatConfig_APCP_regrid @@ -180,6 +180,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -199,6 +200,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_GRIB_lvl_typ_val b/internal/test_unit/config/GridStatConfig_GRIB_lvl_typ_val index 767ef54558..c56f154351 100644 --- a/internal/test_unit/config/GridStatConfig_GRIB_lvl_typ_val +++ b/internal/test_unit/config/GridStatConfig_GRIB_lvl_typ_val @@ -277,6 +277,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; }; // @@ -296,6 +297,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_GRIB_set_attr b/internal/test_unit/config/GridStatConfig_GRIB_set_attr index d4d2b2c387..28a0ab2f19 100644 --- a/internal/test_unit/config/GridStatConfig_GRIB_set_attr +++ b/internal/test_unit/config/GridStatConfig_GRIB_set_attr @@ -209,6 +209,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; }; // @@ -228,6 +229,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_GTG_latlon b/internal/test_unit/config/GridStatConfig_GTG_latlon index afb9706715..115832878f 100644 --- a/internal/test_unit/config/GridStatConfig_GTG_latlon +++ b/internal/test_unit/config/GridStatConfig_GTG_latlon @@ -188,6 +188,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -207,6 +208,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_GTG_lc b/internal/test_unit/config/GridStatConfig_GTG_lc index c8a126fe6d..baed1ceae9 100644 --- a/internal/test_unit/config/GridStatConfig_GTG_lc +++ b/internal/test_unit/config/GridStatConfig_GTG_lc @@ -188,6 +188,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -207,6 +208,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_SEEPS b/internal/test_unit/config/GridStatConfig_SEEPS new file mode 100644 index 0000000000..5bb5d1d514 --- /dev/null +++ b/internal/test_unit/config/GridStatConfig_SEEPS @@ -0,0 +1,221 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// Grid-Stat 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"; + +// +// Output observation type to be written +// +obtype = "ANALYS"; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Verification grid +// +//regrid = { +// to_grid = ${TO_GRID}; +// method = BUDGET; +// width = 2; +// vld_thresh = 0.5; +// shape = SQUARE; +//} + +//////////////////////////////////////////////////////////////////////////////// + +// +// May be set separately in each "field" entry +// +censor_thresh = []; +censor_val = []; +mpr_column = []; +mpr_thresh = []; +cat_thresh = []; +cnt_thresh = [ NA ]; +cnt_logic = UNION; +wind_thresh = [ NA ]; +wind_logic = UNION; +eclv_points = 0.05; +nc_pairs_var_name = ""; +nc_pairs_var_suffix = ""; +hss_ec_value = NA; +rank_corr_flag = FALSE; + +// +// Forecast and observation fields to be verified +// +fcst = { + field = [ + { + name = "precipitation_amount"; + level = [ "(*,*)" ]; + is_precipitation = TRUE; + } + ]; +} +obs = fcst; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Verification masking regions +// +mask = { + grid = [ "FULL" ]; + poly = []; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Confidence interval settings +// +ci_alpha = [ 0.05 ]; + +boot = { + interval = PCTILE; + rep_prop = 1.0; + n_rep = 0; + rng = "mt19937"; + seed = "1"; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Data smoothing methods +// +interp = { + field = BOTH; + vld_thresh = 1.0; + shape = SQUARE; + + type = [ + { + method = NEAREST; + width = 1; + } + ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Neighborhood methods +// +nbrhd = { + width = [ 1 ]; + cov_thresh = [ >=0.5 ]; + vld_thresh = 1.0; + shape = SQUARE; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Fourier decomposition +// +fourier = { + wave_1d_beg = []; + wave_1d_end = []; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Gradient statistics +// May be set separately in each "obs.field" entry +// +gradient = { + dx = [ 1 ]; + dy = [ 1 ]; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Distance Map statistics +// May be set separately in each "obs.field" entry +// +distance_map = { + baddeley_p = 2; + baddeley_max_dist = NA; + fom_alpha = 0.1; + zhu_weight = 0.5; + beta_value(n) = n * n / 2.0; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Statistical output types +// +output_flag = { + fho = NONE; + ctc = NONE; + cts = NONE; + mctc = NONE; + mcts = NONE; + cnt = STAT; + sl1l2 = STAT; + sal1l2 = NONE; + vl1l2 = NONE; + val1l2 = NONE; + vcnt = NONE; + pct = NONE; + pstd = NONE; + pjc = NONE; + prc = NONE; + eclv = NONE; + nbrctc = NONE; + nbrcts = NONE; + nbrcnt = NONE; + grad = NONE; + dmap = NONE; + seeps = ${SEEPS_FLAG}; +} + +// +// NetCDF matched pairs output file +// +nc_pairs_flag = { + latlon = TRUE; + raw = TRUE; + diff = TRUE; + climo = TRUE; + climo_cdp = FALSE; + weight = FALSE; + nbrhd = FALSE; + fourier = FALSE; + gradient = FALSE; + distance_map = FALSE; + apply_mask = TRUE; +} + +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = ${SEEPS_P1_THRESH}; + +//////////////////////////////////////////////////////////////////////////////// + +grid_weight_flag = NONE; +tmp_dir = "/tmp"; +output_prefix = "${OUTPUT_PREFIX}"; +version = "V11.0.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/GridStatConfig_apply_mask b/internal/test_unit/config/GridStatConfig_apply_mask index c3af968573..462a11b57e 100644 --- a/internal/test_unit/config/GridStatConfig_apply_mask +++ b/internal/test_unit/config/GridStatConfig_apply_mask @@ -189,6 +189,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -208,6 +209,11 @@ nc_pairs_flag = { apply_mask = ${APPLY_MASK}; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_climo_WMO b/internal/test_unit/config/GridStatConfig_climo_WMO index c4cedceac9..b991c18ade 100644 --- a/internal/test_unit/config/GridStatConfig_climo_WMO +++ b/internal/test_unit/config/GridStatConfig_climo_WMO @@ -250,6 +250,7 @@ output_flag = { nbrcnt = NONE; grad = STAT; dmap = NONE; + seeps = NONE; } // @@ -269,6 +270,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = COS_LAT; diff --git a/internal/test_unit/config/GridStatConfig_climo_prob b/internal/test_unit/config/GridStatConfig_climo_prob index 15ff20258f..c2646559fa 100644 --- a/internal/test_unit/config/GridStatConfig_climo_prob +++ b/internal/test_unit/config/GridStatConfig_climo_prob @@ -260,6 +260,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -279,6 +280,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = COS_LAT; diff --git a/internal/test_unit/config/GridStatConfig_fbias_perc_thresh b/internal/test_unit/config/GridStatConfig_fbias_perc_thresh index eebddc15fe..eaaa322f39 100644 --- a/internal/test_unit/config/GridStatConfig_fbias_perc_thresh +++ b/internal/test_unit/config/GridStatConfig_fbias_perc_thresh @@ -187,6 +187,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -194,6 +195,11 @@ output_flag = { // nc_pairs_flag = FALSE; +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_fourier b/internal/test_unit/config/GridStatConfig_fourier index f328963f4d..998bd79822 100644 --- a/internal/test_unit/config/GridStatConfig_fourier +++ b/internal/test_unit/config/GridStatConfig_fourier @@ -215,6 +215,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -234,6 +235,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_grid_weight b/internal/test_unit/config/GridStatConfig_grid_weight index f718d258c3..2e3d41fb01 100644 --- a/internal/test_unit/config/GridStatConfig_grid_weight +++ b/internal/test_unit/config/GridStatConfig_grid_weight @@ -200,6 +200,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -219,6 +220,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = ${GRID_WEIGHT}; diff --git a/internal/test_unit/config/GridStatConfig_interp_shape b/internal/test_unit/config/GridStatConfig_interp_shape index 51e0a5ae06..7a9697113e 100644 --- a/internal/test_unit/config/GridStatConfig_interp_shape +++ b/internal/test_unit/config/GridStatConfig_interp_shape @@ -182,6 +182,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -202,6 +203,11 @@ nc_pairs_flag = { } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_mpr_thresh b/internal/test_unit/config/GridStatConfig_mpr_thresh index c266aac664..0bdf102be1 100644 --- a/internal/test_unit/config/GridStatConfig_mpr_thresh +++ b/internal/test_unit/config/GridStatConfig_mpr_thresh @@ -248,6 +248,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -267,6 +268,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = COS_LAT; diff --git a/internal/test_unit/config/GridStatConfig_no_leap b/internal/test_unit/config/GridStatConfig_no_leap index 9e4e39de3b..a2fd732474 100644 --- a/internal/test_unit/config/GridStatConfig_no_leap +++ b/internal/test_unit/config/GridStatConfig_no_leap @@ -189,6 +189,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -208,6 +209,11 @@ nc_pairs_flag = { apply_mask = ${APPLY_MASK}; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_prob_as_scalar b/internal/test_unit/config/GridStatConfig_prob_as_scalar index 29152d64d1..c552888599 100644 --- a/internal/test_unit/config/GridStatConfig_prob_as_scalar +++ b/internal/test_unit/config/GridStatConfig_prob_as_scalar @@ -210,6 +210,7 @@ output_flag = { nbrcnt = STAT; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -229,6 +230,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_python b/internal/test_unit/config/GridStatConfig_python index 933d9ded81..37c1c090e9 100644 --- a/internal/test_unit/config/GridStatConfig_python +++ b/internal/test_unit/config/GridStatConfig_python @@ -186,6 +186,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -205,6 +206,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_python_mixed b/internal/test_unit/config/GridStatConfig_python_mixed index 0968c0f978..352c0ca89c 100644 --- a/internal/test_unit/config/GridStatConfig_python_mixed +++ b/internal/test_unit/config/GridStatConfig_python_mixed @@ -194,6 +194,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -213,6 +214,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_rtma b/internal/test_unit/config/GridStatConfig_rtma index d3ed0cf0ec..2b7b7cafb1 100644 --- a/internal/test_unit/config/GridStatConfig_rtma +++ b/internal/test_unit/config/GridStatConfig_rtma @@ -190,6 +190,7 @@ output_flag = { nbrcnt = BOTH; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -209,6 +210,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_rtma_perc_thresh b/internal/test_unit/config/GridStatConfig_rtma_perc_thresh index a487d3be1f..02caac5087 100644 --- a/internal/test_unit/config/GridStatConfig_rtma_perc_thresh +++ b/internal/test_unit/config/GridStatConfig_rtma_perc_thresh @@ -193,6 +193,7 @@ output_flag = { nbrcnt = BOTH; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -212,6 +213,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_st4 b/internal/test_unit/config/GridStatConfig_st4 index 3b90e8d9da..d4ccf23ea2 100644 --- a/internal/test_unit/config/GridStatConfig_st4 +++ b/internal/test_unit/config/GridStatConfig_st4 @@ -194,6 +194,7 @@ output_flag = { nbrcnt = NONE; grad = BOTH; dmap = NONE; + seeps = NONE; } // @@ -213,6 +214,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/GridStatConfig_st4_censor b/internal/test_unit/config/GridStatConfig_st4_censor index fc023edc06..c07365420a 100644 --- a/internal/test_unit/config/GridStatConfig_st4_censor +++ b/internal/test_unit/config/GridStatConfig_st4_censor @@ -203,6 +203,7 @@ output_flag = { nbrcnt = NONE; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -222,6 +223,11 @@ nc_pairs_flag = { apply_mask = FALSE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_APCP b/internal/test_unit/config/PointStatConfig_APCP index d777ae8583..8da47df1a6 100644 --- a/internal/test_unit/config/PointStatConfig_APCP +++ b/internal/test_unit/config/PointStatConfig_APCP @@ -131,6 +131,11 @@ output_flag = { seeps_mpr = ${SEEPS_FLAG}; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = ${SEEPS_P1_THRESH}; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_APCP_HIRA b/internal/test_unit/config/PointStatConfig_APCP_HIRA index f45252252a..b7755d475f 100644 --- a/internal/test_unit/config/PointStatConfig_APCP_HIRA +++ b/internal/test_unit/config/PointStatConfig_APCP_HIRA @@ -133,6 +133,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_GTG_latlon b/internal/test_unit/config/PointStatConfig_GTG_latlon index 0c7bf07e41..2f7a2ca0c1 100644 --- a/internal/test_unit/config/PointStatConfig_GTG_latlon +++ b/internal/test_unit/config/PointStatConfig_GTG_latlon @@ -152,6 +152,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_GTG_lc b/internal/test_unit/config/PointStatConfig_GTG_lc index 67cb871303..0d29146ca9 100644 --- a/internal/test_unit/config/PointStatConfig_GTG_lc +++ b/internal/test_unit/config/PointStatConfig_GTG_lc @@ -160,6 +160,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_INTERP_OPTS b/internal/test_unit/config/PointStatConfig_INTERP_OPTS index db5d996c24..07013b55b0 100644 --- a/internal/test_unit/config/PointStatConfig_INTERP_OPTS +++ b/internal/test_unit/config/PointStatConfig_INTERP_OPTS @@ -143,6 +143,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_LAND_TOPO_MASK b/internal/test_unit/config/PointStatConfig_LAND_TOPO_MASK index 6b856fcc1e..da96d062a0 100644 --- a/internal/test_unit/config/PointStatConfig_LAND_TOPO_MASK +++ b/internal/test_unit/config/PointStatConfig_LAND_TOPO_MASK @@ -183,6 +183,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_MASK_SID b/internal/test_unit/config/PointStatConfig_MASK_SID index b8c11030cd..5d0aa921ac 100644 --- a/internal/test_unit/config/PointStatConfig_MASK_SID +++ b/internal/test_unit/config/PointStatConfig_MASK_SID @@ -138,6 +138,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_PHYS b/internal/test_unit/config/PointStatConfig_PHYS index d93e5391a7..0c1bd7c2a8 100644 --- a/internal/test_unit/config/PointStatConfig_PHYS +++ b/internal/test_unit/config/PointStatConfig_PHYS @@ -139,6 +139,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_PHYS_pint b/internal/test_unit/config/PointStatConfig_PHYS_pint index d70192a9ae..6ba665fcab 100644 --- a/internal/test_unit/config/PointStatConfig_PHYS_pint +++ b/internal/test_unit/config/PointStatConfig_PHYS_pint @@ -134,6 +134,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_WINDS b/internal/test_unit/config/PointStatConfig_WINDS index 609814fd4f..a990411150 100644 --- a/internal/test_unit/config/PointStatConfig_WINDS +++ b/internal/test_unit/config/PointStatConfig_WINDS @@ -154,6 +154,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_aeronet b/internal/test_unit/config/PointStatConfig_aeronet index bbad04e26c..7da47351fe 100644 --- a/internal/test_unit/config/PointStatConfig_aeronet +++ b/internal/test_unit/config/PointStatConfig_aeronet @@ -203,6 +203,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_airnow b/internal/test_unit/config/PointStatConfig_airnow index 0e4e48d710..0f6844e734 100644 --- a/internal/test_unit/config/PointStatConfig_airnow +++ b/internal/test_unit/config/PointStatConfig_airnow @@ -233,6 +233,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// tmp_dir = "/tmp"; diff --git a/internal/test_unit/config/PointStatConfig_climo b/internal/test_unit/config/PointStatConfig_climo index 1621e61ec5..565463ada5 100644 --- a/internal/test_unit/config/PointStatConfig_climo +++ b/internal/test_unit/config/PointStatConfig_climo @@ -273,6 +273,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry): >=min_p1&&<=max_p1 + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_climo_WMO b/internal/test_unit/config/PointStatConfig_climo_WMO index a703ad8266..cea0081f93 100644 --- a/internal/test_unit/config/PointStatConfig_climo_WMO +++ b/internal/test_unit/config/PointStatConfig_climo_WMO @@ -221,6 +221,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_climo_prob b/internal/test_unit/config/PointStatConfig_climo_prob index 8a924ef106..d8b18fb4c2 100644 --- a/internal/test_unit/config/PointStatConfig_climo_prob +++ b/internal/test_unit/config/PointStatConfig_climo_prob @@ -223,6 +223,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_dup b/internal/test_unit/config/PointStatConfig_dup index efd861241c..4cb4801f87 100644 --- a/internal/test_unit/config/PointStatConfig_dup +++ b/internal/test_unit/config/PointStatConfig_dup @@ -156,6 +156,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = ${DUPLICATE_FLAG}; diff --git a/internal/test_unit/config/PointStatConfig_mpr_thresh b/internal/test_unit/config/PointStatConfig_mpr_thresh index 1a3a33851f..3044babc08 100644 --- a/internal/test_unit/config/PointStatConfig_mpr_thresh +++ b/internal/test_unit/config/PointStatConfig_mpr_thresh @@ -215,6 +215,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_obs_summary b/internal/test_unit/config/PointStatConfig_obs_summary index e2d6bcef3d..394b33c46f 100644 --- a/internal/test_unit/config/PointStatConfig_obs_summary +++ b/internal/test_unit/config/PointStatConfig_obs_summary @@ -145,6 +145,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// rank_corr_flag = FALSE; diff --git a/internal/test_unit/config/PointStatConfig_obs_summary_all b/internal/test_unit/config/PointStatConfig_obs_summary_all index 972fa2dcbe..6d32995f59 100644 --- a/internal/test_unit/config/PointStatConfig_obs_summary_all +++ b/internal/test_unit/config/PointStatConfig_obs_summary_all @@ -214,6 +214,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// rank_corr_flag = FALSE; diff --git a/internal/test_unit/config/PointStatConfig_prob b/internal/test_unit/config/PointStatConfig_prob index 20bd7a217b..33edc81df1 100644 --- a/internal/test_unit/config/PointStatConfig_prob +++ b/internal/test_unit/config/PointStatConfig_prob @@ -141,6 +141,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_python b/internal/test_unit/config/PointStatConfig_python index 8854bf604a..9f9c8c2095 100644 --- a/internal/test_unit/config/PointStatConfig_python +++ b/internal/test_unit/config/PointStatConfig_python @@ -211,6 +211,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// rank_corr_flag = FALSE; diff --git a/internal/test_unit/config/PointStatConfig_qty_inc_exc b/internal/test_unit/config/PointStatConfig_qty_inc_exc index 8dc3f1413a..039008a403 100644 --- a/internal/test_unit/config/PointStatConfig_qty_inc_exc +++ b/internal/test_unit/config/PointStatConfig_qty_inc_exc @@ -200,6 +200,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/PointStatConfig_sid_inc_exc b/internal/test_unit/config/PointStatConfig_sid_inc_exc index e20d929374..2a3ae41828 100644 --- a/internal/test_unit/config/PointStatConfig_sid_inc_exc +++ b/internal/test_unit/config/PointStatConfig_sid_inc_exc @@ -146,6 +146,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// duplicate_flag = NONE; diff --git a/internal/test_unit/config/STATAnalysisConfig_seeps b/internal/test_unit/config/STATAnalysisConfig_seeps new file mode 100644 index 0000000000..f8adae2239 --- /dev/null +++ b/internal/test_unit/config/STATAnalysisConfig_seeps @@ -0,0 +1,118 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// STAT-Analysis configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Filtering input STAT lines by the contents of each column +// +model = []; +desc = []; + +fcst_lead = []; +obs_lead = []; + +fcst_valid_beg = ""; +fcst_valid_end = ""; +fcst_valid_inc = []; +fcst_valid_exc = []; +fcst_valid_hour = []; + +obs_valid_beg = ""; +obs_valid_end = ""; +obs_valid_inc = []; +obs_valid_exc = []; +obs_valid_hour = []; + +fcst_init_beg = ""; +fcst_init_end = ""; +fcst_init_inc = []; +fcst_init_exc = []; +fcst_init_hour = []; + +obs_init_beg = ""; +obs_init_end = ""; +obs_init_inc = []; +obs_init_exc = []; +obs_init_hour = []; + +fcst_var = []; +obs_var = []; + +fcst_lev = []; +obs_lev = []; + +obtype = []; + +vx_mask = []; + +interp_mthd = []; + +interp_pnts = []; + +fcst_thresh = []; +obs_thresh = []; +cov_thresh = []; + +alpha = []; + +line_type = []; + +column = []; + +weight = []; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Array of STAT-Analysis jobs to be performed on the filtered data +// +// Job 1 = filter SEEPS lines +// Job 2 = aggregate SEEPS_MPR lines by interpolation (output equals Job 1) +// Job 3 = aggregate SEEPS lines +// Job 4 = aggregate all SEEPS_MPR lines (output equals Job 3) +// Job 5 = summarize SEEPS_MPR scores +// +jobs = [ + "-job filter -line_type SEEPS -dump_row ${OUTPUT_DIR}/CONFIG_POINT_STAT_filter_seeps.stat", + "-job aggregate_stat -line_type SEEPS_MPR -out_line_type SEEPS -by INTERP_MTHD,INTERP_PNTS", + "-job aggregate -line_type SEEPS", + "-job aggregate_stat -line_type SEEPS_MPR -out_line_type SEEPS", + "-job summary -line_type SEEPS_MPR -column SEEPS -by INTERP_MTHD,INTERP_PNTS" +]; + +//////////////////////////////////////////////////////////////////////////////// + +// +// Confidence interval settings +// +out_alpha = 0.05; + +boot = { + interval = PCTILE; + rep_prop = 1.0; + n_rep = 1000; + rng = "mt19937"; + seed = "1"; +} + +//////////////////////////////////////////////////////////////////////////////// + +// +// Skill score index options +// +ss_index_name = "SS_INDEX"; +ss_index_vld_thresh = 1.0; + +//////////////////////////////////////////////////////////////////////////////// + +hss_ec_value = NA; +rank_corr_flag = TRUE; +vif_flag = FALSE; +tmp_dir = "/tmp"; +version = "V11.0.0"; + +//////////////////////////////////////////////////////////////////////////////// diff --git a/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS b/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS index a2ae20816b..b3018adf16 100644 --- a/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS +++ b/internal/test_unit/config/TCPairsConfig_DIAGNOSTICS @@ -131,7 +131,22 @@ watch_warn = { // // Diagnostics to be extracted // -diag_name = [ ${DIAG_NAME} ]; +diag_info_map = [ + { + diag_source = "CIRA_DIAG_RT"; + track_source = "GFS"; + field_source = "GFS_0p50"; + match_to_track = [ "GFS" ]; + diag_name = [ ${CIRA_RT_DIAG_NAME} ]; + }, + { + diag_source = "SHIPS_DIAG_RT"; + track_source = "OFCL"; + field_source = "GFS_0p50"; + match_to_track = [ "OFCL", "SHIP" ]; + diag_name = [ ${SHIPS_RT_DIAG_NAME} ]; + } +]; // // Unit conversions to be applied based on diagnostic names and units diff --git a/internal/test_unit/hdr/met_11_0.hdr b/internal/test_unit/hdr/met_11_0.hdr index 573457d898..4c496af7f9 100644 --- a/internal/test_unit/hdr/met_11_0.hdr +++ b/internal/test_unit/hdr/met_11_0.hdr @@ -33,6 +33,6 @@ SSIDX : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_L MODE_SOA : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE N_VALID GRID_RES OBJECT_ID OBJECT_CAT CENTROID_X CENTROID_Y CENTROID_LAT CENTROID_LON AXIS_ANG LENGTH WIDTH AREA AREA_THRESH CURVATURE CURVATURE_X CURVATURE_Y COMPLEXITY INTENSITY_10 INTENSITY_25 INTENSITY_50 INTENSITY_75 INTENSITY_90 INTENSITY_50 INTENSITY_SUM MODE_POA : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE N_VALID GRID_RES OBJECT_ID OBJECT_CAT CENTROID_DIST BOUNDARY_DIST CONVEX_HULL_DIST ANGLE_DIFF ASPECT_DIFF AREA_RATIO INTERSECTION_AREA UNION_AREA SYMMETRIC_DIFF INTERSECTION_OVER_AREA CURVATURE_RATIO COMPLEXITY_RATIO PERCENTILE_INTENSITY_RATIO INTEREST MODE_CTS : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE N_VALID GRID_RES FIELD TOTAL FY_OY FY_ON FN_OY FN_ON BASER FMEAN ACC FBIAS PODY PODN POFD FAR CSI GSS HK HSS ODDS -TCST_TCMPR : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE AMODEL BMODEL STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK TOTAL INDEX LEVEL WATCH_WARN INITIALS ALAT ALON BLAT BLON TK_ERR X_ERR Y_ERR ALTK_ERR CRTK_ERR ADLAND BDLAND AMSLP BMSLP AMAX_WIND BMAX_WIND AAL_WIND_34 BAL_WIND_34 ANE_WIND_34 BNE_WIND_34 ASE_WIND_34 BSE_WIND_34 ASW_WIND_34 BSW_WIND_34 ANW_WIND_34 BNW_WIND_34 AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH NUM_MEMBERS TRACK_SPREAD DIST_MEAN MSLP_SPREAD MAX_WIND_SPREAD -TCST_TCDIAG : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE AMODEL BMODEL STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK TOTAL INDEX DIAG_SOURCE N_DIAG _VAR_ +TCST_TCMPR : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE AMODEL BMODEL STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK TOTAL INDEX LEVEL WATCH_WARN INITIALS ALAT ALON BLAT BLON TK_ERR X_ERR Y_ERR ALTK_ERR CRTK_ERR ADLAND BDLAND AMSLP BMSLP AMAX_WIND BMAX_WIND AAL_WIND_34 BAL_WIND_34 ANE_WIND_34 BNE_WIND_34 ASE_WIND_34 BSE_WIND_34 ASW_WIND_34 BSW_WIND_34 ANW_WIND_34 BNW_WIND_34 AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH NUM_MEMBERS TRACK_SPREAD TRACK_STDEV MSLP_STDEV MAX_WIND_STDEV +TCST_TCDIAG : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE AMODEL BMODEL STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK TOTAL INDEX DIAG_SOURCE TRACK_SOURCE FIELD_SOURCE N_DIAG _VAR_ TCST_PROBRIRW : VERSION MODEL DESC FCST_LEAD FCST_VALID_BEG FCST_VALID_END OBS_LEAD OBS_VALID_BEG OBS_VALID_END FCST_VAR FCST_UNITS FCST_LEV OBS_VAR OBS_UNITS OBS_LEV OBTYPE VX_MASK INTERP_MTHD INTERP_PNTS FCST_THRESH OBS_THRESH COV_THRESH ALPHA LINE_TYPE ALAT ALON BLAT BLON INITIALS TK_ERR X_ERR Y_ERR ADLAND BDLAND RI_BEG RI_END RI_WINDOW AWIND_END BWIND_BEG BWIND_END BDELTA BDELTA_MAX BLEVEL_BEG BLEVEL_END N_THRESH _VAR_ diff --git a/internal/test_unit/perl/tcst_conv.pl b/internal/test_unit/perl/tcst_conv.pl index d653e9cb4b..f1b9fd7122 100755 --- a/internal/test_unit/perl/tcst_conv.pl +++ b/internal/test_unit/perl/tcst_conv.pl @@ -30,10 +30,10 @@ () AAL_WIND_50 BAL_WIND_50 ANE_WIND_50 BNE_WIND_50 ASE_WIND_50 BSE_WIND_50 ASW_WIND_50 BSW_WIND_50 ANW_WIND_50 BNW_WIND_50 AAL_WIND_64 BAL_WIND_64 ANE_WIND_64 BNE_WIND_64 ASE_WIND_64 BSE_WIND_64 ASW_WIND_64 BSW_WIND_64 ANW_WIND_64 BNW_WIND_64 ARADP BRADP ARRP BRRP AMRD BMRD AGUSTS BGUSTS AEYE BEYE ADIR BDIR ASPEED BSPEED ADEPTH BDEPTH - NUM_MEMBERS TRACK_SPREAD DIST_MEAN MSLP_SPREAD MAX_WIND_SPREAD); + NUM_MEMBERS TRACK_SPREAD TRACK_STDEV MSLP_STDEV MAX_WIND_STDEV); my @fld_tcdiag = qw(AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK - TOTAL INDEX LEVEL DIAG_SOURCE N_DIAG DIAG_ VALUE_); + TOTAL INDEX LEVEL DIAG_SOURCE TRACK_SOURCE FIELD_SOURCE N_DIAG DIAG_ VALUE_); my @fld_probrirw = qw(AMODEL BMODEL DESC STORM_ID BASIN CYCLONE STORM_NAME INIT_MASK VALID_MASK ALAT ALON BLAT BLON INITIALS TK_ERR X_ERR Y_ERR ADLAND BDLAND RI_BEG RI_END RI_WINDOW @@ -143,9 +143,9 @@ () "%15s" . # BDEPTH "%15s" . # NUM_MEMBERS "%15s" . # TRACK_SPREAD - "%15s" . # DIST_MEAN - "%15s" . # MSLP_SPREAD - "%15s"; # MAX_WIND_SPREAD + "%15s" . # TRACK_STDEV + "%15s" . # MSLP_STDEV + "%15s"; # MAX_WIND_STDEV my $fmt_tcdiag = "%15s" . # AMODEL @@ -159,6 +159,8 @@ () "%15s" . # TOTAL "%15s" . # INDEX "%15s" . # DIAG_SOURCE + "%15s" . # TRACK_SOURCE + "%15s" . # FIELD_SOURCE "%15s"; # N_DIAG my $fmt_probrirw = @@ -353,17 +355,19 @@ () # 79 - BDEPTH # 80 - NUM_MEMBERS # 81 - TRACK_SPREAD -# 82 - DIST_MEAN -# 83 - MSLP_SPREAD -# 84 - MAX_WIND_SPREAD +# 82 - TRACK_STDEV +# 83 - MSLP_STDEV +# 84 - MAX_WIND_STDEV # TCDIAG Line Type # 14 - TOTAL # 15 - INDEX # 16 - DIAG_SOURCE -# 17 - N_DIAG -# 18 - DIAG_i -# 19 - VALUE_i +# 17 - TRACK_SOURCE +# 18 - FIELD_SOURCE +# 19 - N_DIAG +# 20 - DIAG_i +# 21 - VALUE_i # PROBRIRW Line Type # 14 - ALAT diff --git a/internal/test_unit/xml/unit_grid_stat.xml b/internal/test_unit/xml/unit_grid_stat.xml index 9d154a91a6..204e934230 100644 --- a/internal/test_unit/xml/unit_grid_stat.xml +++ b/internal/test_unit/xml/unit_grid_stat.xml @@ -289,4 +289,26 @@ + + &MET_BIN;/grid_stat + + OMP_NUM_THREADS 2 + SEEPS_FLAG BOTH + SEEPS_P1_THRESH NA + OUTPUT_PREFIX SEEPS + MET_SEEPS_GRID_CLIMO_NAME&DATA_DIR_CLIMO;/seeps/PPT24_seepsweights_grid.nc + + \ + &DATA_DIR_MODEL;/seeps/gpm_2021120100_2021120200_trmmgrid.nc \ + &DATA_DIR_MODEL;/seeps/prods_op_gl-mn_2021120100_f024_A24_trmmgrid.nc \ + &CONFIG_DIR;/GridStatConfig_SEEPS \ + -outdir &OUTPUT_DIR;/grid_stat -v 1 + + + &OUTPUT_DIR;/grid_stat/grid_stat_SEEPS_000000L_20211202_000000V.stat + &OUTPUT_DIR;/grid_stat/grid_stat_SEEPS_000000L_20211202_000000V_seeps.txt + &OUTPUT_DIR;/grid_stat/grid_stat_SEEPS_000000L_20211202_000000V_pairs.nc + + + diff --git a/internal/test_unit/xml/unit_point_stat.xml b/internal/test_unit/xml/unit_point_stat.xml index 8eb398fa66..b57d756a14 100644 --- a/internal/test_unit/xml/unit_point_stat.xml +++ b/internal/test_unit/xml/unit_point_stat.xml @@ -167,6 +167,7 @@ FCST_FIELD_LEVEL A3 OBS_DICT fcst SEEPS_FLAG NONE + SEEPS_P1_THRESH NA OUTPUT_PREFIX GRIB1_NAM_TRMM \ @@ -194,6 +195,7 @@ FCST_FIELD_LEVEL A3 OBS_DICT fcst SEEPS_FLAG NONE + SEEPS_P1_THRESH NA OUTPUT_PREFIX GRIB2_SREF_TRMM \ @@ -221,6 +223,7 @@ FCST_FIELD_LEVEL (*,*) OBS_DICT { field = [ { name = "APCP"; level = "A24"; } ]; } SEEPS_FLAG NONE + SEEPS_P1_THRESH NA OUTPUT_PREFIX NCMET_NAM_HMTGAGE \ @@ -248,6 +251,7 @@ FCST_FIELD_LEVEL (*,*) OBS_DICT { field = [ { name = "TP24"; level = "L0"; is_precipitation = TRUE; } ]; } SEEPS_FLAG BOTH + SEEPS_P1_THRESH ge0.1&&le0.85 OUTPUT_PREFIX NCMET_NAM_NDAS_SEEPS \ @@ -277,6 +281,7 @@ FCST_FIELD_LEVEL (0,*,*) OBS_DICT { field = [ { name = "APCP"; level = "A3"; } ]; } SEEPS_FLAG NONE + SEEPS_P1_THRESH NA OUTPUT_PREFIX NCPINT_TRMM \ diff --git a/internal/test_unit/xml/unit_stat_analysis_ps.xml b/internal/test_unit/xml/unit_stat_analysis_ps.xml index dba20346c8..87884061ea 100644 --- a/internal/test_unit/xml/unit_stat_analysis_ps.xml +++ b/internal/test_unit/xml/unit_stat_analysis_ps.xml @@ -100,6 +100,23 @@ + + + OUTPUT_DIR &OUTPUT_DIR;/stat_analysis_ps + + &MET_BIN;/stat_analysis + \ + -lookin &OUTPUT_DIR;/point_stat/point_stat_NCMET_NAM_NDAS_SEEPS_360000L_20120410_120000V.stat \ + -config &CONFIG_DIR;/STATAnalysisConfig_seeps \ + -out &OUTPUT_DIR;/stat_analysis_ps/POINT_STAT_SEEPS.out \ + -v 1 + + + &OUTPUT_DIR;/stat_analysis_ps/CONFIG_POINT_STAT_filter_seeps.stat + &OUTPUT_DIR;/stat_analysis_ps/POINT_STAT_SEEPS.out + + + OUTPUT_DIR &OUTPUT_DIR;/stat_analysis_ps diff --git a/internal/test_unit/xml/unit_tc_pairs.xml b/internal/test_unit/xml/unit_tc_pairs.xml index 0080eb53b1..24bec42c00 100644 --- a/internal/test_unit/xml/unit_tc_pairs.xml +++ b/internal/test_unit/xml/unit_tc_pairs.xml @@ -188,13 +188,14 @@ STORM_ID "AL092022" INIT_BEG "20220926_00" INIT_END "20220926_18" - DIAG_NAME "DTL", "PW01", "SHRD", "TPW", "LAND", "SHR_MAG", "STM_SPD" + CIRA_RT_DIAG_NAME "TPW", "LAND", "SHR_MAG", "STM_SPD" + SHIPS_RT_DIAG_NAME "DTL", "PW01", "SHRD" \ -adeck &DATA_DIR;/adeck/aal092022_OFCL_SHIP_AVNO.dat \ -bdeck &DATA_DIR;/bdeck/bal092022.dat \ - -diag TCDIAG &DATA_DIR;/tcdiag/2022/sal092022_avno_doper_20220926*_diag.dat \ - -diag LSDIAG_RT &DATA_DIR;/lsdiag_rt/2022/220926*AL0922_lsdiag.dat model=OFCL,SHIP \ + -diag CIRA_DIAG_RT &DATA_DIR;/diag/cira_diag_rt/2022/sal092022_avno_doper_20220926*_diag.dat \ + -diag SHIPS_DIAG_RT &DATA_DIR;/diag/ships_diag_rt/2022/220926*AL0922_lsdiag.dat \ -config &CONFIG_DIR;/TCPairsConfig_DIAGNOSTICS \ -out &OUTPUT_DIR;/tc_pairs/al092022_20220926_DIAGNOSTICS \ -log &OUTPUT_DIR;/tc_pairs/tc_pairs_DIAGNOSTICS.log \ diff --git a/scripts/config/GridStatConfig_APCP_12 b/scripts/config/GridStatConfig_APCP_12 index 12213ba271..5779600b7b 100644 --- a/scripts/config/GridStatConfig_APCP_12 +++ b/scripts/config/GridStatConfig_APCP_12 @@ -189,6 +189,7 @@ output_flag = { nbrcnt = BOTH; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -208,6 +209,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/scripts/config/GridStatConfig_APCP_24 b/scripts/config/GridStatConfig_APCP_24 index 7b9b82c7a4..a925c2e872 100644 --- a/scripts/config/GridStatConfig_APCP_24 +++ b/scripts/config/GridStatConfig_APCP_24 @@ -199,6 +199,7 @@ output_flag = { nbrcnt = BOTH; grad = NONE; dmap = NONE; + seeps = NONE; } // @@ -218,6 +219,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/scripts/config/GridStatConfig_POP_12 b/scripts/config/GridStatConfig_POP_12 index 03cdd78e64..86f92a8220 100644 --- a/scripts/config/GridStatConfig_POP_12 +++ b/scripts/config/GridStatConfig_POP_12 @@ -218,6 +218,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/scripts/config/GridStatConfig_all b/scripts/config/GridStatConfig_all index a782bcfb67..ae38f06cd5 100644 --- a/scripts/config/GridStatConfig_all +++ b/scripts/config/GridStatConfig_all @@ -230,6 +230,7 @@ output_flag = { nbrcnt = BOTH; grad = BOTH; dmap = BOTH; + seeps = NONE; } // @@ -249,6 +250,11 @@ nc_pairs_flag = { apply_mask = TRUE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// grid_weight_flag = NONE; diff --git a/scripts/config/PointStatConfig b/scripts/config/PointStatConfig index b4d26b94b7..6b8a88bf75 100644 --- a/scripts/config/PointStatConfig +++ b/scripts/config/PointStatConfig @@ -208,6 +208,11 @@ output_flag = { seeps_mpr = NONE; } +//////////////////////////////////////////////////////////////////////////////// +// Threshold for SEEPS p1 (Probability of being dry) + +seeps_p1_thresh = NA; + //////////////////////////////////////////////////////////////////////////////// rank_corr_flag = TRUE; diff --git a/src/basic/vx_config/config_constants.h b/src/basic/vx_config/config_constants.h index 1141e2c75e..fb49e5cce7 100644 --- a/src/basic/vx_config/config_constants.h +++ b/src/basic/vx_config/config_constants.h @@ -90,14 +90,28 @@ enum TrackType { // enum DiagType { - DiagType_None, // Default - TCDiagType, // Tropical Cyclone Diagnostics - LSDiagRTType, // Realtime Large Scale Diagnostics - LSDiagDevType // Development Large Scale Diagnostics + DiagType_None, // Default + DiagType_CIRA_RT, // Realtime CIRA Tropical Cyclone Diagnostics + DiagType_CIRA_Dev, // Developmental CIRA Tropical Cyclone Diagnostics + DiagType_SHIPS_RT, // Realtime SHIPS Large Scale Diagnostics + DiagType_SHIPS_Dev // Developmental SHIPS Large Scale Diagnostics }; //////////////////////////////////////////////////////////////////////// +// +// Corresponding diagnostic type strings +// + +static const char cira_diag_str[] = "CIRA_DIAG"; +static const char cira_diag_rt_str[] = "CIRA_DIAG_RT"; +static const char cira_diag_dev_str[] = "CIRA_DIAG_DEV"; +static const char ships_diag_str[] = "SHIPS_DIAG"; +static const char ships_diag_rt_str[] = "SHIPS_DIAG_RT"; +static const char ships_diag_dev_str[] = "SHIPS_DIAG_DEV"; + +//////////////////////////////////////////////////////////////////////// + // // Enumeration for 12-hour interpolation logic // @@ -634,6 +648,7 @@ static const char conf_key_raw_flag[] = "raw"; static const char conf_key_diff_flag[] = "diff"; static const char conf_key_climo_flag[] = "climo"; static const char conf_key_climo_cdp_flag[] = "climo_cdp"; +static const char conf_key_seeps_flag[] = "seeps"; static const char conf_key_apply_mask_flag[] = "apply_mask"; static const char conf_key_object_raw_flag[] = "object_raw"; static const char conf_key_object_id_flag[] = "object_id"; @@ -654,6 +669,7 @@ static const char conf_key_obs_to_qc_map[] = "obs_to_qc_map"; static const char conf_key_missing_thresh[] = "missing_thresh"; static const char conf_key_control_id[] = "control_id"; static const char conf_key_ens_member_ids[] = "ens_member_ids"; +static const char conf_key_seeps_p1_thresh[] = "seeps_p1_thresh"; // // Entries to override file metadata @@ -1081,9 +1097,13 @@ static const char conf_key_dland_file[] = "dland_file"; static const char conf_key_basin_file[] = "basin_file"; static const char conf_key_track_watch_warn[] = "track_watch_warn"; static const char conf_key_watch_warn[] = "watch_warn"; +static const char conf_key_diag_info_map[] = "diag_info_map"; +static const char conf_key_diag_source[] = "diag_source"; +static const char conf_key_track_source[] = "track_source"; +static const char conf_key_field_source[] = "field_source"; +static const char conf_key_match_to_track[] = "match_to_track"; static const char conf_key_diag_name[] = "diag_name"; static const char conf_key_diag_convert_map[] = "diag_convert_map"; -static const char conf_key_source[] = "source"; static const char conf_key_basin_map[] = "basin_map"; static const char conf_key_time_offset[] = "time_offset"; static const char conf_key_amodel[] = "amodel"; diff --git a/src/basic/vx_config/config_util.cc b/src/basic/vx_config/config_util.cc index 607844db71..52817f1b6d 100644 --- a/src/basic/vx_config/config_util.cc +++ b/src/basic/vx_config/config_util.cc @@ -1121,7 +1121,7 @@ map parse_conf_key_convert_map( exit(1); } - // Conf: tcdiag_convert_map, lsdiag_convert_map, etc + // Conf: diag_convert_map map_dict = dict->lookup_array(conf_key_map_name); // Loop through the array entries @@ -2757,16 +2757,13 @@ ConcatString tracktype_to_string(TrackType type) { DiagType string_to_diagtype(const char *s) { DiagType t = DiagType_None; - // Convert string to enumerated DiagType - if(strcasecmp(s, conf_val_none) == 0) t = DiagType_None; - else if(strcasecmp(s, "TCDIAG") == 0) t = TCDiagType; - else if(strcasecmp(s, "LSDIAG_RT") == 0) t = LSDiagRTType; - else if(strcasecmp(s, "LSDIAG_DEV") == 0) t = LSDiagDevType; - else { - mlog << Error << "\nstring_to_diagtype() -> " - << "Unexpected DiagType string \"" << s << "\".\n\n"; - exit(1); - } + // Convert string to enumerated DiagType, storing unknown strings + // as the default type + if(strcasecmp(s, cira_diag_rt_str) == 0) t = DiagType_CIRA_RT; + else if(strcasecmp(s, cira_diag_dev_str) == 0) t = DiagType_CIRA_Dev; + else if(strcasecmp(s, ships_diag_rt_str) == 0) t = DiagType_SHIPS_RT; + else if(strcasecmp(s, ships_diag_dev_str) == 0) t = DiagType_SHIPS_Dev; + else t = DiagType_None; return(t); } @@ -2778,10 +2775,11 @@ ConcatString diagtype_to_string(DiagType type) { // Convert enumerated DiagType to string switch(type) { - case(DiagType_None): s = conf_val_none; break; - case(TCDiagType): s = "TCDIAG"; break; - case(LSDiagRTType): s = "LSDIAG_RT"; break; - case(LSDiagDevType): s = "LSDIAG_DEV"; break; + case(DiagType_None): s = conf_val_none; break; + case(DiagType_CIRA_RT): s = cira_diag_rt_str; break; + case(DiagType_CIRA_Dev): s = cira_diag_dev_str; break; + case(DiagType_SHIPS_RT): s = ships_diag_rt_str; break; + case(DiagType_SHIPS_Dev): s = ships_diag_dev_str; break; default: mlog << Error << "\ndiagtype_to_string() -> " << "Unexpected DiagType value of " << type << ".\n\n"; diff --git a/src/libcode/vx_analysis_util/stat_job.cc b/src/libcode/vx_analysis_util/stat_job.cc index f82b0f6151..ef6e5030a4 100644 --- a/src/libcode/vx_analysis_util/stat_job.cc +++ b/src/libcode/vx_analysis_util/stat_job.cc @@ -2223,6 +2223,14 @@ void STATAnalysisJob::dump_stat_line(const STATLine &line) { write_header_row(ssvar_columns, n_ssvar_columns, 1, dump_at, 0, 0); break; + case(stat_seeps): + write_header_row(seeps_columns, n_seeps_columns, 1, dump_at, 0, 0); + break; + + case(stat_seeps_mpr): + write_header_row(seeps_mpr_columns, n_seeps_mpr_columns, 1, dump_at, 0, 0); + break; + // Just write a STAT header line for indeterminant line types case(stat_mctc): case(stat_mcts): diff --git a/src/libcode/vx_seeps/seeps.cc b/src/libcode/vx_seeps/seeps.cc index b0367af7b2..469217251c 100644 --- a/src/libcode/vx_seeps/seeps.cc +++ b/src/libcode/vx_seeps/seeps.cc @@ -27,36 +27,182 @@ using namespace netCDF; bool standalone_debug_seeps = false; -static SeepsClimo *seeps_Climo = 0; +static SeepsClimo *seeps_climo = 0; +static std::map seeps_climo_grid_map_00; static const char *def_seeps_filename = "MET_BASE/climo/seeps/PPT24_seepsweights.nc"; +static const char *def_seeps_grid_filename = "MET_BASE/climo/seeps/PPT24_seepsweights_grid.nc"; static const char *var_name_sid = "sid"; static const char *var_name_lat = "lat"; static const char *var_name_lon = "lon"; static const char *var_name_elv = "elv"; +static const char *dim_name_lat = "latitude"; +static const char *dim_name_lon = "longitude"; + +double weighted_average(double, double, double, double); + //////////////////////////////////////////////////////////////////////// SeepsClimo *get_seeps_climo() { - if (! seeps_Climo) seeps_Climo = new SeepsClimo(); - return seeps_Climo; + if (! seeps_climo) seeps_climo = new SeepsClimo(); + return seeps_climo; } //////////////////////////////////////////////////////////////////////// void release_seeps_climo() { - if (seeps_Climo) { delete seeps_Climo; seeps_Climo = 0; } + if (seeps_climo) { delete seeps_climo; seeps_climo = 0; } +} + +//////////////////////////////////////////////////////////////////////// + +SeepsClimoGrid *get_seeps_climo_grid(int month, int hour) { + bool not_found = true; + SeepsClimoGrid *seeps_climo_grid = NULL; + for (map::iterator it=seeps_climo_grid_map_00.begin(); + it!=seeps_climo_grid_map_00.end(); ++it) { + if (it->first == month) { + not_found = false; + seeps_climo_grid = (SeepsClimoGrid *)it->second; + break; + } + } + + if (not_found) { + seeps_climo_grid = new SeepsClimoGrid(month, hour); + seeps_climo_grid_map_00[month] = seeps_climo_grid; + } + return seeps_climo_grid; +} + +//////////////////////////////////////////////////////////////////////// + +void release_seeps_climo_grid(int month, int hour) { + for (map::iterator it=seeps_climo_grid_map_00.begin(); + it!=seeps_climo_grid_map_00.end(); ++it) { + if (it->first == month) { + delete it->second; + seeps_climo_grid_map_00.erase(it); + break; + } + } } +//////////////////////////////////////////////////////////////////////// + +double weighted_average(double v1, double w1, double v2, double w2) { + return(is_bad_data(v1) || is_bad_data(v2) ? + bad_data_double : + v1 * w1 + v2 * w2); +} + + +//////////////////////////////////////////////////////////////////////// + + +void SeepsAggScore::clear() { + + n_obs = 0; + c12 = c13 = c21 = c23 = c31 = c32 = 0; + s12 = s13 = s21 = s23 = s31 = s32 = 0.; + pv1 = pv2 = pv3 = 0.; + pf1 = pf2 = pf3 = 0.; + mean_fcst = mean_obs = bad_data_double; + weighted_score = score = bad_data_double; + +} + +//////////////////////////////////////////////////////////////////////// + +SeepsAggScore & SeepsAggScore::operator+=(const SeepsAggScore &c) { + + // Check for degenerate case + if(n_obs == 0 && c.n_obs == 0) return(*this); + + // Compute weights + double w1 = (double) n_obs / (n_obs + c.n_obs); + double w2 = (double) c.n_obs / (n_obs + c.n_obs); + + // Increment number of obs + n_obs += c.n_obs; + + // Increment counts + c12 += c.c12; + c13 += c.c13; + c21 += c.c21; + c23 += c.c23; + c31 += c.c31; + c32 += c.c32; + + // Compute weighted averages + s12 = weighted_average(s12, w1, c.s12, w2); + s13 = weighted_average(s13, w1, c.s13, w2); + s21 = weighted_average(s21, w1, c.s21, w2); + s23 = weighted_average(s23, w1, c.s23, w2); + s31 = weighted_average(s31, w1, c.s31, w2); + s32 = weighted_average(s32, w1, c.s32, w2); + + pv1 = weighted_average(pv1, w1, c.pv1, w2); + pv2 = weighted_average(pv2, w1, c.pv2, w2); + pv3 = weighted_average(pv3, w1, c.pv3, w2); + + pf1 = weighted_average(pf1, w1, c.pf1, w2); + pf2 = weighted_average(pf2, w1, c.pf2, w2); + pf3 = weighted_average(pf3, w1, c.pf3, w2); + + mean_fcst = weighted_average(mean_fcst, w1, c.mean_fcst, w2); + mean_obs = weighted_average(mean_obs, w1, c.mean_obs, w2); + + score = weighted_average(score, w1, c.score, w2); + weighted_score = weighted_average(weighted_score, w1, c.weighted_score, w2); + + return(*this); +} + + +//////////////////////////////////////////////////////////////////////// + + +SeepsClimoBase::SeepsClimoBase() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +SeepsClimoBase::~SeepsClimoBase() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +void SeepsClimoBase::clear() { + seeps_ready = false; + filtered_count = 0; + seeps_p1_thresh.clear(); +} + +//////////////////////////////////////////////////////////////////////// + +void SeepsClimoBase::set_p1_thresh(const SingleThresh &p1_thresh) { + seeps_p1_thresh = p1_thresh; +} + + //////////////////////////////////////////////////////////////////////// SeepsClimo::SeepsClimo() { - ConcatString seeps_name= get_seeps_climo_filename(); - if (file_exists(seeps_name.c_str())) { - read_records(seeps_name); + ConcatString seeps_name = get_seeps_climo_filename(); + seeps_ready = file_exists(seeps_name.c_str()); + if (seeps_ready) read_seeps_scores(seeps_name); + else { + mlog << Error << "\nSeepsClimo::SeepsClimo() -> " + << "The SEEPS climo data \"" << seeps_name << "\" is missing!" + << " Turn off SEEPS and SEEPS_MPR to continue\n\n"; + exit(1); } } @@ -69,6 +215,7 @@ SeepsClimo::~SeepsClimo() { //////////////////////////////////////////////////////////////////////// void SeepsClimo::clear() { + SeepsClimoBase::clear(); for (map::iterator it=seeps_score_00_map.begin(); it!=seeps_score_00_map.end(); ++it) { delete it->second; @@ -86,8 +233,8 @@ void SeepsClimo::clear() { //////////////////////////////////////////////////////////////////////// SeepsClimoRecord *SeepsClimo::create_climo_record( - int sid, float lat, float lon, float elv, - float *p1, float *p2, float *t1, float *t2, float *scores) { + int sid, double lat, double lon, double elv, + double *p1, double *p2, double *t1, double *t2, double *scores) { int offset; SeepsClimoRecord *record = new SeepsClimoRecord(); @@ -106,7 +253,8 @@ SeepsClimoRecord *SeepsClimo::create_climo_record( if (standalone_debug_seeps && SAMPLE_STATION_ID == sid) { cout << str_format("\t%2d: %6.3f %6.3f %6.3f %6.3f ", - (idx+1), record->p1[idx], record->p2[idx], record->t1[idx], record->t2[idx]); + (idx+1), record->p1[idx], record->p2[idx], + record->t1[idx], record->t2[idx]); } for (int idx_m=0; idx_msecond; } - if (climo_record) { - - record = new SeepsRecord; - record->sid = climo_record->sid; - record->lat = climo_record->lat; - record->lon = climo_record->lon; - record->elv = climo_record->elv; - record->month = month; - record->p1 = climo_record->p1[month-1]; - record->p2 = climo_record->p2[month-1]; - record->t1 = climo_record->t1[month-1]; - record->t2 = climo_record->t2[month-1]; - for (int idx=0; idxscores[idx] = climo_record->scores[month-1][idx]; + if (NULL != climo_record) { + double p1 = climo_record->p1[month-1]; + if (seeps_p1_thresh.check(p1)) { + record = new SeepsRecord; + record->sid = climo_record->sid; + record->lat = climo_record->lat; + record->lon = climo_record->lon; + record->elv = climo_record->elv; + record->month = month; + record->p1 = climo_record->p1[month-1]; + record->p2 = climo_record->p2[month-1]; + record->t1 = climo_record->t1[month-1]; + record->t2 = climo_record->t2[month-1]; + for (int idx=0; idxscores[idx] = climo_record->scores[month-1][idx]; + } + } + else if (!is_eq(p1, bad_data_double)) { + filtered_count++; + mlog << Debug(7) << method_name << " filtered by threshold p1=" + << climo_record->p1[month-1] <<"\n"; } } } @@ -172,9 +327,7 @@ ConcatString SeepsClimo::get_seeps_climo_filename() { // Use the MET_TMP_DIR environment variable, if set. bool use_env = get_env(MET_ENV_SEEPS_CLIMO_NAME, seeps_filename); - if(use_env) { - seeps_filename = replace_path(seeps_filename); - } + if(use_env) seeps_filename = replace_path(seeps_filename); else seeps_filename = replace_path(def_seeps_filename); if (seeps_ready = file_exists(seeps_filename.c_str())) { @@ -197,9 +350,9 @@ ConcatString SeepsClimo::get_seeps_climo_filename() { //////////////////////////////////////////////////////////////////////// -float SeepsClimo::get_score(int sid, float p_fcst, float p_obs, - int month, int hour) { - float score = (float)bad_data_double; +double SeepsClimo::get_score(int sid, double p_fcst, double p_obs, + int month, int hour) { + double score = bad_data_double; SeepsRecord *record = get_record(sid, month, hour); if (NULL != record) { @@ -216,21 +369,22 @@ float SeepsClimo::get_score(int sid, float p_fcst, float p_obs, //////////////////////////////////////////////////////////////////////// -SeepsScore *SeepsClimo::get_seeps_score(int sid, float p_fcst, - float p_obs, int month, int hour) { +SeepsScore *SeepsClimo::get_seeps_score(int sid, double p_fcst, double p_obs, + int month, int hour) { SeepsScore *score = NULL; SeepsRecord *record = get_record(sid, month, hour); if (NULL != record) { score = new SeepsScore(); score->p1 = record->p1; - score->p2 = record->p1; + score->p2 = record->p2; score->t1 = record->t1; score->t2 = record->t2; score->obs_cat = (p_obs>record->t1)+(p_obs>record->t2); - score->model_cat = (p_fcst>record->t1)+(p_fcst>record->t2); - score->score = record->scores[(score->model_cat*3)+score->obs_cat]; + score->fcst_cat = (p_fcst>record->t1)+(p_fcst>record->t2); + score->s_idx = (score->fcst_cat*3)+score->obs_cat; + score->score = record->scores[score->s_idx]; delete record; } @@ -298,206 +452,504 @@ void SeepsClimo::print_record(SeepsRecord *record, bool with_header) { //////////////////////////////////////////////////////////////////////// -void SeepsClimo::read_records(ConcatString filename) { +void SeepsClimo::read_seeps_scores(ConcatString filename) { clock_t clock_time = clock(); - float p1_00_buf[SEEPS_MONTH]; - float p2_00_buf[SEEPS_MONTH]; - float t1_00_buf[SEEPS_MONTH]; - float t2_00_buf[SEEPS_MONTH]; - float p1_12_buf[SEEPS_MONTH]; - float p2_12_buf[SEEPS_MONTH]; - float t1_12_buf[SEEPS_MONTH]; - float t2_12_buf[SEEPS_MONTH]; - float matrix_00_buf[SEEPS_MONTH*SEEPS_MATRIX_SIZE]; - float matrix_12_buf[SEEPS_MONTH*SEEPS_MATRIX_SIZE]; - NcFile *nc_file = open_ncfile(filename.c_str()); const char *method_name = "SeepsClimo::read_records() -> "; - // dimensions: month = 12 ; nstn = 5293 ; nmatrix = 9 ; - get_dim(nc_file, dim_name_nstn, nstn, true); - mlog << Debug(3) << method_name << "dimension nstn = " << nstn << "\n"; - if (standalone_debug_seeps) cout << "dimension nstn = " << nstn << "\n"; - - int *sid_array = new int[nstn]; - float *lat_array = new float[nstn]; - float *lon_array = new float[nstn]; - float *elv_array = new float[nstn]; - float *p1_00_array = new float[nstn*SEEPS_MONTH]; - float *p2_00_array = new float[nstn*SEEPS_MONTH]; - float *t1_00_array = new float[nstn*SEEPS_MONTH]; - float *t2_00_array = new float[nstn*SEEPS_MONTH]; - float *p1_12_array = new float[nstn*SEEPS_MONTH]; - float *p2_12_array = new float[nstn*SEEPS_MONTH]; - float *t1_12_array = new float[nstn*SEEPS_MONTH]; - float *t2_12_array = new float[nstn*SEEPS_MONTH]; - float *matrix_00_array = new float[nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE]; - float *matrix_12_array = new float[nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE]; - - NcVar var_sid = get_nc_var(nc_file, var_name_sid); - NcVar var_lat = get_nc_var(nc_file, var_name_lat); - NcVar var_lon = get_nc_var(nc_file, var_name_lon); - NcVar var_elv = get_nc_var(nc_file, var_name_elv); - NcVar var_p1_00 = get_nc_var(nc_file, var_name_p1_00); - NcVar var_p2_00 = get_nc_var(nc_file, var_name_p2_00); - NcVar var_t1_00 = get_nc_var(nc_file, var_name_t1_00); - NcVar var_t2_00 = get_nc_var(nc_file, var_name_t2_00); - NcVar var_p1_12 = get_nc_var(nc_file, var_name_p1_12); - NcVar var_p2_12 = get_nc_var(nc_file, var_name_p2_12); - NcVar var_t1_12 = get_nc_var(nc_file, var_name_t1_12); - NcVar var_t2_12 = get_nc_var(nc_file, var_name_t2_12); - NcVar var_matrix_00 = get_nc_var(nc_file, var_name_matrix_00); - NcVar var_matrix_12 = get_nc_var(nc_file, var_name_matrix_12); - - if (IS_INVALID_NC(var_sid) || !get_nc_data(&var_sid, sid_array)) { - mlog << Error << "\n" << method_name - << "Did not get sid\n\n"; - exit(1); - } - if (IS_INVALID_NC(var_lat) || !get_nc_data(&var_lat, lat_array)) { - mlog << Error << "\n" << method_name - << "Did not get lat\n\n"; - exit(1); - } - if (IS_INVALID_NC(var_lon) || !get_nc_data(&var_lon, lon_array)) { - mlog << Error << "\n" << method_name - << "Did not get lon\n\n"; - exit(1); - } - if (IS_INVALID_NC(var_elv) || !get_nc_data(&var_elv, elv_array)) { - mlog << Error << "\n" << method_name - << "Did not get elv\n\n"; - exit(1); - } - if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_00_array)) { + try { + double p1_00_buf[SEEPS_MONTH]; + double p2_00_buf[SEEPS_MONTH]; + double t1_00_buf[SEEPS_MONTH]; + double t2_00_buf[SEEPS_MONTH]; + double p1_12_buf[SEEPS_MONTH]; + double p2_12_buf[SEEPS_MONTH]; + double t1_12_buf[SEEPS_MONTH]; + double t2_12_buf[SEEPS_MONTH]; + double matrix_00_buf[SEEPS_MONTH*SEEPS_MATRIX_SIZE]; + double matrix_12_buf[SEEPS_MONTH*SEEPS_MATRIX_SIZE]; + NcFile *nc_file = open_ncfile(filename.c_str()); + + // dimensions: month = 12 ; nstn = 5293 ; nmatrix = 9 ; + get_dim(nc_file, dim_name_nstn, nstn, true); + mlog << Debug(6) << method_name << "dimensions nstn = " << nstn << "\n"; + if (standalone_debug_seeps) cout << "dimensions nstn = " << nstn << "\n"; + + int *sid_array = new int[nstn]; + double *lat_array = new double[nstn]; + double *lon_array = new double[nstn]; + double *elv_array = new double[nstn]; + double *p1_00_array = new double[nstn*SEEPS_MONTH]; + double *p2_00_array = new double[nstn*SEEPS_MONTH]; + double *t1_00_array = new double[nstn*SEEPS_MONTH]; + double *t2_00_array = new double[nstn*SEEPS_MONTH]; + double *p1_12_array = new double[nstn*SEEPS_MONTH]; + double *p2_12_array = new double[nstn*SEEPS_MONTH]; + double *t1_12_array = new double[nstn*SEEPS_MONTH]; + double *t2_12_array = new double[nstn*SEEPS_MONTH]; + double *matrix_00_array = new double[nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE]; + double *matrix_12_array = new double[nstn*SEEPS_MONTH*SEEPS_MATRIX_SIZE]; + + NcVar var_sid = get_nc_var(nc_file, var_name_sid); + NcVar var_lat = get_nc_var(nc_file, var_name_lat); + NcVar var_lon = get_nc_var(nc_file, var_name_lon); + NcVar var_elv = get_nc_var(nc_file, var_name_elv); + NcVar var_p1_00 = get_nc_var(nc_file, var_name_p1_00); + NcVar var_p2_00 = get_nc_var(nc_file, var_name_p2_00); + NcVar var_t1_00 = get_nc_var(nc_file, var_name_t1_00); + NcVar var_t2_00 = get_nc_var(nc_file, var_name_t2_00); + NcVar var_p1_12 = get_nc_var(nc_file, var_name_p1_12); + NcVar var_p2_12 = get_nc_var(nc_file, var_name_p2_12); + NcVar var_t1_12 = get_nc_var(nc_file, var_name_t1_12); + NcVar var_t2_12 = get_nc_var(nc_file, var_name_t2_12); + NcVar var_matrix_00 = get_nc_var(nc_file, var_name_matrix_00); + NcVar var_matrix_12 = get_nc_var(nc_file, var_name_matrix_12); + + if (IS_INVALID_NC(var_sid) || !get_nc_data(&var_sid, sid_array)) { + mlog << Error << "\n" << method_name + << "Did not get sid\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_lat) || !get_nc_data(&var_lat, lat_array)) { + mlog << Error << "\n" << method_name + << "Did not get lat\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_lon) || !get_nc_data(&var_lon, lon_array)) { + mlog << Error << "\n" << method_name + << "Did not get lon\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_elv) || !get_nc_data(&var_elv, elv_array)) { + mlog << Error << "\n" << method_name + << "Did not get elv\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_00_array)) { + mlog << Error << "\n" << method_name + << "Did not get p1_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_00_array)) { + mlog << Error << "\n" << method_name + << "Did not get p2_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_00_array)) { + mlog << Error << "\n" << method_name + << "Did not get t1_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_00_array)) { + mlog << Error << "\n" << method_name + << "Did not get t2_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_p1_12) || !get_nc_data(&var_p1_12, p1_12_array)) { + mlog << Error << "\n" << method_name + << "Did not get p1_12\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_p2_12) || !get_nc_data(&var_p2_12, p2_12_array)) { + mlog << Error << "\n" << method_name + << "Did not get p2_12\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_t1_12) || !get_nc_data(&var_t1_12, t1_12_array)) { + mlog << Error << "\n" << method_name + << "Did not get t1_12\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_t2_12) || !get_nc_data(&var_t2_12, t2_12_array)) { + mlog << Error << "\n" << method_name + << "Did not get t2_12\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_matrix_00) || !get_nc_data(&var_matrix_00, matrix_00_array)) { + mlog << Error << "\n" << method_name + << "Did not get matrix_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_matrix_12) || !get_nc_data(&var_matrix_12, matrix_12_array)) { + mlog << Error << "\n" << method_name + << "Did not get matrix_12\n\n"; + exit(1); + } + + SeepsClimoRecord *rec_00; + SeepsClimoRecord *rec_12; + for (int idx=0; idxclose(); + + float duration = (float)(clock() - clock_time)/CLOCKS_PER_SEC; + mlog << Debug(6) << method_name << "took " << duration << " seconds\n"; + if (standalone_debug_seeps) cout << method_name << "took " << duration << " seconds\n"; + } + catch(int i_err) { + + seeps_ready = false; mlog << Error << "\n" << method_name - << "Did not get p1_00\n\n"; + << "encountered an error on reading " << filename << ". ecception_code=" << i_err << "\n\n"; + + exit(i_err); + } // end catch block + +} + + + +//////////////////////////////////////////////////////////////////////// + + +SeepsClimoGrid::SeepsClimoGrid(int month, int hour) : month{month}, hour{hour} +{ + + p1_buf = p2_buf = t1_buf = t2_buf = NULL; + s12_buf = s13_buf = s21_buf = s23_buf = s31_buf = s32_buf = NULL; + + ConcatString seeps_name = get_seeps_climo_filename(); + seeps_ready = file_exists(seeps_name.c_str()); + if (seeps_ready) read_seeps_scores(seeps_name); + else { + mlog << Error << "\nSeepsClimoGrid::SeepsClimoGrid -> " + << "The SEEPS climo data \"" << seeps_name << "\" is missing!" + << " Turn off SEEPS to continue\n\n"; exit(1); } - if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_00_array)) { - mlog << Error << "\n" << method_name - << "Did not get p2_00\n\n"; - exit(1); + +} + +//////////////////////////////////////////////////////////////////////// + +SeepsClimoGrid::~SeepsClimoGrid() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +void SeepsClimoGrid::clear() { + SeepsClimoBase::clear(); + if (NULL != p1_buf) { delete [] p1_buf; p1_buf = NULL; } + if (NULL != p2_buf) { delete [] p2_buf; p2_buf = NULL; } + if (NULL != t1_buf) { delete [] t1_buf; t1_buf = NULL; } + if (NULL != t2_buf) { delete [] t2_buf; t2_buf = NULL; } + if (NULL != s12_buf) { delete [] s12_buf; s12_buf = NULL; } + if (NULL != s13_buf) { delete [] s13_buf; s13_buf = NULL; } + if (NULL != s21_buf) { delete [] s21_buf; s21_buf = NULL; } + if (NULL != s23_buf) { delete [] s23_buf; s23_buf = NULL; } + if (NULL != s31_buf) { delete [] s31_buf; s31_buf = NULL; } + if (NULL != s32_buf) { delete [] s32_buf; s32_buf = NULL; } +}; + +//////////////////////////////////////////////////////////////////////// + +SeepsScore *SeepsClimoGrid::get_record(int ix, int iy, + double p_fcst, double p_obs) { + SeepsScore *seeps_record = NULL; + const char *method_name = "SeepsClimoGrid::get_record() -> "; + if (!is_eq(p_fcst, -9999.0) && !is_eq(p_obs, -9999.0)) { + int offset = iy * nx + ix; + double p1 = p1_buf[offset]; + + if (seeps_p1_thresh.check(p1)) { + // Determine location in contingency table + int ic = (p_obs>t1_buf[offset])+(p_obs>t2_buf[offset]); + int jc = (p_fcst>t1_buf[offset])+(p_fcst>t2_buf[offset]); + double score = get_score(offset, ic, jc); + + seeps_record = new SeepsScore(); + seeps_record->obs_cat = ic; + seeps_record->fcst_cat = jc; + seeps_record->s_idx = (jc*3)+ic; + seeps_record->p1 = p1_buf[offset]; + seeps_record->p2 = p2_buf[offset]; + seeps_record->t1 = t1_buf[offset]; + seeps_record->t2 = t2_buf[offset]; + seeps_record->score = score; + } + else if (~is_eq(p1, bad_data_double)) { + filtered_count++; + mlog << Debug(7) << method_name << " filtered by threshold p1=" << p1_buf[offset] <<"\n"; + } } - if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_00_array)) { - mlog << Error << "\n" << method_name - << "Did not get t1_00\n\n"; - exit(1); + + return seeps_record; +} + +//////////////////////////////////////////////////////////////////////// + +double SeepsClimoGrid::get_score(int offset, int obs_cat, int fcst_cat) { + double score = bad_data_double; + + if (offset >= (nx * ny)) { + mlog << Error << "\nSeepsClimoGrid::get_score() --> offset (" << offset + << " is too big (" << (nx*ny) << ")\n"; + return score; } - if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_00_array)) { - mlog << Error << "\n" << method_name - << "Did not get t2_00\n\n"; - exit(1); + + if (obs_cat == 0) { + if (fcst_cat == 1) score = s12_buf[offset]; + else if (fcst_cat == 2) score = s13_buf[offset]; + else score = 0.; } - if (IS_INVALID_NC(var_p1_12) || !get_nc_data(&var_p1_12, p1_12_array)) { - mlog << Error << "\n" << method_name - << "Did not get p1_12\n\n"; - exit(1); + else if (obs_cat == 1) { + if (fcst_cat == 0) score = s21_buf[offset]; + else if (fcst_cat == 2) score = s23_buf[offset]; + else score = 0.; } - if (IS_INVALID_NC(var_p2_12) || !get_nc_data(&var_p2_12, p2_12_array)) { - mlog << Error << "\n" << method_name - << "Did not get p2_12\n\n"; - exit(1); + else { + if (fcst_cat == 0) score = s31_buf[offset]; + else if (fcst_cat == 1) score = s32_buf[offset]; + else score = 0.; } - if (IS_INVALID_NC(var_t1_12) || !get_nc_data(&var_t1_12, t1_12_array)) { - mlog << Error << "\n" << method_name - << "Did not get t1_12\n\n"; - exit(1); + + return score; +} + +//////////////////////////////////////////////////////////////////////// + +double SeepsClimoGrid::get_score(int ix, int iy, double p_fcst, double p_obs) { + double score = bad_data_double; + + if (!is_eq(p_fcst, -9999.0) && !is_eq(p_obs, -9999.0)) { + int offset = iy * nx + ix; + // Determine location in contingency table + int ic = (p_obs>t1_buf[offset])+(p_obs>t2_buf[offset]); + int jc = (p_fcst>t1_buf[offset])+(p_fcst>t2_buf[offset]); + score = get_score(offset, ic, jc); } - if (IS_INVALID_NC(var_t2_12) || !get_nc_data(&var_t2_12, t2_12_array)) { - mlog << Error << "\n" << method_name - << "Did not get t2_12\n\n"; - exit(1); + + return score; +} + +//////////////////////////////////////////////////////////////////////// + +ConcatString SeepsClimoGrid::get_seeps_climo_filename() { + ConcatString seeps_filename; + const char *method_name = "SeepsClimoGrid::get_seeps_climo_filename() -> "; + + // Use the MET_TMP_DIR environment variable, if set. + bool use_env = get_env(MET_ENV_SEEPS_CLIMO_GRID_NAME, seeps_filename); + if(use_env) { + seeps_filename = replace_path(seeps_filename); } - if (IS_INVALID_NC(var_matrix_00) || !get_nc_data(&var_matrix_00, matrix_00_array)) { - mlog << Error << "\n" << method_name - << "Did not get matrix_00\n\n"; - exit(1); + else seeps_filename = replace_path(def_seeps_grid_filename); + + if (seeps_ready = file_exists(seeps_filename.c_str())) { + mlog << Debug(7) << method_name << "SEEPS climo name=\"" + << seeps_filename.c_str() << "\"\n"; } - if (IS_INVALID_NC(var_matrix_12) || !get_nc_data(&var_matrix_12, matrix_12_array)) { - mlog << Error << "\n" << method_name - << "Did not get matrix_12\n\n"; - exit(1); + else { + ConcatString message = " "; + if (use_env) { + message.add("from the env. name "); + message.add(MET_ENV_SEEPS_CLIMO_GRID_NAME); + } + mlog << Warning << "\n" << method_name + << "The SEEPS climo name \"" << seeps_filename.c_str() + << "\"" << message << " does not exist!\n\n"; } - SeepsClimoRecord *rec_00; - SeepsClimoRecord *rec_12; - for (int idx=0; idxclose(); - - float duration = (float)(clock() - clock_time)/CLOCKS_PER_SEC; - mlog << Debug(4) << method_name << "took " << duration << " seconds\n"; - if (standalone_debug_seeps) cout << method_name << "took " << duration << " seconds\n"; - + return seeps_filename; } //////////////////////////////////////////////////////////////////////// -void SeepsAggScore::init() { +void SeepsClimoGrid::read_seeps_scores(ConcatString filename) { + clock_t clock_time = clock(); + const char *method_name = "SeepsClimoGrid::read_seeps_scores() -> "; - n_obs = 0; - c12 = 0; - c13 = 0; - c21 = 0; - c23 = 0; - c31 = 0; - c32 = 0; - s12 = 0.; - s13 = 0.; - s21 = 0.; - s23 = 0.; - s31 = 0.; - s32 = 0.; - pv1 = 0.; - pv2 = 0.; - pv3 = 0.; - pf1 = 0.; - pf2 = 0.; - pf3 = 0.; - mean_fcst = bad_data_double; - mean_obs = bad_data_double; - score = bad_data_double; + try { + NcFile *nc_file = open_ncfile(filename.c_str()); + + // dimensions: month = 12; + if (!has_dim(nc_file, dim_name_lat) || !has_dim(nc_file, dim_name_lon)) { + mlog << Error << "\n" << method_name + << "\"" << filename << "\" is not valid SEEPS climo file\n\n"; + //exit(1); + } + + get_dim(nc_file, dim_name_lat, ny, true); + get_dim(nc_file, dim_name_lon, nx, true); + mlog << Debug(6) << method_name << "dimensions lon = " << nx << " lat = " << ny + << " month=" << month << "\n"; + if (standalone_debug_seeps) cout << "dimensions lon = " << nx << " lat = " << ny + << " month=" << month << "\n";; + + p1_buf = new double[nx*ny]; + p2_buf = new double[nx*ny]; + t1_buf = new double[nx*ny]; + t2_buf = new double[nx*ny]; + s12_buf = new double[nx*ny]; + s13_buf = new double[nx*ny]; + s21_buf = new double[nx*ny]; + s23_buf = new double[nx*ny]; + s31_buf = new double[nx*ny]; + s32_buf = new double[nx*ny]; + + long curs[3] = { month-1, 0, 0 }; + long dims[3] = { 1, ny, nx }; + NcVar var_p1_00 = get_nc_var(nc_file, var_name_p1_00); + NcVar var_p2_00 = get_nc_var(nc_file, var_name_p2_00); + NcVar var_t1_00 = get_nc_var(nc_file, var_name_t1_00); + NcVar var_t2_00 = get_nc_var(nc_file, var_name_t2_00); + NcVar var_s12_00 = get_nc_var(nc_file, var_name_s12_00); + NcVar var_s13_00 = get_nc_var(nc_file, var_name_s13_00); + NcVar var_s21_00 = get_nc_var(nc_file, var_name_s21_00); + NcVar var_s23_00 = get_nc_var(nc_file, var_name_s23_00); + NcVar var_s31_00 = get_nc_var(nc_file, var_name_s31_00); + NcVar var_s32_00 = get_nc_var(nc_file, var_name_s32_00); + + if (IS_INVALID_NC(var_p1_00) || !get_nc_data(&var_p1_00, p1_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get p1_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_p2_00) || !get_nc_data(&var_p2_00, p2_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get p2_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_t1_00) || !get_nc_data(&var_t1_00, t1_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get t1_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_t2_00) || !get_nc_data(&var_t2_00, t2_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get t2_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_s12_00) || !get_nc_data(&var_s12_00, s12_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get s12_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_s13_00) || !get_nc_data(&var_s12_00, s13_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get s13_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_s21_00) || !get_nc_data(&var_s21_00, s21_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get s21_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_s23_00) || !get_nc_data(&var_s23_00, s23_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get s23_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_s31_00) || !get_nc_data(&var_s31_00, s31_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get s31_00\n\n"; + exit(1); + } + if (IS_INVALID_NC(var_s32_00) || !get_nc_data(&var_s32_00, s32_buf, dims, curs)) { + mlog << Error << "\n" << method_name + << "Did not get s32_00\n\n"; + exit(1); + } + nc_file->close(); + + float duration = (float)(clock() - clock_time)/CLOCKS_PER_SEC; + mlog << Debug(6) << method_name << "took " << duration << " seconds\n"; + if (standalone_debug_seeps) cout << method_name << "took " << duration << " seconds\n"; + } + catch(...) { + + seeps_ready = false; + mlog << Error << "\n" << method_name + << "encountered an error on reading " << filename << ".\n\n"; + + exit(-1); + } // end catch block } //////////////////////////////////////////////////////////////////////// + +void SeepsClimoGrid::print_all() { + const char *method_name = "SeepsClimoGrid::print_all() -> "; + if (standalone_debug_seeps) { + int offset; + offset = 0; + cout << method_name << " p1_buf[" << offset << "] = " << p1_buf[offset] << "\n"; + cout << method_name << " p2_buf[" << offset << "] = " << p2_buf[offset] << "\n"; + cout << method_name << " t1_buf[" << offset << "] = " << t1_buf[offset] << "\n"; + cout << method_name << " t2_buf[" << offset << "] = " << t2_buf[offset] << "\n"; + cout << method_name << "s12_buf[" << offset << "] = " << s12_buf[offset] << "\n"; + cout << method_name << "s13_buf[" << offset << "] = " << s13_buf[offset] << "\n"; + cout << method_name << "s21_buf[" << offset << "] = " << s21_buf[offset] << "\n"; + cout << method_name << "s23_buf[" << offset << "] = " << s23_buf[offset] << "\n"; + cout << method_name << "s31_buf[" << offset << "] = " << s31_buf[offset] << "\n"; + cout << method_name << "s32_buf[" << offset << "] = " << s32_buf[offset] << "\n"; + + offset = 400; + cout << method_name << " p1_buf[" << offset << "] = " << p1_buf[offset] << "\n"; + cout << method_name << " p2_buf[" << offset << "] = " << p2_buf[offset] << "\n"; + cout << method_name << " t1_buf[" << offset << "] = " << t1_buf[offset] << "\n"; + cout << method_name << " t2_buf[" << offset << "] = " << t2_buf[offset] << "\n"; + cout << method_name << "s12_buf[" << offset << "] = " << s12_buf[offset] << "\n"; + cout << method_name << "s13_buf[" << offset << "] = " << s13_buf[offset] << "\n"; + cout << method_name << "s21_buf[" << offset << "] = " << s21_buf[offset] << "\n"; + cout << method_name << "s23_buf[" << offset << "] = " << s23_buf[offset] << "\n"; + cout << method_name << "s31_buf[" << offset << "] = " << s31_buf[offset] << "\n"; + cout << method_name << "s32_buf[" << offset << "] = " << s32_buf[offset] << "\n"; + + offset = (nx*ny) - 1; + cout << method_name << " p1_buf[" << offset << "] = " << p1_buf[offset] << "\n"; + cout << method_name << " p2_buf[" << offset << "] = " << p2_buf[offset] << "\n"; + cout << method_name << " t1_buf[" << offset << "] = " << t1_buf[offset] << "\n"; + cout << method_name << " t2_buf[" << offset << "] = " << t2_buf[offset] << "\n"; + cout << method_name << "s12_buf[" << offset << "] = " << s12_buf[offset] << "\n"; + cout << method_name << "s13_buf[" << offset << "] = " << s13_buf[offset] << "\n"; + cout << method_name << "s21_buf[" << offset << "] = " << s21_buf[offset] << "\n"; + cout << method_name << "s23_buf[" << offset << "] = " << s23_buf[offset] << "\n"; + cout << method_name << "s31_buf[" << offset << "] = " << s31_buf[offset] << "\n"; + cout << method_name << "s32_buf[" << offset << "] = " << s32_buf[offset] << "\n"; + } +} + +//////////////////////////////////////////////////////////////////////// + diff --git a/src/libcode/vx_seeps/seeps.h b/src/libcode/vx_seeps/seeps.h index 6fd08ca12e..ae5915298e 100644 --- a/src/libcode/vx_seeps/seeps.h +++ b/src/libcode/vx_seeps/seeps.h @@ -24,7 +24,8 @@ //////////////////////////////////////////////////////////////////////// -static const char *MET_ENV_SEEPS_CLIMO_NAME = "MET_SEEPS_CLIMO_NAME"; +static const char *MET_ENV_SEEPS_CLIMO_NAME = "MET_SEEPS_POINT_CLIMO_NAME"; +static const char *MET_ENV_SEEPS_CLIMO_GRID_NAME = "MET_SEEPS_GRID_CLIMO_NAME"; static const char *dim_name_nstn = "nstn"; @@ -38,24 +39,43 @@ static const char *var_name_t1_12 = "t1_12"; static const char *var_name_t2_12 = "t2_12"; static const char *var_name_matrix_00 = "matrix_00"; static const char *var_name_matrix_12 = "matrix_12"; +static const char *var_name_s12_00 = "s12_00"; +static const char *var_name_s13_00 = "s13_00"; +static const char *var_name_s21_00 = "s21_00"; +static const char *var_name_s23_00 = "s23_00"; +static const char *var_name_s31_00 = "s31_00"; +static const char *var_name_s32_00 = "s32_00"; +static const char *var_name_s12_12 = "s12_12"; +static const char *var_name_s13_12 = "s13_12"; +static const char *var_name_s21_12 = "s21_12"; +static const char *var_name_s23_12 = "s23_12"; +static const char *var_name_s31_12 = "s31_12"; +static const char *var_name_s32_12 = "s32_12"; + +//density_radius = 0.75 degrees (83km; this is described as “the smallest possible +// value that ensures approximately equal representation of all subregions of Europe”.) +static double density_radius = 0.75; +const double density_radius_rad = density_radius * rad_per_deg; //////////////////////////////////////////////////////////////////////// -struct SeepsScore { +struct SeepsScore { // For SEEPS_MPR int obs_cat; // i = obs category 0,1,2 - int model_cat; // j = model category 0,1,2 - float p1; - float p2; - float t1; - float t2; - float score; + int fcst_cat; // j = model category 0,1,2 + int s_idx; // index for 3 by 3 matrix as 1 dimensional (fcst_cat*3)+obs_cat + double p1; + double p2; + double t1; + double t2; + double score; }; //////////////////////////////////////////////////////////////////////// -struct SeepsAggScore { - void init(); - +struct SeepsAggScore { // For SEEPS + void clear(); + SeepsAggScore & operator+=(const SeepsAggScore &); + int n_obs; int c12; int c13; @@ -63,21 +83,22 @@ struct SeepsAggScore { int c23; int c31; int c32; - float s12; - float s13; - float s21; - float s23; - float s31; - float s32; - float pv1; // marginal probabilities of the observed values - float pv2; - float pv3; - float pf1; // marginal probabilities of the forecast values - float pf2; - float pf3; - float mean_fcst; - float mean_obs; - float score; + double s12; + double s13; + double s21; + double s23; + double s31; + double s32; + double pv1; // marginal probabilities of the observed values + double pv2; + double pv3; + double pf1; // marginal probabilities of the forecast values + double pf2; + double pf3; + double mean_fcst; + double mean_obs; + double score; + double weighted_score; }; //////////////////////////////////////////////////////////////////////// @@ -86,57 +107,80 @@ struct SeepsAggScore { struct SeepsRecord { int sid; int month; // 1 to 12, not 0 to 11 - float lat; - float lon; - float elv; - float p1; - float p2; - float t1; - float t2; - float scores[SEEPS_MATRIX_SIZE]; + double lat; + double lon; + double elv; + double p1; + double p2; + double t1; + double t2; + double scores[SEEPS_MATRIX_SIZE]; }; //////////////////////////////////////////////////////////////////////// struct SeepsClimoRecord { int sid; - float lat; - float lon; - float elv; - float p1[SEEPS_MONTH]; - float p2[SEEPS_MONTH]; - float t1[SEEPS_MONTH]; - float t2[SEEPS_MONTH]; - float scores[SEEPS_MONTH][SEEPS_MATRIX_SIZE]; + double lat; + double lon; + double elv; + double p1[SEEPS_MONTH]; + double p2[SEEPS_MONTH]; + double t1[SEEPS_MONTH]; + double t2[SEEPS_MONTH]; + double scores[SEEPS_MONTH][SEEPS_MATRIX_SIZE]; }; //////////////////////////////////////////////////////////////////////// -class SeepsClimo { +class SeepsClimoBase { - private: + protected: bool seeps_ready; + int filtered_count; + SingleThresh seeps_p1_thresh; // Range of SEEPS p1 (probability of being dry)std::map seeps_score_00_map; + + virtual void clear(); + + public: + + SeepsClimoBase(); + ~SeepsClimoBase(); + + void set_p1_thresh(const SingleThresh &p1_thresh); + int get_filtered_count(); + +}; + +//////////////////////////////////////////////////////////////////////// + +class SeepsClimo : public SeepsClimoBase { + + private: + int nstn; std::map seeps_score_00_map; std::map seeps_score_12_map; - SeepsClimoRecord *create_climo_record(int sid, float lat, float lon, float elv, - float *p1, float *p2, float *t1, float *t2, - float *scores); - ConcatString get_seeps_climo_filename(); + SeepsClimoRecord *create_climo_record(int sid, double lat, double lon, double elv, + double *p1, double *p2, double *t1, double *t2, + double *scores); void print_record(SeepsClimoRecord *record, bool with_header=false); void read_records(ConcatString filename); + ConcatString get_seeps_climo_filename(); + void read_seeps_scores(ConcatString filename); + public: SeepsClimo(); ~SeepsClimo(); - void clear(); - SeepsRecord *get_record(int sid, int month, int hour=0); - float get_score(int sid, float p_fcst, float p_obs, int month, int hour=0); - SeepsScore *get_seeps_score(int sid, float p_fcst, float p_obs, int month, int hour=0); + void clear() override; + SeepsRecord *get_record(int sid, int month, int hour); + double get_score(int sid, double p_fcst, double p_obs, int month, int hour); + SeepsScore *get_seeps_score(int sid, double p_fcst, double p_obs, int month, int hour); void print_all(); void print_record(SeepsRecord *record, bool with_header=false); @@ -150,8 +194,57 @@ class SeepsClimo { //////////////////////////////////////////////////////////////////////// +class SeepsClimoGrid : public SeepsClimoBase { + + private: + + int month; + int hour; + int nx; + int ny; + double *p1_buf; + double *p2_buf; + double *t1_buf; + double *t2_buf; + double *s12_buf; + double *s13_buf; + double *s21_buf; + double *s23_buf; + double *s31_buf; + double *s32_buf; + + ConcatString get_seeps_climo_filename(); + void read_seeps_scores(ConcatString filename); + + public: + + SeepsClimoGrid(int month, int hour); + ~SeepsClimoGrid(); + + void clear() override; + SeepsScore *get_record(int ix, int iy, double p_fcst, double p_obs); + double get_score(int offset, int obs_cat, int fcst_cat); + double get_score(int ix, int iy, double p_fcst, double p_obs); + void print_all(); + + // + // + // + +}; + + +//////////////////////////////////////////////////////////////////////// + +inline int SeepsClimoBase::get_filtered_count() { return filtered_count; } + +//////////////////////////////////////////////////////////////////////// + extern SeepsClimo *get_seeps_climo(); +extern SeepsClimoGrid *get_seeps_climo_grid(int month, int hour=0); + extern void release_seeps_climo(); +extern void release_seeps_climo_grid(); //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_stat_out/stat_columns.cc b/src/libcode/vx_stat_out/stat_columns.cc index 7844205758..a106f389a9 100644 --- a/src/libcode/vx_stat_out/stat_columns.cc +++ b/src/libcode/vx_stat_out/stat_columns.cc @@ -24,6 +24,10 @@ using namespace std; //////////////////////////////////////////////////////////////////////// +const bool use_weighted_seeps = false; + +//////////////////////////////////////////////////////////////////////// + void parse_row_col(const char *col_name, int &r, int &c) { const char *ptr = (const char *) 0; @@ -1567,7 +1571,7 @@ void write_mpr_row(StatHdrColumns &shc, const PairDataPoint *pd_ptr, //////////////////////////////////////////////////////////////////////// -void write_seeps_row(StatHdrColumns &shc, const PairDataPoint *pd_ptr, +void write_seeps_row(StatHdrColumns &shc, const SeepsAggScore *seeps, STATOutputType out_type, AsciiTable &stat_at, int &stat_row, AsciiTable &txt_at, int &txt_row, @@ -1589,14 +1593,13 @@ void write_seeps_row(StatHdrColumns &shc, const PairDataPoint *pd_ptr, shc.set_alpha(bad_data_double); // Write a line - if(pd_ptr->seeps.n_obs > 0 && - !is_bad_data(pd_ptr->seeps.score)) { + if(seeps->n_obs > 0 && !is_bad_data(seeps->score)) { // Write the header columns write_header_cols(shc, stat_at, stat_row); // Write the data columns - write_seeps_cols(pd_ptr, stat_at, stat_row, n_header_columns); + write_seeps_cols(seeps, stat_at, stat_row, n_header_columns); // If requested, copy row to the text file if(out_type == STATOutputType_Both) { @@ -4068,8 +4071,8 @@ void write_mpr_cols(const PairDataPoint *pd_ptr, int i, //////////////////////////////////////////////////////////////////////// -void write_seeps_cols(const PairDataPoint *pd_ptr, - AsciiTable &at, int r, int c) { +void write_seeps_cols(const SeepsAggScore *seeps, + AsciiTable &at, int r, int c) { // // Stable Equitable Error in Probability Space (SEEPS) // Dump out the SEEPS line: @@ -4081,27 +4084,27 @@ void write_seeps_cols(const PairDataPoint *pd_ptr, // SEEPS // - at.set_entry(r, c+0, pd_ptr->seeps.n_obs); // Total Number of Pairs + at.set_entry(r, c+0, seeps->n_obs); // Total Number of Pairs - at.set_entry(r, c+1, bad_data_double); // pd_ptr->seeps.s12); // s12 - at.set_entry(r, c+2, bad_data_double); // pd_ptr->seeps.s13); // S13 - at.set_entry(r, c+3, bad_data_double); // pd_ptr->seeps.s21); // s21 - at.set_entry(r, c+4, bad_data_double); // pd_ptr->seeps.s23); // s23 - at.set_entry(r, c+5, bad_data_double); // pd_ptr->seeps.s31); // s31 - at.set_entry(r, c+6, bad_data_double); // pd_ptr->seeps.s32); // s32 + at.set_entry(r, c+1, seeps->s12); // s12 + at.set_entry(r, c+2, seeps->s13); // s13 + at.set_entry(r, c+3, seeps->s21); // s21 + at.set_entry(r, c+4, seeps->s23); // s23 + at.set_entry(r, c+5, seeps->s31); // s31 + at.set_entry(r, c+6, seeps->s32); // s32 - at.set_entry(r, c+7, bad_data_double); // pd_ptr->seeps.pf1); // pf1 - at.set_entry(r, c+8, bad_data_double); // pd_ptr->seeps.pf2); // pf2 - at.set_entry(r, c+9, bad_data_double); // pd_ptr->seeps.pf3); // pf3 + at.set_entry(r, c+7, seeps->pf1); // pf1 + at.set_entry(r, c+8, seeps->pf2); // pf2 + at.set_entry(r, c+9, seeps->pf3); // pf3 - at.set_entry(r, c+10, bad_data_double); // pd_ptr->seeps.pv1); // pv1 - at.set_entry(r, c+11, bad_data_double); // pd_ptr->seeps.pv2); // pv2 - at.set_entry(r, c+12, bad_data_double); // pd_ptr->seeps.pv3); // pv3 + at.set_entry(r, c+10, seeps->pv1); // pv1 + at.set_entry(r, c+11, seeps->pv2); // pv2 + at.set_entry(r, c+12, seeps->pv3); // pv3 - at.set_entry(r, c+13, pd_ptr->seeps.mean_fcst); // mean_fcst - at.set_entry(r, c+14, pd_ptr->seeps.mean_obs); // mean_obs + at.set_entry(r, c+13, seeps->mean_fcst); // mean_fcst + at.set_entry(r, c+14, seeps->mean_obs); // mean_obs - at.set_entry(r, c+15, pd_ptr->seeps.score); // SEEPS score + at.set_entry(r, c+15, (use_weighted_seeps ? seeps->weighted_score : seeps->score)); // SEEPS score/weighted score return; } @@ -4132,7 +4135,7 @@ void write_seeps_mpr_cols(const PairDataPoint *pd_ptr, int i, at.set_entry(r, c+5, (string)pd_ptr->o_qc_sa[i]); // Observation Quality Control - at.set_entry(r, c+6, pd_ptr->seeps_mpr[i]->model_cat); // model category + at.set_entry(r, c+6, pd_ptr->seeps_mpr[i]->fcst_cat); // model category at.set_entry(r, c+7, pd_ptr->seeps_mpr[i]->obs_cat); // observation category at.set_entry(r, c+8, pd_ptr->seeps_mpr[i]->p1); // p1 diff --git a/src/libcode/vx_stat_out/stat_columns.h b/src/libcode/vx_stat_out/stat_columns.h index 27b16aa66b..8bcfddb789 100644 --- a/src/libcode/vx_stat_out/stat_columns.h +++ b/src/libcode/vx_stat_out/stat_columns.h @@ -104,7 +104,7 @@ extern void write_dmap_row (StatHdrColumns &, const DMAPInfo &, STATOutputType, extern void write_mpr_row (StatHdrColumns &, const PairDataPoint *, STATOutputType, AsciiTable &, int &, AsciiTable &, int &, bool update_thresh = true); -extern void write_seeps_row (StatHdrColumns &, const PairDataPoint *, STATOutputType, +extern void write_seeps_row (StatHdrColumns &, const SeepsAggScore *, STATOutputType, AsciiTable &, int &, AsciiTable &, int &, bool update_thresh = true); extern void write_seeps_mpr_row (StatHdrColumns &, const PairDataPoint *, STATOutputType, @@ -179,7 +179,7 @@ extern void write_dmap_cols (const DMAPInfo &, AsciiTable &, int, int); extern void write_mpr_cols (const PairDataPoint *, int, AsciiTable &, int, int); -extern void write_seeps_cols (const PairDataPoint *, +extern void write_seeps_cols (const SeepsAggScore *, AsciiTable &, int, int); extern void write_seeps_mpr_cols (const PairDataPoint *, int, AsciiTable &, int, int); diff --git a/src/libcode/vx_statistics/compute_stats.cc b/src/libcode/vx_statistics/compute_stats.cc index 8ad12fa36e..50bfdb3638 100644 --- a/src/libcode/vx_statistics/compute_stats.cc +++ b/src/libcode/vx_statistics/compute_stats.cc @@ -23,6 +23,11 @@ using namespace std; #include "vx_util.h" #include "vx_log.h" +//////////////////////////////////////////////////////////////////////// + +const int detailed_debug_level = 5; + + //////////////////////////////////////////////////////////////////////// void compute_cntinfo(const SL1L2Info &s, bool aflag, CNTInfo &cnt_info) { @@ -1412,12 +1417,12 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps) { SeepsScore *seeps_mpr; int count, count_diagonal; int c12, c13, c21, c23, c31, c32; - double score, obs_sum, fcst_sum; + double score_sum, obs_sum, fcst_sum; + vector seeps_mprs; - score = obs_sum = fcst_sum = 0.0; + score_sum = obs_sum = fcst_sum = 0.0; count = count_diagonal = c12 = c13 = c21 = c23 = c31 = c32 = 0; for(int i=0; in_obs; i++) { - if (i >= pd->seeps_mpr.size()) break; seeps_mpr = pd->seeps_mpr[i]; if (!seeps_mpr || is_eq(seeps_mpr->score, bad_data_double)) continue; @@ -1425,32 +1430,93 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps) { count++; obs_sum += pd->o_na[i]; // Observation Value fcst_sum += pd->f_na[i]; // Forecast Value - score += seeps_mpr->score; - - if (seeps_mpr->model_cat == 0) { + score_sum += seeps_mpr->score; + if (seeps_mpr->fcst_cat == 0) { if (seeps_mpr->obs_cat == 1) c12++; else if(seeps_mpr->obs_cat == 2) c13++; else count_diagonal++; - // pd.pf1 } - else if (seeps_mpr->model_cat == 1) { + else if (seeps_mpr->fcst_cat == 1) { if (seeps_mpr->obs_cat == 0) c21++; else if(seeps_mpr->obs_cat == 2) c23++; else count_diagonal++; - // pd.pf2 } - else if (seeps_mpr->model_cat == 2) { + else if (seeps_mpr->fcst_cat == 2) { if (seeps_mpr->obs_cat == 0) c31++; else if (seeps_mpr->obs_cat == 1) c32++; else count_diagonal++; - // pd.pf3 } + seeps_mprs.push_back(seeps_mpr); } if (count > 0) { + double *density_vector; + double pvf[SEEPS_MATRIX_SIZE]; + double weighted_score, weight_sum, weight[count]; + seeps->n_obs = count; seeps->mean_fcst = fcst_sum / count; seeps->mean_obs = obs_sum / count; - seeps->score = score / count; + seeps->score = score_sum / count; + density_vector = compute_seeps_density_vector(pd, seeps); + + weighted_score = 0.; + for (int i=0; iscore * weight[i]; + //IDL: svf(cat{i)) = svf(cat{i)) + c(4+cat(i) * w{i) + //IDL: pvf(cat{i)) = pvf(cat{i)) + w{i) + pvf[seeps_mpr->s_idx] += weight[i]; + } + else mlog << Debug(1) << method_name + << "the length of density vector (" << count << ") is less than MPR.\n"; + } + } + + if (density_vector != NULL) delete [] density_vector; + } + seeps_mprs.clear(); + + // The weight for s12 to s32 should come from climo file, but not available yet + seeps->pv1 = pvf[0] + pvf[3] + pvf[6]; // sum by column for obs + seeps->pv2 = pvf[1] + pvf[4] + pvf[7]; // sum by column for obs + seeps->pv3 = pvf[2] + pvf[5] + pvf[8]; // sum by column for obs + seeps->pf1 = pvf[0] + pvf[1] + pvf[2]; // sum by row for forecast + seeps->pf2 = pvf[3] + pvf[4] + pvf[5]; // sum by row for forecast + seeps->pf3 = pvf[6] + pvf[7] + pvf[8]; // sum by row for forecast + seeps->s12 = c12 * seeps->pf1 * seeps->pv2; + seeps->s13 = c13 * seeps->pf1 * seeps->pv3; + seeps->s21 = c21 * seeps->pf2 * seeps->pv1; + seeps->s23 = c23 * seeps->pf2 * seeps->pv3; + seeps->s31 = c31 * seeps->pf3 * seeps->pv1; + seeps->s32 = c32 * seeps->pf3 * seeps->pv2; + seeps->weighted_score = weighted_score; + + mlog << Debug(7) << method_name + << "SEEPS score=" << seeps->score << " weighted_score=" << weighted_score + << " pv1=" << seeps->pv1 << " pv2=" << seeps->pv2 << " pv3=" << seeps->pv3 + << " pf1=" << seeps->pf1 << " pf2=" << seeps->pf2 << " pf3=" << seeps->pf3 << "\n"; + } + else { + mlog << Debug(5) << method_name + << "no SEEPS_MPR available\n"; } seeps->c12 = c12; seeps->c13 = c13; @@ -1460,7 +1526,7 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps) { seeps->c32 = c32; if (count != (c12+c13+c21+c23+c31+c32+count_diagonal)){ - mlog << Debug(2) << method_name + mlog << Debug(6) << method_name << "INFO check count: all=" << count << " s12=" << c12<< " s13=" << c13 << " s21=" << c21 << " s23=" << c23 << " s31=" << c31 << " s32=" << c32 << "\n"; @@ -1470,3 +1536,264 @@ void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps) { } //////////////////////////////////////////////////////////////////////// + +void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &obs_dp, + DataPlane &seeps_dp, DataPlane &seeps_dp_fcat, + DataPlane &seeps_dp_ocat,SeepsAggScore *seeps, + int month, int hour, const SingleThresh &seeps_p1_thresh) { + int fcst_cat, obs_cat; + int seeps_count, count_diagonal, nan_count, bad_count; + int nx = fcst_dp.nx(); + int ny = fcst_dp.ny(); + int dp_size = (nx * ny); + int pvf_cnt[SEEPS_MATRIX_SIZE]; + double pvf[SEEPS_MATRIX_SIZE]; + int c12, c13, c21, c23, c31, c32; + double obs_sum, fcst_sum; + double seeps_score, seeps_score_sum, seeps_score_partial_sum; + SeepsScore *seeps_mpr; + static const char *method_name = "compute_aggregated_seeps_grid() -> "; + + seeps_dp.set_size(nx, ny); + seeps_dp_fcat.set_size(nx, ny); + seeps_dp_ocat.set_size(nx, ny); + obs_sum = fcst_sum = seeps_score_sum = 0.; + seeps_count = count_diagonal = nan_count = bad_count = 0; + c12 = c13 = c21 = c23 = c31 = c32 = 0; + + seeps->clear(); + SeepsClimoGrid *seeps_climo = get_seeps_climo_grid(month); + seeps_climo->set_p1_thresh(seeps_p1_thresh); + for (int i=0; iget_record(ix, iy, fcst_value, obs_value); + if (seeps_mpr != NULL) { + fcst_cat = seeps_mpr->fcst_cat; + obs_cat = seeps_mpr->obs_cat; + if (fcst_cat == 0) { + if (obs_cat == 1) c12++; + else if(obs_cat == 2) c13++; + else count_diagonal++; + } + else if (fcst_cat == 1) { + if (obs_cat == 0) c21++; + else if(obs_cat == 2) c23++; + else count_diagonal++; + } + else if (fcst_cat == 2) { + if (obs_cat == 0) c31++; + else if (obs_cat == 1) c32++; + else count_diagonal++; + } + seeps_score = seeps_mpr->score; + if (isnan(seeps_score)) { + nan_count++; + seeps_score = bad_data_double; + } + else if (is_eq(seeps_score, bad_data_double)) bad_count++; + else { + seeps_count++; + obs_sum += obs_value; + fcst_sum += fcst_value; + seeps_score_partial_sum += seeps_score; + //IDL: s = s + c(4+cat(i) * w{i) <-- this is done by summing scores + //IDL: svf(cat{i)) = svf(cat{i)) + c(4+cat(i) * w{i) + //IDL: pvf(cat{i)) = pvf(cat{i)) + w{i) + //pvf[seeps_mpr->s_idx] += weight; + pvf_cnt[seeps_mpr->s_idx] += 1; + } + + delete seeps_mpr; + } + } + seeps_dp.set(seeps_score, ix, iy); + seeps_dp_fcat.set(fcst_cat, ix, iy); + seeps_dp_ocat.set(obs_cat, ix, iy); + } + seeps_score_sum += seeps_score_partial_sum; + } + int cell_count = dp_size - nan_count - bad_count; + if (cell_count > 0) { + seeps->weighted_score = seeps_score_sum/cell_count; + for (int i=0; in_obs = seeps_count; + seeps->c12 = c12; + seeps->c13 = c13; + seeps->c21 = c21; + seeps->c23 = c23; + seeps->c31 = c31; + seeps->c32 = c32; + if (seeps_count > 0) { + seeps->mean_fcst = fcst_sum / seeps_count; + seeps->mean_obs = obs_sum / seeps_count; + + seeps->pv1 = pvf[0] + pvf[3] + pvf[6]; // sum by column for obs + seeps->pv2 = pvf[1] + pvf[4] + pvf[7]; // sum by column for obs + seeps->pv3 = pvf[2] + pvf[5] + pvf[8]; // sum by column for obs + seeps->pf1 = pvf[0] + pvf[1] + pvf[2]; // sum by row for forecast + seeps->pf2 = pvf[3] + pvf[4] + pvf[5]; // sum by row for forecast + seeps->pf3 = pvf[6] + pvf[7] + pvf[8]; // sum by row for forecast + seeps->s12 = c12 * seeps->pf1 * seeps->pv2; + seeps->s13 = c13 * seeps->pf1 * seeps->pv3; + seeps->s21 = c21 * seeps->pf2 * seeps->pv1; + seeps->s23 = c23 * seeps->pf2 * seeps->pv3; + seeps->s31 = c31 * seeps->pf3 * seeps->pv1; + seeps->s32 = c32 * seeps->pf3 * seeps->pv2; + seeps->score = seeps_score_sum / seeps_count; + } + mlog << Debug(6) << method_name + << "SEEPS score=" << seeps->score << " weighted_score=" << seeps->weighted_score + << " pv1=" << seeps->pv1 << " pv2=" << seeps->pv2 << " pv3=" << seeps->pv3 + << " pf1=" << seeps->pf1 << " pf2=" << seeps->pf2 << " pf3=" << seeps->pf3 << "\n"; + if(mlog.verbosity_level() >= detailed_debug_level) { + char buffer[100]; + ConcatString log_message; + for (int i=0; i 0) { + snprintf ( buffer, 100, " nan: %d ", nan_count); + log_message.add(buffer); + } + mlog << Debug(7) << method_name << "pvf = " << log_message + << " weight=" << (1. / cell_count) << " (1/" << cell_count << ")" << "\n"; + } + +} + +//////////////////////////////////////////////////////////////////////// +// IDL code: +// +// my_pi = 3.1415926 +// r0 = density_radius*my_pi/180 +// lat = reform(float(data_geo(0,*))) +// lon = reform( ((float(data_geo(1,*)+360.) mod 360.))) +// n = n_elements(lat) +// lat = lat * my_pi / 180 +// slat = sin(lat) +// clat = cos(lat) +// lon = lon * my_pi / 180 +// slon = sin(lon) +// clon = cos(lon) +// r=(clat#transpose(clat))*(slon#transpose(slon)) + (clon#transpose(clon))*(slat#transpose(slat)) +// r=r*((r lt 1.) and (r gt -1.)) + (r ge 1.) - (r le -1.) +// r=acos(r) +// if r0 gt 0.0 then r = exp( -(r / r0) ^ 2) * ( r le 4. * r0) else r = (r*0.) + 1 +// r=sum(r,1) +// +// SUM function (not built-in) +// PRINT, array2 +// ; PV-WAVE prints the following: +// ; 0.00000 1.00000 +// ; 2.00000 3.00000 +// PRINT, SUM(array2, 0) +// ; PV-WAVE prints: 1.00000 5.00000 +// PRINT, SUM(array2, 1) +// ; PV-WAVE prints: 2.00000 4.00000 +//////////////////////////////////////////////////////////////////////// + +double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps) { + int seeps_idx; + SeepsScore *seeps_mpr; + int seeps_cnt = seeps->n_obs; + vector rlat(seeps_cnt), rlon(seeps_cnt); + vector clat(seeps_cnt), clon(seeps_cnt); + vector slat(seeps_cnt), slon(seeps_cnt); + vector< vector > clat_m(seeps_cnt, vector (seeps_cnt)); + vector< vector > clon_m(seeps_cnt, vector (seeps_cnt)); + vector< vector > slat_m(seeps_cnt, vector (seeps_cnt)); + vector< vector > slon_m(seeps_cnt, vector (seeps_cnt)); + vector< vector > clon_slat(seeps_cnt, vector (seeps_cnt)); + vector< vector > density_m(seeps_cnt, vector (seeps_cnt)); + static const char *method_name = "compute_seeps_density_vector() -> "; + + if (seeps_cnt == 0) { + mlog << Debug(1) << method_name + << "no SEEPS_MPR available.\n"; + return NULL; + } + + // Get lat/lon & convert them to radian and get sin/cos values + seeps_idx = 0; + double *density_vector = new double[seeps_cnt]; + for(int i=0; in_obs; i++) { + if (i >= pd->seeps_mpr.size()) break; + seeps_mpr = pd->seeps_mpr[i]; + if (!seeps_mpr || is_eq(seeps_mpr->score, bad_data_double)) continue; + + rlat[seeps_idx] = pd->lat_na[i] * rad_per_deg; // radian of lat + rlon[seeps_idx] = fmod((pd->lon_na[i] + 360.), 360.) * rad_per_deg; // radian of long + clat[seeps_idx] = cos(rlat[seeps_idx]); + clon[seeps_idx] = cos(rlon[seeps_idx]); + slat[seeps_idx] = sin(rlat[seeps_idx]); + slon[seeps_idx] = sin(rlon[seeps_idx]); + + seeps_idx++; + } + + // prooducs n by n matrix by multipling transposed vector + int v_count; + double density, mask1, mask2, mask3, temp; + // Initialize + v_count = 0; + if (seeps_idx < seeps_cnt) seeps_cnt = seeps_idx; + for(int i=0; i -1.0) ? 1. : 0.; + mask2 = (density >= 1.0 ) ? 1. : 0.; + mask3 = (density <= -1.0) ? 1. : 0.; + density = density * mask1 + mask2 - mask3; + //IDL: r = acos(r) + density = acos(density); + //IDL: if r0 gt 0.0 then r = exp(-(r/r0)^2) * (r le 4. * r0) else r = (r*0.)+1. + if (density_radius_rad <= 0.) density = 1.0; + else { + mask3 = (density <= 4.0) ? 1. : 0.; + temp = density / density_radius_rad; + density = exp(-(temp * temp)) * mask3 * density_radius_rad; + } + density_vector[j] += density; + } + if (!is_eq(density_vector[j], 0.)) v_count++; + } + if (v_count == 0) { + mlog << Debug(4) << method_name + << "no non-zero values at density_vector\n"; + } + if (seeps_cnt > 0) { + mlog << Debug(10) << method_name + << " non zero count=" << v_count + << " density_vector[0]=" << density_vector[0] + << " density_vector[" << (seeps_cnt-1) << "]=" << density_vector[seeps_cnt-1] << "\n"; + } + + return density_vector; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_statistics/compute_stats.h b/src/libcode/vx_statistics/compute_stats.h index 2d46624ad7..dc443b39b6 100644 --- a/src/libcode/vx_statistics/compute_stats.h +++ b/src/libcode/vx_statistics/compute_stats.h @@ -61,6 +61,11 @@ extern void compute_i_mean_stdev(const NumArray &, CIInfo &, CIInfo &); extern void compute_aggregated_seeps(const PairDataPoint *pd, SeepsAggScore *seeps); +extern double *compute_seeps_density_vector(const PairDataPoint *pd, SeepsAggScore *seeps); +extern void compute_aggregated_seeps_grid(const DataPlane &fcst_dp, const DataPlane &obs_dp, + DataPlane &seeps_dp, DataPlane &seeps_dp_fcat, + DataPlane &seeps_dp_ocat,SeepsAggScore *seeps, + int month, int hour, const SingleThresh &seeps_p1_thresh); //////////////////////////////////////////////////////////////////////// // diff --git a/src/libcode/vx_statistics/pair_data_point.cc b/src/libcode/vx_statistics/pair_data_point.cc index 48d8b6c068..ff7df4262e 100644 --- a/src/libcode/vx_statistics/pair_data_point.cc +++ b/src/libcode/vx_statistics/pair_data_point.cc @@ -69,6 +69,7 @@ PairDataPoint & PairDataPoint::operator=(const PairDataPoint &pd) { void PairDataPoint::init_from_scratch() { seeps_mpr.clear(); + seeps.clear(); clear(); seeps_climo = get_seeps_climo(); @@ -86,7 +87,7 @@ void PairDataPoint::clear() { if (seeps_mpr[idx]) delete seeps_mpr[idx]; } seeps_mpr.clear(); - seeps.init(); + seeps.clear(); return; } @@ -179,6 +180,12 @@ bool PairDataPoint::add_point_pair(const char *sid, double lat, double lon, //////////////////////////////////////////////////////////////////////// +void PairDataPoint::set_seeps_thresh(const SingleThresh &p1_thresh) { + seeps_climo->set_p1_thresh(p1_thresh); +} + +//////////////////////////////////////////////////////////////////////// + void PairDataPoint::set_seeps_score(SeepsScore *seeps, int index) { if (seeps) { int seeps_count = seeps_mpr.size(); @@ -1463,6 +1470,18 @@ void VxPairDataPoint::set_obs_perc_value(int percentile) { //////////////////////////////////////////////////////////////////////// +void VxPairDataPoint::set_seeps_thresh(const SingleThresh &p1_thresh) { + for(int i=0; i < n_msg_typ; i++){ + for(int j=0; j < n_mask; j++){ + for(int k=0; k < n_interp; k++){ + pd[i][j][k].set_seeps_thresh(p1_thresh); + } + } + } +} + +//////////////////////////////////////////////////////////////////////// + void VxPairDataPoint::print_obs_summary() { if(n_msg_typ == 0 || n_mask == 0 || n_interp == 0) { diff --git a/src/libcode/vx_statistics/pair_data_point.h b/src/libcode/vx_statistics/pair_data_point.h index 95576b3104..ab44e87d17 100644 --- a/src/libcode/vx_statistics/pair_data_point.h +++ b/src/libcode/vx_statistics/pair_data_point.h @@ -33,7 +33,7 @@ class PairDataPoint : public PairBase { void init_from_scratch(); void assign(const PairDataPoint &); - SeepsClimo *seeps_climo; + SeepsClimo *seeps_climo; public: PairDataPoint(); @@ -58,6 +58,7 @@ class PairDataPoint : public PairBase { bool add_point_pair(const char *, double, double, double, double, unixtime, double, double, double, double, const char *, double, double, double); + void set_seeps_thresh(const SingleThresh &p1_thresh); void set_seeps_score(SeepsScore *, int index=-1); void set_point_pair(int, const char *, double, double, double, double, @@ -226,7 +227,7 @@ class VxPairDataPoint { void set_mpr_thresh(const StringArray &, const ThreshArray &); - void set_seeps_thresh(const StringArray &, const ThreshArray &); + void set_seeps_thresh(const SingleThresh &p1_thresh); void set_climo_cdf_info_ptr(const ClimoCDFInfo *); diff --git a/src/libcode/vx_tc_util/diag_file.cc b/src/libcode/vx_tc_util/diag_file.cc index 32a37d2f07..3a4bced924 100644 --- a/src/libcode/vx_tc_util/diag_file.cc +++ b/src/libcode/vx_tc_util/diag_file.cc @@ -27,19 +27,17 @@ using namespace std; //////////////////////////////////////////////////////////////////////// -static const char default_lsdiag_technique[] = "OFCL"; - -static const int lsdiag_wdth[] = { +static const int ships_wdth[] = { 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 }; -static int n_lsdiag_wdth = sizeof(lsdiag_wdth)/sizeof(*lsdiag_wdth); +static int n_ships_wdth = sizeof(ships_wdth)/sizeof(*ships_wdth); -static const int tcdiag_fill_value = 9999; -static const int lsdiag_fill_value = 9999; +static const int cira_fill_value = 9999; +static const int ships_fill_value = 9999; -static const char lsdiag_skip_str[] = "TIME,DELV,MTPW,IR00,IRM1,IRM3,PC00,PCM1,PCM3,PSLV,IRXX"; +static const char ships_skip_str[] = "TIME,DELV,MTPW,IR00,IRM1,IRM3,PC00,PCM1,PCM3,PSLV,IRXX"; //////////////////////////////////////////////////////////////////////// // @@ -167,7 +165,9 @@ void DiagFile::init_from_scratch() { void DiagFile::clear() { // Initialize values - SourceType = DiagType_None; + DiagSource = DiagType_None; + TrackSource.clear(); + FieldSource.clear(); StormId.clear(); Basin.clear(); Cyclone.clear(); @@ -209,30 +209,38 @@ void DiagFile::add_technique(const string &str) { //////////////////////////////////////////////////////////////////////// -void DiagFile::read(const DiagType source, - const ConcatString &path, const StringArray &model_names, - const std::map *convert_map) { +void DiagFile::read(const ConcatString &path, const ConcatString &diag_source, + const ConcatString &track_source, const ConcatString &field_source, + const StringArray &match_to_track, + const std::map *convert_map) { + + DiagType type = string_to_diagtype(diag_source.c_str()); - // Read diagnostics baesd on the source type - if(source == TCDiagType) { - read_tcdiag(path, model_names, convert_map); + // Read diagnostics based on the source type + if(type == DiagType_CIRA_RT) { + read_cira_rt(path, convert_map); } - else if(source == LSDiagRTType) { - read_lsdiag_rt(path, model_names, convert_map); + else if(type == DiagType_SHIPS_RT) { + read_ships_rt(path, convert_map); } else { mlog << Error << "\nDiagFile::read() -> " - << "diagnostics of type " << diagtype_to_string(source) - << " are not currently supported!\n\n"; + << "diagnostics of type \"" << diag_source + << "\" are not currently supported!\n\n"; exit(1); } + // Store the metadata + TrackSource = track_source; + FieldSource = field_source; + set_technique(match_to_track); + return; } //////////////////////////////////////////////////////////////////////// -void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_names, +void DiagFile::read_cira_rt(const ConcatString &path, const map *convert_map) { int i; double v_in, v_out; @@ -243,10 +251,7 @@ void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_na clear(); // Store the file type - SourceType = TCDiagType; - - // Store user-specified model names or read it from the file below - if(model_names.n() > 0) set_technique(model_names); + DiagSource = DiagType_CIRA_RT; // Open the file open(path.c_str()); @@ -306,7 +311,7 @@ void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_na // Check for the expected number of items if(NTime != LeadTime.n() || NTime != Lat.n() || NTime != Lon.n()) { - mlog << Error << "\nDiagFile::read_tcdiag() -> " + mlog << Error << "\nDiagFile::read_cira_rt() -> " << "the NTIME value (" << NTime << ") does not match the actual number of times (" << LeadTime.n() << "), latitudes (" << Lat.n() @@ -350,7 +355,7 @@ void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_na v_in = atof(dl[i]); // Check for bad data and apply conversions - if(atoi(dl[i]) == tcdiag_fill_value) v_out = bad_data_double; + if(atoi(dl[i]) == cira_fill_value) v_out = bad_data_double; else if(fx_ptr) v_out = (*fx_ptr)(v_in); else v_out = v_in; @@ -360,7 +365,7 @@ void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_na // Check for the expected number of items if(NTime != data.n()) { - mlog << Error << "\nDiagFile::read_tcdiag() -> " + mlog << Error << "\nDiagFile::read_cira_rt() -> " << "the number of \"" << cs << "\" diagnostic values (" << data.n() << ") does not match the expected number (" << NTime << "): " << path << "\n\n"; @@ -385,8 +390,8 @@ void DiagFile::read_tcdiag(const ConcatString &path, const StringArray &model_na //////////////////////////////////////////////////////////////////////// -void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model_names, - const map *convert_map) { +void DiagFile::read_ships_rt(const ConcatString &path, + const map *convert_map) { int i, v_int; double v_dbl; NumArray data; @@ -396,18 +401,14 @@ void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model clear(); // Store the file type - SourceType = LSDiagRTType; - - // Store user-specified model names or the default one - if(model_names.n() > 0) set_technique(model_names); - else add_technique(default_lsdiag_technique); + DiagSource = DiagType_SHIPS_RT; // Open the file open(path.c_str()); // Diagnostic names to ignore StringArray skip_diag_sa; - skip_diag_sa.parse_css(lsdiag_skip_str); + skip_diag_sa.parse_css(ships_skip_str); // Parse the header information from the first line DataLine dl; @@ -416,7 +417,7 @@ void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model // Check for the expected number of items if(dl.n_items() != 9 || strncasecmp(dl[8], "HEAD", strlen("HEAD") != 0)) { - mlog << Error << "\nDiagFile::read_lsdiag_rt() -> " + mlog << Error << "\nDiagFile::read_ships_rt() -> " << "unexpected header line: " << path << "\n\n"; exit(1); } @@ -438,7 +439,7 @@ void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model int data_start_location = in->tellg(); // Parse time and location info - while(read_fwf_line(dl, lsdiag_wdth, n_lsdiag_wdth)) { + while(read_fwf_line(dl, ships_wdth, n_ships_wdth)) { // Fixed width column 24 has the data name cs = dl[23]; @@ -473,7 +474,7 @@ void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model in->seekg(data_start_location); // Store the diagnostics data - while(read_fwf_line(dl, lsdiag_wdth, n_lsdiag_wdth)) { + while(read_fwf_line(dl, ships_wdth, n_ships_wdth)) { // Skip empty lines if(dl.n_items() == 0) continue; @@ -506,7 +507,7 @@ void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model v_int = atoi(dl[i]); // Check for bad data and apply conversions - if(v_int == lsdiag_fill_value) v_dbl = bad_data_double; + if(v_int == ships_fill_value) v_dbl = bad_data_double; else if(fx_ptr) v_dbl = (*fx_ptr)(v_int); else v_dbl = (double) v_int; @@ -516,7 +517,7 @@ void DiagFile::read_lsdiag_rt(const ConcatString &path, const StringArray &model // Check for the expected number of items if(NTime != data.n()) { - mlog << Error << "\nDiagFile::read_lsdiag_rt() -> " + mlog << Error << "\nDiagFile::read_ships_rt() -> " << "the number of \"" << cs << "\" diagnostic values (" << data.n() << ") does not match the expected number (" << NTime << ")!\n\n"; diff --git a/src/libcode/vx_tc_util/diag_file.h b/src/libcode/vx_tc_util/diag_file.h index 6db6848870..c9109f9475 100644 --- a/src/libcode/vx_tc_util/diag_file.h +++ b/src/libcode/vx_tc_util/diag_file.h @@ -23,7 +23,7 @@ //////////////////////////////////////////////////////////////////////// // -// TCDIAG files: +// Realtime CIRA Diagnostics files: // - Add link to sample data // - Header: // * ATCF_ID YYYYMMDDHH * @@ -33,7 +33,7 @@ // BB is the 2-letter basin name // CC is the 2-digit cyclone number // -// Real-time LSDIAG files: +// Real-time SHIPS Diagnostics files: // - https://ftp.nhc.noaa.gov/atcf/lsdiag // - Header: // BBCC YYMMDD HH WS LAT LON 9999 BBCCYYYY @@ -42,7 +42,7 @@ // YYMMDD HH is the initialization time // WS is the wind speed // -// Developmental LSDIAG files (not currently supported): +// Developmental SHIPS Diagnostics files (not currently supported): // - https://rammb2.cira.colostate.edu/research/tropical-cyclones/ships // //////////////////////////////////////////////////////////////////////// @@ -58,8 +58,10 @@ class DiagFile : public LineDataFile { void init_from_scratch(); - // Diagnostics file type - DiagType SourceType; + // Diagnostics Metadata + DiagType DiagSource; + ConcatString TrackSource; + ConcatString FieldSource; // Storm and model identification ConcatString StormId; @@ -108,7 +110,9 @@ class DiagFile : public LineDataFile { double lat(int) const; double lon(int) const; - DiagType source() const; + DiagType diag_source() const; + const ConcatString & track_source() const; + const ConcatString & field_source() const; int n_diag() const; const StringArray & diag_name() const; bool has_diag(const std::string &) const; @@ -118,27 +122,30 @@ class DiagFile : public LineDataFile { // do stuff // - void read (const DiagType, - const ConcatString &, const StringArray &, - const std::map *); - void read_tcdiag (const ConcatString &, const StringArray &, - const std::map *); - void read_lsdiag_rt (const ConcatString &, const StringArray &, - const std::map *); + void read (const ConcatString &, const ConcatString &, + const ConcatString &, const ConcatString &, + const StringArray &, + const std::map *); + void read_cira_rt (const ConcatString &, + const std::map *); + void read_ships_rt (const ConcatString &, + const std::map *); }; //////////////////////////////////////////////////////////////////////// -inline const ConcatString & DiagFile::storm_id() const { return(StormId); } -inline const ConcatString & DiagFile::basin() const { return(Basin); } -inline const ConcatString & DiagFile::cyclone() const { return(Cyclone); } -inline const StringArray & DiagFile::technique() const { return(Technique); } -inline unixtime DiagFile::init() const { return(InitTime); } -inline int DiagFile::n_time() const { return(NTime); } -inline DiagType DiagFile::source() const { return(SourceType); } -inline int DiagFile::n_diag() const { return(DiagName.n()); } -inline const StringArray & DiagFile::diag_name() const { return(DiagName); } +inline const ConcatString & DiagFile::storm_id() const { return(StormId); } +inline const ConcatString & DiagFile::basin() const { return(Basin); } +inline const ConcatString & DiagFile::cyclone() const { return(Cyclone); } +inline const StringArray & DiagFile::technique() const { return(Technique); } +inline unixtime DiagFile::init() const { return(InitTime); } +inline int DiagFile::n_time() const { return(NTime); } +inline DiagType DiagFile::diag_source() const { return(DiagSource); } +inline const ConcatString & DiagFile::track_source() const { return(TrackSource); } +inline const ConcatString & DiagFile::field_source() const { return(FieldSource); } +inline int DiagFile::n_diag() const { return(DiagName.n()); } +inline const StringArray & DiagFile::diag_name() const { return(DiagName); } //////////////////////////////////////////////////////////////////////// diff --git a/src/libcode/vx_tc_util/tc_columns.cc b/src/libcode/vx_tc_util/tc_columns.cc index 4cb90895f4..b561c448c9 100644 --- a/src/libcode/vx_tc_util/tc_columns.cc +++ b/src/libcode/vx_tc_util/tc_columns.cc @@ -342,8 +342,8 @@ void write_tc_mpr_cols(const TrackPairInfo &p, int i, at.set_entry(r, c++, systemsdepth_to_string(p.bdeck()[i].depth())); at.set_entry(r, c++, p.adeck()[i].num_members()); - at.set_entry(r, c++, p.adeck()[i].spread()); - at.set_entry(r, c++, p.adeck()[i].dist_mean()); + at.set_entry(r, c++, p.adeck()[i].track_spread()); + at.set_entry(r, c++, p.adeck()[i].track_stdev()); at.set_entry(r, c++, p.adeck()[i].mslp_stdev()); at.set_entry(r, c++, p.adeck()[i].v_max_stdev()); @@ -360,6 +360,8 @@ void write_tc_diag_cols(const TrackPairInfo &p, int i, at.set_entry(r, c++, p.n_points()); at.set_entry(r, c++, i+1); at.set_entry(r, c++, diagtype_to_string(p.adeck().diag_source())); + at.set_entry(r, c++, p.adeck().track_source()); + at.set_entry(r, c++, p.adeck().field_source()); at.set_entry(r, c++, p.adeck()[i].n_diag()); // Check the number of names and values match diff --git a/src/libcode/vx_tc_util/tc_columns.h b/src/libcode/vx_tc_util/tc_columns.h index 6910590efc..8236e141df 100644 --- a/src/libcode/vx_tc_util/tc_columns.h +++ b/src/libcode/vx_tc_util/tc_columns.h @@ -70,10 +70,8 @@ static const char * tc_mpr_cols [] = { "ASPEED", "BSPEED", "ADEPTH", "BDEPTH", "NUM_MEMBERS", - "TRACK_SPREAD", - "DIST_MEAN", - "MSLP_SPREAD", - "MAX_WIND_SPREAD" + "TRACK_SPREAD", "TRACK_STDEV", + "MSLP_STDEV", "MAX_WIND_STDEV" }; static const int n_tc_mpr_cols = sizeof(tc_mpr_cols)/sizeof(*tc_mpr_cols); @@ -122,8 +120,9 @@ static const int n_tc_cols_xy = sizeof(tc_cols_xy)/sizeof(*tc_cols_xy); //////////////////////////////////////////////////////////////////////// static const char * tc_diag_cols [] = { - "TOTAL", "INDEX", "DIAG_SOURCE", - "N_DIAG", "DIAG_", "VALUE_" + "TOTAL", "INDEX", + "DIAG_SOURCE", "TRACK_SOURCE", "FIELD_SOURCE", + "N_DIAG", "DIAG_", "VALUE_" }; static const int n_tc_diag_cols = sizeof(tc_diag_cols)/sizeof(*tc_diag_cols); diff --git a/src/libcode/vx_tc_util/track_info.cc b/src/libcode/vx_tc_util/track_info.cc index cddb2b1536..f99d95cf3e 100644 --- a/src/libcode/vx_tc_util/track_info.cc +++ b/src/libcode/vx_tc_util/track_info.cc @@ -92,6 +92,8 @@ void TrackInfo::clear() { MinWarmCore = (unixtime) 0; MaxWarmCore = (unixtime) 0; DiagSource = DiagType_None; + TrackSource.clear(); + FieldSource.clear(); DiagName.clear(); TrackLines.clear(); @@ -133,6 +135,8 @@ void TrackInfo::dump(ostream &out, int indent_depth) const { out << prefix << "MinWarmCore = \"" << (MinWarmCore > 0 ? unix_to_yyyymmdd_hhmmss(MinWarmCore).text() : na_str) << "\n"; out << prefix << "MaxWarmCore = \"" << (MaxWarmCore > 0 ? unix_to_yyyymmdd_hhmmss(MaxWarmCore).text() : na_str) << "\n"; out << prefix << "DiagSource = " << diagtype_to_string(DiagSource) << "\n"; + out << prefix << "TrackSource = " << TrackSource.contents() << "\n"; + out << prefix << "FieldSource = " << FieldSource.contents() << "\n"; out << prefix << "NDiag = " << DiagName.n() << "\n"; out << prefix << "NPoints = " << NPoints << "\n"; out << prefix << "NAlloc = " << NAlloc << "\n"; @@ -172,6 +176,8 @@ ConcatString TrackInfo::serialize() const { << ", MinWarmCore = " << (MinWarmCore > 0 ? unix_to_yyyymmdd_hhmmss(MinWarmCore).text() : na_str) << ", MaxWarmCore = " << (MaxWarmCore > 0 ? unix_to_yyyymmdd_hhmmss(MaxWarmCore).text() : na_str) << ", DiagSource = " << diagtype_to_string(DiagSource) + << ", TrackSource = " << TrackSource.contents() + << ", FieldSource = " << FieldSource.contents() << ", NDiag = " << DiagName.n() << ", NPoints = " << NPoints << ", NAlloc = " << NAlloc @@ -223,6 +229,8 @@ void TrackInfo::assign(const TrackInfo &t) { MinWarmCore = t.MinWarmCore; MaxWarmCore = t.MaxWarmCore; DiagSource = t.DiagSource; + TrackSource = t.TrackSource; + FieldSource = t.FieldSource; DiagName = t.DiagName; TrackLines = t.TrackLines; @@ -534,8 +542,10 @@ bool TrackInfo::add_diag_data(DiagFile &diag_file, const StringArray &req_diag_n InitTime != diag_file.init() || !diag_file.technique().has(Technique)) return(false); - // Store the diagnostic source - DiagSource = diag_file.source(); + // Store the diagnostic metadata + DiagSource = diag_file.diag_source(); + TrackSource = diag_file.track_source(); + FieldSource = diag_file.field_source(); // If empty, store all diagnostics bool store_all_diag = (req_diag_name.n() == 0 ? true : false); @@ -563,12 +573,21 @@ bool TrackInfo::add_diag_data(DiagFile &diag_file, const StringArray &req_diag_n if(i_name == 0) { if(!is_eq(diag_file.lat(i_time), Point[i_pnt].lat()) || !is_eq(diag_file.lon(i_time), Point[i_pnt].lon())) { - mlog << Warning << "\nTrackInfo::add_diag_data() -> " - << "the " << StormId << " " << Technique << " " << unix_to_yyyymmddhh(InitTime) + ConcatString cs; + cs << "The " << StormId << " " << Technique << " " << unix_to_yyyymmddhh(InitTime) << " lead time " << sec_to_timestring(diag_file.lead(i_time)) << " track location (" << Point[i_pnt].lat() << ", " - << Point[i_pnt].lon() << ") does not match the diagnostic location (" - << diag_file.lat(i_time) << ", " << diag_file.lon(i_time) << ")\n"; + << Point[i_pnt].lon() << ") does not match the " + << diagtype_to_string(DiagSource) << " diagnostic location (" + << diag_file.lat(i_time) << ", " << diag_file.lon(i_time) << ")"; + + // Print a warning if the TrackSource and Technique match + if(TrackSource == Technique) { + mlog << Warning << "\nTrackInfo::add_diag_data() -> " << cs << "\n\n"; + } + else { + mlog << Debug(4) << cs << "\n"; + } } } @@ -957,7 +976,7 @@ TrackInfo consensus(const TrackInfoArray &tracks, QuadInfo wavg; NumArray plon, plat, pvmax, pmslp; double lon_range, lon_shift, lon_avg; - double track_spread, vmax_stdev, mslp_stdev; + double track_spread, track_stdev, vmax_stdev, mslp_stdev; // Check for at least one track if(tracks.n() == 0) { @@ -1082,18 +1101,15 @@ TrackInfo consensus(const TrackInfoArray &tracks, // Save the number of members that went into the consensus if(pcnt > 0) pavg.set_num_members(pcnt); - // Compute track spread and distance mean, convert to nautical-miles - double track_spread, dist_mean; - compute_gc_dist_stdev(pavg.lat(), pavg.lon(), plat, plon, track_spread, dist_mean); + // Compute track mean and standard deviation, convert to nautical-miles + compute_gc_dist_stdev(pavg.lat(), pavg.lon(), plat, plon, track_spread, track_stdev); if(!is_bad_data(track_spread)) { - track_spread *= tc_nautical_miles_per_km; - pavg.set_spread(track_spread); + pavg.set_track_spread(track_spread*tc_nautical_miles_per_km); } - if(!is_bad_data(dist_mean)) { - dist_mean *= tc_nautical_miles_per_km; - pavg.set_dist_mean(dist_mean); + if(!is_bad_data(track_stdev)) { + pavg.set_track_stdev(track_stdev*tc_nautical_miles_per_km); } // Compute wind-speed (v_max) and pressure (mslp) standard deviation @@ -1132,28 +1148,20 @@ TrackInfo consensus(const TrackInfoArray &tracks, //////////////////////////////////////////////////////////////////////// -void compute_gc_dist_stdev(const double lat, const double lon, const NumArray &lats, const NumArray &lons, double &spread, double &mean) { - - int i, count; +void compute_gc_dist_stdev(const double lat, const double lon, + const NumArray &lats, const NumArray &lons, + double &mean, double &stdev) { NumArray dist_na; // Loop over member lat/lon track values, calculate great-circle distance between memmber values and consensus track - for(i=0, count=0; i " << "the diagnostic source type has changed (" << diagtype_to_string(ADeck.diag_source()) << " != " - << diagtype_to_string(t) << ")!\n\n"; + << diagtype_to_string(diag_source) << ")!\n\n"; + exit(1); + } + + // Make sure TRACK_SOURCE does not change + ConcatString track_source = l.get_item("TRACK_SOURCE"); + if(ADeck.track_source().length() > 0 && + ADeck.track_source() != track_source) { + mlog << Error << "\nTrackPairInfo::add_tcdiag_line() -> " + << "the diagnostic track source has changed (" + << ADeck.track_source() << " != " + << track_source << ")!\n\n"; + exit(1); + } + + // Make sure FIELD_SOURCE does not change + ConcatString field_source = l.get_item("FIELD_SOURCE"); + if(ADeck.field_source().length() > 0 && + ADeck.field_source() != field_source) { + mlog << Error << "\nTrackPairInfo::add_tcdiag_line() -> " + << "the diagnostic field source has changed (" + << ADeck.field_source() << " != " + << field_source << ")!\n\n"; exit(1); } - // Store the source type - ADeck.set_diag_source(t); + // Store the diagnostic metadata + ADeck.set_diag_source(diag_source); + ADeck.set_track_source(track_source.c_str()); + ADeck.set_field_source(field_source.c_str()); // Number of diagnostics n_diag = atoi(l.get_item("N_DIAG")); diff --git a/src/libcode/vx_tc_util/track_point.cc b/src/libcode/vx_tc_util/track_point.cc index 0badcb827c..541a5ec1cd 100644 --- a/src/libcode/vx_tc_util/track_point.cc +++ b/src/libcode/vx_tc_util/track_point.cc @@ -381,11 +381,11 @@ TrackPoint & TrackPoint::operator+=(const TrackPoint &p) { else Speed += p.speed(); // Set consensus (spread) variables to missing - NumMembers = bad_data_int; - Spread = bad_data_double; - DistMean = bad_data_double; - VmaxStdev = bad_data_double; - MSLPStdev = bad_data_double; + NumMembers = bad_data_int; + TrackSpread = bad_data_double; + TrackStdev = bad_data_double; + VmaxStdev = bad_data_double; + MSLPStdev = bad_data_double; // Increment wind quadrants for(i=0; iis_precipitation() + && conf_info.vx_opt[i].obs_info->is_precipitation()) { + SeepsAggScore seeps; + int month, day, year, hour, minute, second; + + unix_to_mdyhms(fcst_dp.valid(), month, day, year, hour, minute, second); + compute_aggregated_seeps_grid(fcst_dp_smooth, obs_dp_smooth, + seeps_dp, seeps_dp_fcat, seeps_dp_ocat, + &seeps, month, hour, + conf_info.vx_opt[i].seeps_p1_thresh); + + write_nc("SEEPS_MPR_SCORE", seeps_dp, i, mthd, pnts, + conf_info.vx_opt[i].interp_info.field); + write_nc("SEEPS_MPR_FCAT", seeps_dp_fcat, i, mthd, pnts, + conf_info.vx_opt[i].interp_info.field); + write_nc("SEEPS_MPR_OCAT", seeps_dp_ocat, i, mthd, pnts, + conf_info.vx_opt[i].interp_info.field); + write_seeps_row(shc, &seeps, conf_info.output_flag[i_seeps], + stat_at, i_stat_row, txt_at[i_seeps], i_txt_row[i_seeps]); + } + // Compute gradient statistics if requested in the config file if(!conf_info.vx_opt[i].fcst_info->is_prob() && conf_info.vx_opt[i].output_flag[i_grad] != STATOutputType_None) { @@ -2660,6 +2687,26 @@ void write_nc(const ConcatString &field_name, const DataPlane &dp, level_att = shc.get_obs_lev(); units_att = conf_info.vx_opt[i_vx].obs_info->units_attr(); } + else if(strncmp(field_name.c_str(), "SEEPS_MPR", 9) == 0) { + ConcatString seeps_desc; + var_name << cs_erase << field_name << "_" + << obs_name << var_suffix << "_" << mask_str; + if(field_type == FieldType_Obs || + field_type == FieldType_Both) { + var_name << interp_str; + } + if(strncmp(field_name.c_str(), "SEEPS_MPR_SCORE", 15) == 0) + seeps_desc = "score"; + else if(strncmp(field_name.c_str(), "SEEPS_MPR_FCAT", 14) == 0) + seeps_desc = "forecast category"; + else if(strncmp(field_name.c_str(), "SEEPS_MPR_OCAT", 14) == 0) + seeps_desc = "observation category"; + long_att << cs_erase + << "SEEPS MPR " << seeps_desc << " for " + << obs_long_name; + level_att = shc.get_obs_lev(); + units_att = ""; + } else { mlog << Error << "\nwrite_nc() -> " << "unexpected field name \"" << field_name << "\"\n\n"; diff --git a/src/tools/core/grid_stat/grid_stat.h b/src/tools/core/grid_stat/grid_stat.h index 555941e10c..3154d5f834 100644 --- a/src/tools/core/grid_stat/grid_stat.h +++ b/src/tools/core/grid_stat/grid_stat.h @@ -73,7 +73,8 @@ static const char **txt_columns[n_txt] = { val1l2_columns, pct_columns, pstd_columns, pjc_columns, prc_columns, eclv_columns, nbrctc_columns, nbrcts_columns, nbrcnt_columns, - grad_columns, vcnt_columns, dmap_columns + grad_columns, vcnt_columns, dmap_columns, + seeps_columns }; // Length of header columns @@ -84,7 +85,8 @@ static const int n_txt_columns[n_txt] = { n_val1l2_columns, n_pct_columns, n_pstd_columns, n_pjc_columns, n_prc_columns, n_eclv_columns, n_nbrctc_columns, n_nbrcts_columns, n_nbrcnt_columns, - n_grad_columns, n_vcnt_columns, n_dmap_columns + n_grad_columns, n_vcnt_columns, n_dmap_columns, + n_seeps_columns }; // Text file abbreviations @@ -95,7 +97,8 @@ static const char *txt_file_abbr[n_txt] = { "val1l2", "pct", "pstd", "pjc", "prc", "eclv", "nbrctc", "nbrcts", "nbrcnt", - "grad", "vcnt", "dmap" + "grad", "vcnt", "dmap", + "seeps" }; //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/grid_stat/grid_stat_conf_info.cc b/src/tools/core/grid_stat/grid_stat_conf_info.cc index 7cdae0c351..2a5f447972 100644 --- a/src/tools/core/grid_stat/grid_stat_conf_info.cc +++ b/src/tools/core/grid_stat/grid_stat_conf_info.cc @@ -538,6 +538,8 @@ void GridStatVxOpt::clear() { hss_ec_value = bad_data_double; rank_corr_flag = false; + seeps_p1_thresh.clear(); + for(i=0; i " << "unexpected output type index value: " << i_txt_row diff --git a/src/tools/core/grid_stat/grid_stat_conf_info.h b/src/tools/core/grid_stat/grid_stat_conf_info.h index c75693bcde..d8d9281c8a 100644 --- a/src/tools/core/grid_stat/grid_stat_conf_info.h +++ b/src/tools/core/grid_stat/grid_stat_conf_info.h @@ -59,7 +59,9 @@ static const int i_vcnt = 19; static const int i_dmap = 20; -static const int n_txt = 21; +static const int i_seeps = 21; + +static const int n_txt = 22; // Text file type static const STATLineType txt_file_type[n_txt] = { @@ -90,7 +92,8 @@ static const STATLineType txt_file_type[n_txt] = { stat_grad, // 18 stat_vcnt, // 19 - stat_dmap // 20 + stat_dmap, // 20 + stat_seeps // 21 }; @@ -159,6 +162,8 @@ class GridStatVxOpt { ThreshArray owind_ta; // obs wind speed thresholds SetLogic wind_logic; // wind speed field logic + SingleThresh seeps_p1_thresh; // SEESP p1 threshold + StringArray mask_grid; // Masking grid strings StringArray mask_poly; // Masking polyline strings diff --git a/src/tools/core/point_stat/point_stat.cc b/src/tools/core/point_stat/point_stat.cc index 257c06e86f..2da8c28106 100644 --- a/src/tools/core/point_stat/point_stat.cc +++ b/src/tools/core/point_stat/point_stat.cc @@ -166,7 +166,6 @@ static void do_vl1l2 (VL1L2Info *&, int, const PairDataPoint *, const PairDa static void do_pct (const PointStatVxOpt &, const PairDataPoint *); static void do_hira_ens ( int, const PairDataPoint *); static void do_hira_prob ( int, const PairDataPoint *); -static void do_seeps_agg (const PairDataPoint *); static void finish_txt_files(); @@ -1054,7 +1053,7 @@ void process_scores() { // Write out the SEEPS lines if(conf_info.vx_opt[i].output_flag[i_seeps] != STATOutputType_None) { compute_aggregated_seeps(pd_ptr, &pd_ptr->seeps); - write_seeps_row(shc, pd_ptr, + write_seeps_row(shc, &pd_ptr->seeps, conf_info.vx_opt[i].output_flag[i_seeps], stat_at, i_stat_row, txt_at[i_seeps], i_txt_row[i_seeps]); diff --git a/src/tools/core/point_stat/point_stat_conf_info.cc b/src/tools/core/point_stat/point_stat_conf_info.cc index e16bc6bad4..9d7c7fec2e 100644 --- a/src/tools/core/point_stat/point_stat_conf_info.cc +++ b/src/tools/core/point_stat/point_stat_conf_info.cc @@ -670,6 +670,8 @@ void PointStatVxOpt::clear() { msg_typ.clear(); + seeps_p1_thresh.clear(); + duplicate_flag = DuplicateType_None; obs_summary = ObsSummary_None; obs_perc = bad_data_int; @@ -916,6 +918,9 @@ void PointStatVxOpt::process_config(GrdFileType ftype, // Conf: rank_corr_flag rank_corr_flag = odict.lookup_bool(conf_key_rank_corr_flag); + // Conf: threshold for SEEPS p1 + seeps_p1_thresh = odict.lookup_thresh(conf_key_seeps_p1_thresh); + // Conf: message_type msg_typ = parse_conf_message_type(&odict); @@ -1079,7 +1084,7 @@ void PointStatVxOpt::set_vx_pd(PointStatConfInfo *conf_info) { vx_pd.set_duplicate_flag(duplicate_flag); vx_pd.set_obs_summary(obs_summary); vx_pd.set_obs_perc_value(obs_perc); - + vx_pd.set_seeps_thresh(seeps_p1_thresh); return; } diff --git a/src/tools/core/point_stat/point_stat_conf_info.h b/src/tools/core/point_stat/point_stat_conf_info.h index ec27d61cc1..5c4d29df62 100644 --- a/src/tools/core/point_stat/point_stat_conf_info.h +++ b/src/tools/core/point_stat/point_stat_conf_info.h @@ -131,6 +131,7 @@ class PointStatVxOpt { StringArray mpr_sa; // MPR column names ThreshArray mpr_ta; // MPR column thresholds + SingleThresh seeps_p1_thresh; // SEESP p1 threshold // Vector of MaskLatLon objects defining Lat/Lon Point masks vector mask_llpnt; diff --git a/src/tools/core/stat_analysis/aggr_stat_line.cc b/src/tools/core/stat_analysis/aggr_stat_line.cc index fe4a55d2ec..8df6643c46 100644 --- a/src/tools/core/stat_analysis/aggr_stat_line.cc +++ b/src/tools/core/stat_analysis/aggr_stat_line.cc @@ -35,6 +35,8 @@ // 015 01/24/20 Halley Gotway Add aggregate RPS lines. // 016 04/12/21 Halley Gotway MET #1735 Support multiple // -out_thresh and -out_line_type options. +// 017 11/10/22 Halley Gotway MET #2339 Add SEEPS and SEEPS_MPR +// line types. // //////////////////////////////////////////////////////////////////////// @@ -3356,6 +3358,184 @@ void aggr_ssvar_lines(LineDataFile &f, STATAnalysisJob &job, //////////////////////////////////////////////////////////////////////// +void aggr_seeps_lines(LineDataFile &f, STATAnalysisJob &job, + map &m, + int &n_in, int &n_out) { + STATLine line; + AggrSEEPSInfo aggr; + SeepsAggScore cur; + ConcatString key, fcst_var, obs_var; + + // + // Process the STAT lines + // + while(f >> line) { + + if(line.is_header()) continue; + + n_in++; + + if(job.is_keeper(line)) { + + job.dump_stat_line(line); + + if(line.type() != stat_seeps) { + mlog << Error << "\naggr_seeps_lines() -> " + << "should only encounter SEEPS line types.\n" + << "ERROR occurred on STAT line:\n" << line << "\n\n"; + throw(1); + } + + // + // Only aggregate consistent variable names + // + if(fcst_var.empty()) fcst_var = line.fcst_var(); + if(obs_var.empty()) obs_var = line.obs_var(); + + if(fcst_var != line.fcst_var() || + obs_var != line.obs_var()) { + mlog << Error << "\naggr_seeps_lines() -> " + << "both the forecast and observation variable types must " + << "remain constant. Try setting \"-fcst_var\" and/or " + << "\"-obs_var\".\n" + << "ERROR occurred on STAT line:\n" << line << "\n\n"; + throw(1); + } + + // + // Parse the current SSVAR line + // + parse_seeps_line(line, cur); + + // + // Build the map key for the current line + // + key = job.get_case_info(line); + + // + // Add a new map entry, if necessary + // + if(m.count(key) == 0) { + aggr.agg_score = cur; + aggr.hdr.clear(); + m[key] = aggr; + } + // + // Increment counts in the existing map entry + // + else { + m[key].agg_score += cur; + } + + // + // Keep track of the unique header column entries + // + m[key].hdr.add(line); + + n_out++; + } + } // end while + + return; +} + + +//////////////////////////////////////////////////////////////////////// + +void aggr_seeps_mpr_lines(LineDataFile &f, STATAnalysisJob &job, + map &m, + int &n_in, int &n_out) { + STATLine line; + AggrSEEPSMPRInfo aggr; + SEEPSMPRData cur; + ConcatString key, fcst_var, obs_var; + + // + // Process the STAT lines + // + while(f >> line) { + + if(line.is_header()) continue; + + n_in++; + + if(job.is_keeper(line)) { + + job.dump_stat_line(line); + + if(line.type() != stat_seeps_mpr) { + mlog << Error << "\naggr_seeps_mpr_lines() -> " + << "should only encounter SEEPS_MPR line types.\n" + << "ERROR occurred on STAT line:\n" << line << "\n\n"; + throw(1); + } + + // + // Only aggregate consistent variable names + // + if(fcst_var.empty()) fcst_var = line.fcst_var(); + if(obs_var.empty()) obs_var = line.obs_var(); + + if(fcst_var != line.fcst_var() || + obs_var != line.obs_var()) { + mlog << Error << "\naggr_seeps_mpr_lines() -> " + << "both the forecast and observation variable types must " + << "remain constant. Try setting \"-fcst_var\" and/or " + << "\"-obs_var\".\n" + << "ERROR occurred on STAT line:\n" << line << "\n\n"; + throw(1); + } + + // + // Parse the current SSVAR line + // + parse_seeps_mpr_line(line, cur); + + // + // Build the map key for the current line + // + key = job.get_case_info(line); + + // + // Add a new map entry, if necessary + // + if(m.count(key) == 0) { + aggr.pd.n_obs = 0; + aggr.pd.f_na.clear(); + aggr.pd.o_na.clear(); + aggr.pd.seeps_mpr.clear(); + aggr.hdr.clear(); + m[key] = aggr; + } + + // + // Increment counts + // + m[key].pd.n_obs += 1; + m[key].pd.f_na.add(cur.fcst); + m[key].pd.o_na.add(cur.obs); + m[key].pd.lat_na.add(cur.obs_lat); + m[key].pd.lon_na.add(cur.obs_lon); + + // Allocated here but deallocated by PairDataPoint + SeepsScore *score = new SeepsScore; + *score = cur.seeps_mpr; + m[key].pd.seeps_mpr.push_back(score); + + // + // Keep track of the unique header column entries + // + m[key].hdr.add(line); + + n_out++; + } + } // end while + + return; +} + +//////////////////////////////////////////////////////////////////////// + void aggr_time_series_lines(LineDataFile &f, STATAnalysisJob &job, map &m, int &n_in, int &n_out) { diff --git a/src/tools/core/stat_analysis/aggr_stat_line.h b/src/tools/core/stat_analysis/aggr_stat_line.h index 0297334f3d..f5e9616e16 100644 --- a/src/tools/core/stat_analysis/aggr_stat_line.h +++ b/src/tools/core/stat_analysis/aggr_stat_line.h @@ -192,6 +192,16 @@ struct AggrSSVARInfo { std::map ssvar_bins; }; +struct AggrSEEPSInfo { + StatHdrInfo hdr; + SeepsAggScore agg_score; +}; + +struct AggrSEEPSMPRInfo { + StatHdrInfo hdr; + PairDataPoint pd; +}; + struct AggrTimeSeriesInfo { StatHdrInfo hdr; ConcatString fcst_var, obs_var; @@ -290,6 +300,16 @@ extern void aggr_ssvar_lines( std::map &, int &, int &); +extern void aggr_seeps_lines( + LineDataFile &, STATAnalysisJob &, + std::map &, + int &, int &); + +extern void aggr_seeps_mpr_lines( + LineDataFile &, STATAnalysisJob &, + std::map &, + int &, int &); + extern void aggr_time_series_lines( LineDataFile &, STATAnalysisJob &, std::map &, diff --git a/src/tools/core/stat_analysis/parse_stat_line.cc b/src/tools/core/stat_analysis/parse_stat_line.cc index a5410e2a53..8fe2b63621 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.cc +++ b/src/tools/core/stat_analysis/parse_stat_line.cc @@ -29,6 +29,8 @@ // 009 10/09/17 Halley Gotway Add GRAD line type. // 010 04/25/18 Halley Gotway Add ECNT line type. // 011 01/24/20 Halley Gotway Add RPS line type. +// 012 11/10/22 Halley Gotway MET #2339 Add SEEPS and SEEPS_MPR +// line types. // //////////////////////////////////////////////////////////////////////// @@ -549,3 +551,57 @@ void parse_ssvar_line(STATLine &l, SSVARInfo &ssvar_info) { } //////////////////////////////////////////////////////////////////////// + +void parse_seeps_line(STATLine &l, SeepsAggScore &agg_score) { + + agg_score.clear(); + + agg_score.n_obs = atoi(l.get_item("TOTAL")); + + agg_score.s12 = atoi(l.get_item("S12")); + agg_score.s13 = atoi(l.get_item("S13")); + agg_score.s21 = atof(l.get_item("S21")); + agg_score.s23 = atof(l.get_item("S23")); + agg_score.s31 = atof(l.get_item("S31")); + agg_score.s32 = atof(l.get_item("S32")); + + agg_score.pf1 = atof(l.get_item("PF1")); + agg_score.pf2 = atof(l.get_item("PF2")); + agg_score.pf3 = atof(l.get_item("PF3")); + + agg_score.pv1 = atof(l.get_item("PV1")); + agg_score.pv2 = atof(l.get_item("PV2")); + agg_score.pv3 = atof(l.get_item("PV3")); + + agg_score.mean_fcst = atof(l.get_item("MEAN_FCST")); + agg_score.mean_obs = atof(l.get_item("MEAN_OBS")); + + agg_score.score = atof(l.get_item("SEEPS")); + agg_score.weighted_score = agg_score.score; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void parse_seeps_mpr_line(STATLine &l, SEEPSMPRData &s_data) { + + s_data.obs_sid = l.get_item("OBS_SID"); + s_data.obs_lat = atof(l.get_item("OBS_LAT")); + s_data.obs_lon = atof(l.get_item("OBS_LON")); + s_data.fcst = atof(l.get_item("FCST")); + s_data.obs = atof(l.get_item("OBS")); + s_data.obs_qc = l.get_item("OBS_QC"); + + s_data.seeps_mpr.fcst_cat = atoi(l.get_item("FCST_CAT")); + s_data.seeps_mpr.obs_cat = atoi(l.get_item("OBS_CAT")); + s_data.seeps_mpr.p1 = atof(l.get_item("P1")); + s_data.seeps_mpr.p2 = atof(l.get_item("P2")); + s_data.seeps_mpr.t1 = atof(l.get_item("T1")); + s_data.seeps_mpr.t2 = atof(l.get_item("T2")); + s_data.seeps_mpr.score = atof(l.get_item("SEEPS")); + + return; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/stat_analysis/parse_stat_line.h b/src/tools/core/stat_analysis/parse_stat_line.h index 62718981b2..ba08d05216 100644 --- a/src/tools/core/stat_analysis/parse_stat_line.h +++ b/src/tools/core/stat_analysis/parse_stat_line.h @@ -26,6 +26,8 @@ // 009 04/25/18 Halley Gotway Add ECNT line type. // 010 01/24/20 Halley Gotway Add RPS line type. // 011 09/28/22 Prestopnik MET #2227 Remove namespace std +// 012 11/10/22 Halley Gotway MET #2339 Add SEEPS and SEEPS_MPR +// line types. // //////////////////////////////////////////////////////////////////////// @@ -45,6 +47,7 @@ #include "vx_analysis_util.h" #include "vx_util.h" #include "vx_statistics.h" +#include "seeps.h" //////////////////////////////////////////////////////////////////////// @@ -102,6 +105,13 @@ struct ORANKData { NumArray ens_na; }; +// SEEPS Matched Pair (SEEPS_MPR) data structure +struct SEEPSMPRData { + ConcatString obs_sid, obs_qc; + double obs_lat, obs_lon, fcst, obs; + SeepsScore seeps_mpr; +}; + //////////////////////////////////////////////////////////////////////// extern void parse_fho_ctable (STATLine &, TTContingencyTable &); @@ -128,6 +138,9 @@ extern void parse_relp_line (STATLine &, RELPData &); extern void parse_orank_line (STATLine &, ORANKData &); extern void parse_ssvar_line (STATLine &, SSVARInfo &); +extern void parse_seeps_line (STATLine &, SeepsAggScore &); +extern void parse_seeps_mpr_line(STATLine &, SEEPSMPRData &); + //////////////////////////////////////////////////////////////////////// #endif // __PARSE_STAT_LINE_H__ diff --git a/src/tools/core/stat_analysis/stat_analysis_job.cc b/src/tools/core/stat_analysis/stat_analysis_job.cc index bceb2af367..c87c87ffaa 100644 --- a/src/tools/core/stat_analysis/stat_analysis_job.cc +++ b/src/tools/core/stat_analysis/stat_analysis_job.cc @@ -43,6 +43,8 @@ // 021 04/12/21 Halley Gotway MET #1735 Support multiple // -out_thresh and -out_line_type options. // 022 05/28/21 Halley Gotway Add MCTS HSS_EC output. +// 023 11/10/22 Halley Gotway MET #2339 Add SEEPS and SEEPS_MPR +// line types. // //////////////////////////////////////////////////////////////////////// @@ -421,6 +423,7 @@ void do_job_aggr(const ConcatString &jobstring, LineDataFile &f, map ens_map; map rps_map; map ssvar_map; + map seeps_map; // // Check that the -line_type option has been supplied exactly once @@ -449,13 +452,14 @@ void do_job_aggr(const ConcatString &jobstring, LineDataFile &f, lt != stat_grad && lt != stat_ecnt && lt != stat_rps && lt != stat_rhist && lt != stat_phist && lt != stat_relp && - lt != stat_ssvar && lt != stat_isc) { + lt != stat_ssvar && lt != stat_isc && + lt != stat_seeps) { mlog << Error << "\ndo_job_aggr() -> " << "the \"-line_type\" option must be set to one of:\n" << "\tFHO, CTC, MCTC,\n" << "\tSL1L2, SAL1L2, VL1L2, VAL1L2,\n" << "\tPCT, NBRCTC, NBRCNT, GRAD, ISC,\n" - << "\tECNT, RPS, RHIST, PHIST, RELP, SSVAR\n\n"; + << "\tECNT, RPS, RHIST, PHIST, RELP, SSVAR, SEEPS\n\n"; throw(1); } @@ -571,6 +575,14 @@ void do_job_aggr(const ConcatString &jobstring, LineDataFile &f, write_job_aggr_ssvar(job, lt, ssvar_map, out_at); } + // + // Sum the SEEPS line types + // + else if(lt == stat_seeps) { + aggr_seeps_lines(f, job, seeps_map, n_in, n_out); + write_job_aggr_seeps(job, lt, seeps_map, out_at); + } + // // Check for no matching STAT lines // @@ -609,13 +621,14 @@ void do_job_aggr_stat(const ConcatString &jobstring, LineDataFile &f, AsciiTable out_at; int i, n; - map ctc_map; - map mctc_map; - map pct_map; - map psum_map; - map wind_map; - map orank_map; - map mpr_map; + map ctc_map; + map mctc_map; + map pct_map; + map psum_map; + map wind_map; + map orank_map; + map mpr_map; + map seeps_mpr_map; // // Check that the -line_type and -out_line_type options have been @@ -658,6 +671,7 @@ void do_job_aggr_stat(const ConcatString &jobstring, LineDataFile &f, // PCT, PSTD, PJC, PRC, ECLV, // WDIR (wind direction) // -line_type ORANK, -out_line_type RHIST, PHIST, RELP, SSVAR + // -line_type SEEPS_MPR, -out_line_type SEEPS // // @@ -810,6 +824,21 @@ void do_job_aggr_stat(const ConcatString &jobstring, LineDataFile &f, } } + // + // Sum the SEEPS_MPR lines: + // SEEPS_MPR -> SEEPS + // + else if(in_lt == stat_seeps_mpr && + has_line_type(out_lt, stat_seeps)) { + + aggr_seeps_mpr_lines(f, job, seeps_mpr_map, n_in, n_out); + for(it=out_lt.begin(); it!=out_lt.end(); it++) { + write_job_aggr_seeps_mpr(job, *it, seeps_mpr_map, out_at); + if(!job.stat_out) write_table(out_at, sa_out); + } + } + + // // Read the matched pair lines: // MPR -> FHO, CTC, CTS, MCTC, MCTS, CNT, @@ -2677,6 +2706,140 @@ void write_job_aggr_ssvar(STATAnalysisJob &job, STATLineType lt, //////////////////////////////////////////////////////////////////////// +void write_job_aggr_seeps(STATAnalysisJob &job, STATLineType lt, + map &m, + AsciiTable &at) { + map::iterator it; + int n, n_row, n_col, r, c; + StatHdrColumns shc; + + // + // Setup the output table + // + n_row = 1 + m.size(); + n_col = 1 + job.by_column.n() + n_seeps_columns; + write_job_aggr_hdr(job, n_row, n_col, at); + + // + // Write the rest of the header row + // + c = 1 + job.by_column.n(); + write_header_row(seeps_columns, n_seeps_columns, 0, at, 0, c); + + // + // Setup the output STAT file + // + job.setup_stat_file(n_row, n); + + mlog << Debug(2) << "Computing output for " + << (int) m.size() << " case(s).\n"; + + // + // Loop through the map + // + for(it = m.begin(), r=1; it != m.end(); it++) { + + // + // Write the output STAT header columns + // + shc = it->second.hdr.get_shc(it->first, job.by_column, + job.hdr_name, job.hdr_value, lt); + + // + // Initialize + // + c = 0; + + // + // SEEPS output line + // + if(job.stat_out) { + write_header_cols(shc, job.stat_at, job.stat_row); + write_seeps_cols(&it->second.agg_score, job.stat_at, + job.stat_row++, n_header_columns); + } + else { + at.set_entry(r, c++, "SEEPS:"); + write_case_cols(it->first, at, r, c); + write_seeps_cols(&it->second.agg_score, at, r++, c); + } + } // end for it + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void write_job_aggr_seeps_mpr(STATAnalysisJob &job, STATLineType lt, + map &m, + AsciiTable &at) { + map::iterator it; + int n, n_row, n_col, r, c; + StatHdrColumns shc; + SeepsAggScore agg_score; + + // + // Setup the output table + // + n_row = 1 + m.size(); + n_col = 1 + job.by_column.n() + n_seeps_columns; + write_job_aggr_hdr(job, n_row, n_col, at); + + // + // Write the rest of the header row + // + c = 1 + job.by_column.n(); + write_header_row(seeps_columns, n_seeps_columns, 0, at, 0, c); + + // + // Setup the output STAT file + // + job.setup_stat_file(n_row, n); + + mlog << Debug(2) << "Computing output for " + << (int) m.size() << " case(s).\n"; + + // + // Loop through the map + // + for(it = m.begin(), r=1; it != m.end(); it++) { + + // + // Write the output STAT header columns + // + shc = it->second.hdr.get_shc(it->first, job.by_column, + job.hdr_name, job.hdr_value, lt); + + // + // Initialize + // + c = 0; + + // + // Compute the aggregated SEEPS score + // + compute_aggregated_seeps(&it->second.pd, &agg_score); + + // + // SEEPS output line + // + if(job.stat_out) { + write_header_cols(shc, job.stat_at, job.stat_row); + write_seeps_cols(&agg_score, job.stat_at, + job.stat_row++, n_header_columns); + } + else { + at.set_entry(r, c++, "SEEPS:"); + write_case_cols(it->first, at, r, c); + write_seeps_cols(&agg_score, at, r++, c); + } + } // end for it + + return; +} + +//////////////////////////////////////////////////////////////////////// + void write_job_aggr_orank(STATAnalysisJob &job, STATLineType lt, map &m, AsciiTable &at, gsl_rng *rng_ptr) { diff --git a/src/tools/core/stat_analysis/stat_analysis_job.h b/src/tools/core/stat_analysis/stat_analysis_job.h index a963fd1944..6340a0549a 100644 --- a/src/tools/core/stat_analysis/stat_analysis_job.h +++ b/src/tools/core/stat_analysis/stat_analysis_job.h @@ -118,6 +118,12 @@ extern void write_job_aggr_relp(STATAnalysisJob &, STATLineType, extern void write_job_aggr_ssvar(STATAnalysisJob &, STATLineType, std::map &, AsciiTable &); +extern void write_job_aggr_seeps(STATAnalysisJob &, STATLineType, + std::map &, AsciiTable &); + +extern void write_job_aggr_seeps_mpr(STATAnalysisJob &, STATLineType, + std::map &, AsciiTable &); + extern void write_job_aggr_orank(STATAnalysisJob &, STATLineType, std::map &, AsciiTable &, gsl_rng *); diff --git a/src/tools/tc_utils/tc_pairs/tc_pairs.cc b/src/tools/tc_utils/tc_pairs/tc_pairs.cc index 1861e16eb0..734f83d256 100644 --- a/src/tools/tc_utils/tc_pairs/tc_pairs.cc +++ b/src/tools/tc_utils/tc_pairs/tc_pairs.cc @@ -247,9 +247,8 @@ void process_command_line(int argc, char **argv) { for(i=0; i n_diag_map; int i, j, n; - map n_diag_map; - // Process the diagnostic inputs - for(i=0; i " + << "the diagnostic path and source arrays must have equal length!\n\n"; + exit(1); + } - // Process the current source - cur_path.clear(); - cur_path.add(diag_path[i]); - cur_model_name.clear(); - cur_model_name.add(diag_model_name[i]); + // Process the diagnostic inputs + for(i=0; i " + << "no \"" << conf_key_diag_info_map + << "\" entry found for diagnostic source \"" + << diag_source[i] << "\"!\n\n"; + exit(1); + } + info_ptr = &conf_info.DiagInfoMap.at(diag_source[i]); - // Pointer to the conversion map for this source + // Check for diagnostic conversion map entry map * convert_ptr = 0; - if(conf_info.DiagConvertMap.count(diag_source[i]) > 0) { - convert_ptr = &conf_info.DiagConvertMap.at(diag_source[i]); + map< ConcatString,map >::iterator it; + for(it = conf_info.DiagConvertMap.begin(); + it != conf_info.DiagConvertMap.end(); it++) { + + // Partial string matching for diag_source + if(check_reg_exp(it->first.c_str(), diag_source[i].c_str())) { + convert_ptr = &it->second; + break; + } } // Process the diagnostic files - diag_file.read(diag_source[i], files[j], model_names, convert_ptr); + diag_file.read(cur_files[j], diag_source[i], + info_ptr->TrackSource, info_ptr->FieldSource, + info_ptr->MatchToTrack, convert_ptr); // Add diagnostics data to existing tracks - n += tracks.add_diag_data(diag_file, conf_info.DiagName); + n += tracks.add_diag_data(diag_file, info_ptr->DiagName); } // end for j @@ -1058,11 +1072,11 @@ void process_diags(TrackInfoArray &tracks) { } // end for i // Print the diagnostic track counts - for(map::iterator it = n_diag_map.begin(); + for(map::iterator it = n_diag_map.begin(); it != n_diag_map.end(); it++) { mlog << Debug(3) - << "Added " << diagtype_to_string(it->first) - << " diagnostics to " << it->second << " tracks.\n"; + << "Added " << it->first << " diagnostics to " + << it->second << " tracks.\n"; } return; @@ -2233,8 +2247,8 @@ void usage() { << "\t\t\"-diag source path\" is used one or more times to " << "specify a file or top-level directory containing tropical " << "cyclone diagnostics \"" << atcf_suffix - << "\" data to process. The supported formats are TCDIAG, " - << "LSDIAG_RT, LSDIAG_DEV (optional).\n" + << "\" data to process. The supported sources are CIRA_DIAG_RT " + << "and SHIPS_DIAG_RT (optional).\n" << "\t\t\"-out base\" overrides the default output file base " << "(" << out_base << ") (optional).\n" @@ -2247,11 +2261,7 @@ void usage() { << "\tNote: The \"-adeck\", \"-edeck\", and \"-bdeck\" options " << "may include \"suffix=string\" to modify the model names " - << "from that path.\n\n" - - << "\tNote: The \"-diag\" option may include \"model=string\" to " - << "override the model name of the tracks to which those diagnostics " - << "correspond.\n\n"; + << "from that path.\n\n"; exit(1); } @@ -2313,10 +2323,6 @@ void set_atcf_path(const StringArray &a, //////////////////////////////////////////////////////////////////////// void set_diag(const StringArray &a) { - int i; - StringArray sa; - DiagType source; - ConcatString cs, model_name; // Should have length of at least 2 if(a.n() < 2) { @@ -2326,33 +2332,10 @@ void set_diag(const StringArray &a) { exit(1); } - // Parse the DiagType source - source = string_to_diagtype(a[0].c_str()); - - // Check for the model sub-argument - for(i=1; i " - << "the model name must be specified as " - << "\"model=string\".\n\n"; - usage(); - } - else { - model_name = sa[1]; - } - } - } - // Parse the remaining paths - for(i=1; i diag_source; -static StringArray diag_path, diag_model_name; +static StringArray diag_source, diag_path; static ConcatString config_file; static TCPairsConfInfo conf_info; diff --git a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc index fa6c9b318c..b4faab8aef 100644 --- a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc +++ b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.cc @@ -25,8 +25,11 @@ using namespace std; //////////////////////////////////////////////////////////////////////// +void parse_conf_diag_info_map(Dictionary *, + map &); + void parse_conf_diag_convert_map(Dictionary *, - map > &); + map< ConcatString,map > &); //////////////////////////////////////////////////////////////////////// // @@ -101,7 +104,7 @@ void TCPairsConfInfo::clear() { DLandFile.clear(); WatchWarnFile.clear(); WatchWarnOffset = bad_data_int; - DiagName.clear(); + DiagInfoMap.clear(); DiagConvertMap.clear(); BasinMap.clear(); Version.clear(); @@ -312,8 +315,8 @@ void TCPairsConfInfo::process_config() { // Conf: WatchWarnOffset WatchWarnOffset = dict->lookup_int(conf_key_time_offset); - // Conf: DiagName - DiagName = dict->lookup_string_array(conf_key_diag_name); + // Conf: DiagInfoMap + parse_conf_diag_info_map(dict, DiagInfoMap); // Conf: DiagConvertMap parse_conf_diag_convert_map(dict, DiagConvertMap); @@ -330,14 +333,74 @@ void TCPairsConfInfo::process_config() { // //////////////////////////////////////////////////////////////////////// +void DiagInfo::clear() { + + TrackSource.clear(); + FieldSource.clear(); + MatchToTrack.clear(); + DiagName.clear(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void parse_conf_diag_info_map(Dictionary *dict, map &source_map) { + int i; + Dictionary *map_dict = (Dictionary *) 0; + ConcatString diag_source; + DiagInfo cur_info; + + const char *method_name = "parse_conf_diag_info_map() -> "; + + if(!dict) { + mlog << Error << "\n" << method_name + << "empty dictionary!\n\n"; + exit(1); + } + + // Conf: diag_info_map + map_dict = dict->lookup_array(conf_key_diag_info_map); + + // Loop through the array entries + for(i=0; in_entries(); i++) { + + // Initialize the current map + cur_info.clear(); + + // Lookup the metadata entries + diag_source = (*map_dict)[i]->dict_value()->lookup_string(conf_key_diag_source); + cur_info.TrackSource = (*map_dict)[i]->dict_value()->lookup_string(conf_key_track_source); + cur_info.FieldSource = (*map_dict)[i]->dict_value()->lookup_string(conf_key_field_source); + cur_info.MatchToTrack = (*map_dict)[i]->dict_value()->lookup_string_array(conf_key_match_to_track); + cur_info.DiagName = (*map_dict)[i]->dict_value()->lookup_string_array(conf_key_diag_name); + + // diag_source entries must be unique + if(source_map.count(diag_source) > 0) { + mlog << Error << "\n" << method_name + << "multiple entries found for diag_source \"" << diag_source + << "\" in \"" << conf_key_diag_info_map << "\"!\n\n"; + exit(1); + } + // Add a new source entry + else { + source_map[diag_source] = cur_info; + } + + } // end for i + + return; +} + +//////////////////////////////////////////////////////////////////////// + void parse_conf_diag_convert_map(Dictionary *dict, - map > &source_map) { + map< ConcatString,map > &source_map) { int i, j; Dictionary *map_dict = (Dictionary *) 0; map cur_map; - DiagType source; + ConcatString diag_source, key; StringArray sa; - ConcatString key; UserFunc_1Arg fx; const char *method_name = "parse_conf_diag_convert_map() -> "; @@ -358,9 +421,8 @@ void parse_conf_diag_convert_map(Dictionary *dict, cur_map.clear(); // Lookup the source, key, and convert function - source = string_to_diagtype( - (*map_dict)[i]->dict_value()->lookup_string(conf_key_source).c_str()); - sa = (*map_dict)[i]->dict_value()->lookup_string_array(conf_key_key); + diag_source = (*map_dict)[i]->dict_value()->lookup_string(conf_key_diag_source); + sa = (*map_dict)[i]->dict_value()->lookup_string_array(conf_key_key); fx.clear(); fx.set((*map_dict)[i]->dict_value()->lookup(conf_key_convert)); @@ -378,15 +440,15 @@ void parse_conf_diag_convert_map(Dictionary *dict, } // Append to the existing source entry - if(source_map.count(source) > 0) { + if(source_map.count(diag_source) > 0) { for(map::iterator it = cur_map.begin(); it != cur_map.end(); it++) { - source_map.at(source).insert(pair(it->first, it->second)); + source_map.at(diag_source).insert(pair(it->first, it->second)); } } // Add a new source entry else { - source_map.insert(pair >(source, cur_map)); + source_map.insert(pair< ConcatString,map >(diag_source, cur_map)); } } // end for i diff --git a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h index 12dc7eea14..bd65d31d45 100644 --- a/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h +++ b/src/tools/tc_utils/tc_pairs/tc_pairs_conf_info.h @@ -33,6 +33,16 @@ struct ConsensusInfo { //////////////////////////////////////////////////////////////////////// +struct DiagInfo { + ConcatString TrackSource; + ConcatString FieldSource; + StringArray MatchToTrack; + StringArray DiagName; + void clear(); +}; + +//////////////////////////////////////////////////////////////////////// + class TCPairsConfInfo { private: @@ -111,11 +121,11 @@ class TCPairsConfInfo { // Watch/warnings time offset int WatchWarnOffset; - // Diagnostics to be extracted - StringArray DiagName; + // Diagnostics metadata + std::map DiagInfoMap; // Diagnostic conversions - std::map< DiagType, std::map > DiagConvertMap; + std::map< ConcatString,std::map > DiagConvertMap; // Basin Map std::map BasinMap;