diff --git a/met/docs/Users_Guide/tc-gen.rst b/met/docs/Users_Guide/tc-gen.rst index b67f1a26b9..74271ec9e1 100644 --- a/met/docs/Users_Guide/tc-gen.rst +++ b/met/docs/Users_Guide/tc-gen.rst @@ -6,19 +6,25 @@ TC-Gen Tool Introduction ____________ -The TC-Gen tool provides verification of deterministic and probabilistic tropical cyclone genesis forecasts in the ATCF file format. Producing reliable tropical cyclone genesis forecasts is an important metric for global numerical weather prediction models. This tool ingests deterministic model output post-processed by a genesis tracking software (e.g. GFDL vortex tracker), ATCF edeck files containing probability of genesis forecasts, and ATCF reference track dataset(s) (e.g. Best Track analysis and CARQ operational tracks). It writes categorical counts and statistics. The capability to modify the spatial and temporal tolerances when matching forecasts to reference genesis events, as well as scoring those matched pairs, gives users the ability to condition the criteria based on model performance and/or conduct sensitivity analyses. Statistical aspects are outlined in :numref:`tc-gen_stat_aspects` and practical aspects of the TC-Gen tool are described in :numref:`tc-gen_practical_info`. +The TC-Gen tool provides verification of deterministic and probabilistic tropical cyclone genesis forecasts in the ATCF file and shapefile formats. Producing reliable tropical cyclone genesis forecasts is an important metric for global numerical weather prediction models. This tool ingests deterministic model output post-processed by genesis tracking software (e.g. GFDL vortex tracker), ATCF edeck files containing probability of genesis forecasts, operational shapefile warning areas, and ATCF reference track dataset(s) (e.g. Best Track analysis and CARQ operational tracks). It writes categorical counts and statistics. The capability to modify the spatial and temporal tolerances when matching forecasts to reference genesis events, as well as scoring those matched pairs, gives users the ability to condition the criteria based on model performance and/or conduct sensitivity analyses. Statistical aspects are outlined in :numref:`tc-gen_stat_aspects` and practical aspects of the TC-Gen tool are described in :numref:`tc-gen_practical_info`. .. _tc-gen_stat_aspects: Statistical aspects ___________________ -The TC-Gen tool processes both deterministic and probabilistic forecasts. For deterministic forecasts specified using the **-track** command line option, it identifies genesis events in both the forecasts and reference datasets, typically Best tracks. It applies user-specified configuration options to pair up the forecast and reference genesis events and categorize each pair as a hit, miss, or false alarm. +The TC-Gen tool processes both deterministic and probabilistic forecasts. + +For deterministic forecasts specified using the **-track** command line option, it identifies genesis events in both the forecasts and reference datasets, typically Best tracks. It applies user-specified configuration options to pair up the forecast and reference genesis events and categorize each pair as a hit, miss, or false alarm. As with other extreme events (where the event occurs much less frequently than the non-event), the correct negative category is not computed since the non-events would dominate the contingency table. Therefore, only statistics that do not include correct negatives should be considered for this tool. The following CTS statistics are relevant: Base rate (BASER), Mean forecast (FMEAN), Frequency Bias (FBIAS), Probability of Detection (PODY), False Alarm Ratio (FAR), Critical Success Index (CSI), Gilbert Skill Score (GSS), Extreme Dependency Score (EDS), Symmetric Extreme Dependency Score (SEDS), Bias Adjusted Gilbert Skill Score (BAGSS). For probabilistic forecasts specified using the **-edeck** command line option, it identifies genesis events in the reference dataset. It applies user-specified configuration options to pair the forecast probabilities to the reference genesis events. These pairs are added to an Nx2 probabilistic contingency table. If the reference genesis event occurs within in the predicted time window, the pair is counted in the observation-yes column. Otherwise, it is added to the observation-no column. +For warning area shapefiles specified using the **-shape** command line option, it processes metadata from the corresponding database files. The database file is assumed to exist at exactly the same path as the shapefile, but with a ".dbf" suffix instead of ".shp". Note that only shapefiles exactly following the NOAA National Hurricane Center's (NHC) "gtwo_areas_YYYYMMDDHHMM.shp" file naming and corresonding metadata conventions are supported. For each shapefile record, the database file defines up to three corresponding probability values. The first percentage is interpreted as the probability of genesis inside the shape within 48 hours. The second and, if provided, third percentages are interpreted as the 120-hour and 168-hour probabilities, respectively. Care is taken to identify and either ignore or update duplicate shapes found in the input. + +The shapes are then subset based on the filtering criteria in the configuration file. For each probability and shape, the reference genesis events are searched for a match within the defined time window. These pairs are added to an Nx2 probabilistic contingency table. The probabilistic contingeny tables and statistics are computed and reported separately for filter defined and lead hour encountered in the input. + Other considerations for interpreting the output of the TC-Gen tool involve the size of the contingency table output. The size of the contingency table will change depending on the number of matches. Additionally, the number of misses is based on the forecast duration and interval (specified in the configuration file). This change is due to the number of model opportunities to forecast the event, which is determined by the specified duration/interval. Care should be taken when interpreting the statistics for filtered data. In some cases, variables (e.g. storm name) are only available in either the forecast or reference datasets, rather than both. When filtering on a field that is only present in one dataset, the contingency table counts will be impacted. Similarly, the initialization field only impacts the model forecast data. If the valid time (which will impact the reference dataset) isn't also specified, the forecasts will be filtered and matched such that the number of misses will erroneously increase. See :numref:`tc-gen_practical_info` for more detail. @@ -28,7 +34,7 @@ Care should be taken when interpreting the statistics for filtered data. In some Practical information _____________________ -This section describes how to configure and run the TC-Gen tool. The TC-Gen tool identifies tropical cyclone genesis events in both genesis forecasts and ATCF track datasets. It applies configurable logic to process the forecast and observed genesis events, classify them, and populate a contingency table with hits, misses, and false alarms. It writes the categorical counts and statistics to the output file(s). The tool can be configured to apply one or more sets of filtering criteria in a single run. The following sections describe the usage statement, required arguments, and optional arguments for tc_gen. +This section describes how to configure and run the TC-Gen tool. The following sections describe the usage statement, required arguments, and optional arguments for tc_gen. tc_gen usage ~~~~~~~~~~~~ @@ -38,7 +44,9 @@ The usage statement for tc_gen is shown below: .. code-block:: none Usage: tc_gen - -genesis source and/or -edeck source + -genesis source + -edeck source + -shape source -track source -config file [-out base] @@ -52,20 +60,24 @@ Required arguments for tc_gen 1. The **-genesis source** argument is the path to one or more ATCF or fort.66 (see documentation listed below) files generated by the Geophysical Fluid Dynamics Laboratory (GFDL) Vortex Tracker when run in tcgen mode or an ASCII file list or a top-level directory containing them. The required file format is described in the "Output formats" section of the `GFDL Vortex Tracker users guide. `_ -2. The **-edeck source** argument is the path to one or more ATCF edeck files, an ASCII file list containing them, or a top-level directory with files matching the regular expression ".dat". The probability of genesis are read from each edeck input file and verified against at the **-track** data. The **-genesis** or **-edeck** option must be used at least once. +2. The **-edeck source** argument is the path to one or more ATCF edeck files, an ASCII file list containing them, or a top-level directory with files matching the regular expression ".dat". The probability of genesis are read from each edeck input file and verified against at the **-track** data. -3. The **-track source** argument is one or more ATCF reference track files or an ASCII file list or top-level directory containing them, with files ending in ".dat". This tool processes either Best track data from bdeck files, or operational track data (e.g. CARQ) from adeck files, or both. Providing both bdeck and adeck files will result in a richer dataset to match with the **-genesis** files. Both adeck and bdeck data should be provided using the **-track** option. The **-track** option must be used at least once. +3. The **-shape source** argument is the path to one or more NHC genesis warning area shapefiles, an ASCII file list containing them, or a top-level directory with files matching the regular expression "gtwo_areas.*.shp". The genesis warning areas and corresponding 2, 5, and 7 day probability values area verified against the **-track** data. -4. The **-config** file argument indicates the name of the configuration file to be used. The contents of the configuration file are discussed below. +Note: The **-genesis**, **-edeck**, or **-shape** options must be used at least once. + +4. The **-track source** argument is one or more ATCF reference track files or an ASCII file list or top-level directory containing them, with files ending in ".dat". This tool processes either Best track data from bdeck files, or operational track data (e.g. CARQ) from adeck files, or both. Providing both bdeck and adeck files will result in a richer dataset to match with the **-genesis** files. Both adeck and bdeck data should be provided using the **-track** option. The **-track** option must be used at least once. + +5. The **-config** file argument indicates the name of the configuration file to be used. The contents of the configuration file are discussed below. Optional arguments for tc_gen ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -5. The **-out base** argument indicates the path of the output file base. This argument overrides the default output file base (./tc_gen) +6. The **-out base** argument indicates the path of the output file base. This argument overrides the default output file base (./tc_gen) -6. The **-log file** option directs output and errors to the specified log file. All messages will be written to that file as well as standard out and error. Thus, users can save the messages without having to redirect the output on the command line. The default behavior is no log file. +7. The **-log file** option directs output and errors to the specified log file. All messages will be written to that file as well as standard out and error. Thus, users can save the messages without having to redirect the output on the command line. The default behavior is no log file. -7. The **-v level** option indicates the desired level of verbosity. The contents of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity above 1 will increase the amount of logging. +8. The **-v level** option indicates the desired level of verbosity. The contents of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity above 1 will increase the amount of logging. Scoring Logic ^^^^^^^^^^^^^ @@ -120,6 +132,29 @@ The TC-Gen tool implements the following logic: * Report the Nx2 probabilistic contingency table counts and statistics for each forecast model, lead time, and configuration file filter. These counts and statistics are identified in the output files as *PROB_GENESIS*. +* For **-shape** inputs: + + * For each input shapefile, parse the timestamp from the "gtwo_areas_YYYYMMDDHHMM.shp" naming convention, and error out otherwise. Round the timestamp to the nearest synoptic time (e.g. 00, 06, 12, 18) and store that as the issuance time. + + * Open the shapefile and corresponding database file. Process each record. + + * For each record, extract the shape and metadata which defines the basin and 2, 5, and 7 day probabilities. + + * Check if this shape is a duplicate that has already been processed. If it is an exact duplicate, with the same basin, file timestamp, issue time, and min/max lat/lon values, ignore it. If the file timestamp is older than the existing shape, also ignore it. If the file timestamp is newer than the existing shape, replace the existing shape with the new one. + + * Loop over the filters defined in the configuration file and apply the following logic for each. + + * Subset the list of genesis shapes based on the current filter criteria. + + * Search the Best track genesis events to see if any occurred inside the shape within 7 days of the issuance time. If multiple genesis events occurred, choose the one closest to the issuance time. + + * If not found, score each probability as a miss. + + * If found, further check the 2 and 5 day time windows to classify each probability as a hit or miss. + + * Add each probability pair to an Nx2 probabilistic contingency table, tracking results separately for each lead time. + + * Report the Nx2 probabilistic contingency table counts and statistics for each lead time. These counts and statistics are identified in the output files as *GENESIS_SHAPE*. tc_gen configuration file ~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -239,6 +274,8 @@ ______________________ The **init_beg**, **init_end**, **init_inc**, and **init_exc** entries define strings in YYYYMMDD[_HH[MMSS]] format which defines which forecast and operational tracks initializations to be processed. If left empty, all tracks will be used. Otherwise, only those tracks whose initialization time meets all the criteria will be processed. The initialization time must fall between **init_beg**, and **init_end**, must appear in **init_inc** inclusion list, and must not appear in the **init_exc** exclusion list. Note that these settings only apply to the forecast and operational tracks, not the Best tracks, for which the initialization time is undefined. Care should be given when interpreting the contingency table results for filtered data. +For genesis shapes, these options are used to filter the warning issuance time. + ______________________ .. code-block:: none @@ -257,6 +294,8 @@ ______________________ The **init_hour** and **lead** entries are arrays of strings in HH[MMSS] format defining which forecast tracks should be included. If left empty, all tracks will be used. Otherwise, only those forecast tracks whose initialization hour and lead times appear in the list will be used. Note that these settings only apply to the forecast tracks, not the Best tracks, for which the initialization time is undefined. Care should be given when interpreting the contingency table results for filtered data. +For genesis shapes, the **init_hour** option is used to filter the warning issuance hour. + ______________________ .. code-block:: none @@ -265,6 +304,8 @@ ______________________ The **vx_mask** entry is a string defining the path to a Lat/Lon polyline file or a gridded data file that MET can read to subset the results spatially. If specified, only those genesis events whose Lat/Lon location falls within the specified area will be included. +If specified for genesis shapes, the lat/lon of the central location of the shape will be checked. The central location is computed as the average of the min/max lat/lon values of the shape points. + ______________________ .. code-block:: none @@ -275,6 +316,8 @@ The **basin_mask** entry is an array of strings listing tropical cycline basin a The **vx_mask** and **basin_mask** names are concatenated and written to the **VX_MASK** output column. +If **vx_mask** is not specified for genesis shapes and **basin_mask** is, the basin name is extracted from the shapefile metadata and compared to the **basin_mask** list. + ______________________ .. code-block:: none @@ -415,7 +458,7 @@ ______________________ prob_genesis_thresh = ==0.25; -The **prob_genesis_thresh** entry defines the probability thresholds used to create the output Nx2 contingency table when verifying edeck probability of genesis forecasts. The default is probability bins of width 0.25. These probabilities may be specified as a list (>0.00,>0.25,>0.50,>0.75,>1.00) or using shorthand notation (==0.25) for bins of equal width. +The **prob_genesis_thresh** entry defines the probability thresholds used to create the output Nx2 contingency table when verifying edeck probability of genesis forecasts and probabilistic shapefile warning areas. The default is probability bins of width 0.25. These probabilities may be specified as a list (>0.00,>0.25,>0.50,>0.75,>1.00) or using shorthand notation (==0.25) for bins of equal width. ______________________ @@ -435,12 +478,12 @@ ______________________ dland_file = "MET_BASE/tc_data/dland_global_tenth_degree.nc"; version = "VN.N"; -The configuration options listed above are common to many MET tools and are described in :numref:`config_options`. TC-Gen writes output for 2x2 contingency tables to the **FHO**, **CTC**, and **CTS** line types when verifying deterministic genesis forecasts specified using the **-track** command line option. TC-Gen writes output for Nx2 probabilistic contingency tables to the **PCT**, **PSTD**, **PJC**, and **PRC** line types when verifying the probability of genesis forecasts specified using the **-edeck** command line option. Note that the **genmpr** line type is specific to TC-Gen and describes individual genesis matched pairs. +The configuration options listed above are common to many MET tools and are described in :numref:`config_options`. TC-Gen writes output for 2x2 contingency tables to the **FHO**, **CTC**, and **CTS** line types when verifying deterministic genesis forecasts specified using the **-track** command line option. TC-Gen writes output for Nx2 probabilistic contingency tables to the **PCT**, **PSTD**, **PJC**, and **PRC** line types when verifying the probability of genesis forecasts specified using the **-edeck** command line option and probabilistic shapefiles using the **-shape** command line option. Note that the **genmpr** line type is specific to TC-Gen and describes individual genesis matched pairs. tc_gen output ~~~~~~~~~~~~~ -TC-Gen produces output in STAT and, optionally, ASCII and NetCDF formats. The ASCII output duplicates the STAT output but has the data organized by line type. The output files are created based on the **-out** command line argument. The default output base name, **./tc_gen** writes output files in the current working directory named **tc_gen.stat** and, optionally, **tc_gen_fho.txt, tc_gen_ctc.txt**, **tc_gen_cts.txt**, **tc_gen_genmpr.txt**, and **tc_gen_pairs.nc**. The format of the STAT and ASCII output of the TC-Gen tool matches the output of other MET tools with the exception of the genesis matched pair line type. Please refer to the tables in :numref:`point_stat-output` for a description of the common output line types. The genesis matched pair line type and NetCDF output file are described below. +TC-Gen produces output in STAT and, optionally, ASCII and NetCDF formats. The ASCII output duplicates the STAT output but has the data organized by line type. The output files are created based on the **-out** command line argument. The default output base name, **./tc_gen** writes output files in the current working directory named **tc_gen.stat** and, optionally, **tc_gen_pairs.nc** and **tc_gen_{TYPE}.txt** for each of the supported output line types. These output files can easily be redirected to another location using the **-out** command line option. The format of the STAT and ASCII output of the TC-Gen tool matches the output of other MET tools with the exception of the genesis matched pair line type. Please refer to the tables in :numref:`point_stat-output` for a description of the common output line types. The genesis matched pair line type and NetCDF output file are described below. .. _table_TG_header_info_tg_outputs: @@ -483,7 +526,7 @@ TC-Gen produces output in STAT and, optionally, ASCII and NetCDF formats. The AS - Maximum Best track valid time in YYYYMMDD_HHMMSS format * - 10 - FCST_VAR - - Genesis methodology (GENESIS_DEV, GENESIS_OPS, or PROB_GENESIS) + - Genesis methodology (GENESIS_DEV, GENESIS_OPS, PROB_GENESIS, or GENESIS_SHAPE) * - 11 - FCST_UNITS - Does not apply and is set to NA @@ -492,7 +535,7 @@ TC-Gen produces output in STAT and, optionally, ASCII and NetCDF formats. The AS - Does not apply and is set to NA * - 13 - OBS_VAR - - Genesis methodology (GENESIS_DEV, GENESIS_OPS, or PROB_GENESIS) + - Genesis methodology (GENESIS_DEV, GENESIS_OPS, PROB_GENESIS, or GENESIS_SHAPE) * - 14 - OBS_UNITS - Does not apply and is set to NA diff --git a/met/src/basic/vx_config/Makefile.am b/met/src/basic/vx_config/Makefile.am index bce95fa29c..b4d2d52ce6 100644 --- a/met/src/basic/vx_config/Makefile.am +++ b/met/src/basic/vx_config/Makefile.am @@ -11,7 +11,6 @@ include ${top_srcdir}/Make-include # Yacc/lex flags AM_YFLAGS = --defines=config.tab.h -p config -#AM_LFLAGS = --prefix=config --outfile=lex.yy.c # The library diff --git a/met/src/libcode/vx_color/Makefile.am b/met/src/libcode/vx_color/Makefile.am index 0035825c40..8c147693d7 100644 --- a/met/src/libcode/vx_color/Makefile.am +++ b/met/src/libcode/vx_color/Makefile.am @@ -11,7 +11,6 @@ include ${top_srcdir}/Make-include # YACC/LEX flags AM_YFLAGS = --defines=color_parser_yacc.h -p color -#AM_LFLAGS = --prefix=color --outfile=lex.yy.c # The library diff --git a/met/src/libcode/vx_gis/Makefile.am b/met/src/libcode/vx_gis/Makefile.am index fcc506f49f..f214481071 100644 --- a/met/src/libcode/vx_gis/Makefile.am +++ b/met/src/libcode/vx_gis/Makefile.am @@ -12,24 +12,16 @@ include ${top_srcdir}/Make-include noinst_LIBRARIES = libvx_gis.a -libvx_gis_a_SOURCES = shapetype_to_string.h \ - shapetype_to_string.cc \ - dbf_file.cc \ - dbf_file.h \ - shp_array.h \ - shp_file.cc \ - shp_file.h \ - shp_point_record.cc \ - shp_point_record.h \ - shp_poly_record.cc \ - shp_poly_record.h \ - shp_types.h \ - shx_file.cc \ - shx_file.h - - - - +libvx_gis_a_SOURCES = \ + shapetype_to_string.cc shapetype_to_string.h \ + dbf_file.cc dbf_file.h \ + shx_file.cc shx_file.h \ + shp_file.cc shp_file.h \ + shp_types.h \ + shp_array.h \ + shp_point_record.cc shp_point_record.h \ + shp_poly_record.cc shp_poly_record.h \ + vx_gis.h libvx_gis_a_CPPFLAGS = ${MET_CPPFLAGS} # If we are in development mode, generate the "to_string" files and diff --git a/met/src/libcode/vx_gis/dbf_file.cc b/met/src/libcode/vx_gis/dbf_file.cc index eb0fdc921f..0a5049d150 100644 --- a/met/src/libcode/vx_gis/dbf_file.cc +++ b/met/src/libcode/vx_gis/dbf_file.cc @@ -1,4 +1,3 @@ - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) @@ -9,7 +8,6 @@ //////////////////////////////////////////////////////////////////////// - using namespace std; #include @@ -311,15 +309,12 @@ unsigned char buf[32]; // seek to the first record // -// pos = pos_first_record; - pos = 32; -if ( lseek(fd, pos, SEEK_SET) < 0 ) { +if ( ::lseek(fd, pos, SEEK_SET) < 0 ) { - mlog << Error - << "\n\n DbfHeader::set_subrecords(int) -> lseek error ... " - << strerror(errno) << "\n\n"; + mlog << Error << "\nDbfHeader::set_subrecords(int) -> " + << "lseek error ... " << strerror(errno) << "\n\n"; exit ( 1 ); @@ -339,8 +334,8 @@ for (j=0; j read error ... n_read = " << n_read << "\n\n"; + mlog << Error << "\nDbfHeader::set_subrecords(int) -> " + << "read error ... n_read = " << n_read << "\n\n"; exit ( 1 ); @@ -372,8 +367,8 @@ DbfSubRecord * DbfHeader::lookup_subrec(const char * text) const if ( empty(text) ) { - mlog << Error - << "\n\n DbfHeader::set_subrecords(const char *) -> empty string!\n\n"; + mlog << Error << "\nDbfHeader::set_subrecords(const char *) -> " + << "empty string!\n\n"; exit ( 1 ); @@ -552,7 +547,6 @@ out << prefix << "field_flags = " << ((int) field_flags) << '\n'; out << prefix << "autoinc_next_value = " << autoinc_next_value << '\n'; out << prefix << "autoinc_step_value = " << autoinc_step_value << '\n'; - // // done // @@ -589,9 +583,6 @@ memcpy(&autoinc_next_value, buf + 19, 4); autoinc_step_value = (int) (buf[23]); - - - // // done // @@ -651,7 +642,6 @@ for (j=0; j<(h.n_subrecs); ++j) { } - // for (k=0; k " + << "should never be called!\n\n"; + +exit ( 1 ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +DbfFile & DbfFile::operator=(const DbfFile &) + +{ + +mlog << Error << "\nDbfFile::operator=(const DbfFile &) -> " + << "should never be called!\n\n"; + +exit ( 1 ); + +return ( * this ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +void DbfFile::init_from_scratch() + +{ + +fd = -1; + +close(); + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + +void DbfFile::close() + +{ + +if ( fd >= 0 ) ::close(fd); + +fd = -1; + +memset(&Header, 0, sizeof(Header)); + +Filename.clear(); + +At_Eof = false; + + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + +const char * DbfFile::filename() const + +{ + +if ( Filename.empty() ) return ( (const char *) 0 ); + +return ( Filename.text() ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool DbfFile::open(const char * path) + +{ + +size_t n_read, bytes; +int j; + +close(); + +if ( (fd = met_open(path, O_RDONLY)) < 0 ) { + + fd = -1; + + return ( false ); + +} + +Filename = path; + + // + // read the file header + // + +const size_t buf_size = 65536; +unsigned char buf[buf_size]; + +bytes = 32; + +if((n_read = ::read(fd, buf, bytes)) != bytes) { + mlog << Warning << "\nDbfFile::open(const char * path) -> " + << "trouble reading file header from dbf file \"" + << path << "\"\n\n"; + close(); + + return(false); +} + +Header.set_header(buf); + + // + // subrecords + // + +Header.set_subrecords(fd); + + // + // done + // + +return ( true ); + +} + +//////////////////////////////////////////////////////////////////////// + + +int DbfFile::position() const + +{ + +if ( ! is_open() ) { + + mlog << Error << "\nDbfFile::position() -> " + << "no file open!\n\n"; + + exit ( 1 ); + +} + +const int pos = ::lseek(fd, 0, SEEK_CUR); + +if ( pos < 0 ) { + + mlog << Error << "\nDbfFile::position() const -> " + << "lseek error ... " << strerror(errno) << "\n\n"; + + exit ( 1 ); + +} + + +return ( pos ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +void DbfFile::lseek(int offset, int whence) + +{ + + +if ( ! is_open() ) { + + mlog << Error << "\nDbfFile::lseek() -> " + << "no file open!\n\n"; + + exit ( 1 ); + +} + + +if ( ::lseek(fd, offset, whence) < 0 ) { + + mlog << Error << "\nDbfFile::lseek() -> " + << "lseek(2) failed! ... " << strerror(errno) << "\n\n"; + + exit ( 1 ); + +} + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool DbfFile::read(unsigned char * buf, int n_bytes) + +{ + +if ( ! is_open() ) { + + mlog << Error << "\nDbfFile::read() -> " + << "no file open!\n\n"; + + exit ( 1 ); + +} + +int n_read; + +n_read = ::read(fd, buf, n_bytes); + +if ( n_read < 0 ) { // some kind of error + + mlog << Error << "\nDbfFile::read() -> " + << "read error on dbf file \"" << Filename << "\"\n\n"; + + exit ( 1 ); + +} + +if ( n_read == 0 ) return ( false ); + +if ( n_read == n_bytes ) return ( true ); + +return ( false ); // gotta return something + +} + + +//////////////////////////////////////////////////////////////////////// + + +StringArray DbfFile::subrecord_names() + +{ + +StringArray sa; + + // + // subrecords (if any) + // + +if ( Header.subrec ) { + + for (int j=0; j (Header.n_records) ) { + + mlog << Error << "\nDbfFile::subrecord_values(int i_rec) -> " + << "requested record index (" << i_rec << ") out of range (0," + << Header.n_records << ")!\n\n"; + + exit ( 1 ); +} + +bytes = Header.record_length; + +if ( ::lseek(fd, Header.pos_first_record + i_rec * bytes, SEEK_SET) < 0 ) { + + mlog << Error << "\nDbfFile::subrecord_values(int i_rec) -> " + << "lseek error\n\n"; + + exit ( 1 ); + +} + + // + // read the requested record into the buffer + // + +n_read = ::read(fd, buf, bytes < buf_size ? bytes : buf_size); + +if ( n_read != bytes ) { + + mlog << Error << "\nDbfFile::subrecord_values(int i_rec) -> " + << "read error ... n_read = " << n_read << "\n\n"; + + exit ( 1 ); + +} + +if ( Header.record_length < buf_size) buf[Header.record_length] = 0; + +std::string s = (const char *) buf+1; // skip first byte + + // + // parse each subrecord value + // + +for (j=0,pos=0; j<(Header.n_subrecs); ++j) { + + cs = s.substr(pos, Header.subrec[j].field_length); + cs.ws_strip(); + sa.add(cs); + + pos += Header.subrec[j].field_length; + +} + +return ( sa ); + +} + + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_gis/dbf_file.h b/met/src/libcode/vx_gis/dbf_file.h index 6bc2a354f4..7e60e681e5 100644 --- a/met/src/libcode/vx_gis/dbf_file.h +++ b/met/src/libcode/vx_gis/dbf_file.h @@ -6,17 +6,13 @@ // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - //////////////////////////////////////////////////////////////////////// - #ifndef __VX_SHAPEFILES_DBF_FILE_H__ #define __VX_SHAPEFILES_DBF_FILE_H__ - //////////////////////////////////////////////////////////////////////// - // // Got info on the dbf file format at // @@ -27,16 +23,12 @@ // http://web.archive.org/web/20150323061445/http://ulisse.elettra.trieste.it/services/doc/dbase/DBFstruct.htm // - //////////////////////////////////////////////////////////////////////// - struct DbfSubRecord; // forward reference - //////////////////////////////////////////////////////////////////////// - class DbfHeader { private: @@ -85,10 +77,8 @@ class DbfHeader { }; - //////////////////////////////////////////////////////////////////////// - class DbfSubRecord { private: @@ -133,19 +123,83 @@ class DbfSubRecord { }; +//////////////////////////////////////////////////////////////////////// + +void dump_record(ostream &, const int depth, const unsigned char * buf, const DbfHeader &); //////////////////////////////////////////////////////////////////////// +class DbfFile { -void dump_record(ostream &, const int depth, const unsigned char * buf, const DbfHeader &); + private: + + DbfFile(const DbfFile &); + DbfFile & operator=(const DbfFile &); + + protected: + + void init_from_scratch(); + + int fd; + + bool At_Eof; + + DbfHeader Header; + + ConcatString Filename; + + public: + + DbfFile(); + ~DbfFile(); + + bool open(const char * path); + + void close(); + + // + // set stuff + // + + // + // get stuff + // + + const DbfHeader * header() const; + + const char * filename() const; + + bool at_eof() const; + + int position() const; // offset in bytes from beginning of file + + bool is_open() const; + + // + // do stuff + // + + void lseek(int offset, int whence = SEEK_SET); // just like lseek(2) + + bool read(unsigned char * buf, int nbytes); + + StringArray subrecord_names(); + + StringArray subrecord_values(int i_rec); + +}; //////////////////////////////////////////////////////////////////////// +inline const DbfHeader * DbfFile::header() const { return ( &Header ); } -#endif /* __VX_SHAPEFILES_DBF_FILE_H__ */ +inline bool DbfFile::at_eof() const { return ( At_Eof ); } +inline bool DbfFile::is_open() const { return ( fd >= 0 ); } //////////////////////////////////////////////////////////////////////// +#endif /* __VX_SHAPEFILES_DBF_FILE_H__ */ +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_gis/shp_array.h b/met/src/libcode/vx_gis/shp_array.h index 419d1ca5bf..bfeeba578b 100644 --- a/met/src/libcode/vx_gis/shp_array.h +++ b/met/src/libcode/vx_gis/shp_array.h @@ -97,6 +97,8 @@ class Shp_Array { int n() const { return ( Nelements ); } + int n_elements() const { return ( Nelements ); } + // // do stuff // diff --git a/met/src/libcode/vx_gis/shp_file.cc b/met/src/libcode/vx_gis/shp_file.cc index 14417c0048..1f311fce2a 100644 --- a/met/src/libcode/vx_gis/shp_file.cc +++ b/met/src/libcode/vx_gis/shp_file.cc @@ -1,8 +1,3 @@ - - -//////////////////////////////////////////////////////////////////////// - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) diff --git a/met/src/libcode/vx_gis/shp_file.h b/met/src/libcode/vx_gis/shp_file.h index 77b2c48418..7a362d9e52 100644 --- a/met/src/libcode/vx_gis/shp_file.h +++ b/met/src/libcode/vx_gis/shp_file.h @@ -1,8 +1,3 @@ - - -//////////////////////////////////////////////////////////////////////// - - // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2021 // ** University Corporation for Atmospheric Research (UCAR) diff --git a/met/src/libcode/vx_gis/shp_poly_record.cc b/met/src/libcode/vx_gis/shp_poly_record.cc index b540b78779..b5876253ab 100644 --- a/met/src/libcode/vx_gis/shp_poly_record.cc +++ b/met/src/libcode/vx_gis/shp_poly_record.cc @@ -50,6 +50,27 @@ static SmartBuffer Buf; //////////////////////////////////////////////////////////////////////// +void ShpPolyRecord::clear() + +{ + + memset(&rh, 0, sizeof(rh)); + memset(&bbox, 0, sizeof(bbox)); + + shape_type = bad_data_int; + + n_parts = 0; + n_points = 0; + + parts.clear(); + points.clear(); + +} + + +//////////////////////////////////////////////////////////////////////// + + void ShpPolyRecord::set(unsigned char * buf) { diff --git a/met/src/libcode/vx_gis/shp_poly_record.h b/met/src/libcode/vx_gis/shp_poly_record.h index faad26478d..18131cadaa 100644 --- a/met/src/libcode/vx_gis/shp_poly_record.h +++ b/met/src/libcode/vx_gis/shp_poly_record.h @@ -30,13 +30,6 @@ #include "shp_array.h" -//////////////////////////////////////////////////////////////////////// - - -// static const int max_shp_poly_parts = 600; -// static const int max_shp_poly_points = (1 << 17); - - //////////////////////////////////////////////////////////////////////// @@ -63,6 +56,8 @@ struct ShpPolyRecord { // this should really be a class, not a struct // + void clear(); + double x_min() const; double x_max() const; @@ -81,7 +76,7 @@ struct ShpPolyRecord { // this should really be a class, not a struct bool is_closed() const; - void toggle_longitudes(); + void toggle_longitudes(); }; diff --git a/met/src/libcode/vx_gis/vx_gis.h b/met/src/libcode/vx_gis/vx_gis.h new file mode 100644 index 0000000000..40ade039ec --- /dev/null +++ b/met/src/libcode/vx_gis/vx_gis.h @@ -0,0 +1,29 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + +#ifndef __VX_GIS_H__ +#define __VX_GIS_H__ + +//////////////////////////////////////////////////////////////////////// + +#include "dbf_file.h" +#include "shx_file.h" +#include "shp_file.h" +#include "shp_types.h" +#include "shp_array.h" +#include "shapetype_to_string.h" +#include "shp_point_record.h" +#include "shp_poly_record.h" + +//////////////////////////////////////////////////////////////////////// + +#endif // __VX_GIS_H__ + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_statistics/Makefile.am b/met/src/libcode/vx_statistics/Makefile.am index 8d3a183cbd..bf1cd72c74 100644 --- a/met/src/libcode/vx_statistics/Makefile.am +++ b/met/src/libcode/vx_statistics/Makefile.am @@ -13,6 +13,7 @@ include ${top_srcdir}/Make-include noinst_LIBRARIES = libvx_statistics.a libvx_statistics_a_SOURCES = \ apply_mask.cc apply_mask.h \ + grid_closed_poly.cc grid_closed_poly.h \ compute_ci.cc compute_ci.h \ contable.cc contable.h \ contable_stats.cc \ diff --git a/met/src/tools/other/gen_vx_mask/grid_closed_poly.cc b/met/src/libcode/vx_statistics/grid_closed_poly.cc similarity index 53% rename from met/src/tools/other/gen_vx_mask/grid_closed_poly.cc rename to met/src/libcode/vx_statistics/grid_closed_poly.cc index a39c63457a..e010cca11b 100644 --- a/met/src/tools/other/gen_vx_mask/grid_closed_poly.cc +++ b/met/src/libcode/vx_statistics/grid_closed_poly.cc @@ -8,7 +8,7 @@ /////////////////////////////////////////////////////////////////////////////// // -// Filename: polyline.cc +// Filename: grid_closed_poly.cc // // Description: // @@ -34,257 +34,187 @@ using namespace std; #include "vx_util.h" #include "grid_closed_poly.h" - /////////////////////////////////////////////////////////////////////////////// // // Code for class GridClosedPoly // /////////////////////////////////////////////////////////////////////////////// - -GridClosedPoly::GridClosedPoly() - -{ - -gcp_init_from_scratch(); - +GridClosedPoly::GridClosedPoly() { + gcp_init_from_scratch(); } - /////////////////////////////////////////////////////////////////////////////// - -GridClosedPoly::~GridClosedPoly() - -{ - -} - +GridClosedPoly::~GridClosedPoly() { } /////////////////////////////////////////////////////////////////////////////// +void GridClosedPoly::clear() { -void GridClosedPoly::clear() - -{ + Polyline::clear(); -Polyline::clear(); - -u_min = u_max = v_min = v_max = 0.0; - -return; + u_min = u_max = v_min = v_max = 0.0; + return; } - /////////////////////////////////////////////////////////////////////////////// +GridClosedPoly::GridClosedPoly(const GridClosedPoly & g) { -GridClosedPoly::GridClosedPoly(const GridClosedPoly & g) - -{ - -gcp_init_from_scratch(); - -gcp_assign(g); + gcp_init_from_scratch(); + gcp_assign(g); } - /////////////////////////////////////////////////////////////////////////////// +GridClosedPoly & GridClosedPoly::operator=(const GridClosedPoly & g) { -GridClosedPoly & GridClosedPoly::operator=(const GridClosedPoly & g) - -{ + if(this == &g) return(*this); -if ( this == &g ) return ( * this ); - -gcp_assign(g); - - -return ( * this ); + gcp_assign(g); + return(*this); } - /////////////////////////////////////////////////////////////////////////////// +void GridClosedPoly::gcp_init_from_scratch() { -void GridClosedPoly::gcp_init_from_scratch() - -{ - -u_min = u_max = v_min = v_max = 0.0; - -return; + u_min = u_max = v_min = v_max = 0.0; + return; } - /////////////////////////////////////////////////////////////////////////////// +void GridClosedPoly::gcp_assign(const GridClosedPoly & g) { -void GridClosedPoly::gcp_assign(const GridClosedPoly & g) + Polyline::assign(g); -{ + u_min = g.u_min; + u_max = g.u_max; -Polyline::assign(g); - -u_min = g.u_min; -u_max = g.u_max; - -v_min = g.v_min; -v_max = g.v_max; - -return; + v_min = g.v_min; + v_max = g.v_max; + return; } - /////////////////////////////////////////////////////////////////////////////// +void GridClosedPoly::add_point(double uu, double vv) { -void GridClosedPoly::add_point(double uu, double vv) + if(n_points == 0) { + u_min = u_max = uu; + v_min = v_max = vv; + } + else { + if(uu < u_min) u_min = uu; + if(uu > u_max) u_max = uu; -{ + if(vv < v_min) v_min = vv; + if(vv > v_max) v_max = vv; + } -if ( n_points == 0 ) { - - u_min = u_max = uu; v_min = v_max = vv; - -} else { - - if ( uu < u_min ) u_min = uu; - if ( uu > u_max ) u_max = uu; - - if ( vv < v_min ) v_min = vv; - if ( vv > v_max ) v_max = vv; - -} - -Polyline::add_point(uu, vv); - - -return; + Polyline::add_point(uu, vv); + return; } - /////////////////////////////////////////////////////////////////////////////// -int GridClosedPoly::is_inside(double u_test, double v_test) const - -{ +int GridClosedPoly::is_inside(double u_test, double v_test) const { // // test bounding box // -if ( u_test < u_min ) return ( 0 ); -if ( u_test > u_max ) return ( 0 ); + if(u_test < u_min) return(0); + if(u_test > u_max) return(0); -if ( v_test < v_min ) return ( 0 ); -if ( v_test > v_max ) return ( 0 ); + if(v_test < v_min) return(0); + if(v_test > v_max) return(0); // // test polyline // -const int status = Polyline::is_inside(u_test, v_test); - -return ( status ); + const int status = Polyline::is_inside(u_test, v_test); + return(status); } - /////////////////////////////////////////////////////////////////////////////// // // Code for class GridClosedPolyArray // /////////////////////////////////////////////////////////////////////////////// +bool GridClosedPolyArray::is_inside(double u_test, double v_test) const { -bool GridClosedPolyArray::is_inside(double u_test, double v_test) const - -{ + if(Nelements == 0) return(false); -if ( Nelements == 0 ) return ( false ); - -int j, status; + int j, status; // - // if the test point is inside **ANY** of the element polylines, return true + // return true if the test point is inside **ANY** of the polylines // + for(j=0; jis_inside(u_test, v_test); + status = e[j]->is_inside(u_test, v_test); - if ( status != 0 ) return ( true ); - -} + if(status != 0) return(true); + } - -return ( false ); + return(false); } - /////////////////////////////////////////////////////////////////////////////// +void GridClosedPolyArray::set(const ShpPolyRecord &r, const Grid &grid) { -void GridClosedPolyArray::set(const ShpPolyRecord & r, const Grid & grid) - -{ - -clear(); - -int i, j, k, m; -int start, stop; -double lat, lon; -double x, y; -GridClosedPoly g; - + clear(); -if ( r.n_parts == 0 ) return; + int i, j, k, m; + int start, stop; + double lat, lon; + double x, y; + GridClosedPoly p; + for(j=0; j<(r.n_parts); j++) { -for (j=0; j<(r.n_parts); ++j) { + p.clear(); - g.clear(); + start = r.start_index(j); + stop = r.stop_index(j); - start = r.start_index(j); + i = stop - start + 1; - stop = r.stop_index(j); + for(k=0; k #include "polyline.h" @@ -33,13 +32,9 @@ #include "shp_poly_record.h" #include "vx_grid.h" - /////////////////////////////////////////////////////////////////////////////// - -class GridClosedPoly : public Polyline - -{ +class GridClosedPoly : public Polyline { protected: @@ -47,8 +42,7 @@ class GridClosedPoly : public Polyline void gcp_assign(const GridClosedPoly &); - // bounding box info - + // bounding box limits double u_min, u_max; double v_min, v_max; @@ -61,36 +55,28 @@ class GridClosedPoly : public Polyline void clear(); - int is_inside(double u_test, double v_test) const; // test bounding box first + // tests bounding box first for efficiency + int is_inside(double u_test, double v_test) const; - void add_point(double, double); // updates bounding box + // updates bounding box for each new point + void add_point(double, double); }; - /////////////////////////////////////////////////////////////////////////////// - -class GridClosedPolyArray : public NCRR_Array - -{ +class GridClosedPolyArray : public NCRR_Array { public: bool is_inside(double u_test, double v_test) const; - void set (const ShpPolyRecord &, const Grid &); + void set(const ShpPolyRecord &, const Grid &); }; - /////////////////////////////////////////////////////////////////////////////// - -#endif // __DATA2D_UTIL_GCP_H__ - +#endif // __GRID_CLOSED_POLY_H__ /////////////////////////////////////////////////////////////////////////////// - - - diff --git a/met/src/libcode/vx_statistics/vx_statistics.h b/met/src/libcode/vx_statistics/vx_statistics.h index db1d6685d1..11c06f68e8 100644 --- a/met/src/libcode/vx_statistics/vx_statistics.h +++ b/met/src/libcode/vx_statistics/vx_statistics.h @@ -19,6 +19,7 @@ #include "apply_mask.h" +#include "grid_closed_poly.h" #include "compute_ci.h" #include "contable.h" #include "met_stats.h" diff --git a/met/src/libcode/vx_tc_util/Makefile.am b/met/src/libcode/vx_tc_util/Makefile.am index 733ef9cd6b..390eb1e5ae 100644 --- a/met/src/libcode/vx_tc_util/Makefile.am +++ b/met/src/libcode/vx_tc_util/Makefile.am @@ -24,6 +24,7 @@ libvx_tc_util_a_SOURCES = \ prob_gen_info.h prob_gen_info.cc \ prob_rirw_pair_info.h prob_rirw_pair_info.cc \ genesis_info.cc genesis_info.h \ + gen_shape_info.cc gen_shape_info.h \ pair_data_genesis.cc pair_data_genesis.h \ tc_hdr_columns.cc tc_hdr_columns.h \ tc_columns.cc tc_columns.h \ diff --git a/met/src/libcode/vx_tc_util/gen_shape_info.cc b/met/src/libcode/vx_tc_util/gen_shape_info.cc new file mode 100644 index 0000000000..f51b92c4e7 --- /dev/null +++ b/met/src/libcode/vx_tc_util/gen_shape_info.cc @@ -0,0 +1,377 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + +using namespace std; + +#include +#include +#include +#include +#include +#include + +#include "gen_shape_info.h" + +//////////////////////////////////////////////////////////////////////// +// +// Code for class GenShapeInfo +// +//////////////////////////////////////////////////////////////////////// + +GenShapeInfo::GenShapeInfo() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +GenShapeInfo::~GenShapeInfo() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +GenShapeInfo::GenShapeInfo(const GenShapeInfo &g) { + clear(); + assign(g); +} + +//////////////////////////////////////////////////////////////////////// + +GenShapeInfo & GenShapeInfo::operator=(const GenShapeInfo &g) { + + if(this == &g) return(*this); + + assign(g); + + return(*this); +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfo::clear() { + + Basin.clear(); + + FileTime = (unixtime) 0; + IssueTime = (unixtime) 0; + + Poly.clear(); + LeadSec.clear(); + ProbVal.clear(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +ConcatString GenShapeInfo::serialize() const { + ConcatString s; + + s << "GenShapeInfo: " + << "Basin = \"" << Basin << "\"" + << ", FileTime = \"" << (FileTime > 0 ? + unix_to_yyyymmdd_hhmmss(FileTime).text() : na_str) << "\"" + << ", IssueTime = \"" << (IssueTime > 0 ? + unix_to_yyyymmdd_hhmmss(IssueTime).text() : na_str) << "\"" + << ", NPoints = " << Poly.n_points + << ", Lat = " << Poly.y_min() << " to " << Poly.y_max() + << ", Lon = " << Poly.x_min() << " to " << Poly.x_max(); + + return(s); + +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfo::assign(const GenShapeInfo &gsi) { + + clear(); + + Basin = gsi.Basin; + + FileTime = gsi.FileTime; + IssueTime = gsi.IssueTime; + + Poly = gsi.Poly; + LeadSec = gsi.LeadSec; + ProbVal = gsi.ProbVal; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfo::set_basin(const char *s) { + Basin = s; + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfo::set_time(unixtime file_ut) { + + FileTime = file_ut; + + // Round the file timestamp to the nearest synoptic time (00, 06, 12, 18) + int h3 = 3*sec_per_hour; + int h6 = 6*sec_per_hour; + IssueTime = ( (int) (file_ut + h3) / h6 ) * h6; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfo::set_poly(const ShpPolyRecord &rec) { + Poly = rec; + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfo::add_prob(int sec, double prob) { + + // Check range of values + if(sec < 0 || prob < 0.0 || prob > 1.0) { + mlog << Error << "\nGenShapeInfo::add_prob() -> " + << "unexpected lead time (" << sec + << ") or probability value (" << prob + << ")!\n\n"; + exit(1); + } + + LeadSec.add(sec); + ProbVal.add(prob); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +const ShpPolyRecord &GenShapeInfo::poly() const { + return(Poly); +} + +//////////////////////////////////////////////////////////////////////// + +double GenShapeInfo::center_lat() const { + return((Poly.y_min() + Poly.y_max())/2.0); +} + +//////////////////////////////////////////////////////////////////////// + +double GenShapeInfo::center_lon() const { + return((Poly.x_min() + Poly.x_max())/2.0); +} + +//////////////////////////////////////////////////////////////////////// + +int GenShapeInfo::n_prob() const { + return(ProbVal.n()); +} + +//////////////////////////////////////////////////////////////////////// + +int GenShapeInfo::lead_sec(int i) const { + + // Check range of values + if(i < 0 || i > LeadSec.n()) { + mlog << Error << "\nGenShapeInfo::lead_sec(int) -> " + << "range check error (" << i << ")!\n\n"; + exit(1); + } + + return(LeadSec[i]); +} + +//////////////////////////////////////////////////////////////////////// + +double GenShapeInfo::prob_val(int i) const { + + // Check range of values + if(i < 0 || i > ProbVal.n()) { + mlog << Error << "\nGenShapeInfo::prob_val(int) -> " + << "range check error (" << i << ")!\n\n"; + exit(1); + } + + return(ProbVal[i]); +} + +//////////////////////////////////////////////////////////////////////// +// +// Code for class GenShapeInfoArray +// +//////////////////////////////////////////////////////////////////////// + +GenShapeInfoArray::GenShapeInfoArray() { + + init_from_scratch(); +} + +//////////////////////////////////////////////////////////////////////// + +GenShapeInfoArray::~GenShapeInfoArray() { + + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +GenShapeInfoArray::GenShapeInfoArray(const GenShapeInfoArray & t) { + + init_from_scratch(); + + assign(t); +} + +//////////////////////////////////////////////////////////////////////// + +GenShapeInfoArray & GenShapeInfoArray::operator=(const GenShapeInfoArray & t) { + + if(this == &t) return(*this); + + assign(t); + + return(*this); +} + +//////////////////////////////////////////////////////////////////////// +// +// Define equality as duplicates with the same file time +// +//////////////////////////////////////////////////////////////////////// + +bool GenShapeInfo::operator==(const GenShapeInfo & gsi) const { + + return(is_duplicate(gsi) && FileTime == gsi.FileTime); +} + + +//////////////////////////////////////////////////////////////////////// +// +// Duplicates have the same basin, issue time, number of points, +// and range of lat/lons +// +//////////////////////////////////////////////////////////////////////// + +bool GenShapeInfo::is_duplicate(const GenShapeInfo & gsi) const { + + return(Basin == gsi.Basin && + IssueTime == gsi.IssueTime && + Poly.n_points == gsi.Poly.n_points && + Poly.x_min() == gsi.Poly.x_min() && + Poly.x_max() == gsi.Poly.x_max() && + Poly.y_min() == gsi.Poly.y_min() && + Poly.y_max() == gsi.Poly.y_max()); +} + + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfoArray::init_from_scratch() { + + clear(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfoArray::clear() { + + GenShape.clear(); + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void GenShapeInfoArray::assign(const GenShapeInfoArray &ga) { + + clear(); + + for(int i=0; i= GenShape.size())) { + mlog << Error << "\nGenShapeInfoArray::operator[](int) -> " + << "range check error for index value " << n << "\n\n"; + exit(1); + } + + return(GenShape[n]); +} + +//////////////////////////////////////////////////////////////////////// + +bool GenShapeInfoArray::add(const GenShapeInfo &gsi, bool check_dup) { + + // Check for duplicates: + // - Skip exact duplicates + // - Replace older duplicates with newer ones + // - Skip older duplicates + if(check_dup) { + + // Loop over existing shapes + for(int i=0; i +#include + +#include "vx_cal.h" +#include "vx_math.h" +#include "vx_util.h" +#include "vx_gis.h" + +//////////////////////////////////////////////////////////////////////// +// +// GenShapeInfo class stores information about genesis shapefiles. +// +//////////////////////////////////////////////////////////////////////// + +class GenShapeInfo { + + private: + + void assign(const GenShapeInfo &); + + // Shapefile information + ConcatString Basin; + unixtime FileTime; + unixtime IssueTime; + + ShpPolyRecord Poly; + + IntArray LeadSec; + NumArray ProbVal; + + public: + + GenShapeInfo(); + ~GenShapeInfo(); + GenShapeInfo(const GenShapeInfo &); + + GenShapeInfo & operator=(const GenShapeInfo &); + bool operator==(const GenShapeInfo &) const; + bool is_duplicate(const GenShapeInfo &) const; + + void clear(); + ConcatString serialize() const; + + // + // set stuff + // + + void set_basin(const char *); + void set_time(unixtime); + void set_poly(const ShpPolyRecord &); + void add_prob(int, double); + + // + // get stuff + // + + const ConcatString & basin() const; + + unixtime file_time() const; + unixtime issue_time() const; + int issue_hour() const; + + const ShpPolyRecord & poly() const; + + double center_lat() const; + double center_lon() const; + + int n_prob() const; + int lead_sec(int) const; + double prob_val(int) const; + + // + // do stuff + // +}; + +//////////////////////////////////////////////////////////////////////// + +inline const ConcatString & GenShapeInfo::basin() const { return(Basin); } + +inline unixtime GenShapeInfo::file_time() const { return(FileTime); } +inline unixtime GenShapeInfo::issue_time() const { return(IssueTime); } +inline int GenShapeInfo::issue_hour() const { return(unix_to_sec_of_day(IssueTime)); } + +//////////////////////////////////////////////////////////////////////// +// +// GenShapeInfoArray class stores an array of GenShapeInfo objects. +// +//////////////////////////////////////////////////////////////////////// + +class GenShapeInfoArray { + + private: + + void init_from_scratch(); + void assign(const GenShapeInfoArray &); + + vector GenShape; + + public: + + GenShapeInfoArray(); + ~GenShapeInfoArray(); + GenShapeInfoArray(const GenShapeInfoArray &); + GenShapeInfoArray & operator=(const GenShapeInfoArray &); + + void clear(); + + ConcatString serialize() const; + + // + // set stuff + // + + bool add(const GenShapeInfo &, bool check_dup = false); + bool has(const GenShapeInfo &) const; + + // + // get stuff + // + + const GenShapeInfo & operator[](int) const; + int n() const; +}; + +//////////////////////////////////////////////////////////////////////// + +inline int GenShapeInfoArray::n() const { return(GenShape.size()); } + +//////////////////////////////////////////////////////////////////////// + +#endif /* __VX_GEN_SHAPE_INFO_H__ */ + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_tc_util/vx_tc_util.h b/met/src/libcode/vx_tc_util/vx_tc_util.h index 76e5448429..89a4434432 100644 --- a/met/src/libcode/vx_tc_util/vx_tc_util.h +++ b/met/src/libcode/vx_tc_util/vx_tc_util.h @@ -24,6 +24,7 @@ #include "prob_rirw_info.h" #include "prob_rirw_pair_info.h" #include "genesis_info.h" +#include "gen_shape_info.h" #include "tc_columns.h" #include "tc_hdr_columns.h" #include "tc_stat_line.h" diff --git a/met/src/tools/core/grid_stat/grid_stat.cc b/met/src/tools/core/grid_stat/grid_stat.cc index 9b2447818a..f420e568c4 100644 --- a/met/src/tools/core/grid_stat/grid_stat.cc +++ b/met/src/tools/core/grid_stat/grid_stat.cc @@ -85,7 +85,7 @@ // output line types. // 034 05/10/16 Halley Gotway Add grid weighting. // 035 05/20/16 Prestopnik J Removed -version (now in command_line.cc) -// 036 02/14/17 Win MET#621 Add nc_pairs_flag.apply_mask +// 036 02/14/17 Win MET #621 Add nc_pairs_flag.apply_mask // 037 05/15/17 Prestopnik P Add shape for regrid, nbrhd and interp // 038 06/26/17 Halley Gotway Add ECLV line type. // 039 08/18/17 Halley Gotway Add fourier decomposition. diff --git a/met/src/tools/core/point_stat/point_stat.cc b/met/src/tools/core/point_stat/point_stat.cc index 7e78611c59..f62f3ae511 100644 --- a/met/src/tools/core/point_stat/point_stat.cc +++ b/met/src/tools/core/point_stat/point_stat.cc @@ -98,7 +98,7 @@ // 047 08/23/21 Seth Linden Add ORANK line type for HiRA. // 048 09/13/21 Seth Linden Changed obs_qty to obs_qty_inc. // Added code for obs_qty_exc. -// 049 12/11/21 Halley Gotway MET#1991 Fix VCNT output. +// 049 12/11/21 Halley Gotway MET #1991 Fix VCNT output. // //////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/other/gen_vx_mask/Makefile.am b/met/src/tools/other/gen_vx_mask/Makefile.am index f8b773183b..5eb4938ea7 100644 --- a/met/src/tools/other/gen_vx_mask/Makefile.am +++ b/met/src/tools/other/gen_vx_mask/Makefile.am @@ -12,8 +12,7 @@ include ${top_srcdir}/Make-include bin_PROGRAMS = gen_vx_mask -gen_vx_mask_SOURCES = gen_vx_mask.cc \ - grid_closed_poly.cc +gen_vx_mask_SOURCES = gen_vx_mask.cc gen_vx_mask_CPPFLAGS = ${MET_CPPFLAGS} gen_vx_mask_LDFLAGS = ${MET_LDFLAGS} gen_vx_mask_LDADD = -lvx_shapedata \ @@ -51,5 +50,4 @@ gen_vx_mask_LDADD = -lvx_shapedata \ -lvx_log \ -lm -lnetcdf_c++4 -lnetcdf -lgsl -lgslcblas -EXTRA_DIST = gen_vx_mask.h \ - grid_closed_poly.h +EXTRA_DIST = gen_vx_mask.h diff --git a/met/src/tools/other/gen_vx_mask/gen_vx_mask.cc b/met/src/tools/other/gen_vx_mask/gen_vx_mask.cc index 1fd7bb9bc2..ed5ac52802 100644 --- a/met/src/tools/other/gen_vx_mask/gen_vx_mask.cc +++ b/met/src/tools/other/gen_vx_mask/gen_vx_mask.cc @@ -25,8 +25,8 @@ // 007 04/08/19 Halley Gotway Add percentile thresholds. // 008 04/06/20 Halley Gotway Generalize input_grid option. // 009 06/01/21 Seth Linden Change -type from optional to required. -// 010 08/30/21 Halley Gotway MET#1891 Fix input and mask fields. -// 011 12/13/21 Halley Gotway MET#1993 Fix -type grid. +// 010 08/30/21 Halley Gotway MET #1891 Fix input and mask fields. +// 011 12/13/21 Halley Gotway MET #1993 Fix -type grid. // //////////////////////////////////////////////////////////////////////// @@ -549,10 +549,8 @@ void apply_poly_mask(DataPlane & dp) { bool inside; double lat, lon; - n_in = 0; - // Check each grid point being inside the polyline - for(x=0; x 0) { + mlog << Debug(2) + << "Scoring track genesis forecasts.\n"; score_track_genesis(best_ga, oper_ta); } // Score EDECK genesis probabilities and write output if(edeck_source.n() > 0) { + mlog << Debug(2) + << "Scoring probability of genesis forecasts.\n"; score_genesis_prob(best_ga, oper_ta); } + // Score genesis shapefiles and write output + if(shape_source.n() > 0) { + mlog << Debug(2) + << "Scoring genesis shapefiles.\n"; + score_genesis_shape(best_ga); + } + // Finish output files finish_txt_files(); @@ -205,6 +233,7 @@ void process_command_line(int argc, char **argv) { // Add function calls for the arguments cline.add(set_genesis, "-genesis", -1); cline.add(set_edeck, "-edeck", -1); + cline.add(set_shape, "-shape", -1); cline.add(set_track, "-track", -1); cline.add(set_config, "-config", 1); cline.add(set_out, "-out", 1); @@ -224,10 +253,10 @@ void process_command_line(int argc, char **argv) { } // Check for the minimum number of arguments - if(genesis_source.n() == 0 && edeck_source.n() == 0) { + if((genesis_source.n() + edeck_source.n() + shape_source.n()) == 0) { mlog << Error << "\nprocess_command_line(int argc, char **argv) -> " - << "at least one of the \"-genesis\" or \"-edeck\" command " - << "line options are required\n\n"; + << "at least one of the \"-genesis\", \"-edeck\", or \"-shape\" " + << "command line options are required\n\n"; usage(); } @@ -255,6 +284,13 @@ void process_command_line(int argc, char **argv) { << ", Model Suffix: " << edeck_model_suffix[i] << "\n"; } + // List the input shapefiles + for(i=0; igenesis_time() - fcst_sa[i].issue_time()) < + (fcst_sa[i].lead_sec(j)); + + if(is_event) { + case_cs << "MATCHES " + << unix_to_yyyymmdd_hhmmss(bgi->genesis_time()) + << " BEST track " << bgi->storm_id() + << " genesis at (" << bgi->lat() << ", " + << bgi->lon() << ").\n"; + } + else { + case_cs << "has NO MATCH in the BEST track.\n"; + } + } + // No Best track match + else { + case_cs << "has NO MATCH in the BEST track.\n"; + } + + mlog << Debug(5) << case_cs; + + // Store pair info + pgi.add_genshape(fcst_sa[i], j, bgi, is_event); + + } // end for j + } // end for i + + return; +} //////////////////////////////////////////////////////////////////////// @@ -966,6 +1131,48 @@ int find_probgen_match(const ProbGenInfo &prob_gi, return(i_best); } +//////////////////////////////////////////////////////////////////////// + +int find_genshape_match(const GenShapeInfo &gsi, + const GenesisInfoArray &bga) { + int i, i_bga; + double x, y; + unixtime min_ut; + + // Time window + unixtime beg_ut = gsi.issue_time(); + unixtime end_ut = beg_ut + shape_prob_search_sec; + + // Load the shape + GridClosedPolyArray p; + p.set(gsi.poly(), conf_info.DLandGrid); + + // Search Best track genesis events for a match + for(i=0,i_bga=bad_data_int; i end_ut) continue; + + // Second, check the polyline + conf_info.DLandGrid.latlon_to_xy(bga[i].lat(), -1.0*bga[i].lon(), x, y); + if(p.is_inside(x, y)) { + + // First match found + if(is_bad_data(i_bga)) { + i_bga = i; + min_ut = bga[i].genesis_time(); + } + // Better match found + else if(bga[i].genesis_time() < min_ut) { + i_bga = i; + min_ut = bga[i].genesis_time(); + } + } + } // end for i + + return(i_bga); +} //////////////////////////////////////////////////////////////////////// @@ -1358,6 +1565,140 @@ void process_edecks(const StringArray &files, return; } +//////////////////////////////////////////////////////////////////////// + +void process_shapes(const StringArray &files, + GenShapeInfoArray &shapes) { + int i, j, k, n_rec, total_probs; + unixtime file_ut; + ConcatString shp_file_name, dbf_file_name; + StringArray sa; + ConcatString cs; + GenShapeInfo gsi; + DbfFile dbf_file; + ShpFile shp_file; + ShpPolyRecord poly_rec; + + // Initialize + total_probs = 0; + + // Process each shapefile + for(i=0; i " + << "the specified shapefile (" << shp_file_name + << ") and corresponding database file (" << dbf_file_name + << ") must both exists!\n\n"; + } + + // Extract the issue time from the filename: + // gtwo_areas_YYYYMMDDHHMM.shp + sa = shp_file_name.basename().split("_."); + if(sa.n() >= 3 && sa[2].length() >= 12) { + cs << cs_erase + << sa[2].substr(0, 8).c_str() << "_" + << sa[2].substr(8, 4).c_str() << "00"; + file_ut = yyyymmdd_hhmmss_to_unix(cs.c_str()); + } + else { + mlog << Error << "\nprocess_shapes() -> " + << "can't extract the issue time from \"" << shp_file_name + << "\" since it does not match the regular expression \"" + << gen_shp_reg_exp << "\"!\n\n"; + exit(1); + } + + // Open the dbf and shp files + dbf_file.open(dbf_file_name.c_str()); + shp_file.open(shp_file_name.c_str()); + + // Store the number of records + n_rec = dbf_file.header()->n_records; + + // Check expected shape types + const ShapeType shape_type = (ShapeType) (shp_file.header()->shape_type); + if(shape_type != shape_type_polygon) { + mlog << Error << "\nprocess_shapes() -> " + << "shapefile type \"" << shapetype_to_string(shape_type) + << "\" in file \"" << shp_file_name + << "\" is not supported!\n\n"; + exit(1); + } + + mlog << Debug(4) << "[File " << i+1 << " of " << files.n() << "]: " + << "Found " << n_rec << " records in file \"" << shp_file_name << "\".\n"; + + // Process each shape + for(j=0; j " + << "hit shp file EOF before reading all records!\n\n"; + exit(1); + } + + // Read the current shape and metadata + shp_file >> poly_rec; + poly_rec.toggle_longitudes(); + sa = dbf_file.subrecord_values(j); + + // Initialize GenShapeInfo + gsi.clear(); + gsi.set_time(file_ut); + gsi.set_basin(string_to_basin_abbr(sa[0]).c_str()); + gsi.set_poly(poly_rec); + + // Parse probabilities from the subrecord values + for(k=0; k= max_n_shape_prob) { + mlog << Warning << "\nprocess_shapes() -> " + << "unexpected number of shapefile probabilities (" + << gsi.n_prob() << ") in record " << j+1 + << " of file \"" << dbf_file_name + << "\"!\n\n"; + continue; + } + + // Store the probability info + gsi.add_prob(shape_prob_lead_hr[gsi.n_prob()]*sec_per_hour, + atoi(sa[k].c_str())/100.0); + } + } // end for k + + // Store this shape, check for duplicates + if(shapes.add(gsi, true)) { + mlog << Debug(5) << "Add new " << gsi.serialize() << "\n"; + total_probs += gsi.n_prob(); + } + + } // end for j + + // Close files + dbf_file.close(); + shp_file.close(); + + } // end for i + + mlog << Debug(3) << "Found a total of " << total_probs + << " probabilities in " << shapes.n() << " genesis shapes from " + << files.n() << " input files.\n"; + + return; +} + //////////////////////////////////////////////////////////////////////// // // Setup the output ASCII files @@ -1391,7 +1732,7 @@ void setup_txt_files(int n_model, int max_n_prob, int n_pair) { break; // Nx2 probabilistic contingency table output: - // 1 header + 1 vx method * # models * #probs * # filters + // 1 header + 1 vx method * # models * # probs * # filters case(i_pct): case(i_pstd): case(i_pjc): @@ -1818,16 +2159,19 @@ void write_ctc_genmpr_row(StatHdrColumns &shc, //////////////////////////////////////////////////////////////////////// void write_pct_stats(ProbGenPCTInfo &pgi) { + + // Check for no data to write + if(pgi.PCTMap.size() == 0) return; + int i, lead_hr, lead_sec; // Setup header columns shc.set_model(pgi.Model.c_str()); shc.set_desc(pgi.VxOpt->Desc.c_str()); shc.set_obtype(conf_info.BestEventInfo.Technique.c_str()); - shc.set_mask(pgi.VxOpt->VxMaskName.empty() ? - na_str : pgi.VxOpt->VxMaskName.c_str()); - shc.set_fcst_var(prob_genesis_name); - shc.set_obs_var (prob_genesis_name); + shc.set_mask(pgi.VxMask.c_str()); + shc.set_fcst_var(pgi.VarName.c_str()); + shc.set_obs_var (pgi.VarName.c_str()); // Write results for each lead time for(i=0; ioutput_map(stat_pjc) != STATOutputType_None) { - write_pct_row(shc, pgi.PCTMap[lead_hr], + write_pjc_row(shc, pgi.PCTMap[lead_hr], pgi.VxOpt->output_map(stat_pjc), 1, 1, stat_at, i_stat_row, txt_at[i_pjc], i_txt_row[i_pjc]); @@ -1871,22 +2215,19 @@ void write_pct_stats(ProbGenPCTInfo &pgi) { // Write PRC output if(pgi.VxOpt->output_map(stat_pjc) != STATOutputType_None) { - write_pct_row(shc, pgi.PCTMap[lead_hr], - pgi.VxOpt->output_map(stat_pjc), + write_prc_row(shc, pgi.PCTMap[lead_hr], + pgi.VxOpt->output_map(stat_prc), 1, 1, stat_at, i_stat_row, txt_at[i_prc], i_txt_row[i_prc]); } // Write out GENMPR if(pgi.VxOpt->output_map(stat_genmpr) != STATOutputType_None) { - shc.set_fcst_var(prob_genesis_name); - shc.set_obs_var (prob_genesis_name); write_pct_genmpr_row(shc, pgi, lead_hr, pgi.VxOpt->output_map(stat_genmpr), stat_at, i_stat_row, txt_at[i_genmpr], i_txt_row[i_genmpr]); } - } // end for i return; @@ -1912,10 +2253,10 @@ void write_pct_genmpr_row(StatHdrColumns &shc, shc.set_alpha(bad_data_double); // Write a line for each matched pair - for(i=0; i 50 are for testing or invests static const int max_best_cyclone_number = 50; +// 2, 5, and 7 days shapefile probabilities +static const int max_n_shape_prob = 3; +static const int shape_prob_lead_hr[max_n_shape_prob] = { + 48, 120, 168 +}; +static const int shape_prob_search_sec = 168*sec_per_hour; + //////////////////////////////////////////////////////////////////////// // // Variables for Command Line Arguments @@ -95,6 +100,7 @@ static const int max_best_cyclone_number = 50; // Input files static StringArray genesis_source, genesis_model_suffix; static StringArray edeck_source, edeck_model_suffix; +static StringArray shape_source; static StringArray track_source, track_model_suffix; static ConcatString config_file; static TCGenConfInfo conf_info; diff --git a/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.cc b/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.cc index 0592c667b4..a39a7fc06c 100644 --- a/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.cc +++ b/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.cc @@ -150,6 +150,7 @@ void TCGenVxOpt::clear() { ValidBeg = ValidEnd = (unixtime) 0; InitHour.clear(); Lead.clear(); + VxMaskConf.clear(); VxMaskName.clear(); VxPolyMask.clear(); VxGridMask.clear(); @@ -228,8 +229,9 @@ void TCGenVxOpt::process_config(Dictionary &dict) { } // Conf: vx_mask - if(nonempty(dict.lookup_string(conf_key_vx_mask).c_str())) { - file_name = replace_path(dict.lookup_string(conf_key_vx_mask)); + VxMaskConf = dict.lookup_string(conf_key_vx_mask); + if(VxMaskConf.nonempty()) { + file_name = replace_path(VxMaskConf); mlog << Debug(2) << "Masking File: " << file_name << "\n"; parse_poly_mask(file_name, VxPolyMask, VxGridMask, VxAreaMask, VxMaskName); @@ -359,7 +361,8 @@ void TCGenVxOpt::process_basin_mask(const Grid &basin_grid, if(!basin_abbr.has(VxBasinMask[i], j)) { mlog << Error << "\nTCGenConfInfo::process_basin_mask() -> " << "\"" << VxBasinMask[i] - << "\" is not a valid basin name!\n\n"; + << "\" is not a valid basin name (" + << write_css(basin_abbr) << ")!\n\n"; exit(1); } @@ -570,6 +573,60 @@ bool TCGenVxOpt::is_keeper(const ProbGenInfo &gi) const { //////////////////////////////////////////////////////////////////////// +bool TCGenVxOpt::is_keeper(const GenShapeInfo &gsi) const { + bool keep = true; + + // ATCF ID and storm name do not apply + + // Initialization time + if((InitBeg > 0 && InitBeg > gsi.issue_time()) || + (InitEnd > 0 && InitEnd < gsi.issue_time()) || + (InitInc.n() > 0 && !InitInc.has(gsi.issue_time())) || + (InitExc.n() > 0 && InitExc.has(gsi.issue_time()))) + keep = false; + + // Initialization hours + if(InitHour.n() > 0 && !InitHour.has(gsi.issue_hour())) + keep = false; + + // Lead and valid times: + // GenShapeInfo objects contain multiple lead/valid times. + // Do not filter by them here. + + // If VxMaskConf set, filter spatially by the center of the shape + if(keep && VxMaskConf.nonempty()) { + + // Poly masking: use center lat/lon + if(VxPolyMask.n_points() > 0 && + !VxPolyMask.latlon_is_inside(gsi.center_lat(), gsi.center_lon())) + keep = false; + + // Area masking: use center lat/lon + if(!VxAreaMask.is_empty()) { + double x, y; + VxGridMask.latlon_to_xy(gsi.center_lat(), -1.0*gsi.center_lon(), x, y); + if(x < 0 || x >= VxGridMask.nx() || + y < 0 || y >= VxGridMask.ny()) { + keep = false; + } + else { + keep = VxAreaMask(nint(x), nint(y)); + } + } + } + // Otherwise, if VxBasinMask set, filter by the GenShapeArea basin + else if(keep && VxBasinMask.n() > 0) { + keep = VxBasinMask.has(gsi.basin()); + } + + // Distance to land does not apply + + // Return the keep status + return(keep); +} + +//////////////////////////////////////////////////////////////////////// + STATOutputType TCGenVxOpt::output_map(STATLineType t) const { return(OutputMap.at(t)); } @@ -1242,12 +1299,16 @@ void ProbGenPCTInfo::clear() { DefaultPCT.clear(); Model.clear(); + VarName.clear(); + VxMask.clear(); + InitBeg = InitEnd = (unixtime) 0; BestBeg = BestEnd = (unixtime) 0; PCTMap.clear(); - FcstGenMap.clear(); - FcstIdxMap.clear(); + ProbGenMap.clear(); + GenShapeMap.clear(); + ProbIdxMap.clear(); BestGenMap.clear(); BestEvtMap.clear(); @@ -1276,16 +1337,18 @@ void ProbGenPCTInfo::set_vx_opt(const TCGenVxOpt *vx_opt) { //////////////////////////////////////////////////////////////////////// -void ProbGenPCTInfo::add(const ProbGenInfo &fgi, int index, - const GenesisInfo *bgi, bool is_event) { - int i; +void ProbGenPCTInfo::add_probgen(const ProbGenInfo &pgi, int index, + const GenesisInfo *bgi, bool is_event) { unixtime ut; - // Store the model name - if(Model.empty()) Model = fgi.technique(); + // Store the model and variable names + if(Model.empty()) Model = pgi.technique(); + if(VarName.empty()) VarName = prob_genesis_name; + if(VxMask.empty()) VxMask = (VxOpt->VxMaskName.empty() ? + na_str : VxOpt->VxMaskName); // Track the range of forecast initalization times - ut = fgi.init(); + ut = pgi.init(); if(InitBeg == 0 || InitBeg > ut) InitBeg = ut; if(InitEnd == 0 || InitEnd < ut) InitEnd = ut; @@ -1297,28 +1360,28 @@ void ProbGenPCTInfo::add(const ProbGenInfo &fgi, int index, } // Current lead time and probability value - int lead_hr = nint(fgi.prob_item(index)); - double prob = fgi.prob(index) / 100.0; + int lead_hr = nint(pgi.prob_item(index)); + double prob = pgi.prob(index) / 100.0; // Add new map entries, if needed if(!LeadTimes.has(lead_hr)) { LeadTimes.add(lead_hr); - vector empty_fgi; + vector empty_pgi; vector empty_idx; vector empty_bgi; vector empty_evt; PCTMap [lead_hr] = DefaultPCT; - FcstGenMap[lead_hr] = empty_fgi; - FcstIdxMap[lead_hr] = empty_idx; + ProbGenMap[lead_hr] = empty_pgi; + ProbIdxMap[lead_hr] = empty_idx; BestGenMap[lead_hr] = empty_bgi; BestEvtMap[lead_hr] = empty_evt; } // Update map entries - FcstGenMap[lead_hr].push_back(&fgi); - FcstIdxMap[lead_hr].push_back(index); + ProbGenMap[lead_hr].push_back(&pgi); + ProbIdxMap[lead_hr].push_back(index); BestGenMap[lead_hr].push_back(bgi); BestEvtMap[lead_hr].push_back(is_event); @@ -1330,3 +1393,58 @@ void ProbGenPCTInfo::add(const ProbGenInfo &fgi, int index, } //////////////////////////////////////////////////////////////////////// + +void ProbGenPCTInfo::add_genshape(const GenShapeInfo &gsi, int index, + const GenesisInfo *bgi, bool is_event) { + unixtime ut; + int lead_hr; + + // Store the model and variable names + if(Model.empty()) Model = "OPER"; + if(VarName.empty()) VarName = genesis_shape_name; + if(VxMask.empty()) VxMask = (VxOpt->VxMaskName.empty() ? + na_str : VxOpt->VxMaskName); + + // Track the range of forecast issue times + ut = gsi.issue_time(); + if(InitBeg == 0 || InitBeg > ut) InitBeg = ut; + if(InitEnd == 0 || InitEnd < ut) InitEnd = ut; + + // Track the range of verifying BEST genesis events + if(bgi) { + ut = bgi->genesis_time(); + if(BestBeg == 0 || BestBeg > ut) BestBeg = ut; + if(BestEnd == 0 || BestEnd < ut) BestEnd = ut; + } + + // Add new map entries, if needed + lead_hr = gsi.lead_sec(index)/sec_per_hour; + if(!LeadTimes.has(lead_hr)) { + + LeadTimes.add(lead_hr); + vector empty_gsi; + vector empty_idx; + vector empty_bgi; + vector empty_evt; + + PCTMap [lead_hr] = DefaultPCT; + GenShapeMap[lead_hr] = empty_gsi; + ProbIdxMap[lead_hr] = empty_idx; + BestGenMap[lead_hr] = empty_bgi; + BestEvtMap[lead_hr] = empty_evt; + } + + // Update map entries + GenShapeMap[lead_hr].push_back(&gsi); + ProbIdxMap[lead_hr].push_back(index); + BestGenMap[lead_hr].push_back(bgi); + BestEvtMap[lead_hr].push_back(is_event); + + // Increment counts + if(is_event) PCTMap[lead_hr].pct.inc_event (gsi.prob_val(index)); + else PCTMap[lead_hr].pct.inc_nonevent(gsi.prob_val(index)); + + return; +} + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.h b/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.h index b3512c0d8c..487ce9245b 100644 --- a/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.h +++ b/met/src/tools/tc_utils/tc_gen/tc_gen_conf_info.h @@ -49,6 +49,14 @@ static const STATLineType txt_file_type[n_txt] = { stat_genmpr // 7 }; +// Output data type names + +static const string genesis_name ("GENESIS"); +static const string genesis_dev_name ("GENESIS_DEV"); +static const string genesis_ops_name ("GENESIS_OPS"); +static const string prob_genesis_name ("PROB_GENESIS"); +static const string genesis_shape_name("GENESIS_SHAPE"); + // Names for output data plane types static const string fgen_str = "fcst_genesis"; @@ -129,6 +137,7 @@ class TCGenVxOpt { NumArray Lead; // Spatial masking information + ConcatString VxMaskConf; ConcatString VxMaskName; MaskPoly VxPolyMask; Grid VxGridMask; @@ -168,8 +177,9 @@ class TCGenVxOpt { const StringArray &); void parse_nc_info(Dictionary &); - bool is_keeper(const GenesisInfo &) const; - bool is_keeper(const ProbGenInfo &) const; + bool is_keeper(const GenesisInfo &) const; + bool is_keeper(const ProbGenInfo &) const; + bool is_keeper(const GenShapeInfo &) const; STATOutputType output_map(STATLineType) const; }; @@ -331,6 +341,9 @@ class ProbGenPCTInfo { ////////////////////////////////////////////////////////////////// ConcatString Model; + ConcatString VarName; + ConcatString VxMask; + unixtime InitBeg, InitEnd; unixtime BestBeg, BestEnd; const TCGenVxOpt* VxOpt; @@ -340,10 +353,11 @@ class ProbGenPCTInfo { map PCTMap; // Map of lead times to vectors of pair info - map > FcstGenMap; - map > FcstIdxMap; - map > BestGenMap; - map > BestEvtMap; + map > ProbGenMap; + map > GenShapeMap; + map > ProbIdxMap; + map > BestGenMap; + map > BestEvtMap; ////////////////////////////////////////////////////////////////// @@ -351,8 +365,11 @@ class ProbGenPCTInfo { void set_vx_opt(const TCGenVxOpt *); - void add(const ProbGenInfo &, int, - const GenesisInfo *, bool); + void add_probgen(const ProbGenInfo &, int, + const GenesisInfo *, bool); + + void add_genshape(const GenShapeInfo &, int, + const GenesisInfo *, bool); }; diff --git a/test/config/TCGenConfig_shape b/test/config/TCGenConfig_shape new file mode 100644 index 0000000000..767349741e --- /dev/null +++ b/test/config/TCGenConfig_shape @@ -0,0 +1,292 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// TC-Gen configuration file. +// +// For additional information, please see the MET User's Guide. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// ATCF file format reference: +// http://www.nrlmry.navy.mil/atcf_web/docs/database/new/abrdeck.html +// + +//////////////////////////////////////////////////////////////////////////////// +// +// Genesis event definition criteria +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Model initialization frequency in hours, starting at 0 +// +init_freq = 6; + +// +// Valid hour frequency to be analyzed in hours, starting at 0 +// +valid_freq = 6; + +// +// Forecast hours to be searched for genesis events +// +fcst_hr_window = { + beg = 6; + end = 120; +} + +// +// Minimum track duration for genesis event in hours +// +min_duration = 12; + +// +// Forecast genesis event criteria. Defined as tracks reaching the specified +// intensity category, maximum wind speed threshold, and minimum sea-level +// pressure threshold. The forecast genesis time is the valid time of the first +// track point where all of these criteria are met. +// +fcst_genesis = { + vmax_thresh = NA; + mslp_thresh = NA; +} + +// +// BEST track genesis event criteria. Defined as tracks reaching the specified +// intensity category, maximum wind speed threshold, and minimum sea-level +// pressure threshold. The BEST track genesis time is the valid time of the +// first track point where all of these criteria are met. +// +best_genesis = { + technique = "BEST"; + category = [ "TD", "TS" ]; + vmax_thresh = NA; + mslp_thresh = NA; +} + +// +// Operational track technique name +// +oper_technique = "CARQ"; + +//////////////////////////////////////////////////////////////////////////////// +// +// Track filtering options +// May be specified separately in each filter array entry. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Array of dictionaries containing the track filtering options +// If empty, a single filter is defined using the top-level settings. +// +filter = [ + { desc = "SUMMER_INIT"; init_beg = "20210601"; init_end = "20210901"; }, + { desc = "00Z_12Z_INIT"; init_hour = [ "00", "12" ]; }, + { desc = "06Z_18Z_INIT"; init_hour = [ "06", "18" ]; }, + { desc = "AL_EP_BASIN"; basin_mask = ["AL", "EP" ]; }, + { desc = "AL_BASIN"; basin_mask = "AL"; }, + { desc = "EP_BASIN"; basin_mask = "EP"; } +]; + +// +// Description written to output DESC column +// +desc = "ALL"; + +// +// Forecast ATCF ID's +// If empty, all ATCF ID's found will be processed. +// Statistics will be generated separately for each ATCF ID. +// +model = []; + +// +// BEST and operational track storm identifiers +// +storm_id = []; + +// +// BEST and operational track storm names +// +storm_name = []; + +// +// Forecast and operational initialization times to include or exclude +// +init_beg = ""; +init_end = ""; +init_inc = []; +init_exc = []; + +// +// Forecast, BEST, and operational valid time window +// +valid_beg = ""; +valid_end = ""; + +// +// Forecast and operational initialization hours +// +init_hour = []; + +// +// Forecast and operational lead times in hours +// +lead = []; + +// +// Spatial masking region (path to gridded data file or polyline file) +// +vx_mask = ""; + +// +// Spatial masking of hurricane basin names from the basin_file +// +basin_mask = []; + +// +// Distance to land threshold +// +dland_thresh = NA; + +//////////////////////////////////////////////////////////////////////////////// +// +// Matching and scoring options +// May be specified separately in each filter array entry. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Genesis matching logic. Compare the forecast genesis point to all points in +// the Best track (TRUE) or the single Best track genesis point (FALSE). +// +genesis_match_point_to_track = TRUE; + +// +// Radius in km to search for a matching genesis event +// +genesis_match_radius = 500; + +// +// Time window in hours, relative to the model genesis time, to search for a +// matching Best track point +// +genesis_match_window = { + beg = 0; + end = 0; +} + +// +// Radius in km for a development scoring method hit +// +dev_hit_radius = 500; + +// +// Time window in hours, relative to the model genesis time, for a development +// scoring method hit +// +dev_hit_window = { + beg = -24; + end = 24; +} + +// +// Time window in hours for the Best track genesis minus model initialization +// time difference for an operational scoring method hit +// +ops_hit_window = { + beg = 0; + end = 48; +} + +// +// Discard genesis forecasts for initializations at or after the matching +// BEST track genesis time +// +discard_init_post_genesis_flag = TRUE; + +// +// Scoring methods to be applied +// +dev_method_flag = TRUE; +ops_method_flag = TRUE; + +//////////////////////////////////////////////////////////////////////////////// +// +// Output options +// May be specified separately in each filter array entry. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Probability of genesis thresholds +// +prob_genesis_thresh = ==0.25; + +// +// Confidence interval alpha value +// +ci_alpha = 0.05; + +// +// Statistical output types +// +output_flag = { + fho = NONE; + ctc = NONE; + cts = NONE; + pct = BOTH; + pstd = STAT; + pjc = STAT; + prc = STAT; + genmpr = NONE; +} + + +// +// NetCDF genesis pair counts +// +nc_pairs_flag = FALSE; + +// +// Specify which track points should be counted by thresholding the track point +// valid time minus genesis time difference. +// +valid_minus_genesis_diff_thresh = NA; + +// +// Count unique BEST track genesis event locations (TRUE) versus counting the +// location for all pairs (FALSE). +// +best_unique_flag = TRUE; + +//////////////////////////////////////////////////////////////////////////////// +// +// Global settings +// May only be specified once. +// +//////////////////////////////////////////////////////////////////////////////// + +// +// Specify the NetCDF output of the gen_dland tool containing a gridded +// representation of the minimum distance to land. +// +dland_file = "MET_BASE/tc_data/dland_global_tenth_degree.nc"; + +// +// Specify the NetCDF file containing a gridded representation of the +// global basins. +// +basin_file = "MET_BASE/tc_data/basin_global_tenth_degree.nc"; + +// +// NetCDF genesis pairs grid +// +nc_pairs_grid = "G003"; + +// +// Indicate a version number for the contents of this configuration file. +// The value should generally not be modified. +// +version = "V10.1.0"; diff --git a/test/xml/unit_tc_gen.xml b/test/xml/unit_tc_gen.xml index 86245745d8..6146f1d3ec 100644 --- a/test/xml/unit_tc_gen.xml +++ b/test/xml/unit_tc_gen.xml @@ -53,4 +53,23 @@ + + + + &MET_BIN;/tc_gen + \ + -shape &DATA_DIR;/genesis/shape/atl \ + -shape &DATA_DIR;/genesis/shape/epac \ + -track &DATA_DIR;/genesis/atcf/2021 \ + -config &CONFIG_DIR;/TCGenConfig_shape \ + -out &OUTPUT_DIR;/tc_gen/tc_gen_2021_shape \ + -log &OUTPUT_DIR;/tc_gen/tc_gen_2021_shape.log \ + -v 3 + + + &OUTPUT_DIR;/tc_gen/tc_gen_2021_shape.stat + &OUTPUT_DIR;/tc_gen/tc_gen_2021_shape_pct.txt + + +