diff --git a/CMakeLists.txt b/CMakeLists.txt index 4b3831e3..bdc4c7ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,7 +6,7 @@ ENABLE_LANGUAGE(CXX) ENABLE_TESTING() SET(GFAST_VERSION_MAJOR 1) SET(GFAST_VERSION_MINOR 2) -SET(GFAST_VERSION_PATCH 3) +SET(GFAST_VERSION_PATCH 4) SET(GFAST_VERSION ${GFAST_VERSION_MAJOR}.${GFAST_VERSION_MINOR}.${GFAST_VERSION_PATCH}) IF ("${GFAST_INSTANCE}" STREQUAL "") SET(GFAST_INSTANCE "DEFAULT") @@ -114,7 +114,8 @@ SET(SRCS_CORE #ADD_SUBDIRECTORY(src/eewUtils) SET(SRCS_EEW src/eewUtils/driveCMT.c src/eewUtils/driveFF.c src/eewUtils/driveGFAST.c src/eewUtils/drivePGD.c src/eewUtils/makeXML.c src/eewUtils/parseCoreXML.c - src/eewUtils/setLogFileNames.c) + src/eewUtils/setLogFileNames.c src/eewUtils/fillCoreEventInfo.c + src/eewUtils/sendXMLFilter.c) #ADD_SUBDIRECTORY(src/traceBuffer) SET(SRCS_H5TB src/traceBuffer/h5/copyTraceBufferToGFAST.c src/traceBuffer/h5/finalize.c src/traceBuffer/h5/getData.c src/traceBuffer/h5/getDoubleArray.c diff --git a/CMakeModules/FindXML2.cmake b/CMakeModules/FindXML2.cmake index 4923848b..690aba65 100644 --- a/CMakeModules/FindXML2.cmake +++ b/CMakeModules/FindXML2.cmake @@ -19,7 +19,7 @@ find_library(LIBXML2_LIBRARY NAMES xml2 libxml2 include(FindPackageHandleStandardArgs) # handle the QUIETLY and REQUIRED arguments and set LIBXML2_FOUND to TRUE # if all listed variables are TRUE -find_package_handle_standard_args(LibXml2 DEFAULT_MSG +find_package_handle_standard_args(XML2 DEFAULT_MSG LIBXML2_LIBRARY LIBXML2_INCLUDE_DIR) mark_as_advanced(LIBXML2_INCLUDE_DIR LIBXML2_LIBRARY ) diff --git a/ChangeLog b/ChangeLog new file mode 100644 index 00000000..7de780cb --- /dev/null +++ b/ChangeLog @@ -0,0 +1,83 @@ +version 1.2.4 2023-08-03 Carl Ulberg + * Summary: + - Add internal/external message versions in driveGFAST. The internal versions start at v-1, + then decrement by 1 on every single iteration. The external version numbers start at v0, + and increment by 1 on pgd xml messages that are actually sent to the Solution Aggregator. + Every message is still saved in the run/out directory - if the message was sent to the SA, + it will have the positive, external version number. If it was not sent, it will have the + negative, internal version number. + - Change how pgd xml timestamp is defined. Previously it represented the time at the start + of the current iteration. Now it represents the time closest to when the message is + actually sent (really from when the message is generated, but then it is sent in the next + step). + - Add publish change thresholds. Among the other send filter logic, the change thresholds + ensure that GFAST won't send identical (or nearly identical) messages to the SA. + - Refactor some code out of driveGFAST (easier to unit test), fillCoreEventInfo and + sendXMLFilter + - New unit tests +version 1.2.3 2023-04-10 Carl Ulberg + * Summary: + - Fully use Q channel with "goodness/use" flag. If configuration is set, GFAST will ignore + data points and PGD observations with a use flag of 0, as determined by the CWU/Fastlane + processing. + - Use "sleep" instead of "while true" loop for each 1s iteration. Saves processor time, + although may lead to slight creep in iteration start time. +version 1.2.2 2022-12-28 Carl Ulberg + * Summary: + - Bring Q channel through processing. Not entirely set up, but the Q channel is mostly + available when calculating PGD, along with the 6 other channels: ENZ/ENZ uncertainty. +version 1.2.1 2022-11-28 Carl Ulberg + * Summary: + - Add hashmap to tb2 processing to speed up get/unpack messages. + - Add raw observation uncertainty thresholds in PGD processing. This is intended to handle + clock correction steps in the data. +version 1.2.0 2022-10-04 Carl Ulberg + * Summary: + - Add GFAST to ShakeAlert repo + - Use DMMessageReceiver for reading AMQ messages. + - Add GFAST to ShakeAlert unit test coverage report. + - Add config options for turning on/off pgd/cmt/ff processing. + - Reduce hdf5 archive size +### Versions before this were in the PNSN github repo at https://github.com/pnsn/GFAST ### +version 1.1.10 2022-07-27 Carl Ulberg + * Summary: + - logging updates, etc before adding to ShakeAlert repo +version 1.1.9 2022-06-29 Carl Ulberg + * Summary: + - pgd_mag_sigma throttling and updated instance header +version 1.1.8 2022-05-27 Carl Ulberg + * Summary: + - add 'assoc' attribute to pgd xml +version 1.1.7 2022-04-28 Carl Ulberg + * Summary: + - min/max pgd thresholds +version 1.1.6 2022-04-01 Carl Ulberg + * Summary: + - update cmt quakeml +version 1.1.5 2022-03-01 Carl Ulberg + * Summary: + - add time/mag based magnitude uncertainty +version 1.1.4b 2022-02-18 Carl Ulberg + * Summary: + - throttle based on time-varying pgd_obs +version 1.1.4 2022-01-14 Carl Ulberg + * Summary: + - throttle based on pgd_obs +version 1.1.3b 2021-12-22 Carl Ulberg + * Summary: + - use nan uncert +version 1.1.3 2021-12-02 Carl Ulberg + * Summary: + - Throttle xml based on SA mag +version 1.1.2 2021-11-17 Carl Ulberg + * Summary: + - Add strike, dip to ff xml message (using old xml writer) +version 1.1.1 2021-11-09 Carl Ulberg + * Summary: + - bug fix from 1.1.0 to only add pgd observations that were used in pgd inversion +version 1.1.0 2021 Carl Ulberg + * Summary: + - Use dmlib to encode pgd xml, with pgd observations +version 1.0.0 2021 Carl Ulberg + * Summary: + - Initial merge of SAdev and 2020 branches from GFAST PNSN github repo \ No newline at end of file diff --git a/README.md b/README.md index fd2309c5..2eb857fb 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,7 @@ This is the source code for Geodetic First Approximation of Size and Timing (GFA ## Recent updates +2023/12/20 CWU: Sync most recent changes from ShakeAlert GFAST version to the PNSN repo (GFAST v1.2.4) 2023/03/20 CWU: Sync most recent changes from ShakeAlert GFAST version to the PNSN repo (GFAST v1.2.3-beta) ~2021/10 CWU: merge 2020 and SAdev branches into one branch to use going forward, called "development" 2020/05/27 MTH: Working on new branch (=2020) with cleaned up dependencies and build instructions. @@ -14,7 +15,7 @@ This is the source code for Geodetic First Approximation of Size and Timing (GFA 2. include contains the GFAST C include files. 3. legacy is the original Python source code. 4. src contains the source code. - + src/activeMQ contains C++ readers/writers and C interfaes for using ActiveMQ. + + src/activeMQ contains C++ readers/writers and C interfaces for using ActiveMQ. + src/core contains the core GFAST computations. + src/dmlib contains functionality making use of some ShakeAlert libraries, primarily related to ActiveMQ connections and xml encoding/decoding. Access to the ShakeAlert code repository would be necessary to build using these codes. + src/eewUtils contains application specific functions for performing the earthquake early warning tasks. @@ -22,6 +23,7 @@ This is the source code for Geodetic First Approximation of Size and Timing (GFA + src/traceBuffer contains routines for reading an Earthworm ring and converting to a GFAST specific buffer. + src/uw contains functions specific to the University of Washington and Amazon project. + src/xml contains functions for certainly writing and potentially reading QuakeML and ShakeAlert specific XML. + + src/tests contains googletest unit tests, primarily tuned to work within a ShakeAlert environment/repo 5. unit\_tests contains some simple regression tests for the core modules. # Building GFAST @@ -186,7 +188,7 @@ But here's how to do it in stand-alone >make // will build lib/libiscl_shared.so and lib/libiscl_static.a -7. [compearth/mtbeach](https://github.com/carltape/mtbeach.git) for the moment tensor decompositions. Note, this was in Ben Baker's [compearth repository](https://github.com/bakerb845/compearth), but was merged back into Carl Tape's repository, and the moment tensor part was renames mtbeach. Also, you'll need to descend into c_src. +7. [compearth/mtbeach](https://github.com/carltape/mtbeach.git) for the moment tensor decompositions. Note, this was in Ben Baker's [compearth repository](https://github.com/bakerb845/compearth), but was merged back into Carl Tape's repository, and the moment tensor part was renamed mtbeach. Also, you'll need to descend into c_src. Clone and build >git clone https://github.com/carltape/mtbeach.git @@ -243,7 +245,7 @@ But here's how to do it in stand-alone 9a. [ActiveMQ](http://activemq.apache.org/) The C++ portion. These will require other things that you likely already have like libssl, libcrypto, and the Apache runtime library. - Note: [MTH] I haven't found the APR to be necessary + Note: [MTH] I haven't found the APR to be necessary Note: For building GFAST, the instructions here for downloading and installing activemq-cpp are sufficient. However, as soon as gfast_eew starts, it will try to connect diff --git a/include/gfast_activeMQ.h b/include/gfast_activeMQ.h index 48505bbe..573fe715 100644 --- a/include/gfast_activeMQ.h +++ b/include/gfast_activeMQ.h @@ -50,10 +50,6 @@ void activeMQ_producer_finalize(void *producerIn); /* Send a text message */ int activeMQ_producer_sendMessage(void *producerIn, const char *message); -/* Convenience function to set the tcp URI request */ -char *activeMQ_setTcpURIRequest(const char *url, - const int msReconnect, - const int maxAttempts); /* Initialize the ActiveMQ producer */ void *activeMQ_producer_initialize(const char AMQuser[], const char AMQpassword[], @@ -63,6 +59,11 @@ void *activeMQ_producer_initialize(const char AMQuser[], const bool clientAck, const int verbose, int *ierr); +/* Convenience function to set the tcp URI request */ +char *activeMQ_setTcpURIRequest(const char *url, + const int msReconnect, + const int maxAttempts); + #ifndef __cplusplus #define GFAST_activeMQ_consumer_initialize(...) \ diff --git a/include/gfast_config.h b/include/gfast_config.h index 1099cd7f..1a006a61 100644 --- a/include/gfast_config.h +++ b/include/gfast_config.h @@ -1,7 +1,7 @@ #ifndef _gfast_config__h_ #define _gfast_config__h_ 1 -#define GFAST_VERSION "gfast-1.2.3-beta" -#define GFAST_MAXMSG_LEN 4096 +#define GFAST_VERSION "gfast-1.2.4-2023-10-12" +#define GFAST_MAXMSG_LEN 8192 #define MAX_OUTPUT_INTERVALS 16 #define MAX_THROTTLING_THRESHOLDS 100 #define MAX_SIGMA_LOOKUP_VALUES 100 diff --git a/include/gfast_core.h b/include/gfast_core.h index 7ae99062..f16a7792 100644 --- a/include/gfast_core.h +++ b/include/gfast_core.h @@ -153,10 +153,10 @@ int core_data_readSiteMaskFile(const char *siteMaskFile, //----------------------------------------------------------------------------// /* Frees memory on an event structure */ void core_events_freeEvents(struct GFAST_activeEvents_struct *events); -/* Convenience function to find min origin time in event list */ -double core_events_getMinOriginTime(struct GFAST_props_struct props, - struct GFAST_activeEvents_struct events, - bool *lnoEvents); +// /* Convenience function to find min origin time in event list */ +// double core_events_getMinOriginTime(struct GFAST_props_struct props, +// struct GFAST_activeEvents_struct events, +// bool *lnoEvents); /* Adds a new event to the event list */ bool core_events_newEvent(struct GFAST_shakeAlert_struct SA, struct GFAST_activeEvents_struct *events, @@ -166,12 +166,12 @@ bool core_events_syncXMLStatusWithEvents(struct GFAST_activeEvents_struct *event /* Print the events in the event list */ void core_events_printEvents(struct GFAST_shakeAlert_struct SA); -/* Remove a cancelled event from the events list */ -bool core_events_removeCancelledEvent(const char *evid, - const double currentTime, - const int verbose, - struct GFAST_shakeAlert_struct SA, - struct GFAST_activeEvents_struct *events); +// /* Remove a cancelled event from the events list */ +// bool core_events_removeCancelledEvent(const char *evid, +// const double currentTime, +// const int verbose, +// struct GFAST_shakeAlert_struct SA, +// struct GFAST_activeEvents_struct *events); /* Look through the events list and remove expired events */ int core_events_removeExpiredEvents(const double maxTime, const double currenTime, @@ -183,10 +183,10 @@ bool core_events_removeExpiredEvent(const double maxtime, const int verbose, struct GFAST_shakeAlert_struct SA, struct GFAST_activeEvents_struct *events); -/* Potentially add the shakeAlert event to the event list */ -bool core_events_updateEvent(struct GFAST_shakeAlert_struct SA, - struct GFAST_activeEvents_struct *events, - int *ierr); +// /* Potentially add the shakeAlert event to the event list */ +// bool core_events_updateEvent(struct GFAST_shakeAlert_struct SA, +// struct GFAST_activeEvents_struct *events, +// int *ierr); //----------------------------------------------------------------------------// @@ -626,9 +626,12 @@ int core_scaling_readRawSigmaThresholdLookupFile(const char *rawSigmaThresholdLo /* Parse the packed quality channel, return chi^2 value */ double core_waveformProcessor_parseQChannelChi2CWU( const double q_value); -/* Parse the packed quality channel, return mapped chi^2 value*/ +/* Parse the packed quality channel, return mapped chi^2 value */ int core_waveformProcessor_parseQChannelChi2CWUmap( const double q_value); +/* Parse the packed quality channel, return goodness (0 or 1) */ +int core_waveformProcessor_parseQChannelGoodness( + const double q_value); /* Compute the offset in north, east, and up */ int core_waveformProcessor_offset(const int utm_zone, const double svel_window, @@ -705,24 +708,23 @@ double core_waveformProcessor_peakDisplacementHelper( core_data_readMetaDataFile(__VA_ARGS__) #define GFAST_core_data_readSiteMaskFile(...) \ core_data_readSiteMaskFile(__VA_ARGS__) - + #define GFAST_core_events_freeEvents(...) \ core_events_freeEvents(__VA_ARGS__) -#define GFAST_core_events_getMinOriginTime(...) \ - core_events_getMinOriginTime(__VA_ARGS__) #define GFAST_core_events_newEvent(...) \ core_events_newEvent(__VA_ARGS__) #define GFAST_core_events_syncXMLStatusWithEvents(...) \ core_events_syncXMLStatusWithEvents(__VA_ARGS__) - #define GFAST_core_events_printEvents(...) \ core_events_printEvents(__VA_ARGS__) -#define GFAST_core_events_removeCancelledEvent(...) \ - core_events_removeCancelledEvent(__VA_ARGS__) #define GFAST_core_events_removeExpiredEvent(...) \ core_events_removeExpiredEvent(__VA_ARGS__) -#define GFAST_core_events_updateEvent(...) \ - core_events_updateEvent(__VA_ARGS__) +// #define GFAST_core_events_getMinOriginTime(...) +// core_events_getMinOriginTime(__VA_ARGS__) +// #define GFAST_core_events_removeCancelledEvent(...) +// core_events_removeCancelledEvent(__VA_ARGS__) +// #define GFAST_core_events_updateEvent(...) +// core_events_updateEvent(__VA_ARGS__) #define GFAST_core_ff_faultPlaneGridSearch(...) \ core_ff_faultPlaneGridSearch(__VA_ARGS__) diff --git a/include/gfast_eewUtils.h b/include/gfast_eewUtils.h index 7eced1c4..6303f7c4 100644 --- a/include/gfast_eewUtils.h +++ b/include/gfast_eewUtils.h @@ -96,6 +96,31 @@ void eewUtils_setLogFileNames(const char *eventid, char infoLogFileName[PATH_MAX], char debugLogFileName[PATH_MAX], char warnLogFileName[PATH_MAX]); +/* Fill a coreInfo struct with the given values */ +int eewUtils_fillCoreEventInfo( + const char *evid, + const int version, + const double SA_lat, + const double SA_lon, + const double SA_depth, + const double SA_mag, + const double SA_time, + const int num_stations, + struct coreInfo_struct *core); +/* Helper function for eewUtils_sendXMLFilter, check change thresholds */ +bool eewUtils_changeThresholdsExceeded( + const struct GFAST_props_struct *props, + const struct coreInfo_struct *core, + const struct coreInfo_struct *last_sent_core); +/* Determine if PGD message should be sent */ +bool eewUtils_sendXMLFilter( + const struct GFAST_props_struct *props, + const struct GFAST_shakeAlert_struct *SA, + const struct GFAST_pgdResults_struct *pgd, + const struct GFAST_peakDisplacementData_struct *pgd_data, + const struct coreInfo_struct *core, + const struct coreInfo_struct *last_sent_core, + const double age_of_event); #define GFAST_eewUtils_driveCMT(...) \ eewUtils_driveCMT(__VA_ARGS__) diff --git a/include/gfast_struct.h b/include/gfast_struct.h index 51b8b201..d0792f26 100644 --- a/include/gfast_struct.h +++ b/include/gfast_struct.h @@ -56,11 +56,18 @@ struct GFAST_pgd_props_struct double u_raw_sigma_threshold; /*!< Maximum raw sigma (up) to allow in pd calculations (cm) */ double n_raw_sigma_threshold; /*!< Maximum raw sigma (east) to allow in pd calculations (cm) */ double e_raw_sigma_threshold; /*!< Maximum raw sigma (north) to allow in pd calculations (cm) */ + int q_value_threshold; /*!< Don't use pd observations for a station if q value is below threshold */ double pgd_sigma_throttle; /*!< Maximum pgd magnitude uncertainty to allow in sending xml */ double SA_mag_threshold; /*!< Magnitude threshold above which to send xml messages */ double minimum_pgd_cm; /*!< Minimum value to include a pgd value in inversion (cm) */ double maximum_pgd_cm; /*!< Maximum value to include a pgd value in inversion (cm) */ int max_assoc_stations; /*!< Maximum stations to add the 'assoc' tag to in xml */ + double change_threshold_mag; /*!< Minimum magnitude change threshold */ + double change_threshold_mag_uncer; /*!< Minimum magnitude uncertainty change threshold */ + double change_threshold_lat; /*!< Minimum latitude change threshold */ + double change_threshold_lon; /*!< Minimum longitude change threshold */ + double change_threshold_orig_time; /*!< Minimum origin time change threshold */ + int change_threshold_num_stations; /*!< Minimum number of stations change threshold */ }; struct GFAST_cmt_props_struct @@ -493,20 +500,6 @@ struct GFAST_activeEvents_struct char pad1[4]; }; -struct GFAST_xml_output_status -{ - char eventid[128]; /*!< Event ID. */ - int version; /*!< GFAST iteration version for this event ID. */ - bool interval_complete[3][MAX_OUTPUT_INTERVALS]; /*!< intervals_complete[] for pgd + cmt + ff = 3 */ -}; - -struct GFAST_activeEvents_xml_status -{ - struct GFAST_xml_output_status *SA_status; - int nev; - //char pad1[4]; /* MTH: seems to be for byte alignment but **check this** */ -}; - struct coreInfo_struct { char id[128]; /*!< Event ID */ @@ -582,6 +575,22 @@ struct coreInfo_struct int numStations; /*!< Number of stations used in event. */ }; +struct GFAST_xml_output_status +{ + char eventid[128]; /*!< Event ID. */ + int internal_version; /*!< GFAST internal version for this event ID. Always < 0 */ + int external_version; /*!< GFAST external version for this event ID. Always >= 0 */ + struct coreInfo_struct last_sent_core; /*!< Core info of the last sent message */ + bool interval_complete[3][MAX_OUTPUT_INTERVALS]; /*!< intervals_complete[] for pgd + cmt + ff = 3 */ +}; + +struct GFAST_activeEvents_xml_status +{ + struct GFAST_xml_output_status *SA_status; + int nev; + //char pad1[4]; /* MTH: seems to be for byte alignment but **check this** */ +}; + struct GFAST_xmlMessages_struct { char **evids; /*!< Event IDs. */ diff --git a/run/bin/gfast_eew_env.sh b/run/bin/gfast_eew_env.sh index 070ecc43..0f52b24f 100644 --- a/run/bin/gfast_eew_env.sh +++ b/run/bin/gfast_eew_env.sh @@ -4,7 +4,7 @@ echo -e 'Version: $Id: gfast_env.sh 000 2023-01-11 18:00:00Z ulbergc $\n' # With "LEAPSECONDS=", times will be handled consistently internally # pending addressing iscl <-> dmlib issues. -export LEAPSECONDS= +export LEAPSECONDS=no-correction-for-gfast # export LEAPSECONDS=/app/share/etc/leapseconds #env | grep LEAPSECONDS diff --git a/run/params/gfast.props.SA b/run/params/gfast.props.SA index 3f3015fd..df390f4a 100644 --- a/run/params/gfast.props.SA +++ b/run/params/gfast.props.SA @@ -104,16 +104,25 @@ pgd_min_sites = 4 # Lookup table for non-density corrected value to generate magnitude # uncertainty based on time and magnitude sigmaLookupFile = /app/gfast/run/params/M99.txt + # Lookup table for time-dependant PGD thresholds # columns: time (s) | threshold (cm) | n stations pgdThresholdLookupFile = /app/gfast/run/params/pgd_threshold_PW.txt + # Positional uncertainty thresholds for including peak displacement observations # If not specified, allow any or NAN uncertainties rawSigmaThresholdLookupFile = /app/gfast/run/params/raw_sigma_threshold_PW.txt + +# Q channel threshold (as of 3/22/23, using the last digit, which is the +# encoded goodness value). If goodness is below this threshold (<), ignore the observation +q_value_threshold = -1 + # If defined, don't send xml if pgd_sigma is greater than this value pgd_sigma_throttle = 0.5 + # Message throttling logic SA_mag_threshold = 6.0 + # Minimum and maximum thresholds for using pgd values in inversion (cm). Any # values equal or outside these bounds will be ignored in the inversion. minimum_pgd_cm = 0.0 @@ -122,6 +131,17 @@ maximum_pgd_cm = 3500.0 # by the ShakeAlert Solution Aggregator in associating algorithm messages. max_assoc_stations = 6 +# Change thresholds to decide when to send a message, assuming other throttles +# pass. If any of these are met or exceeded, a new message would be sent. If +# none are exceeded, message would not be sent. +# If not present, threshold will be ignored +change_threshold_mag = 0.05 +change_threshold_mag_uncer = 0.05 +change_threshold_lat = 0.05 +change_threshold_lon = 0.05 +change_threshold_orig_time = 1 +change_threshold_num_stations = 1 + ################################################################################ # CMT/Finite Fault Properties # ################################################################################ diff --git a/src/Make.include.Linux b/src/Make.include.Linux index 49c2e98b..3e247c0a 100644 --- a/src/Make.include.Linux +++ b/src/Make.include.Linux @@ -38,8 +38,8 @@ LAPACKE_INCL=-I /usr/include/lapacke LAPACKE_LIB= $(LIB_DIR)/liblapack.so.3 $(LIB_DIR)/liblapacke.so.3 #LAPACKE_LIB= $(LIB_DIR)/liblapack_atlas.a #LAPACKE_LIB=-llapack -llapacke -XML_INCL=-I /usr/include/libxml2 -XML_LIB=-lxml2 +XML2_INCL=-I /usr/include/libxml2 +XML2_LIB=-lxml2 FFTW_INCL= FFTW_LIB= SSL_LIB=/usr/lib64/libssl.so.10 diff --git a/src/Makefile b/src/Makefile index ca38c521..ceb6c25d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -20,7 +20,7 @@ traceBuffer/traceBuffer.a hdf5/hdf5.a activeMQ/gfast_amq.a ifdef GFAST_USE_DMLIB OBJS += dmlib/dmlibWrapper.cpp.o endif -LIBS = $(EW_LIB) $(ISCL_LIB) $(INIPARSER_LIB) $(XML_LIB) $(HDF5_LIB) $(COMPEARTH_LIB) \ +LIBS = $(EW_LIB) $(ISCL_LIB) $(INIPARSER_LIB) $(XML2_LIB) $(HDF5_LIB) $(COMPEARTH_LIB) \ $(ACTIVEMQ_LIB) $(DM_LIB) $(APR_LIB) $(LAPACKE_LIB) $(CBLAS_LIB) $(FFTW_LIB) $(SSL_LIB) \ $(XERCES_LIB) $(QLIB2_LIB) $(CRYPTO_LIB) SYSLIBS = -ldl -lm -lnsl -lrt -pthread -lz -lsz -lstdc++ @@ -69,20 +69,26 @@ install: make -C $$dir install ; \ done +TESTDIR = ./tests +ut: + # make -C $(TESTDIR) ut +test: + # make -C $(TESTDIR) test + clean: - -rm -f *.o $(EXE) $(OBJS) + -rm -f *.o $(EXE) for dir in $(SUBDIRS) ; do \ make -C $$dir clean ; \ done + # make -C $(TESTDIR) clean veryclean: cleandocs for dir in $(SUBDIRS) ; do \ make -C $$dir veryclean ; \ done + # make -C $(TESTDIR) veryclean depend: for dir in $(SUBDIRS) ; do \ make -C $$dir depend ; \ done - -test: # no-op diff --git a/src/activeMQ/Makefile b/src/activeMQ/Makefile index e0dbabc4..6cd92f8c 100644 --- a/src/activeMQ/Makefile +++ b/src/activeMQ/Makefile @@ -6,9 +6,7 @@ MAKEHOST=$(shell uname) endif include $(GFASTROOT)/Make.include.$(MAKEHOST) -OBJS = activeMQ.c.o activeMQ.cpp.o \ -ShakeAlertConsumer.cpp.o ShakeAlertProducer.cpp.o \ -consumer.cpp.o producer.cpp.o readIni.c.o +OBJS = activeMQ.c.o activeMQ.cpp.o readIni.c.o LIB = gfast_amq.a DEBUG = -g #CCFLAGS += -std=c++0x -O0 -D_REENTRANT -Dstatic_config diff --git a/src/core/Makefile b/src/core/Makefile index d554634c..88cca0db 100644 --- a/src/core/Makefile +++ b/src/core/Makefile @@ -44,6 +44,7 @@ clean: for dir in $(SDIRS) ; do \ make -C $$dir clean ; \ done + -rm -f $(LIB) veryclean: cleandocs for dir in $(SDIRS) ; do \ diff --git a/src/core/events/Makefile b/src/core/events/Makefile index 42c5aa63..fbb37c9f 100644 --- a/src/core/events/Makefile +++ b/src/core/events/Makefile @@ -7,8 +7,12 @@ MAKEHOST=$(shell uname) endif include $(GFASTROOT)/Make.include.$(MAKEHOST) -OBJS = freeEvents.o newEvent.o removeCancelledEvent.o removeExpiredEvents.o\ - getMinOriginTime.o printEvent.o removeExpiredEvent.o updateEvent.o syncXMLStatusWithEvents.o +OBJS = freeEvents.o newEvent.o removeExpiredEvents.o\ + printEvent.o removeExpiredEvent.o syncXMLStatusWithEvents.o + +# Remove some files from compile for now since they are not used - CWU 8/18/23 +# OBJS = freeEvents.o newEvent.o removeCancelledEvent.o removeExpiredEvents.o\ +# getMinOriginTime.o printEvent.o removeExpiredEvent.o updateEvent.o syncXMLStatusWithEvents.o DEBUG = -g CCFLAGS += -O0 -D_REENTRANT -Dstatic_config diff --git a/src/core/events/newEvent.c b/src/core/events/newEvent.c index dd12e4b0..30c2c6d0 100644 --- a/src/core/events/newEvent.c +++ b/src/core/events/newEvent.c @@ -97,7 +97,10 @@ bool core_events_newEvent(struct GFAST_shakeAlert_struct SA, // MTH: Same thing for xml output status structs // Create a status output_status from incoming SA eventid strcpy(output_status.eventid, SA.eventid); - output_status.version=-1; /*will increment to 0 on first xml update*/ + output_status.internal_version = 0; /*will decrement to -1 on first xml update*/ + output_status.external_version = -1; /*will increment to 0 on first xml update*/ + memset(&output_status.last_sent_core, 0, sizeof(struct coreInfo_struct)); + // Append this output_status to the incoming xml_status list of records Xtemp.nev = nev0 + 1; Xtemp.SA_status = (struct GFAST_xml_output_status *) calloc((size_t) Xtemp.nev, sizeof(struct GFAST_xml_output_status)); diff --git a/src/core/log/Makefile b/src/core/log/Makefile index 170c3654..524ca472 100644 --- a/src/core/log/Makefile +++ b/src/core/log/Makefile @@ -13,10 +13,10 @@ OBJS = log.c.o endif DEBUG = -g -# CCFLAGS += -O0 -D_REENTRANT -Dstatic_config +CCFLAGS += -O0 -D_REENTRANT -Dstatic_config INCL = -I ../../../include $(ISCL_INCL) \ - -I$(EEWDIR)/third_party \ + -I$(THIRD_PARTY) \ -I$(EEWDIR)/libs/utils SYSLIBS = -lpthread -lnsl -lrt -lz diff --git a/src/core/properties/print.c b/src/core/properties/print.c index e72648e7..caab9dd3 100644 --- a/src/core/properties/print.c +++ b/src/core/properties/print.c @@ -231,6 +231,12 @@ void core_properties_print(struct GFAST_props_struct props) } else { LOG_DEBUGMSG("%s Will allow any E raw sigma values in pd calculations", lspace); } + if (props.pgd_props.q_value_threshold >= 0) { + LOG_DEBUGMSG("%s Will ignore pd values if Q value < %d", + lspace, props.pgd_props.q_value_threshold); + } else { + LOG_DEBUGMSG("%s Will allow any Q values in pd calculations", lspace); + } LOG_DEBUGMSG("%s Will throttle messages below SA mag %f", lspace, props.pgd_props.SA_mag_threshold); LOG_DEBUGMSG("%s Will throttle messages above pgd mag sigma %f", @@ -242,6 +248,42 @@ void core_properties_print(struct GFAST_props_struct props) lspace, props.pgd_props.maximum_pgd_cm); LOG_DEBUGMSG("%s GFAST Maximum stations to include assoc tag: %d", lspace, props.pgd_props.max_assoc_stations); + + // Change thresholds to prevent sending the same (or very similar) solution + bool change_thresholds_not_set = true; + if (props.pgd_props.change_threshold_mag >= 0) { + change_thresholds_not_set = false; + LOG_DEBUGMSG("%s Change threshold, magnitude: %f", + lspace, props.pgd_props.change_threshold_mag); + } + if (props.pgd_props.change_threshold_mag_uncer >= 0) { + change_thresholds_not_set = false; + LOG_DEBUGMSG("%s Change threshold, magnitude uncertainty: %f", + lspace, props.pgd_props.change_threshold_mag_uncer); + } + if (props.pgd_props.change_threshold_lat >= 0) { + change_thresholds_not_set = false; + LOG_DEBUGMSG("%s Change threshold, latitude: %f", + lspace, props.pgd_props.change_threshold_lat); + } + if (props.pgd_props.change_threshold_lon >= 0) { + change_thresholds_not_set = false; + LOG_DEBUGMSG("%s Change threshold, longitude: %f", + lspace, props.pgd_props.change_threshold_lon); + } + if (props.pgd_props.change_threshold_orig_time >= 0) { + change_thresholds_not_set = false; + LOG_DEBUGMSG("%s Change threshold, origin time: %f (s)", + lspace, props.pgd_props.change_threshold_orig_time); + } + if (props.pgd_props.change_threshold_num_stations >= 0) { + change_thresholds_not_set = false; + LOG_DEBUGMSG("%s Change threshold, number of stations: %d", + lspace, props.pgd_props.change_threshold_num_stations); + } + if (change_thresholds_not_set) { + LOG_DEBUGMSG("%s Change thresholds not set!", lspace); + } } else { diff --git a/src/core/scaling/pgd_readIni.c b/src/core/scaling/pgd_readIni.c index 71372e16..e459c78f 100644 --- a/src/core/scaling/pgd_readIni.c +++ b/src/core/scaling/pgd_readIni.c @@ -166,11 +166,13 @@ int core_scaling_pgd_readIni(const char *propfilename, ierr = core_scaling_readSigmaLookupFile(sigmaLookupFile, pgd_props); + // only use PGD observations if Q value >= this threshold + setVarName(group, "q_value_threshold\0", var); + pgd_props->q_value_threshold = iniparser_getint(ini, var, -1); // only send XML for pgd magnitude sigma below this threshold setVarName(group, "pgd_sigma_throttle\0", var); pgd_props->pgd_sigma_throttle = iniparser_getdouble(ini, var, 10); - // props->pgd_sigma_throttle = iniparser_getdouble(ini, "general:pgd_sigma_throttle\0", 10); if (pgd_props->pgd_sigma_throttle <= 0) { LOG_ERRMSG("Error pgd_sigma_throttle must be positive: %f", @@ -181,7 +183,6 @@ int core_scaling_pgd_readIni(const char *propfilename, // only send XML for SA magnitude above this threshold setVarName(group, "SA_mag_threshold\0", var); pgd_props->SA_mag_threshold = iniparser_getdouble(ini, var, -10.0); - // props->SA_mag_threshold = iniparser_getdouble(ini, "general:SA_mag_threshold\0", -10.0); setVarName(group, "minimum_pgd_cm\0", var); pgd_props->minimum_pgd_cm = iniparser_getdouble(ini, var, -1.0); @@ -201,6 +202,20 @@ int core_scaling_pgd_readIni(const char *propfilename, goto ERROR; } + // Change thresholds to prevent sending the same (or very similar) message + setVarName(group, "change_threshold_mag\0", var); + pgd_props->change_threshold_mag = iniparser_getdouble(ini, var, -1.0); + setVarName(group, "change_threshold_mag_uncer\0", var); + pgd_props->change_threshold_mag_uncer = iniparser_getdouble(ini, var, -1.0); + setVarName(group, "change_threshold_lat\0", var); + pgd_props->change_threshold_lat = iniparser_getdouble(ini, var, -1.0); + setVarName(group, "change_threshold_lon\0", var); + pgd_props->change_threshold_lon = iniparser_getdouble(ini, var, -1.0); + setVarName(group, "change_threshold_orig_time\0", var); + pgd_props->change_threshold_orig_time = iniparser_getdouble(ini, var, -1.0); + setVarName(group, "change_threshold_num_stations\0", var); + pgd_props->change_threshold_num_stations = iniparser_getint(ini, var, -1); + ERROR:; iniparser_freedict(ini); return ierr; diff --git a/src/core/waveformProcessor/parseQChannel.c b/src/core/waveformProcessor/parseQChannel.c index 94e90a1d..083ac4d3 100644 --- a/src/core/waveformProcessor/parseQChannel.c +++ b/src/core/waveformProcessor/parseQChannel.c @@ -65,4 +65,36 @@ int core_waveformProcessor_parseQChannelChi2CWUmap( } return chimap; +} + + +/*! + * @brief Unpack the quality channel. The most recent definition as of 3/22/23 was: + * Q = 100000*int(10*PDOP) + 1000*nsat + 10*fixtyp + goodness + * + * @param[in] q_value Q channel value. It is a double to stay consistent + * with other data buffers. But treated as integer here. + * + * @result the goodness value (0 or 1) - a combination of chi^2, nsats, clocks, etc. + * + * @author Brendan Crowell (PNSN) and Carl Ulberg (PNSN) + * + * @date Mar 2023 + * + */ +int core_waveformProcessor_parseQChannelGoodness( + const double q_value) +{ + int goodness; + + // Extract the last 1 digit, should only be 0 or 1, so modulo 2 should work + if (q_value < 0) { + LOG_WARNMSG("parseQChannelChi2CWU - q_value is negative (%lf), changing to positive!", + q_value); + goodness = (int)(q_value * -1) % 2; + } else { + goodness = (int)(q_value) % 2; + } + + return goodness; } \ No newline at end of file diff --git a/src/core/waveformProcessor/peakDisplacement.c b/src/core/waveformProcessor/peakDisplacement.c index 718bd7b4..0774a946 100644 --- a/src/core/waveformProcessor/peakDisplacement.c +++ b/src/core/waveformProcessor/peakDisplacement.c @@ -47,7 +47,8 @@ int core_waveformProcessor_peakDisplacement( { double currentTime, distance, effectiveHypoDist, epoch, peakDisp, x1, x2, y1, y2, tmin, tmax, obsTime, - uMaxUncertainty, nMaxUncertainty, eMaxUncertainty, qMax, qRef, qPeak; + uMaxUncertainty, nMaxUncertainty, eMaxUncertainty; + int qMin, qRef, qPeak; int k, nsites, zone_loc, iRef, iPeak; //unused int i; bool lnorthp; @@ -78,6 +79,20 @@ int core_waveformProcessor_peakDisplacement( const double u_raw_sigma_threshold = pgd_props->u_raw_sigma_threshold; const double n_raw_sigma_threshold = pgd_props->n_raw_sigma_threshold; const double e_raw_sigma_threshold = pgd_props->e_raw_sigma_threshold; + // Threshold value for Q channel. If the observed Q value for the reference + // or peak displacement is <= this threshold, ignore the pd observation. + // If q_value_threshold is < 0, allow any Q value + const int q_value_threshold = pgd_props->q_value_threshold; + + LOG_MSG("Calculating PD, svel=(%.1f/%.1f), pgd_bounds=(%.1f/%.1f), sigma_thresh=(%.1f/%.1f/%.1f), q_thresh=%d", + svel_window, + min_svel_window, + min_pgd_cm, + max_pgd_cm, + u_raw_sigma_threshold, + n_raw_sigma_threshold, + e_raw_sigma_threshold, + q_value_threshold); //------------------------------------------------------------------------// // @@ -90,7 +105,9 @@ int core_waveformProcessor_peakDisplacement( uMaxUncertainty = 0.0; eMaxUncertainty = 0.0; nMaxUncertainty = 0.0; - qMax = 0.0; + qMin = 0; + qRef = 0; + qPeak = 0; if (gps_data.stream_length != pgd_data->nsites) { LOG_ERRMSG("Inconsistent structure sizes %d %d", @@ -125,7 +142,9 @@ int core_waveformProcessor_peakDisplacement( uMaxUncertainty = 0.0; eMaxUncertainty = 0.0; nMaxUncertainty = 0.0; - qMax = 0.0; + qMin = 0; + qRef = 0; + qPeak = 0; // Make sure I have the latest/greatest site location pgd_data->sta_lat[k] = gps_data.data[k].sta_lat; pgd_data->sta_lon[k] = gps_data.data[k].sta_lon; @@ -187,7 +206,6 @@ int core_waveformProcessor_peakDisplacement( gps_data.data[k].netw, gps_data.data[k].loc); */ - //if (distance < effectiveHypoDist) if (distance < effectiveHypoDist) { // Compute the peak displacement max(norm(u + n + e, 2)) @@ -254,12 +272,12 @@ M 9 at 100km: -6.687 + 150 - 21.4*2 = 100 cm(?) eMaxUncertainty = fmax(gps_data.data[k].esigmabuff[iRef], gps_data.data[k].esigmabuff[iPeak]); } if (!isnan(gps_data.data[k].qbuff[iRef]) && !isnan(gps_data.data[k].qbuff[iPeak])) { - qRef = core_waveformProcessor_parseQChannelChi2CWU(gps_data.data[k].qbuff[iRef]); - qPeak = core_waveformProcessor_parseQChannelChi2CWU(gps_data.data[k].qbuff[iPeak]); - qMax = fmax(qRef, qPeak); + qRef = core_waveformProcessor_parseQChannelGoodness(gps_data.data[k].qbuff[iRef]); + qPeak = core_waveformProcessor_parseQChannelGoodness(gps_data.data[k].qbuff[iPeak]); + qMin = fmin(qRef, qPeak); } - LOG_MSG("%s.%s.%s.%s peakDisp=%f dist=%.2f, peakSigmas=(%.4f,%.4f,%.4f), peakQ=%.4f", + LOG_MSG("%s.%s.%s.%s peakDisp=%f dist=%.2f, peakSigmas=(%.4f,%.4f,%.4f), minQ=%d", gps_data.data[k].stnm, gps_data.data[k].chan[0], gps_data.data[k].netw, gps_data.data[k].loc, peakDisp, @@ -267,7 +285,7 @@ M 9 at 100km: -6.687 + 150 - 21.4*2 = 100 cm(?) uMaxUncertainty, nMaxUncertainty, eMaxUncertainty, - qMax); + qMin); } // Is the observation above the defined minimum? @@ -313,6 +331,18 @@ M 9 at 100km: -6.687 + 150 - 21.4*2 = 100 cm(?) ); } + // Is the observation under the Q value threshold? + if ((q_value_threshold >= 0) && (qMin < q_value_threshold)) { + l_use_observation = false; + LOG_DEBUGMSG("CCC PeakDisp, %s.%s.%s.%s ignoring observation, qRef/qPeak: (%d/%d), threshold: %d", + gps_data.data[k].stnm, gps_data.data[k].chan[0], + gps_data.data[k].netw, gps_data.data[k].loc, + qRef, + qPeak, + q_value_threshold + ); + } + // If it isn't a NaN and within the sanity bounds then retain it for processing if (l_use_observation) { pgd_data->pd_time[k] = obsTime; // epoch diff --git a/src/core/waveformProcessor/peakDisplacementHelper.c b/src/core/waveformProcessor/peakDisplacementHelper.c index 0633c77a..6b675686 100644 --- a/src/core/waveformProcessor/peakDisplacementHelper.c +++ b/src/core/waveformProcessor/peakDisplacementHelper.c @@ -14,7 +14,7 @@ /*! * @brief Waveform processor to estimate the peak displacement observed * on a 3 channel GPS stream where the peak displacement at any - * sample is Euclidean norm of it's displacement. + * sample is Euclidean norm of its displacement. * * @param[in] npts number of points in time series * @param[in] dt sampling period (s) of GPS buffers @@ -23,21 +23,14 @@ * @param[in] ubuff vertical position [npts] * @param[in] nbuff north position [npts] * @param[in] ebuff east position [npts] - * @param[in] usigmabuff vertical position uncertainty [npts] - * @param[in] nsigmabuff north position uncertainty [npts] - * @param[in] esigmabuff east position uncertainty [npts] * @param[in] nMaxLeader latest time index to start measurement in * case of nan's * @param[in] tmin distance / svel_window (s) * @param[in] tmax distance / min_svel_window (s) * * @param[out] obsTime epochal time of peak displacment (s) - * @param[out] uMaxUncertainty peak vertical uncertainty at the reference - * or peak displacement times (m) - * @param[out] nMaxUncertainty peak north uncertainty at the reference - * or peak displacement times (m) - * @param[out] eMaxUncertainty peak east uncertainty at the reference - * or peak displacement times (m) + * @param[out] iRef index of reference value + * @param[out] iPeak index of peak value * * * @result the peak displacement observed on a trace. this has the same diff --git a/src/dmlib/dmlibWrapper.cpp b/src/dmlib/dmlibWrapper.cpp index 09c39bc7..c6911885 100644 --- a/src/dmlib/dmlibWrapper.cpp +++ b/src/dmlib/dmlibWrapper.cpp @@ -470,9 +470,8 @@ char *dmlibWrapper_createPGDXML(const double currentTime, // required to create FiniteFaultMessage enum FaultSegment::FaultSegmentShape shape = FaultSegment::UNKNOWN_SEGMENT; - // Get time stamp for when message is sent. Corresponds to currentTime from - // the outer driveGFAST loop to be consistent with the last time we have - // data for. + // Get time stamp for when the message is sent, based on the timestamp sent from the outer + // driveGFAST loop. Could be based on the actual current time, or the last time we have data for. int rc; char cnow[128]; rc = xml_epoch2string(currentTime, cnow); @@ -500,8 +499,6 @@ char *dmlibWrapper_createPGDXML(const double currentTime, LOG_DEBUGMSG("old origTime: %lf, origTimeChar: %s, origTimeStr: %s, new origTime: %lf", core->origTime, origTimeChar, origTimeStr.c_str(), algMessage.getOriginTime()); - // LOG_MSG("%s", "createEventXML - created algMessage"); - // Now add pgd observations to algMessage int i, j; int scnl_n = 8; @@ -544,7 +541,6 @@ char *dmlibWrapper_createPGDXML(const double currentTime, // skip site if it wasn't used if (!pgd->lsiteUsed[i]) { continue; } - // LOG_MSG("createEventXML - i = %d, setting chars", i); // see core/data/readMetaDataFile for similar SNCL parsing memset(obs_sta, 0, scnl_n*sizeof(char)); memset(obs_net, 0, scnl_n*sizeof(char)); @@ -554,7 +550,6 @@ char *dmlibWrapper_createPGDXML(const double currentTime, work = (char *)calloc(strlen(pgd_data->stnm[i])+1, sizeof(char)); strcpy(work, pgd_data->stnm[i]); - // LOG_MSG("%s", "createEventXML - starting NSCL tokenizing"); token = strtok(work, "."); int i_tok = 0; while (token) @@ -566,10 +561,8 @@ char *dmlibWrapper_createPGDXML(const double currentTime, i_tok++; token = strtok(NULL, "."); } - // LOG_MSG("%s", "createEventXML - done NSCL tokenizing, freeing work"); delete work; work = NULL; - // LOG_MSG("%s", "createEventXML - freed work, adding gmobs"); assoc_flag = (n_assoc >= max_assoc_stations) ? false : true; // Make sure longitude follows ShakeAlert convention @@ -594,9 +587,6 @@ char *dmlibWrapper_createPGDXML(const double currentTime, assoc_flag); n_assoc++; - - // LOG_DEBUGMSG("createEventXML - added obs=%f, sta=%s, i=%d, j=%d, dist=%f, assoc=%d, n_assoc=%d", - // pgd_data->pd[i] * 100., pgd_data->stnm[i], i, j, vals[j].dist, assoc_flag, n_assoc); } LOG_MSG("%s", "createEventXML - finished adding gmobs, encoding message"); diff --git a/src/doxygen.conf b/src/doxygen.conf index 40cc57b5..36f42eb2 100644 --- a/src/doxygen.conf +++ b/src/doxygen.conf @@ -14,7 +14,7 @@ LAYOUT_FILE = ../../doxygen.layout PROJECT_NAME = "GFAST" # version number -PROJECT_NUMBER = "GFAST: gfast-1.2.1-2022-11-28" +PROJECT_NUMBER = "GFAST: gfast-1.2.4-2023-10-12" # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/src/eewUtils/Makefile b/src/eewUtils/Makefile index 9f09e25b..a318a2ae 100644 --- a/src/eewUtils/Makefile +++ b/src/eewUtils/Makefile @@ -1,21 +1,23 @@ # Makefile for GFAST/src/eewUtils # # + GFASTROOT=.. ifndef MAKEHOST MAKEHOST=$(shell uname) endif include $(GFASTROOT)/Make.include.$(MAKEHOST) -OBJS = driveCMT.o driveFF.o drivePGD.o makeXML.o\ -parseCoreXML.o setLogFileNames.o driveGFAST.o +OBJS = driveCMT.o driveFF.o drivePGD.o makeXML.o \ +parseCoreXML.o setLogFileNames.o driveGFAST.o fillCoreEventInfo.o \ +sendXMLFilter.o LIB = gfast_utils.a DEBUG = -g -CCFLAGS += -O0 -D_REENTRANT -Dstatic_config ${GFAST_VERSION} +CFLAGS += -O0 -D_REENTRANT -Dstatic_config ${GFAST_VERSION} #iniparser will be deleted -INCL = -I ../../include $(ISCL_INCL) $(XML_INCL) $(COMPEARTH_INCL) $(EW_INCL) $(HDF5_INCL) +INCL = -I ../../include $(ISCL_INCL) $(XML2_INCL) $(COMPEARTH_INCL) $(EW_INCL) $(HDF5_INCL) SYSLIBS = -lpthread -lnsl -lrt -lz @@ -46,7 +48,7 @@ cleandepend: rm -rf $(DEPENDFILES) clean: - -rm -f *.o $(BIN) + -rm -f *.o $(BIN) $(LIB) veryclean: clean rm-ids cleandepend diff --git a/src/eewUtils/driveGFAST.c b/src/eewUtils/driveGFAST.c index 1b68f334..c1287d5d 100644 --- a/src/eewUtils/driveGFAST.c +++ b/src/eewUtils/driveGFAST.c @@ -24,7 +24,6 @@ * @param[in] interval_in_mins is interval in minutes? * @return status code. */ - int eewUtils_writeXML(const char *dirname, const char *eventid, const char *msg_type, const char *message, int interval, bool interval_in_mins); bool check_mins_against_intervals(struct GFAST_props_struct props, @@ -35,36 +34,6 @@ bool check_mins_against_intervals(struct GFAST_props_struct props, bool * interval_complete, double age); -/*! - * @brief Fills a given coreInfo_struct with the appropriate information - * @param[in] evid Event ID - * @param[in] version Event version number - * @param[in] SA_lat Event latitude - * @param[in] SA_lon Event longitude - * @param[in] SA_depth Event depth - * @param[in] SA_mag Event magnitude - * @param[in] SA_time Event origin time (UTC) - * @param[in] num_stations Number of stations contributing - * @param[out] core struct to fill with information - * @return status code. - */ -int fill_core_event_info(const char *evid, - const int version, - const double SA_lat, - const double SA_lon, - const double SA_depth, - const double SA_mag, - const double SA_time, - const int num_stations, - struct coreInfo_struct *core); - -bool send_xml_filter(const struct GFAST_props_struct *props, - const struct GFAST_shakeAlert_struct *SA, - const struct GFAST_pgdResults_struct *pgd, - const struct GFAST_peakDisplacementData_struct *pgd_data, - const struct coreInfo_struct *core, - const double age_of_event); - /*! * @brief Expert earthquake early warning GFAST driver. * @@ -158,11 +127,17 @@ int eewUtils_driveGFAST(const double currentTime, t_event = time_timeStamp(); // Get the streams for this event memcpy(&SA, &events->SA[iev], sizeof(struct GFAST_shakeAlert_struct)); + + // internal_version always decrements by 1 + xml_status->SA_status[iev].internal_version -= 1; + t1 = SA.time; // Origin time t2 = currentTime; age_of_event = (t2 - t1); - LOG_MSG("%s: Starting event time:%lf evid:%s [age_of_event=%f]", - fcnm, t2, SA.eventid, age_of_event); + LOG_MSG("%s: Starting event, current time:%lf evid:%s internal version: %d [age_of_event=%lf]", + fcnm, t2, SA.eventid, xml_status->SA_status[iev].internal_version, age_of_event); + LOG_MSG("%s: SA event info OT:%lf lat:%lf lon:%lf mag:%lf SA version:%d", + fcnm, SA.time, SA.lat, SA.lon, SA.mag, SA.version); // Skip event if the times doen't make sense if (t1 > t2) { @@ -188,9 +163,8 @@ int eewUtils_driveGFAST(const double currentTime, continue; } - // Set the log file names. Comment out in advance of providing compile - // option to use these (for NOAA) or plog (for ShakeAlert) - CWU #ifndef ENABLE_PLOG + // Set the log file names. eewUtils_setLogFileNames(SA.eventid,props.SAoutputDir, errorLogFileName, infoLogFileName, debugLogFileName, warnLogFileName); @@ -396,24 +370,29 @@ int eewUtils_driveGFAST(const double currentTime, // Make xml messages /////////////////////////////////////////////////////////////////////////// t_time0 = time_timeStamp(); + LOG_MSG("driveGFAST: make XML msgs: lpgdSuccess=%d lcmtSuccess=%d lffSuccess=%d", lpgdSuccess, lcmtSuccess, lffSuccess); - xml_status->SA_status[iev].version += 1; - char *message_type = (xml_status->SA_status[iev].version==0)?"new\0":"update\0"; + char sversion[6]; - snprintf(sversion,6,"%d",xml_status->SA_status[iev].version); + char message_type[16]; + snprintf(sversion, 6, "%d", xml_status->SA_status[iev].internal_version); + sprintf(message_type, "%s", "update\0"); // Fill coreInfo_struct to pass to makeXML for pgd and ff struct coreInfo_struct core; memset(&core, 0, sizeof(struct coreInfo_struct)); - ierr = fill_core_event_info(SA.eventid, xml_status->SA_status[iev].version, SA.lat, - SA.lon, SA.dep, SA.mag, SA.time, 0, &core); + ierr = eewUtils_fillCoreEventInfo(SA.eventid, xml_status->SA_status[iev].internal_version, + SA.lat, SA.lon, SA.dep, SA.mag, SA.time, 0, &core); - // Make the PGD xml + // Handle the PGD xml if (props.pgd_props.do_pgd && lpgdSuccess) { if (props.verbose > 2) { LOG_DEBUGMSG("%s", "Generating pgd XML"); } + + bool send_pgd = false; + // Change depth, mag to match optimal pgd (by variance reduction) pgdOpt = array_argmax64f(pgd->ndeps, pgd->dep_vr_pgd, &ierr); core.depth = pgd->srcDepths[pgdOpt]; @@ -421,11 +400,42 @@ int eewUtils_driveGFAST(const double currentTime, core.magUncer = pgd->mpgd_sigma[pgdOpt]; core.numStations = nsites_pgd; +#ifdef GFAST_USE_AMQ + // Determine if message should be sent via ActiveMQ + if (!eewUtils_sendXMLFilter( + &props, + &SA, + pgd, + pgd_data, + &core, + &(xml_status->SA_status[iev].last_sent_core), + age_of_event)) + { + // Update version, message type if message should be sent + send_pgd = true; + xml_status->SA_status[iev].external_version += 1; + core.version = xml_status->SA_status[iev].external_version; + snprintf(sversion, 6, "%d", xml_status->SA_status[iev].external_version); + if (xml_status->SA_status[iev].external_version == 0) { + sprintf(message_type, "%s", "new\0"); + } + LOG_MSG("driveGFAST: PGD message will be sent! Internal version: %d, external version: %d", + xml_status->SA_status[iev].internal_version, + xml_status->SA_status[iev].external_version) + } +#endif /* GFAST_USE_AMQ */ + + // Make PGD xml #ifdef GFAST_USE_DMLIB - // Encode xml with dmlib + // Encode xml with dmlib LOG_MSG("%s", "driveGFAST: CWU_TEST dmlib encoding"); + // Use time closest to when xml is actually sent to be consistent with ShakeAlert + // algorithms. + // Alternative is to use the "currentTime" variable to correspond to the start of the + // processing loop in gfast_eew. + double timestamp = time_timeStamp(); pgdXML = dmlibWrapper_createPGDXML( - currentTime, + timestamp, props.opmode, GFAST_VERSION, program_instance, @@ -461,16 +471,19 @@ int eewUtils_driveGFAST(const double currentTime, SA.eventid, iev, xml_status->SA_status[iev].eventid); } -#if defined GFAST_USE_AMQ && defined GFAST_USE_DMLIB - // Send message via ActiveMQ if appropriate - if (!send_xml_filter(&props, &SA, pgd, pgd_data, &core, age_of_event)) { +#ifdef GFAST_USE_AMQ + if (send_pgd) { if (pgdXML != NULL) { sendEventXML(pgdXML); + memcpy( + &(xml_status->SA_status[iev].last_sent_core), + &core, + sizeof(core)); } - LOG_MSG("== Sending xml, [GFAST t0:%f] evid:%s pgdXML=[%s]\n", - currentTime, SA.eventid, pgdXML); + LOG_MSG("== Sending xml, [GFAST t0:%f] evid:%s version:%d pgdXML=[%s]\n", + currentTime, SA.eventid, core.version, pgdXML); } -#endif /* GFAST_USE_AMQ && GFAST_USE_DMLIB */ +#endif /* GFAST_USE_AMQ */ if (props.output_interval_mins[0] == 0) { // Output at every iteration int index = core.version; @@ -483,10 +496,15 @@ int eewUtils_driveGFAST(const double currentTime, check_mins_against_intervals(props, mins, SA.eventid, "pgd", pgdXML, xml_status->SA_status[iev].interval_complete[0], age_of_event); } + + // reset version, message type to reflect internal state for cmt, ff + core.version = xml_status->SA_status[iev].internal_version; + snprintf(sversion, 6, "%d", xml_status->SA_status[iev].internal_version); + sprintf(message_type, "%s", "update\0"); LOG_MSG("Leaving PGD writeXML ierr=%d", ierr); } //if props.pgd_props.do_pgd && lpgdSuccess - // Make the CMT quakeML + // Handle the CMT quakeML if (props.cmt_props.do_cmt && lcmtSuccess) { if (props.verbose > 2) { LOG_DEBUGMSG("%s", "Generating CMT QuakeML"); @@ -522,7 +540,7 @@ int eewUtils_driveGFAST(const double currentTime, xml_status->SA_status[iev].interval_complete[1], age_of_event); } } // if props.cmt_props.do_cmt && lcmtSuccess - // Make the finite fault XML + // Handle the finite fault XML if (props.ff_props.do_ff && lffSuccess) { if (props.verbose > 2) { LOG_DEBUGMSG("Generating FF XML; preferred plane=%d", @@ -761,160 +779,3 @@ bool check_mins_against_intervals(struct GFAST_props_struct props, return false; } - -int fill_core_event_info(const char *evid, - const int version, - const double SA_lat, - const double SA_lon, - const double SA_depth, - const double SA_mag, - const double SA_time, - const int num_stations, - struct coreInfo_struct *core) -{ - strcpy(core->id, evid); - core->version = version; - core->mag = SA_mag; - core->lhaveMag = true; - core->magUnits = MOMENT_MAGNITUDE; - core->lhaveMagUnits = true; - core->magUncer = 0.5; - core->lhaveMagUncer = true; - core->magUncerUnits = MOMENT_MAGNITUDE; - core->lhaveMagUncerUnits = true; - core->lat = SA_lat; - core->lhaveLat = true; - core->latUnits = DEGREES; - core->lhaveLatUnits = true; - core->latUncer = (double) NAN; - core->lhaveLatUncer = true; - core->latUncerUnits = DEGREES; - core->lhaveLatUncerUnits = true; - // GFAST would call lon -120 as 240 by default. Change this to be - // consistent with ShakeAlert seismic algorithms - core->lon = (SA_lon > 180) ? SA_lon - 360: SA_lon; - core->lhaveLon = true; - core->lonUnits = DEGREES; - core->lhaveLonUnits = true; - core->lonUncer = (double) NAN; - core->lhaveLonUncer = true; - core->lonUncerUnits = DEGREES; - core->lhaveLonUncerUnits = true; - core->depth = SA_depth; - core->lhaveDepth = true; - core->depthUnits = KILOMETERS; - core->lhaveDepthUnits = true; - core->depthUncer = (double) NAN; - core->lhaveDepthUncer = true; - core->depthUncerUnits = KILOMETERS; - core->lhaveDepthUncerUnits = true; - core->origTime = SA_time; - core->lhaveOrigTime = true; - core->origTimeUnits = UTC; - core->lhaveOrigTimeUnits = true; - core->origTimeUncer = (double) NAN; - core->lhaveOrigTimeUncer = true; - core->origTimeUncerUnits = SECONDS; - core->lhaveOrigTimeUncerUnits = true; - core->likelihood = 0.8; - core->lhaveLikelihood = true; - core->numStations = num_stations; - return 0; -} - -/* - * Return true if this message should not be sent (false if it should be sent) - */ -bool send_xml_filter(const struct GFAST_props_struct *props, - const struct GFAST_shakeAlert_struct *SA, - const struct GFAST_pgdResults_struct *pgd, - const struct GFAST_peakDisplacementData_struct *pgd_data, - const struct coreInfo_struct *core, - const double age_of_event) -{ - - // Conditions can be met or exceeded. - // Exceeded is used if a value being more than a threshold means no throttling - // Met is used if a value being less than a threshold means no throttling - bool pgd_exceeded = false; - bool mag_exceeded = false; - bool mag_sigma_met = false; - // Determine if pgd threshold is exceeded n times - int num_pgd_exceeded = 0, i, i_throttle; - - // Find the correct throttle criteria based on the time after origin - i_throttle = -1; - for (i = 0; i < props->pgd_props.n_throttle; i++) { - if (props->pgd_props.throttle_time_threshold[i] > age_of_event) break; - i_throttle++; - } - i_throttle = (i_throttle < 0) ? 0: i_throttle; - if (props->verbose > 2) { - LOG_DEBUGMSG("%s: For age_of_event %.2f, using i_throttle=%d, thresholds for time=%.1f, pgd=%.1f, nsta=%d", - __func__, - age_of_event, - i_throttle, - props->pgd_props.throttle_time_threshold[i_throttle], - props->pgd_props.throttle_pgd_threshold[i_throttle], - props->pgd_props.throttle_num_stations[i_throttle]) - } - - // Assumes pgd->nsites = pgd_data->nsites and indices correspond - if (pgd->nsites != pgd_data->nsites) { - LOG_ERRMSG("%s: nsites don't match for pgd, pgd_data! %d, %d\n", - __func__, pgd->nsites, pgd_data->nsites); - return false; - } - for (i = 0; i < pgd->nsites; i++) { - // skip site if it wasn't used - if (!pgd->lsiteUsed[i]) { continue; } - // pd is in meters, so convert to cm before comparing to threshold - if (pgd_data->pd[i] * 100. > props->pgd_props.throttle_pgd_threshold[i_throttle]) { - num_pgd_exceeded++; - } - } - - if (props->verbose > 2) { - LOG_DEBUGMSG("%s: PGD threshold of %.1f cm exceeded at %d stations (threshold num: %d)", - __func__, props->pgd_props.throttle_pgd_threshold[i_throttle], num_pgd_exceeded, - props->pgd_props.throttle_num_stations[i_throttle]); - } - if (num_pgd_exceeded >= props->pgd_props.throttle_num_stations[i_throttle]) { - LOG_MSG("%s: PGD threshold of %.1f cm exceeded at %d stations (threshold num: %d)", - __func__, props->pgd_props.throttle_pgd_threshold[i_throttle], num_pgd_exceeded, - props->pgd_props.throttle_num_stations[i_throttle]); - pgd_exceeded = true; - } - - // Determine if SA mag threshold is exceeded - if (props->verbose > 2) { - LOG_DEBUGMSG("%s: SA mag: %f, threshold mag: %f", - __func__, SA->mag, props->pgd_props.SA_mag_threshold); - } - if (SA->mag >= props->pgd_props.SA_mag_threshold) { - mag_exceeded = true; - LOG_MSG("%s: SA magnitude exceeded! SA mag: %f, threshold mag: %f", - __func__, SA->mag, props->pgd_props.SA_mag_threshold); - } - - // Determine if pgd magnitude sigma threshold is met - if (props->verbose > 2) { - LOG_DEBUGMSG("%s: PGD mag sigma: %f, threshold mag sigma: %f, PGD mag: %f", - __func__, core->magUncer, props->pgd_props.pgd_sigma_throttle, core->mag); - } - if (core->magUncer < props->pgd_props.pgd_sigma_throttle) { - mag_sigma_met = true; - LOG_MSG("%s: PGD mag sigma met! PGD mag sigma: %f, threshold mag sigma: %f, PGD mag: %f", - __func__, core->magUncer, props->pgd_props.pgd_sigma_throttle, core->mag); - } - - if (pgd_exceeded && mag_exceeded && mag_sigma_met) { - LOG_MSG("%s: Message not throttled (%s v%d)! pgd_exceeded: %d, mag_exceeded: %d, mag_sigma_met: %d", - __func__, core->id, core->version, pgd_exceeded, mag_exceeded, mag_sigma_met); - return false; - } - - LOG_MSG("%s: Message throttled (%s v%d)! pgd_exceeded: %d, mag_exceeded: %d, mag_sigma_met: %d", - __func__, core->id, core->version, pgd_exceeded, mag_exceeded, mag_sigma_met); - return true; -} diff --git a/src/eewUtils/fillCoreEventInfo.c b/src/eewUtils/fillCoreEventInfo.c new file mode 100644 index 00000000..a9918a4d --- /dev/null +++ b/src/eewUtils/fillCoreEventInfo.c @@ -0,0 +1,77 @@ +#include +#include +#include "gfast_struct.h" + +/*! + * @brief Fills a given coreInfo_struct with the appropriate information + * @param[in] evid Event ID + * @param[in] version Event version number + * @param[in] SA_lat Event latitude + * @param[in] SA_lon Event longitude + * @param[in] SA_depth Event depth + * @param[in] SA_mag Event magnitude + * @param[in] SA_time Event origin time (UTC) + * @param[in] num_stations Number of stations contributing + * @param[out] core struct to fill with information + * @return status code. + */ +int eewUtils_fillCoreEventInfo( + const char *evid, + const int version, + const double SA_lat, + const double SA_lon, + const double SA_depth, + const double SA_mag, + const double SA_time, + const int num_stations, + struct coreInfo_struct *core) +{ + strcpy(core->id, evid); + core->version = version; + core->mag = SA_mag; + core->lhaveMag = true; + core->magUnits = MOMENT_MAGNITUDE; + core->lhaveMagUnits = true; + core->magUncer = 0.5; + core->lhaveMagUncer = true; + core->magUncerUnits = MOMENT_MAGNITUDE; + core->lhaveMagUncerUnits = true; + core->lat = SA_lat; + core->lhaveLat = true; + core->latUnits = DEGREES; + core->lhaveLatUnits = true; + core->latUncer = (double) NAN; + core->lhaveLatUncer = true; + core->latUncerUnits = DEGREES; + core->lhaveLatUncerUnits = true; + // GFAST would call lon -120 as 240 by default. Change this to be + // consistent with ShakeAlert seismic algorithms + core->lon = (SA_lon > 180) ? SA_lon - 360: SA_lon; + core->lhaveLon = true; + core->lonUnits = DEGREES; + core->lhaveLonUnits = true; + core->lonUncer = (double) NAN; + core->lhaveLonUncer = true; + core->lonUncerUnits = DEGREES; + core->lhaveLonUncerUnits = true; + core->depth = SA_depth; + core->lhaveDepth = true; + core->depthUnits = KILOMETERS; + core->lhaveDepthUnits = true; + core->depthUncer = (double) NAN; + core->lhaveDepthUncer = true; + core->depthUncerUnits = KILOMETERS; + core->lhaveDepthUncerUnits = true; + core->origTime = SA_time; + core->lhaveOrigTime = true; + core->origTimeUnits = UTC; + core->lhaveOrigTimeUnits = true; + core->origTimeUncer = (double) NAN; + core->lhaveOrigTimeUncer = true; + core->origTimeUncerUnits = SECONDS; + core->lhaveOrigTimeUncerUnits = true; + core->likelihood = 0.8; + core->lhaveLikelihood = true; + core->numStations = num_stations; + return 0; +} \ No newline at end of file diff --git a/src/eewUtils/sendXMLFilter.c b/src/eewUtils/sendXMLFilter.c new file mode 100644 index 00000000..769cd9c1 --- /dev/null +++ b/src/eewUtils/sendXMLFilter.c @@ -0,0 +1,220 @@ +#include +#include +#include "gfast_struct.h" +#include "gfast_core.h" + + +/* + * Return true if the change thresholds are met or exceeded. + * For each change threshold, the threshold value must be >= 0 to be evaluated + * In other words, if the threshold is negative, it won't be used + * If any single threshold is met or exceeded, return true + */ +bool eewUtils_changeThresholdsExceeded( + const struct GFAST_props_struct *props, + const struct coreInfo_struct *core, + const struct coreInfo_struct *last_sent_core) +{ + // If all of the thresholds are unset, just ignore the whole check and return true. + if ((props->pgd_props.change_threshold_mag < 0) && + (props->pgd_props.change_threshold_mag_uncer < 0) && + (props->pgd_props.change_threshold_lat < 0) && + (props->pgd_props.change_threshold_lon < 0) && + (props->pgd_props.change_threshold_orig_time < 0) && + (props->pgd_props.change_threshold_num_stations < 0)) + { + LOG_DEBUGMSG("%s: All change thresholds unset, not checking!", __func__); + return true; + } + + // Check magnitude + if ((props->pgd_props.change_threshold_mag >= 0) && + (fabs(core->mag - last_sent_core->mag) >= props->pgd_props.change_threshold_mag)) + { + LOG_DEBUGMSG("%s: Returning true for mag: %f >= %f", + __func__, + fabs(core->mag - last_sent_core->mag), + props->pgd_props.change_threshold_mag); + return true; + } + // Check magnitude uncertainty + if ((props->pgd_props.change_threshold_mag_uncer >= 0) && + (fabs(core->magUncer - last_sent_core->magUncer) >= props->pgd_props.change_threshold_mag_uncer)) + { + LOG_DEBUGMSG("%s: Returning true for mag uncer: %f >= %f", + __func__, + fabs(core->magUncer - last_sent_core->magUncer), + props->pgd_props.change_threshold_mag_uncer); + return true; + } + // Check latitude + if ((props->pgd_props.change_threshold_lat >= 0) && + (fabs(core->lat - last_sent_core->lat) >= props->pgd_props.change_threshold_lat)) + { + LOG_DEBUGMSG("%s: Returning true for lat: %f >= %f", + __func__, + fabs(core->lat - last_sent_core->lat), + props->pgd_props.change_threshold_lat); + return true; + } + // Check longitude + if ((props->pgd_props.change_threshold_lon >= 0) && + (fabs(core->lon - last_sent_core->lon) >= props->pgd_props.change_threshold_lon)) + { + LOG_DEBUGMSG("%s: Returning true for lon: %f >= %f", + __func__, + fabs(core->lon - last_sent_core->lon), + props->pgd_props.change_threshold_lon); + return true; + } + // Check origin time + if ((props->pgd_props.change_threshold_orig_time >= 0) && + (fabs(core->origTime - last_sent_core->origTime) >= props->pgd_props.change_threshold_orig_time)) + { + LOG_DEBUGMSG("%s: Returning true for ot: %f >= %f", + __func__, + fabs(core->origTime - last_sent_core->origTime), + props->pgd_props.change_threshold_orig_time); + return true; + } + // Check number of stations + if ((props->pgd_props.change_threshold_num_stations >= 0) && + (abs(core->numStations - last_sent_core->numStations) >= props->pgd_props.change_threshold_num_stations)) + { + LOG_DEBUGMSG("%s: Returning true for num stations: %d >= %d", + __func__, + abs(core->numStations - last_sent_core->numStations), + props->pgd_props.change_threshold_num_stations); + return true; + } + + // If none of the thresholds are exceeded, return false + LOG_DEBUGMSG("%s: Returning false, no thresholds met", __func__); + return false; +} + +/* + * Return true if this message should not be sent (false if it should be sent) + */ +bool eewUtils_sendXMLFilter( + const struct GFAST_props_struct *props, + const struct GFAST_shakeAlert_struct *SA, + const struct GFAST_pgdResults_struct *pgd, + const struct GFAST_peakDisplacementData_struct *pgd_data, + const struct coreInfo_struct *core, + const struct coreInfo_struct *last_sent_core, + const double age_of_event) +{ + + // Conditions can be met or exceeded. + // Exceeded is used if a value being more than a threshold means no throttling + // Met is used if a value being less than a threshold means no throttling + bool pgd_exceeded = false; + bool mag_exceeded = false; + bool mag_sigma_met = false; + bool change_thresholds_exceeded = false; + + // Determine if pgd threshold is exceeded n times + int num_pgd_exceeded = 0, i, i_throttle; + + // Find the correct throttle criteria based on the time after origin + i_throttle = -1; + for (i = 0; i < props->pgd_props.n_throttle; i++) { + if (props->pgd_props.throttle_time_threshold[i] > age_of_event) break; + i_throttle++; + } + i_throttle = (i_throttle < 0) ? 0: i_throttle; + if (props->verbose > 2) { + LOG_DEBUGMSG("%s: For age_of_event %.2f, using i_throttle=%d, thresholds for time=%.1f, pgd=%.1f, nsta=%d", + __func__, + age_of_event, + i_throttle, + props->pgd_props.throttle_time_threshold[i_throttle], + props->pgd_props.throttle_pgd_threshold[i_throttle], + props->pgd_props.throttle_num_stations[i_throttle]) + } + + // Assumes pgd->nsites = pgd_data->nsites and indices correspond + if (pgd->nsites != pgd_data->nsites) { + LOG_ERRMSG("%s: nsites don't match for pgd, pgd_data! %d, %d\n", + __func__, pgd->nsites, pgd_data->nsites); + return false; + } + for (i = 0; i < pgd->nsites; i++) { + // skip site if it wasn't used + if (!pgd->lsiteUsed[i]) { continue; } + // pd is in meters, so convert to cm before comparing to threshold + if (pgd_data->pd[i] * 100. >= props->pgd_props.throttle_pgd_threshold[i_throttle]) { + num_pgd_exceeded++; + } + } + + if (props->verbose > 2) { + LOG_DEBUGMSG("%s: PGD threshold of %.1f cm exceeded at %d stations (threshold num: %d)", + __func__, props->pgd_props.throttle_pgd_threshold[i_throttle], num_pgd_exceeded, + props->pgd_props.throttle_num_stations[i_throttle]); + } + if (num_pgd_exceeded >= props->pgd_props.throttle_num_stations[i_throttle]) { + LOG_MSG("%s: PGD threshold of %.1f cm exceeded at %d stations (threshold num: %d)", + __func__, props->pgd_props.throttle_pgd_threshold[i_throttle], num_pgd_exceeded, + props->pgd_props.throttle_num_stations[i_throttle]); + pgd_exceeded = true; + } + + // Determine if SA mag threshold is exceeded + if (props->verbose > 2) { + LOG_DEBUGMSG("%s: SA mag: %f, threshold mag: %f", + __func__, SA->mag, props->pgd_props.SA_mag_threshold); + } + if (SA->mag >= props->pgd_props.SA_mag_threshold) { + mag_exceeded = true; + LOG_MSG("%s: SA magnitude exceeded! SA mag: %f, threshold mag: %f", + __func__, SA->mag, props->pgd_props.SA_mag_threshold); + } + + // Determine if pgd magnitude sigma threshold is met + if (props->verbose > 2) { + LOG_DEBUGMSG("%s: PGD mag sigma: %f, threshold mag sigma: %f, PGD mag: %f", + __func__, core->magUncer, props->pgd_props.pgd_sigma_throttle, core->mag); + } + if (core->magUncer <= props->pgd_props.pgd_sigma_throttle) { + mag_sigma_met = true; + LOG_MSG("%s: PGD mag sigma met! PGD mag sigma: %f, threshold mag sigma: %f, PGD mag: %f", + __func__, core->magUncer, props->pgd_props.pgd_sigma_throttle, core->mag); + } + + // Check change thresholds exceeded from the last sent message (first one automatically counts) + if (props->verbose > 2) { + LOG_DEBUGMSG("%s: PGD change thresholds -> dmag:%f, dmagU:%f, dlat:%f, dlon:%f, dot:%f, dnsta:%d", + __func__, + core->mag - last_sent_core->mag, + core->magUncer - last_sent_core->magUncer, + core->lat - last_sent_core->lat, + core->lon - last_sent_core->lon, + core->origTime - last_sent_core->origTime, + core->numStations - last_sent_core->numStations + ); + } + if (eewUtils_changeThresholdsExceeded(props, core, last_sent_core)) { + change_thresholds_exceeded = true; + LOG_MSG("%s: PGD change thresholds exceeded! dmag:%f, dmagU:%f, dlat:%f, dlon:%f, dot:%f, dnsta:%d", + __func__, + core->mag - last_sent_core->mag, + core->magUncer - last_sent_core->magUncer, + core->lat - last_sent_core->lat, + core->lon - last_sent_core->lon, + core->origTime - last_sent_core->origTime, + core->numStations - last_sent_core->numStations + ); + } + + if (pgd_exceeded && mag_exceeded && mag_sigma_met && change_thresholds_exceeded) { + LOG_MSG("%s: Message not throttled (%s v%d)! pgd_exceeded: %d, mag_exceeded: %d, mag_sigma_met: %d, change_thresholds_exceeded: %d", + __func__, core->id, core->version, pgd_exceeded, mag_exceeded, mag_sigma_met, change_thresholds_exceeded); + return false; + } + + LOG_MSG("%s: Message throttled (%s v%d)! pgd_exceeded: %d, mag_exceeded: %d, mag_sigma_met: %d, change_thresholds_exceeded: %d", + __func__, core->id, core->version, pgd_exceeded, mag_exceeded, mag_sigma_met, change_thresholds_exceeded); + return true; +} \ No newline at end of file diff --git a/src/gfast_eew.c b/src/gfast_eew.c index 4f075d6d..75540fa1 100644 --- a/src/gfast_eew.c +++ b/src/gfast_eew.c @@ -9,7 +9,6 @@ #include "iscl/memory/memory.h" #include "iscl/os/os.h" #include "iscl/time/time.h" -#include "H5public.h" #include #include "dmlibWrapper.h" @@ -56,6 +55,7 @@ void printProgramInfo(const char *fcnm) { #ifdef ENABLE_PLOG LOG_MSG("%s", message); #else + // Need to use printf until GFAST_core_properties_initialize is called printf("%s\n", message); #endif free(message); @@ -218,7 +218,6 @@ int main(int argc, char **argv) if (USE_AMQ) { activeMQ_start(); // start library /* dmlib startup */ - if (USE_DMLIB) { /* Start connection, to be used by DMMessageSender. @@ -387,6 +386,8 @@ int main(int argc, char **argv) if (tloop < props.waitTime) { usleep((props.waitTime - tloop) * 1000 * 1000); + // Put in a continue here since t1 is reset at the top of the loop + // Want to make sure it is representing the actual start continue; } else if ((props.waitTime > 0.0) && (tloop >= 2 * props.waitTime)) { diff --git a/src/tests/CoreDataUT.cc b/src/tests/CoreDataUT.cc index 55e0d8d0..d5cc3673 100644 --- a/src/tests/CoreDataUT.cc +++ b/src/tests/CoreDataUT.cc @@ -76,6 +76,43 @@ TEST(CoreData, testReadMetaDataFileWithMetaDataNetworks) { free(metaDataNetworks); } +TEST(CoreData, testReadSiteMaskFile) { + const char *metaDataFile; + const char *siteMaskFile; + char **metaDataNetworks; + int ierr; + struct GFAST_data_struct gps_data; + + // initialize + metaDataFile = "data/merged_chanfile_coord.dat\0"; + siteMaskFile = "data/site_mask_file.dat\0"; + memset(&gps_data, 0, sizeof(struct GFAST_data_struct)); + metaDataNetworks = NULL; + + ierr = core_data_readMetaDataFile(metaDataFile, + metaDataNetworks, + 0, + &gps_data); + EXPECT_EQ(0, ierr) << "Error reading sites file"; + EXPECT_EQ(492, gps_data.stream_length); + + ierr = core_data_readSiteMaskFile(siteMaskFile, 1, &gps_data); + + // Verify the masks are correct + int npgd = 0, ncmt = 0, nff = 0; + for (int i = 0; i < gps_data.stream_length; i++) { + npgd += !gps_data.data[i].lskip_pgd; + ncmt += !gps_data.data[i].lskip_cmt; + nff += !gps_data.data[i].lskip_ff; + } + EXPECT_EQ(489, npgd); + EXPECT_EQ(490, ncmt); + EXPECT_EQ(491, nff); + + // finalize + GFAST_core_data_finalize(&gps_data); +} + TEST(CoreData, testInitialize) { char propfilename[1024]; int ierr, i; diff --git a/src/tests/CoreEventsUT.cc b/src/tests/CoreEventsUT.cc index 3b7a5ed4..e7c4427b 100644 --- a/src/tests/CoreEventsUT.cc +++ b/src/tests/CoreEventsUT.cc @@ -10,8 +10,8 @@ #include "gfast.h" -// test freeEvents, getMinOriginTime, newEvent, printEvent, removeCancelledEvent, -// removeExpiredEvent, removeExpiredEvents, syncXMLStatusWithEvents, updateEvent +// test freeEvents, newEvent, printEvent, +// removeExpiredEvent, removeExpiredEvents, syncXMLStatusWithEvents void get_SA(GFAST_shakeAlert_struct *SA) { // Initialize new event strcpy(SA->eventid, "12345\0"); @@ -149,3 +149,95 @@ TEST_F(CoreEventsFixture, testNewEventUpdatePreviousEvents) { EXPECT_STREQ(SA2.eventid, xml_status.SA_status[0].eventid); } + +TEST_F(CoreEventsFixture, testPrintEvent) { + EXPECT_NO_THROW(GFAST_core_events_printEvents(SA);); +} + +TEST_F(CoreEventsFixture, testRemoveExpiredEventsNone) { + int ret; + ret = core_events_removeExpiredEvents(100, SA.time, 1, &events); + EXPECT_EQ(0, ret); + EXPECT_EQ(0, events.nev); + EXPECT_EQ(0, xml_status.nev); +} + +TEST_F(CoreEventsFixture, testRemoveExpiredEventsOneSync) { + lnewEvent = GFAST_core_events_newEvent(SA, &events, &xml_status); + EXPECT_EQ(1, events.nev); + EXPECT_EQ(1, xml_status.nev); + int ret; + ret = core_events_removeExpiredEvents(100, SA.time + 200, 1, &events); + EXPECT_EQ(1, ret); + EXPECT_EQ(0, events.nev); + EXPECT_EQ(1, xml_status.nev); + + ret = core_events_syncXMLStatusWithEvents(&events, &xml_status); + EXPECT_EQ(0, xml_status.nev); +} + +TEST_F(CoreEventsFixture, testRemoveExpiredEventsMultipleSync) { + // Make the second event + struct GFAST_shakeAlert_struct SA2; + memset(&SA2, 0, sizeof(struct GFAST_shakeAlert_struct)); + get_SA(&SA2); + strcpy(SA2.eventid, "9876\0"); + SA2.version = 0; + SA2.mag = 3.9; + SA2.time += 300; + + // Add both events, verify + lnewEvent = GFAST_core_events_newEvent(SA, &events, &xml_status); + EXPECT_TRUE(lnewEvent); + + lnewEvent = GFAST_core_events_newEvent(SA2, &events, &xml_status); + EXPECT_TRUE(lnewEvent); + + EXPECT_EQ(2, events.nev); + EXPECT_EQ(2, xml_status.nev); + + // Now remove one + int ret; + ret = core_events_removeExpiredEvents(100, SA.time + 200, 1, &events); + EXPECT_EQ(1, ret); + EXPECT_EQ(1, events.nev); + EXPECT_EQ(2, xml_status.nev); + + ret = core_events_syncXMLStatusWithEvents(&events, &xml_status); + EXPECT_EQ(1, xml_status.nev); +} + +TEST_F(CoreEventsFixture, testRemoveExpiredEventExists) { + lnewEvent = GFAST_core_events_newEvent(SA, &events, &xml_status); + EXPECT_EQ(1, events.nev); + EXPECT_EQ(1, xml_status.nev); + + struct GFAST_shakeAlert_struct SArem; + memset(&SArem, 0, sizeof(struct GFAST_shakeAlert_struct)); + memcpy(&SArem, &SA, sizeof(struct GFAST_shakeAlert_struct)); + + int ret; + ret = core_events_removeExpiredEvent(100, SA.time + 200, 1, SArem, &events); + EXPECT_EQ(1, ret); + EXPECT_EQ(0, events.nev); + EXPECT_EQ(1, xml_status.nev); +} + +TEST_F(CoreEventsFixture, testRemoveExpiredEventNoExists) { + lnewEvent = GFAST_core_events_newEvent(SA, &events, &xml_status); + EXPECT_EQ(1, events.nev); + EXPECT_EQ(1, xml_status.nev); + + struct GFAST_shakeAlert_struct SA2; + memset(&SA2, 0, sizeof(struct GFAST_shakeAlert_struct)); + get_SA(&SA2); + strcpy(SA2.eventid, "9876\0"); + SA2.version = 0; + SA2.mag = 3.9; + + int ret; + ret = core_events_removeExpiredEvent(100, SA.time + 200, 1, SA2, &events); + EXPECT_EQ(0, ret); + EXPECT_EQ(1, events.nev); + EXPECT_EQ(1, xml_status.nev); +} diff --git a/src/tests/CoreWaveformProcessorUT.cc b/src/tests/CoreWaveformProcessorUT.cc index eec6ae82..d5475c5d 100644 --- a/src/tests/CoreWaveformProcessorUT.cc +++ b/src/tests/CoreWaveformProcessorUT.cc @@ -11,6 +11,8 @@ #include "iscl/memory/memory.h" #include "gfast.h" +#include "gfast_ut_utils.h" + /* * Test the peakDisplacementHelper function. With input data and uncertainty @@ -117,33 +119,336 @@ TEST(CoreWaveformProcessor, testPeakDisplacementHelperEarlyNansMiddlePeak) { expPeakDisp, 1, 2); } +TEST(CoreWaveformProcessor, testPeakDisplacementHelperMiddleNans) { + const int NPTS = 4; + double ubuff[NPTS] = {1, 2, NAN, 4}; + double nbuff[NPTS] = {1, 2, NAN, 4}; + double ebuff[NPTS] = {1, 2, NAN, 4}; + int npts = NPTS; + + double expPeakDisp = pow(pow(4 - 1, 2) + pow(4 - 1, 2) + pow(4 - 1, 2), 0.5); + + testPDHelper_helper(npts, ubuff, nbuff, ebuff, + expPeakDisp, 0, 3); +} + TEST(CoreWaveformProcessor, testParseQChannel) { double value; // Temp code to print all values - // int chimap; + // int chimap, goodness; // double chi2; // for (int i = 0; i < 100; i++) { // value = 12300 + i; // chi2 = core_waveformProcessor_parseQChannelChi2CWU(value); // chimap = core_waveformProcessor_parseQChannelChi2CWUmap(value); - // std::cout << value << ": " << chimap << ", " << chi2 << std::endl; + // goodness = core_waveformProcessor_parseQChannelGoodness(value); + // std::cout << value << ": " << chimap << ", " << chi2 << ", " << goodness << std::endl; // } // Test a few values value = 12300; EXPECT_EQ(0, core_waveformProcessor_parseQChannelChi2CWUmap(value)); EXPECT_DOUBLE_EQ(0.0001, core_waveformProcessor_parseQChannelChi2CWU(value)); + EXPECT_EQ(0, core_waveformProcessor_parseQChannelGoodness(value)); value = 612325; EXPECT_EQ(25, core_waveformProcessor_parseQChannelChi2CWUmap(value)); EXPECT_DOUBLE_EQ(0.001, core_waveformProcessor_parseQChannelChi2CWU(value)); + EXPECT_EQ(1, core_waveformProcessor_parseQChannelGoodness(value)); value = 765450; EXPECT_EQ(50, core_waveformProcessor_parseQChannelChi2CWUmap(value)); EXPECT_DOUBLE_EQ(0.01, core_waveformProcessor_parseQChannelChi2CWU(value)); + EXPECT_EQ(0, core_waveformProcessor_parseQChannelGoodness(value)); value = 234575; EXPECT_EQ(75, core_waveformProcessor_parseQChannelChi2CWUmap(value)); EXPECT_DOUBLE_EQ(0.1, core_waveformProcessor_parseQChannelChi2CWU(value)); + EXPECT_EQ(1, core_waveformProcessor_parseQChannelGoodness(value)); +} + +/* + * Fixture for testing peakDisplacement() function + * Has two "normal" stations that can be put out of spec in various ways + * in order to not return peak displacement measurements + */ +class CoreWaveformProcessorFixture : public::testing::Test { + protected: + struct GFAST_pgd_props_struct pgd_props; + struct GFAST_data_struct gps_data; + struct GFAST_peakDisplacementData_struct pgd_data; + struct GFAST_pgdResults_struct pgd; + struct GFAST_offsetData_struct offset_data; + + double lat, lon, dep, time; + int ierr, nsites_pgd, k, j; + + void SetUp() { + const int NPTS = 10; + + // Initialize necessary pgd properties + memset(&pgd_props, 0, sizeof(struct GFAST_pgd_props_struct)); + pgd_props.utm_zone = -12345; + pgd_props.window_vel = 3; + pgd_props.min_window_vel = 0.01; + pgd_props.minimum_pgd_cm = 0.0; + pgd_props.maximum_pgd_cm = 3500; + pgd_props.u_raw_sigma_threshold = 35; + pgd_props.n_raw_sigma_threshold = 17; + pgd_props.e_raw_sigma_threshold = 14; + pgd_props.q_value_threshold = 1; + pgd_props.ngridSearch_deps = 10; + pgd_props.ngridSearch_lats = 1; + pgd_props.ngridSearch_lons = 1; + pgd_props.min_sites = 2; + + // Initialize origin + lat = 46; + lon = -122; + dep = 8; + time = 0; + + // Initialize gps_data + memset(&gps_data, 0, sizeof(struct GFAST_data_struct)); + + gps_data.stream_length = 2; + gps_data.data = (struct GFAST_waveform3CData_struct *) + calloc((size_t) gps_data.stream_length, + sizeof(struct GFAST_waveform3CData_struct)); + fill_gps_data(&gps_data, 0, "UW\0", "TEST1\0", "LYZ\0", "00\0", 46.02, -122.02, 1, 1, 1, NPTS); + fill_gps_data(&gps_data, 1, "UW\0", "TEST2\0", "LYZ\0", "00\0", 46.03, -122.03, 1, 1, 1, NPTS); + + // Initialize pgd_data, pgd + memset(&pgd_data, 0, sizeof( struct GFAST_peakDisplacementData_struct)); + memset(&pgd, 0, sizeof(struct GFAST_pgdResults_struct)); + + ierr = core_scaling_pgd_initialize(pgd_props, gps_data, &pgd, &pgd_data); + + // Also initialize offset data (for cmt, ff) + memset(&offset_data, 0, sizeof(struct GFAST_offsetData_struct)); + initialize_offset_data(&gps_data, &offset_data); + + // Now actually put some data in gps_data + double tdata[NPTS] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; + double udata[NPTS] = {0, 3, 5, 6, 4, 3, -1, 0, 3, 4}; + double ndata[NPTS] = {0, 2, 3, 4, 5, 3, 2, 0, -1, -4}; + double edata[NPTS] = {0, -3, -5, -4, -4, -3, -1, 0, 2, 3}; + double usigmadata[NPTS] = {.1, .1, .1, .1, .1, .1, .1, .1, .1, .1}; + double nsigmadata[NPTS] = {.1, .1, .1, .1, .1, .1, .1, .1, .1, .1}; + double esigmadata[NPTS] = {.1, .1, .1, .1, .1, .1, .1, .1, .1, .1}; + double qdata[NPTS] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + + for (k = 0; k < gps_data.stream_length; k++) { + for (j = 0; j < NPTS; j++) { + gps_data.data[k].ubuff[j] = udata[j]; + gps_data.data[k].nbuff[j] = ndata[j]; + gps_data.data[k].ebuff[j] = edata[j]; + gps_data.data[k].usigmabuff[j] = usigmadata[j]; + gps_data.data[k].nsigmabuff[j] = nsigmadata[j]; + gps_data.data[k].esigmabuff[j] = esigmadata[j]; + gps_data.data[k].qbuff[j] = qdata[j]; + gps_data.data[k].tbuff[j] = tdata[j]; + } + gps_data.data[k].npts = NPTS; + } + } + + void TearDown() {} + + void initialize_offset_data(struct GFAST_data_struct *gps_data, + struct GFAST_offsetData_struct *offset_data) + { + offset_data->stnm = (char **)calloc((size_t) gps_data->stream_length, + sizeof(char *)); + offset_data->ubuff = memory_calloc64f(gps_data->stream_length); + offset_data->nbuff = memory_calloc64f(gps_data->stream_length); + offset_data->ebuff = memory_calloc64f(gps_data->stream_length); + offset_data->wtu = memory_calloc64f(gps_data->stream_length); + offset_data->wtn = memory_calloc64f(gps_data->stream_length); + offset_data->wte = memory_calloc64f(gps_data->stream_length); + offset_data->sta_lat = memory_calloc64f(gps_data->stream_length); + offset_data->sta_lon = memory_calloc64f(gps_data->stream_length); + offset_data->sta_alt = memory_calloc64f(gps_data->stream_length); + offset_data->lmask = memory_calloc8l(gps_data->stream_length); + offset_data->lactive = memory_calloc8l(gps_data->stream_length); + offset_data->nsites = gps_data->stream_length; + + for (int i=0; insites; i++) + { + offset_data->sta_lat[i] = gps_data->data[i].sta_lat; + offset_data->sta_lon[i] = gps_data->data[i].sta_lon; + offset_data->sta_alt[i] = gps_data->data[i].sta_alt; + offset_data->stnm[i] = (char *)calloc(64, sizeof(char)); + strcpy(offset_data->stnm[i], gps_data->data[i].netw); + strcat(offset_data->stnm[i], ".\0"); + strcat(offset_data->stnm[i], gps_data->data[i].stnm); + strcat(offset_data->stnm[i], ".\0"); + strncpy(offset_data->stnm[i], gps_data->data[i].chan[0], 2); + strcat(offset_data->stnm[i], "?.\0"); + if (strlen(gps_data->data[i].loc) > 0) + { + strcat(offset_data->stnm[i], gps_data->data[i].loc); + } + offset_data->lmask[i] = false; + } + } +}; + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementNormal) { + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(2, nsites_pgd); + EXPECT_EQ(1, pgd_data.lactive[0]); + EXPECT_EQ(1, pgd_data.lactive[1]); + + EXPECT_TRUE(lequal(8.24621, pgd_data.pd[0], .00001)) << "Error pgd is wrong"; +} + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementTooDistant) { + // If a station is not within the s-wave mask, it should be ignored + gps_data.data[0].sta_lat = 40; + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(1, nsites_pgd); + EXPECT_EQ(0, pgd_data.lactive[0]); + EXPECT_EQ(1, pgd_data.lactive[1]); +} + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementMasked) { + // If a site is masked (via siteMaskFile) it should be ignored + pgd_data.lmask[1] = true; + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(1, nsites_pgd); + EXPECT_EQ(1, pgd_data.lactive[0]); + EXPECT_EQ(0, pgd_data.lactive[1]); +} + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementMinPGD) { + // A displacement that is 0.0 (flat-lined) should be ignored + for (j = 0; j < 10; j++) { + gps_data.data[0].ubuff[j] = 0; + gps_data.data[0].nbuff[j] = 0; + gps_data.data[0].ebuff[j] = 0; + } + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(1, nsites_pgd); + EXPECT_EQ(0, pgd_data.lactive[0]); + EXPECT_EQ(1, pgd_data.lactive[1]); +} + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementMaxPGD) { + // A displacement that is too large should be ignored + gps_data.data[1].ebuff[3] = 1000; + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(1, nsites_pgd); + EXPECT_EQ(1, pgd_data.lactive[0]); + EXPECT_EQ(0, pgd_data.lactive[1]); +} + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementSigma) { + // An uncertainty that is too large means the pgd should be ignored + gps_data.data[0].esigmabuff[0] = 10; + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(1, nsites_pgd); + EXPECT_EQ(0, pgd_data.lactive[0]); + EXPECT_EQ(1, pgd_data.lactive[1]); +} + +TEST_F(CoreWaveformProcessorFixture, testPeakDisplacementQThreshold) { + // A q value that doesn't meet the criteria means the pgd should be ignored + gps_data.data[1].qbuff[0] = 0; + nsites_pgd = GFAST_core_waveformProcessor_peakDisplacement( + &pgd_props, + lat, + lon, + dep, + time, + gps_data, + &pgd_data, + &ierr); + + EXPECT_EQ(1, nsites_pgd); + EXPECT_EQ(1, pgd_data.lactive[0]); + EXPECT_EQ(0, pgd_data.lactive[1]); +} + +TEST_F(CoreWaveformProcessorFixture, testOffset) { + int nsites, ierr; + // An svel_window of 3.5 will put the swave_times at 2.41 and 2.56 + double svel_window = 3.5; + + nsites = core_waveformProcessor_offset( + -12345, + svel_window, + lat, + lon, + dep, + time, + gps_data, + &offset_data, + &ierr); + + EXPECT_EQ(0, ierr); + EXPECT_EQ(2, nsites); + + double toler = 1e-5; + EXPECT_NEAR(3, offset_data.ubuff[0], toler); + EXPECT_NEAR(1.5, offset_data.nbuff[0], toler); + EXPECT_NEAR(-1.5, offset_data.ebuff[0], toler); + EXPECT_NEAR(2.71429, offset_data.ubuff[1], toler); + EXPECT_NEAR(1.28571, offset_data.nbuff[1], toler); + EXPECT_NEAR(-1, offset_data.ebuff[1], toler); } diff --git a/src/tests/DmLibWrapperUT.cc b/src/tests/DmLibWrapperUT.cc new file mode 100644 index 00000000..e045fe1d --- /dev/null +++ b/src/tests/DmLibWrapperUT.cc @@ -0,0 +1,112 @@ +/** + * @file DmLibWrapperUT.cc + * @author Carl Ulberg, University of Washington (ulbergc@uw.edu) + * @brief This tests files in the dmlib/ directory. + */ + +#include "gtest/gtest.h" +#include "ut_log_init.h" +#include "ut_main.h" +#include "ut_xml_utils.h" + +#include "gfast.h" +#include "gfast_eewUtils.h" +#include "dmlibWrapper.h" + +#include "gfast_ut_utils.h" + +#include "iscl/memory/memory.h" +#include "iscl/array/array.h" +#include "iscl/iscl/iscl_enum.h" + + +TEST(DmLibWrapper, testCreatePGDXml) { + char *pgdXML; + const char *filenm = "data/final_pgd.maule.txt\0"; + const char *program_instance = "gfast@eew-uw-dev1\0"; + const char *message_type = "update\0"; + struct GFAST_pgd_props_struct pgd_props; + struct GFAST_peakDisplacementData_struct pgd_data; + struct GFAST_pgdResults_struct pgd_ref, pgd; + double SA_lat, SA_lon, SA_dep, age_of_event; + int i, ierr, max_assoc_stations, pgdOpt; + double currentTime = 0; + enum isclError_enum ierr_iscl; + + pgdXML = NULL; + memset(&pgd_props, 0, sizeof(pgd_props)); + memset(&pgd_data, 0, sizeof(pgd_data)); + memset(&pgd_ref, 0, sizeof(pgd_ref)); + memset(&pgd, 0, sizeof(pgd)); + age_of_event = 0; + max_assoc_stations = 1; + ierr = read_pgd_results(filenm, + &pgd_props, + &pgd_data, + &pgd_ref, + &SA_lat, &SA_lon, &SA_dep); + + EXPECT_EQ(0, ierr) << "Error reading input file"; + + // Set space + pgd.nsites = pgd_ref.nsites; + pgd.ndeps = pgd_ref.ndeps; + pgd.nlats = 1; + pgd.nlons = 1; + pgd.mpgd = ISCL_memory_calloc__double(pgd.ndeps); + pgd.mpgd_sigma = ISCL_memory_calloc__double(pgd.ndeps); + pgd.mpgd_vr = ISCL_memory_calloc__double(pgd.ndeps); + pgd.dep_vr_pgd = ISCL_memory_calloc__double(pgd.ndeps); + pgd.iqr = ISCL_memory_calloc__double(pgd.ndeps); + pgd.UP = ISCL_memory_calloc__double(pgd.ndeps*pgd.nsites); + pgd.UPinp = ISCL_memory_calloc__double(pgd.nsites); + pgd.srcDepths = ISCL_memory_calloc__double(pgd.ndeps); + pgd.srdist = ISCL_memory_calloc__double(pgd.ndeps*pgd.nsites); + pgd.lsiteUsed = ISCL_memory_calloc__bool(pgd.nsites); + for (i=0; i + + +class Hdf5Fixture : public::testing::Test { + protected: + + const char *dir_tmp = "data/tmp\0"; + const char *evid = "1234\0"; + // Careful with expected_path!! It will be deleted on TearDown + const char *expected_path = "data/tmp/1234_archive.h5\0"; + const char *propfilename = "data/gfast.props\0"; + int ierr; + hid_t file; + + void SetUp() { + ierr = hdf5_initialize(dir_tmp, evid, propfilename); + + EXPECT_EQ(0, ierr); + htri_t is_file; + is_file = H5Fis_hdf5(expected_path); + EXPECT_EQ(1, is_file); + } + + void TearDown() { + if (H5Fis_hdf5(expected_path)) { + file = H5Fopen(expected_path, H5F_ACC_RDWR, H5P_DEFAULT); + H5Fflush(file, H5F_SCOPE_GLOBAL); + H5Fclose(file); + remove(expected_path); + } + } + +}; + +TEST_F(Hdf5Fixture, testSetFileName) { + char fname[PATH_MAX]; + + ierr = hdf5_setFileName(dir_tmp, evid, fname); + + EXPECT_EQ(0, ierr); + EXPECT_STREQ(expected_path, fname); +} + +TEST_F(Hdf5Fixture, testSetFileNameNullEvid) { + char fname[PATH_MAX]; + evid = NULL; + + ierr = hdf5_setFileName(dir_tmp, evid, fname); + + EXPECT_EQ(-1, ierr); +} + +TEST_F(Hdf5Fixture, testSetFileNameEmptyEvid) { + char fname[PATH_MAX]; + evid = ""; + + ierr = hdf5_setFileName(dir_tmp, evid, fname); + + EXPECT_EQ(-1, ierr); +} + +TEST_F(Hdf5Fixture, testSetFileNameNullAdir) { + char fname[PATH_MAX]; + dir_tmp = NULL; + + ierr = hdf5_setFileName(dir_tmp, evid, fname); + + EXPECT_EQ(0, ierr); + EXPECT_STREQ("./1234_archive.h5", fname); +} + +TEST_F(Hdf5Fixture, testSetFileNameEmptyAdir) { + char fname[PATH_MAX]; + dir_tmp = ""; + + ierr = hdf5_setFileName(dir_tmp, evid, fname); + + EXPECT_EQ(0, ierr); + EXPECT_STREQ("./1234_archive.h5", fname); +} + +TEST_F(Hdf5Fixture, testCinterItemExists) { + file = H5Fopen(expected_path, H5F_ACC_RDONLY, H5P_DEFAULT); + ASSERT_EQ(1, h5_item_exists(file, "/DataStructures")); + EXPECT_EQ(1, h5_item_exists(file, "/DataStructures/peakDisplacementDataStructure")); + ASSERT_EQ(1, h5_item_exists(file, "/InitializationFile")); + EXPECT_EQ(1, h5_item_exists(file, "/InitializationFile/IniFile")); + EXPECT_EQ(1, h5_item_exists(file, "/GFAST_History")); + EXPECT_EQ(1, h5_item_exists(file, "/Summary")); + EXPECT_EQ(0, h5_item_exists(file, "/Fake")); + H5Fclose(file); +} diff --git a/src/tests/Makefile b/src/tests/Makefile index 31cf999b..e123eaa2 100644 --- a/src/tests/Makefile +++ b/src/tests/Makefile @@ -1,5 +1,6 @@ # Makefile for gfast unit tests -# Modeled after libs/utils/tests/Makefile, may need additional changes +# Currently setup to run within a ShakeAlert directory. Will need work +# to run in a standalone GFAST installation EEWDIR = ../../.. include $(EEWDIR)/Make.include.$(shell uname) @@ -23,20 +24,23 @@ ifdef ENABLE_PLOG CFLAGS+=$(ENABLE_PLOG) endif -OBJS = ../activeMQ/gfast_amq.a ../eewUtils/gfast_utils.a ../xml/gfast_xml.a ../core/gfast_core.a \ - ../traceBuffer/traceBuffer.a ../hdf5/hdf5.a ../activeMQ/gfast_amq.a +OBJS = ../dmlib/dmlibWrapper.cpp.o ../activeMQ/gfast_amq.a ../eewUtils/gfast_utils.a \ + ../xml/gfast_xml.a ../core/gfast_core.a ../traceBuffer/traceBuffer.a ../hdf5/hdf5.a \ + ../activeMQ/gfast_amq.a DATA = UNITTESTS = \ - CoreWaveformProcessorUT \ - CoreActiveMqUT \ - EewUtilsUT \ - CoreScalingUT \ - CoreCoordtoolsUT \ - CoreDataUT \ - CoreEventsUT \ - TraceBufferEwrrUT + CoreActiveMqUT \ + CoreCoordtoolsUT \ + CoreDataUT \ + CoreEventsUT \ + CoreScalingUT \ + CoreWaveformProcessorUT \ + DmLibWrapperUT \ + EewUtilsUT \ + Hdf5UT \ + TraceBufferEwrrUT INCL = -I$(THIRD_PARTY) \ $(AQMS_INCL) $(OTL_INCL) $(APR_INCL) $(ACTIVEMQ_INCL) $(QLIB2_INCL) \ diff --git a/src/tests/data/eewUtils_quakeML_example.cmt b/src/tests/data/eewUtils_quakeML_example.cmt new file mode 100644 index 00000000..16dad56e --- /dev/null +++ b/src/tests/data/eewUtils_quakeML_example.cmt @@ -0,0 +1,119 @@ + + + + + + + + 9.922327e+00 + + + + 3.000000e+00 + + + 1.000000e+00 + + + 2.000000e+00 + + + 5.000000e+00 + + + -6.000000e+00 + + + -4.000000e+00 + + + + 0.091126 + + + 0.707308 + + + + + + 90.000000 + + + 7.460045 + + + -144.735610 + + + + + 324.964874 + + + 85.701068 + + + -83.897425 + + + + + + + 49.299988 + + + 40.399763 + + + 1.017597e+01 + + + + + 241.537565 + + + 48.949095 + + + -5.668683e+00 + + + + + 144.505696 + + + 6.085340 + + + -4.507288e+00 + + + + + + + -5.425927e+00 + + Mw_gps + + + + + 46.000000 + + + -122.000000 + + + 8.000000 + + + + + diff --git a/src/tests/data/eewUtils_xml_example.ff b/src/tests/data/eewUtils_xml_example.ff new file mode 100644 index 00000000..80de6402 --- /dev/null +++ b/src/tests/data/eewUtils_xml_example.ff @@ -0,0 +1,2 @@ + +8.7531610.500000-35.909000nan-72.733000nan35.000000nan1970-01-01T00:00:00.000Znan0.800000501-38.091856285.13434519.488937-37.598740285.35754119.488937-37.719242285.78487725.693362-38.213361285.56399425.693362-0.0004700.8084060.0003250.80840616.2063408.794004-37.598740285.35754119.488937-37.104947285.57742919.488937-37.224473286.00249125.693362-37.719242285.78487725.693362-1.49836323.2643161.69135123.89850116.2063408.794004-37.104947285.57742919.488937-36.610500285.79411319.488937-36.729075286.21694125.693362-37.224473286.00249125.693362-2.79490823.8625793.88004722.66786416.2063408.794004-36.610500285.79411319.488937-36.115419286.00769519.488937-36.233069286.42832725.693362-36.729075286.21694125.693362-3.49189324.1993635.30090120.60648316.2063408.794004-36.115419286.00769519.488937-35.619727286.21827319.488937-35.736475286.63674825.693362-36.233069286.42832725.693362-3.39320624.3572196.44557223.50557116.2063408.794004-35.619727286.21827319.488937-35.123443286.42594219.488937-35.239314286.84229825.693362-35.736475286.63674825.693362-2.96157623.6246318.84914921.78123016.2063408.794004-35.123443286.42594219.488937-34.626587286.63079419.488937-34.741605287.04506925.693362-35.239314286.84229825.693362-2.30902823.8226919.91075319.09695716.2063408.794004-34.626587286.63079419.488937-34.129178286.83291819.488937-34.243367287.24515125.693362-34.741605287.04506925.693362-1.30373323.0551287.13783523.06638716.2063408.794004-34.129178286.83291819.488937-33.631236287.03240319.488937-33.744618287.44263025.693362-34.243367287.24515125.693362-0.49071120.8933383.84405720.97300716.2063408.794004-33.631236287.03240319.488937-33.132778287.22933319.488937-33.245377287.63759125.693362-33.744618287.44263025.6933620.0002720.8083630.0020460.80822616.2063408.794004-38.213361285.56399425.693362-37.719242285.78487725.693362-37.838282286.21375831.897787-38.333384285.99523931.897787-0.0005430.8082980.0002110.80830316.2063408.794004-37.719242285.78487725.693362-37.224473286.00249125.693362-37.342557286.42904931.897787-37.838282286.21375831.897787-2.07789725.3603552.49285928.05429516.2063408.794004-37.224473286.00249125.693362-36.729075286.21694125.693362-36.846228286.64121731.897787-37.342557286.42904931.897787-3.94980421.9858375.84181323.81722616.2063408.794004-36.729075286.21694125.693362-36.233069286.42832725.693362-36.349316286.85036231.897787-36.846228286.64121731.897787-5.20271323.2114777.79186718.26401216.2063408.794004-36.233069286.42832725.693362-35.736475286.63674825.693362-35.851841287.05658031.897787-36.349316286.85036231.897787-4.47632824.6312088.72546625.25795516.2063408.794004-35.736475286.63674825.693362-35.239314286.84229825.693362-35.353823287.25996831.897787-35.851841287.05658031.897787-3.65311121.39681411.76334421.82075316.2063408.794004-35.239314286.84229825.693362-34.741605287.04506925.693362-34.855281287.46061631.897787-35.353823287.25996831.897787-4.19538022.80240713.74492616.72109016.2063408.794004-34.741605287.04506925.693362-34.243367287.24515125.693362-34.356234287.65861431.897787-34.855281287.46061631.897787-2.64267223.0509779.79572625.23431716.2063408.794004-34.243367287.24515125.693362-33.744618287.44263025.693362-33.856698287.85404831.897787-34.356234287.65861431.897787-0.66503819.5165385.22089522.48429816.2063408.794004-33.744618287.44263025.693362-33.245377287.63759125.693362-33.356693288.04700131.897787-33.856698287.85404831.8977870.0000880.8082260.0030730.80706016.2063408.794004-38.333384285.99523931.897787-37.838282286.21375831.897787-37.955840286.64416638.102213-38.451903286.42806238.102213-0.0002830.808293-0.0000060.80830316.2063408.794004-37.838282286.21375831.897787-37.342557286.42904931.897787-37.459178286.85708638.102213-37.955840286.64416638.102213-1.49880923.3020352.18864727.51705516.2063408.794004-37.342557286.42904931.897787-36.846228286.64121731.897787-36.961939287.06692338.102213-37.459178286.85708638.102213-2.81394814.0033505.35355124.29526216.2063408.794004-36.846228286.64121731.897787-36.349316286.85036231.897787-36.464141287.27378038.102213-36.961939287.06692338.102213-4.53835516.8123307.09157019.88782916.2063408.794004-36.349316286.85036231.897787-35.851841287.05658031.897787-35.965806287.47775238.102213-36.464141287.27378038.102213-2.45685220.2482666.38555623.94707716.2063408.794004-35.851841287.05658031.897787-35.353823287.25996831.897787-35.466951287.67893438.102213-35.965806287.47775238.102213-0.74988313.2587167.45758818.15491116.2063408.794004-35.353823287.25996831.897787-34.855281287.46061631.897787-34.967596287.87741738.102213-35.466951287.67893438.102213-5.43032016.3631939.94584617.24664316.2063408.794004-34.855281287.46061631.897787-34.356234287.65861431.897787-34.467759288.07329038.102213-34.967596287.87741738.102213-3.63657519.3033117.40308223.57533016.2063408.794004-34.356234287.65861431.897787-33.856698287.85404831.897787-33.967458288.26663838.102213-34.467759288.07329038.1022130.02520313.6439094.01266123.15282216.2063408.794004-33.856698287.85404831.897787-33.356693288.04700131.897787-33.466709288.45754638.102213-33.967458288.26663838.1022130.0013570.808060-0.0004960.80656316.2063408.794004-38.451903286.42806238.102213-37.955840286.64416638.102213-38.071893287.07608144.306638-38.568897286.86244144.3066380.0000450.808291-0.0001610.80829516.2063408.794004-37.955840286.64416638.102213-37.459178286.85708638.102213-37.574317287.28657944.306638-38.071893287.07608144.306638-0.49340620.1389281.18203221.85892816.2063408.794004-37.459178286.85708638.102213-36.961939287.06692338.102213-37.076188287.49404044.306638-37.574317287.28657944.306638-0.83308717.4331233.25304919.42794516.2063408.794004-36.961939287.06692338.102213-36.464141287.27378038.102213-36.577525287.69856244.306638-37.076188287.49404044.306638-1.82422418.5749804.41879216.10049116.2063408.794004-36.464141287.27378038.102213-35.965806287.47775238.102213-36.078350287.90024344.306638-36.577525287.69856244.3066380.54735318.2868672.50644619.96443816.2063408.794004-35.965806287.47775238.102213-35.466951287.67893438.102213-35.578679288.09917744.306638-36.078350287.90024344.3066382.68012314.7983901.46955913.55755116.2063408.794004-35.466951287.67893438.102213-34.967596287.87741738.102213-35.078532288.29545344.306638-35.578679288.09917744.306638-4.36011715.5256993.61704710.92210916.2063408.794004-34.967596287.87741738.102213-34.467759288.07329038.102213-34.577926288.48916044.306638-35.078532288.29545344.306638-2.94914716.8009223.13486519.14083416.2063408.794004-34.467759288.07329038.102213-33.967458288.26663838.102213-34.076878288.68038444.306638-34.577926288.48916044.3066380.73496312.9112442.12617619.74075616.2063408.794004-33.967458288.26663838.102213-33.466709288.45754638.102213-33.575406288.86920744.306638-34.076878288.68038444.3066380.0032130.8077730.0109880.80435716.2063408.794004-38.568897286.86244144.306638-38.071893287.07608144.306638-38.186422287.50948150.511063-38.684345287.29835650.5110630.0001250.808321-0.0003880.80832316.2063408.794004-38.071893287.07608144.306638-37.574317287.28657944.306638-37.687952287.71751050.511063-38.186422287.50948150.5110630.0000730.8086080.0003470.80863116.2063408.794004-37.574317287.28657944.306638-37.076188287.49404044.306638-37.188954287.92254650.511063-37.687952287.71751050.5110630.0005500.8085280.0022200.80854116.2063408.794004-37.076188287.49404044.306638-36.577525287.69856244.306638-36.689448288.12468950.511063-37.188954287.92254650.511063-0.0004920.8085520.0036590.80842716.2063408.794004-36.577525287.69856244.306638-36.078350287.90024344.306638-36.189454288.32403450.511063-36.689448288.12468950.5110630.0050140.8083610.0010920.80850116.2063408.794004-36.078350287.90024344.306638-35.578679288.09917744.306638-35.688988288.52067650.511063-36.189454288.32403450.5110630.0176860.8071120.0042710.80793916.2063408.794004-35.578679288.09917744.306638-35.078532288.29545344.306638-35.188069288.71470350.511063-35.688988288.52067650.511063-0.0146150.8074110.0067330.80783516.2063408.794004-35.078532288.29545344.306638-34.577926288.48916044.306638-34.686715288.90620550.511063-35.188069288.71470350.511063-0.0058960.8082690.0011520.80843416.2063408.794004-34.577926288.48916044.306638-34.076878288.68038444.306638-34.184942289.09526550.511063-34.686715288.90620550.5110630.0021670.8082570.0025250.80842916.2063408.794004-34.076878288.68038444.306638-33.575406288.86920744.306638-33.682768289.28196650.511063-34.184942289.09526550.5110630.0009820.8080740.0055430.80757116.2063408.794004 diff --git a/src/tests/data/eewUtils_xml_example.pgd b/src/tests/data/eewUtils_xml_example.pgd new file mode 100644 index 00000000..e3a0b881 --- /dev/null +++ b/src/tests/data/eewUtils_xml_example.pgd @@ -0,0 +1,2 @@ + +6.5000000.40000046.000000nan-122.000000nan8.000000nan1970-01-01T00:00:00.000Znan0.800000 diff --git a/src/tests/data/gfast.props b/src/tests/data/gfast.props index e73ffdcd..ee27a995 100644 --- a/src/tests/data/gfast.props +++ b/src/tests/data/gfast.props @@ -10,7 +10,7 @@ metaDataFile=data/merged_chanfile_coord.dat #SA_events_dir=xxxxx SA_output_dir=data #output_interval_mins=1,2,3,4,5,10 -#siteMaskFile=xxxxx +siteMaskFile=data/site_mask_file.dat # File with station names and sampling periods #dtfile=GFAST_streams_dt.txt # ElarmS message file @@ -97,18 +97,37 @@ pgd_window_vel = 3.0 pgd_min_sites = 4 # Message throttling logic and related +# Lookup table for non-density corrected value to generate magnitude +# uncertainty based on time and magnitude sigmaLookupFile = data/M99.txt +# Lookup table for time-dependant PGD thresholds +# columns: time (s) | threshold (cm) | n stations +pgdThresholdLookupFile = data/pgd_threshold_PW.txt +# Positional uncertainty thresholds for including peak displacement observations +# If not specified, allow any or NAN uncertainties +rawSigmaThresholdLookupFile = data/raw_sigma_threshold_PW.txt # If defined, don't send xml if pgd_sigma is greater than this value pgd_sigma_throttle = 0.5 # Message throttling logic -SA_mag_threshold = -5.0 +SA_mag_threshold = 6.0 # Minimum and maximum thresholds for using pgd values in inversion, in cm minimum_pgd_cm = 0.0 maximum_pgd_cm = 3500.0 max_assoc_stations = 6 +# Change thresholds to decide when to send a message, assuming other throttles +# pass. If any of these are met or exceeded, a new message would be sent. If +# none are exceeded, message would not be sent. +# If not present, threshold will be ignored +change_threshold_mag = 0.05 +change_threshold_mag_uncer = 0.05 +change_threshold_lat = 0.05 +change_threshold_lon = 0.05 +change_threshold_orig_time = 1 +change_threshold_num_stations = 1 + ################################################################################ # CMT/Finite Fault Properties # ################################################################################ diff --git a/src/tests/data/site_mask_file.dat b/src/tests/data/site_mask_file.dat new file mode 100644 index 00000000..9f061254 --- /dev/null +++ b/src/tests/data/site_mask_file.dat @@ -0,0 +1,4 @@ +# net sta loc cha usePGD useCMT useFF +CI 7ODM 20 LYE 0 0 1 +CI ALPP 20 LYE 0 1 1 +CI BCWR 20 LYE 0 0 0 \ No newline at end of file diff --git a/src/tests/gfast_ut_utils.h b/src/tests/gfast_ut_utils.h index 5b24cf34..1c1979db 100644 --- a/src/tests/gfast_ut_utils.h +++ b/src/tests/gfast_ut_utils.h @@ -56,13 +56,16 @@ int read_pgd_results(const char *filenm, if (fgets(cline, sizeof(cline), infl) == NULL){goto ERROR;} sscanf(cline, "%lf %lf %lf\n", SA_lat, SA_lon, SA_dep); // pgd data - pgd_data->pd = ISCL_memory_calloc__double(pgd_data->nsites); - pgd_data->wt = ISCL_memory_calloc__double(pgd_data->nsites); - pgd_data->sta_lat = ISCL_memory_calloc__double(pgd_data->nsites); - pgd_data->sta_lon = ISCL_memory_calloc__double(pgd_data->nsites); - pgd_data->sta_alt = ISCL_memory_calloc__double(pgd_data->nsites); - pgd_data->lmask = ISCL_memory_calloc__bool(pgd_data->nsites); - pgd_data->lactive = ISCL_memory_calloc__bool(pgd_data->nsites); + pgd_data->stnm = (char **)calloc((size_t) pgd_data->nsites, + sizeof(char *)); + pgd_data->pd = memory_calloc64f(pgd_data->nsites); + pgd_data->wt = memory_calloc64f(pgd_data->nsites); + pgd_data->sta_lat = memory_calloc64f(pgd_data->nsites); + pgd_data->sta_lon = memory_calloc64f(pgd_data->nsites); + pgd_data->sta_alt = memory_calloc64f(pgd_data->nsites); + pgd_data->pd_time = memory_calloc64f(pgd_data->nsites); + pgd_data->lmask = memory_calloc8l(pgd_data->nsites); + pgd_data->lactive = memory_calloc8l(pgd_data->nsites); for (i=0; insites; i++) { memset(cline, 0, sizeof(cline)); @@ -70,8 +73,12 @@ int read_pgd_results(const char *filenm, sscanf(cline, "%lf %lf %lf %lf\n", &pgd_data->sta_lat[i], &pgd_data->sta_lon[i], &pgd_data->sta_alt[i], &pgd_data->pd[i]); + pgd_data->pd_time[i] = 0.0; pgd_data->wt[i] = 1.0; pgd_data->lactive[i] = true; + + pgd_data->stnm[i] = (char *)calloc(64, sizeof(char)); + sprintf(pgd_data->stnm[i],"XX%02d", i); } // Results + depths in grid search pgd->nsites = pgd_data->nsites; @@ -311,3 +318,55 @@ int read_ff_results(const char *fname, ERROR:; return ierr; } + +void fill_gps_data(struct GFAST_data_struct *gps_data, + const int k, + const char *netw, + const char *stat, + const char *chan, + const char *loc, + const double lat, + const double lon, + const double elev, + const double dt, + const double gain, + const int mpts) { + strcpy(gps_data->data[k].netw, netw); + strcpy(gps_data->data[k].stnm, stat); + strncpy(gps_data->data[k].chan[0], chan, 2); + strcat( gps_data->data[k].chan[0], "Z\0"); + strncpy(gps_data->data[k].chan[1], chan, 2); + strcat( gps_data->data[k].chan[1], "N\0"); + strncpy(gps_data->data[k].chan[2], chan, 2); + strcat( gps_data->data[k].chan[2], "E\0"); + strncpy(gps_data->data[k].chan[3], chan, 2); + strcat( gps_data->data[k].chan[3], "3\0"); + strncpy(gps_data->data[k].chan[4], chan, 2); + strcat( gps_data->data[k].chan[4], "2\0"); + strncpy(gps_data->data[k].chan[5], chan, 2); + strcat( gps_data->data[k].chan[5], "1\0"); + strncpy(gps_data->data[k].chan[6], chan, 2); + strcat( gps_data->data[k].chan[6], "Q\0"); + strcpy(gps_data->data[k].loc, loc); + gps_data->data[k].sta_lat = lat; + gps_data->data[k].sta_lon = lon; + gps_data->data[k].sta_alt = elev; + gps_data->data[k].dt = dt; + gps_data->data[k].gain[0] = gain; + gps_data->data[k].gain[1] = gain; + gps_data->data[k].gain[2] = gain; + gps_data->data[k].gain[3] = gain; + gps_data->data[k].gain[4] = gain; + gps_data->data[k].gain[5] = gain; + gps_data->data[k].gain[6] = 1; // Quality channel has no gain + + gps_data->data[k].maxpts = mpts; + gps_data->data[k].ubuff = memory_calloc64f(mpts); + gps_data->data[k].nbuff = memory_calloc64f(mpts); + gps_data->data[k].ebuff = memory_calloc64f(mpts); + gps_data->data[k].usigmabuff = memory_calloc64f(mpts); + gps_data->data[k].nsigmabuff = memory_calloc64f(mpts); + gps_data->data[k].esigmabuff = memory_calloc64f(mpts); + gps_data->data[k].qbuff = memory_calloc64f(mpts); + gps_data->data[k].tbuff = memory_calloc64f(mpts); +} diff --git a/src/traceBuffer/ewrr/unpackTraceBuf2Messages.c b/src/traceBuffer/ewrr/unpackTraceBuf2Messages.c index 990a1729..40e523f6 100644 --- a/src/traceBuffer/ewrr/unpackTraceBuf2Messages.c +++ b/src/traceBuffer/ewrr/unpackTraceBuf2Messages.c @@ -94,7 +94,7 @@ int traceBuffer_ewrr_unpackTraceBuf2Messages( //char **msg_logos = (char **)malloc(sizeof(char *) * nRead); //char msg_logos[nRead][15]; - char buf[15]; + char buf[256]; char *logo, *nscl; int debug = 0; diff --git a/src/xml/Makefile b/src/xml/Makefile index 6dcbbbff..63996d0a 100644 --- a/src/xml/Makefile +++ b/src/xml/Makefile @@ -43,13 +43,14 @@ clean: for dir in $(SDIRS) ; do \ make -C $$dir clean ; \ done - rm *.a + -rm -f $(LIB) veryclean: cleandocs for dir in $(SDIRS) ; do \ make -C $$dir veryclean ; \ done - rm *.a + -rm -f $(LIB) + depend: for dir in $(SDIRS) ; do \ make -C $$dir depend ; \ diff --git a/src/xml/quakeML/Makefile b/src/xml/quakeML/Makefile index 2264b0b0..0ed247a1 100644 --- a/src/xml/quakeML/Makefile +++ b/src/xml/quakeML/Makefile @@ -13,7 +13,7 @@ principalAxes.o time.o DEBUG = -g CFLAGS += -O0 -D_REENTRANT -Dstatic_config -INCL = -I ../../../include $(XML_INCL) $(ISCL_INCL) $(COMPEARTH_INCL) +INCL = -I ../../../include $(XML2_INCL) $(ISCL_INCL) $(COMPEARTH_INCL) SYSLIBS = -lpthread -lnsl -lrt -lz diff --git a/src/xml/quakeML/epoch2string.c b/src/xml/quakeML/epoch2string.c index bcd3fcf0..9c4d82bf 100644 --- a/src/xml/quakeML/epoch2string.c +++ b/src/xml/quakeML/epoch2string.c @@ -33,7 +33,7 @@ int xml_epoch2string(double epoch, char cepoch[128]) } sec = (double) nzsec + (double) nzmusec*1.e-6; sprintf(cepoch, - "%04d-%02d-%02dT%02d:%02d:%04.3fZ", + "%04d-%02d-%02dT%02d:%02d:%06.3fZ", nzyear, month, mday, nzhour, nzmin, sec); return 0; } diff --git a/src/xml/shakeAlert/Makefile b/src/xml/shakeAlert/Makefile index e99f8848..bf213b05 100644 --- a/src/xml/shakeAlert/Makefile +++ b/src/xml/shakeAlert/Makefile @@ -11,7 +11,7 @@ OBJS = coreInfo.o segment.o slip.o vertex.o vertices.o geometry.o DEBUG = -g CFLAGS += -O0 -D_REENTRANT -Dstatic_config -INCL = -I ../../../include $(XML_INCL) $(ISCL_INCL) +INCL = -I ../../../include $(XML2_INCL) $(ISCL_INCL) SYSLIBS = -lpthread -lnsl -lrt -lz