diff --git a/DESCRIPTION b/DESCRIPTION index 347aaf9..eafbb5f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,23 +1,19 @@ Type: Package Package: IRISSeismic -Version: 1.0.7-1 +Version: 1.3.8 Title: Classes and Methods for Seismic Data Analysis Authors@R: c( - person("Jonathan", "Callahan", email="jonathan.s.callahan@gmail.com", role=c("aut","cre")), + person("Jonathan", "Callahan", email="jonathan@mazamascience.com", role=c("aut")), person("Rob", "Casey", email="rob@iris.washington.edu", role=c("aut")), + person("Gillian", "Sharer", email="gillian@iris.washington.edu", role=c("aut","cre")), person("Mary", "Templeton", email="met@iris.washington.edu", role=c("aut")) ) -Author: Jonathan Callahan [aut, cre], - Rob Casey [aut], - Mary Templeton [aut] -Maintainer: Jonathan Callahan Depends: R (>= 3.1.0) -Imports: methods, pracma, RCurl, seismicRoll (>= 1.0.1), signal, - stringr, XML -Suggests: mapdata, maps +Imports: methods, pracma, RCurl, seismicRoll (>= 1.1.0), signal, + stringr, XML, stats Description: Provides classes and methods for seismic data analysis. The base classes and methods are inspired by the python code found in - the 'ObsPy' python toolbox . Additional classes and + the 'ObsPy' python toolbox . Additional classes and methods support data returned by web services provided by the 'IRIS DMC' . Collate: Class-Trace.R Class-Stream.R Class-IrisClient.R @@ -25,5 +21,10 @@ Collate: Class-Trace.R Class-Stream.R Class-IrisClient.R License: GPL (>= 2) Repository: CRAN NeedsCompilation: yes -Packaged: 2016-04-20 11:32:51 UTC; ligges -Date/Publication: 2016-04-20 13:33:46 +Packaged: 2016-11-30 15:10:48 UTC; mustang +Author: Jonathan Callahan [aut], + Rob Casey [aut], + Gillian Sharer [aut, cre], + Mary Templeton [aut] +Maintainer: Gillian Sharer +Date/Publication: 2016-11-30 19:51:35 diff --git a/MD5 b/MD5 index 76809d5..039d14d 100644 --- a/MD5 +++ b/MD5 @@ -1,67 +1,68 @@ -0a2db36605c1bd4aec5ab930f09275f2 *DESCRIPTION -f3db61db24e4acafaecd0df30671bc35 *NAMESPACE -d05c6156cc3359a372fdb00b04c1df09 *R/Class-IrisClient.R -96888cfa7f9fb5374092099e05e40f84 *R/Class-Stream.R -b5a0919c1419a7f0322f6d17f09ca8a9 *R/Class-Trace.R -2e1576ddf7cdbdf4654c8e853be8ff97 *R/Utils.R +6589b0833be626e370aec6d01d76ad79 *DESCRIPTION +e6c5c7a1e1c83a4e0a74f1f333266f9e *NAMESPACE +32603fc83505c1c8734db9416a957d5c *R/Class-IrisClient.R +60e3778ead5f65ae20fccaa14cd60c2d *R/Class-Stream.R +c6d6dfc9a5d2b2f5f6409afee552c75c *R/Class-Trace.R +7cff554222174b3f96084700d035051c *R/Utils.R 17258addfe8d9ef274ac38a1964572af *R/mseedWrappers.R -ebe934c163b2950cc4edb70cc24f809d *R/spectralUtils.R -2bdaabd8652ec7d78efb9f3fcfcba6db *build/vignette.rds -c6aca87b5dcb817893f218635ced5d15 *inst/doc/IRISSeismic-intro.R -aded8268940353e7a2cd77ee5e71172f *inst/doc/IRISSeismic-intro.Rnw -c50b5c164b371ffcad59706c8186413c *inst/doc/IRISSeismic-intro.pdf -5d00e8132cc6e4aaa72f875f337975b8 *man/DDT.Rd -00ef5f4b3570ebfb78313db01963e617 *man/IRISSeismic-package.Rd -1a1756bb8840dc2a89e8b9992c19e33a *man/IrisClient-class.Rd -ef4aba79e6f8de037c8b93727566a150 *man/McNamaraBins.Rd -f8991a39ebd0471b35ca79eb85af1264 *man/McNamaraPSD.Rd -8dc030d9ce358fc0564179cf5527d9b8 *man/STALTA.Rd -b5c399462c1be9adf153ed7475143de6 *man/Stream-class.Rd -75f09c997d62d94fbdc87e4859e88a92 *man/Trace-class.Rd -a7bb3ac4f64dc1127202ca98484186bd *man/TraceHeader-class.Rd -ecc8d043395040ab335ee12989341e48 *man/basicStats.Rd -d070d66cfa363b2bb855631b6acd56a2 *man/butterworth.Rd -d01fee1bb1ab8d9d23038a6ba6948bbd *man/crossSpectrum.Rd -4c25aabc4b61662e0011d94be83d7467 *man/envelope.Rd -7d3826005d2c05186152aab9b207a8ab *man/eventWindow.Rd -5bb75f67b71bba12f64711ec06dab439 *man/getAvailability.Rd -b3b2367ea5c2a03fed9e5845091829bb *man/getChannel.Rd -b01c2163afea03230ac4f93bd7507f73 *man/getDataselect.Rd -895e48bb9d6b86fe55588b95f70810fa *man/getDistaz.Rd -4180e61f4a12d9f3104e7fcf26d06d40 *man/getEvalresp.Rd -960292543e332a672fbadc682fb37cb8 *man/getEvent.Rd -8b80afe9584d46a5cab0c36c92625b3a *man/getGaps.Rd -a1ff53c4bd11eb424efcb065972db1c6 *man/getNetwork.Rd -d3af24d58782978111bc844957272315 *man/getRotation.Rd -2ebe8c3ea8dfeb2c607f947e4dddc470 *man/getSNCL.Rd -0be75e77ab8c2619b36e5e1eafec7d8d *man/getStation.Rd -7b82ade8e16a123fff724d294bc310ca *man/getTraveltime.Rd -c49d536638f40f877614035395fb3b69 *man/getUnavailability.Rd -ac25b9f53329d2a6db66fb771d6d7352 *man/getUpDownTimes.Rd -a0c0930071c1c5d9d4207c8996c5bb0f *man/hilbert.Rd -930d02f5ad846632fed3d75253cb7ffe *man/hilbertFFT.Rd -f4e851e97b855504ee4809f05acd71dc *man/mergeTraces.Rd -9697fa479d6041fc5116a2ff8639ff18 *man/mergeUpDownTimes.Rd -d4e59f42c37b6e67814aacf4394f4fab *man/miniseed2Stream.Rd -6208cd695ab95772b25c979494c14547 *man/multiplyBy.Rd -8e9e7bcbac66ccb1dbe92efe38a15900 *man/noiseMatrix2PdfMatrix.Rd -bb7ccd0c30cabf046d7cdf1a225b6b4d *man/noiseModels.Rd -3494ed5914822b126cd0da14d1a223f8 *man/psdDF2NoiseMatrix.Rd -6baeadc79229e727c296e06ee18c059f *man/psdList.Rd -4cee794305058d7396b9c1e3417d87a8 *man/psdList2NoiseMatrix.Rd -30d28350d8cbabc7e4440ecbf0252d44 *man/psdPlot.Rd -12479d5f50714f37746cae158b5858c2 *man/psdStatistics.Rd -a6b1a396b2aad51f650d1b009ae730b0 *man/readMiniseedFile.Rd -91bb9b1f1e8002a6d63ae6a1a79e523f *man/rmsVariance.Rd -22dbad60d2341d359f70c894423983e5 *man/rotate2D.Rd -e0f59f43cef71ca2a338a29a11deebf1 *man/slice.Rd -6e021a2ef50306a6ef0e8b6b0ab6fb5d *man/surfaceDistance.Rd -7030e70919e3392f3f6031226bb33100 *man/triggerOnset.Rd +fef312fbfd82d028c65c2f873fde8b81 *R/spectralUtils.R +b5ee3e84ae0cdb78e3951e6a251c9c14 *build/vignette.rds +57ff67647c9abdd110b756ff92021d4d *inst/doc/IRISSeismic-intro.R +d1b5310214d5df335c6ac77349c09135 *inst/doc/IRISSeismic-intro.Rnw +bbe054d4cfe80a0c4549143ad76e6799 *inst/doc/IRISSeismic-intro.pdf +1c05b5d2bdb4830a01b97aad40abbb13 *man/DDT.Rd +15c83d7a02d07937d74cc1c8e4ecda3d *man/IRISSeismic-package.Rd +41b6577888b6e2318255d2dcbf6683ed *man/IrisClient-class.Rd +f80c9b055055da33fcac654b6c6be8d1 *man/McNamaraBins.Rd +5d32e7c88c0c4d4867d0d9108951d448 *man/McNamaraPSD.Rd +62d597b9242ed9a6a359b3db00f5551e *man/STALTA.Rd +7de125372df93a2e2d5ea33aeb79c60b *man/Stream-class.Rd +c8127fde0f4220d49c4ca6e84c99fb4c *man/Trace-class.Rd +202b3c4e04533a05ead9ffcdc51ee396 *man/TraceHeader-class.Rd +99af400523644812ac913f6454270a79 *man/basicStats.Rd +094dce92b83243166a125d85af9f5255 *man/butterworth.Rd +7bd7183fd51f46c477c9473bc51c54b1 *man/crossSpectrum.Rd +7128cf90d133551be89e3d12cf11fdd9 *man/envelope.Rd +7e2985d7766c1933087f4edfd4a82702 *man/eventWindow.Rd +0bf642e2447ac375fc1a4e4faa82660c *man/getAvailability.Rd +e1a9e08c30fe90c620ba4367044627d7 *man/getChannel.Rd +1874bd9b32008c8316b466fa4cb9a63b *man/getDataselect.Rd +cfdc186a50801d563b51e46a0b0ad5cc *man/getDistaz.Rd +399e0c9fa7ff1f4c2f5b9c2b8ca2f9f8 *man/getEvalresp.Rd +9ebda1bdc7a58639e70b8e50d3f6d9f7 *man/getEvent.Rd +c87c91e732533d6af8661c804604104e *man/getGaps.Rd +859f8de27ead8470236766db5fd15635 *man/getNetwork.Rd +3019f5cbdf76f0df99ee4ede70456cc6 *man/getRotation.Rd +3ab2ce7891d2fd5f9c43531fb656e068 *man/getSNCL.Rd +ffabdf622a814c4d02b47e6cb3441c85 *man/getStation.Rd +35a1c0c650f6b53525a7ad4966bd96cf *man/getTraveltime.Rd +5ea2962e4058cd547b92fe43428e966c *man/getUnavailability.Rd +a064b8a19d7c8e2fe818840f0dbf3981 *man/getUpDownTimes.Rd +c4e314f24b03ab7b427669566949f8cf *man/hilbert.Rd +6746e30f4a7e52ef37e27ab71ab12321 *man/hilbertFFT.Rd +2759015a1931e35ddd07e98d8d4ef220 *man/mergeTraces.Rd +aed38b25f4d1f39932c7a24dfdfb710b *man/mergeUpDownTimes.Rd +88b5cd03e54adc3bafe57d727e5128e1 *man/miniseed2Stream.Rd +7ec09f620a8b108bd4510f654b22d383 *man/multiplyBy.Rd +c55e685f4b033015419f6c70eced55b6 *man/noiseMatrix2PdfMatrix.Rd +45ae151d41695392687d7305b18f2636 *man/noiseModels.Rd +2c60f9b878e0bbc538c6f77aa3362430 *man/psdDF2NoiseMatrix.Rd +aa9fa1bf06babcc0f1f3b9b927d7e920 *man/psdList.Rd +8a362080f1c3931419dcfffc9b744b1f *man/psdList2NoiseMatrix.Rd +69c3c8909192d25c5b1df9d37e6b629e *man/psdPlot.Rd +68677df379b36438c6ee42cd6c18dc2d *man/psdStatistics.Rd +c2fdc12a1e3ba48a79ed9cb06b92aab8 *man/readMiniseedFile.Rd +c3b7520282dd1194f278f33f58da4980 *man/rmsVariance.Rd +c946730b94c676ab23a043a0f7529574 *man/rotate2D.Rd +11f38bcbdf2de4f8172be36948e812a4 *man/slice.Rd +9e988d5a3e78b79320e1998fa32ecadf *man/surfaceDistance.Rd +ecf3932e9a72599a45948143579e5469 *man/triggerOnset.Rd +5502755d9d3bef1e349826e6234b92e8 *man/unHistogram.Rd a9e87b98c746ba1a4dc4347fcc234e9c *src/Makevars -269521d42ebdf71be3abd6e5894a5570 *src/Makevars.win +b77488fe6e1ee245a26626ed67a6422e *src/Makevars.win 1b54014c515969ebbe4b908d2140b259 *src/libmseed/ChangeLog 740d6e91c899fa745cd636e5c913ff98 *src/libmseed/INSTALL -6f47978e3e363d5c9f54ba9dff56c83f *src/libmseed/Makefile +94fe255bf645ad896bce2f33e96d5f5d *src/libmseed/Makefile 67786c1ff8b8e840c6a2a1df221fa99e *src/libmseed/Makefile.wat c57ac4c8e7802072770d822557c3accb *src/libmseed/Makefile.win 6d0adb0d969ecd85d4c88c1236eaa86f *src/libmseed/README @@ -177,11 +178,11 @@ ecd4a381f9527662da568974ff46cb55 *src/libmseed/doc/mstl_init.3 baf046a5e6d6120d7e6a363f43dade7e *src/libmseed/doc/mstl_printgaplist.3 baf046a5e6d6120d7e6a363f43dade7e *src/libmseed/doc/mstl_printsynclist.3 baf046a5e6d6120d7e6a363f43dade7e *src/libmseed/doc/mstl_printtracelist.3 -8b684469667f089fc11144eaa83a4952 *src/libmseed/fileutils.c -1ea371f13659740afe88c2669fba0953 *src/libmseed/genutils.c +62f4d747ff89c8416458e36c43c76395 *src/libmseed/fileutils.c +34e41cd8c7243005f630a087d9685a57 *src/libmseed/genutils.c d838f8fc0f2ccfb59d7a439f81bbbc66 *src/libmseed/gswap.c 9e6a0b2e684c4431320a44992471cd70 *src/libmseed/libmseed.def -c7a3491c8076fae8a6ef9f7e1f534a65 *src/libmseed/libmseed.h +c2acf58d2daa1dd082d5fb9a924ec134 *src/libmseed/libmseed.h 408dd2ed162444f22f318447bdf169af *src/libmseed/lmplatform.c c69c2e5230d846272039211b1589ce5d *src/libmseed/lmplatform.h 4ed474573168d42285bb568c8d2d80db *src/libmseed/logging.c @@ -197,12 +198,12 @@ d010af6f2b9848481039b46e048896c2 *src/libmseed/matlab/mexMsReadTracesNative.c 3fa405182f0d326c30ba8703452ac29d *src/libmseed/parseutils.c 8dd84c5cd97ccd86b758ba62fee03333 *src/libmseed/selection.c fec2e84c442061b782be4234324c682d *src/libmseed/steimdata.h -3e8b1a516abc139284317b19da4766e2 *src/libmseed/tracelist.c +9ba1b0d7ac9d19aefbcd525ff02b0f8e *src/libmseed/tracelist.c 30bb50fb0ea36284f31ed337f399552b *src/libmseed/traceutils.c 894b2846cff735b74b81668008dfe65d *src/libmseed/unpack.c -54eb1cb506a369169c2e5edb9bde0fcb *src/libmseed/unpackdata.c +e7e89691f803198f1a4ff71e314cd5aa *src/libmseed/unpackdata.c 781c74c9e80b88a5cad3f8748f8fa3fb *src/libmseed/unpackdata.h -9f69e19ec3626e72924f20c34a1f4c2c *src/parseMiniSEED.c -aded8268940353e7a2cd77ee5e71172f *vignettes/IRISSeismic-intro.Rnw +154373f0480d13d7ba6cfdb70c739821 *src/parseMiniSEED.c +d1b5310214d5df335c6ac77349c09135 *vignettes/IRISSeismic-intro.Rnw 7c30695653c4728536070dbf0099f699 *vignettes/rstudio-IRISSeismic.png 236b60893fcdcc65318ace39d042c44f *vignettes/rstudio-windows.png diff --git a/NAMESPACE b/NAMESPACE index f6dea4c..efc7a53 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ ############################################################ import(methods) # for methods used in defining S4 classes +import(stringr) ############################################################ # Load dynamic libraries (shared object files) @@ -110,7 +111,8 @@ export("crossSpectrum", "psdList", "psdList2NoiseMatrix", "psdPlot", - "psdStatistics" + "psdStatistics", + "unHistogram" ) diff --git a/R/Class-IrisClient.R b/R/Class-IrisClient.R index be6dc67..0d04879 100644 --- a/R/Class-IrisClient.R +++ b/R/Class-IrisClient.R @@ -19,7 +19,7 @@ ## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ################################################################################ -# R classes for an IRIS web servivces request client. +# R classes for an IRIS web services request client. # # Inspiration for the functionality in this code comes from obspy.iris.client.Client # @@ -62,26 +62,6 @@ # ################################################################################ -#setClass("IrisClient", -# # typed slots (aka attributes) for class IrisClient -# representation(site = "character", -# user = "character", -# password = "character", -# timeout = "integer", -# debug = "logical", -# useragent = "character"), -# # default values for slots -# prototype(site = "http://service.iris.edu", -# user = "", -# password = "", -# timeout = as.integer(10), -# debug = FALSE, -# useragent = paste("R package IRISSeismic 0.1-3,", R.version.string, ",", version$platform) -# ) -#) - -# For now we'll just include the basics: - # REC - check for a user R profile for the IRIS Client site URL # otherwise, use the default URL irisSite <- "http://service.iris.edu" @@ -112,15 +92,15 @@ setClass("IrisClient", useragent = "character"), prototype(site = irisSite, debug = FALSE, - useragent = paste0("IRISSeismic ", + useragent = paste0("IRISSeismic/", ifelse("IRISSeismic" %in% rownames(installed.packages()), installed.packages()["IRISSeismic","Version"], "local code"), - ", RCurl ", + " RCurl/", ifelse("RCurl" %in% rownames(installed.packages()), installed.packages()["RCurl","Version"], "local code"), - ", ", R.version.string, - ", OS ", version$platform, - ", (",irisUserAgent,")")) + " R/", R.version$major,".",R.version$minor, + " ", version$platform, + " (",irisUserAgent,")")) ) @@ -135,12 +115,12 @@ setClass("IrisClient", ################################################################################ if (!isGeneric("getDataselect")) { - setGeneric("getDataselect", function(obj, network, station, location, channel, starttime, endtime, quality, ignoreEpoch) { + setGeneric("getDataselect", function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) { standardGeneric("getDataselect") }) } -getDataselect.IrisClient <- function(obj, network, station, location, channel, starttime, endtime, quality, ignoreEpoch) { +getDataselect.IrisClient <- function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) { # allow authenticated access. # if we have a netrc definition, then use queryauth to access data @@ -159,37 +139,67 @@ getDataselect.IrisClient <- function(obj, network, station, location, channel, s url <- paste(url,"&loc=", stringr::str_replace(location," ","--"),sep="") url <- paste(url,"&cha=", channel,sep="") url <- paste(url,"&start=", format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") - url <- paste(url,"&end=", format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") - url <- paste(url,"&quality=", ifelse(missing(quality),"B",quality),sep="") - + if (! inclusiveEnd) { + endtime <- endtime-0.000001 + url <- paste(url,"&end=", format(endtime,"%Y-%m-%dT%H:%M:%OS6", tz="GMT"),sep="") + } else { + url <- paste(url,"&end=", format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") + } + if (!is.null(quality)){ + url <- paste(url,"&quality=",quality,sep="") + } + if (obj@debug) { write(paste("URL =",url), stdout()) } # Make authenticated request using a netrc file if (! is.null(irisNetrc)) { - result <- try( dataselectResponse <- RCurl::getBinaryURL(url, useragent=obj@useragent, - netrc=1, netrc.file=irisNetrc), - silent=TRUE) - } else { - result <- try( dataselectResponse <- RCurl::getURLContent(url, useragent=obj@useragent), + h <- RCurl::basicTextGatherer() + result <- try( dataselectResponse <- RCurl::getBinaryURL(url, useragent=obj@useragent, + netrc=1, netrc.file=irisNetrc, .opts = list(headerfunction = h$update, followlocation=TRUE)), silent=TRUE) + # Handle error response + if (class(result) == "try-error" ) { + + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("Bad Request",ignore_case=TRUE))) { + stop(paste("getDataselect.IrisClient: Bad Request:",url)) + } else if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getDataselect.IrisClient: Not Found:",url)) + } else { + stop(paste("getDataselect.IrisClient:",err_msg, url)) + } + } + + header <- RCurl::parseHTTPHeader(h$value()) + } else { + result <- try( dataselectResponse <- RCurl::getURLContent(url, useragent=obj@useragent, header=TRUE), silent=TRUE) - } - - # Handle error response - if (class(result) == "try-error" ) { - - err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"Bad Request")) { - stop(paste("getDataselect.IrisClient: Bad Request:",url)) - } else if (stringr::str_detect(err_msg,"Not Found")) { - stop(paste("getDataselect.IrisClient: Not Found:",url)) - } else { - stop(paste("getDataselect.IrisClient:",err_msg)) - } - - } - + # Handle error response + if (class(result) == "try-error" ) { + + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("Bad Request",ignore_case=TRUE))) { + stop(paste("getDataselect.IrisClient: Bad Request:",url)) + } else if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getDataselect.IrisClient: URL Not Found:",url)) + } else { + stop(paste("getDataselect.IrisClient:",err_msg,url)) + } + } + + header <- dataselectResponse$header + dataselectResponse <- dataselectResponse$body + } + + if (header["status"] == "204") { #fdsnws dataselect returned nothing + stop(paste("getDataselect.IrisClient: No Data:",header["status"],url)) + } + + if (header["status"] != "200") { #fdsnws dataselect returned something unexpected, not caught by error + stop(paste("getDataselect.IrisClient: Unexpected http status code:",header["status"],url)) + } + # No errors so proceed # Channel metadata is required to properly apply InstrumentSensitivity corrections @@ -200,10 +210,10 @@ getDataselect.IrisClient <- function(obj, network, station, location, channel, s if (class(result) == "try-error" ) { err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"Bad Request")) { + if (stringr::str_detect(err_msg, regex("Bad Request",ignore_case=TRUE))) { stop(paste("getDataselect.IrisClient: Bad Request:",url)) } else { - stop(paste("getDataselect.IrisClient:",err_msg)) + stop(paste("getDataselect.IrisClient:",err_msg, url)) } } @@ -212,7 +222,7 @@ getDataselect.IrisClient <- function(obj, network, station, location, channel, s # NOTE: each with a different scale or starttime. This still represents a single SNCL. # NOTE: What to do about multiple scales in the following getChannel request? # Solution 1: if ignoreEpoch==TRUE, then just take the first epoch presented - # http://service.iris.edu/fdsnws/station/1/query?net=H2&sta=H2O&loc=00&cha=LHZ&starttime=2001-02-28T18:29:44&endtime=2001-02-28T19:29:44&includerestricted=false&output=text&level=channel + # http://service.iris.edu/fdsnws/station/1/query?net=H2&sta=H2O&loc=00&cha=LHZ&starttime=2001-02-28T18:29:44&endtime=2001-02-28T19:29:44&includerestricted=false&format=text&level=channel sncls <- paste(channels$network,channels$station,channels$location,channels$channel) if (nrow(channels) > 1 && !ignoreEpoch) { stop(paste("getDataselect.IrisClient: Multiple epochs: getChannel returned",length(sncls),"records")) @@ -223,7 +233,8 @@ getDataselect.IrisClient <- function(obj, network, station, location, channel, s # No errors so proceed stream <- miniseed2Stream(dataselectResponse,url,starttime,endtime, - channelInfo$instrument,channelInfo$scale,channelInfo$scaleunits) + channelInfo$instrument,channelInfo$scale,channelInfo$scalefreq,channelInfo$scaleunits,channelInfo$latitude, + channelInfo$longitude,channelInfo$elevation,channelInfo$depth,channelInfo$azimuth,channelInfo$dip) return(stream) } @@ -232,64 +243,145 @@ getDataselect.IrisClient <- function(obj, network, station, location, channel, s setMethod("getDataselect", signature(obj="IrisClient", network="character", station="character", location="character", channel="character", starttime="POSIXct", endtime="POSIXct", - quality="character", ignoreEpoch="logical"), - function(obj, network, station, location, channel, starttime, endtime, quality, ignoreEpoch) + quality="character", inclusiveEnd="logical", ignoreEpoch="logical"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, - quality, ignoreEpoch)) -# quality="missing", use "B" for "Best", ignoreEpoch="missing", default to FALSE -setMethod("getDataselect", signature(obj="IrisClient", - network="character", station="character", location="character", + quality, inclusiveEnd, ignoreEpoch)) + +# quality="missing", default to NULL, ignoreEpoch="missing", default to FALSE, inclusiveEnd="missing", default to TRUE +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", channel="character", starttime="POSIXct", endtime="POSIXct", - quality="missing", ignoreEpoch="missing"), - function(obj, network, station, location, channel, starttime, endtime, quality, ignoreEpoch) + quality="missing", inclusiveEnd="missing", ignoreEpoch="missing"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, - quality="B", ignoreEpoch=FALSE)) + quality=NULL, inclusiveEnd=TRUE, ignoreEpoch=FALSE)) + # ignoreEpoch="missing", default to FALSE -setMethod("getDataselect", signature(obj="IrisClient", - network="character", station="character", location="character", +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", channel="character", starttime="POSIXct", endtime="POSIXct", - quality="character", ignoreEpoch="missing"), - function(obj, network, station, location, channel, starttime, endtime, quality, ignoreEpoch) + quality="character", inclusiveEnd="logical", ignoreEpoch="missing"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, - quality, ignoreEpoch=FALSE)) -# quality="missing", use "B" for "Best" -setMethod("getDataselect", signature(obj="IrisClient", - network="character", station="character", location="character", + quality, inclusiveEnd, ignoreEpoch=FALSE)) + +# inclusiveEnd="missing", default to TRUE +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", channel="character", starttime="POSIXct", endtime="POSIXct", - quality="missing", ignoreEpoch="logical"), - function(obj, network, station, location, channel, starttime, endtime, quality, ignoreEpoch) + quality="character", inclusiveEnd="missing", ignoreEpoch="logical"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, - quality="B", ignoreEpoch)) + quality, inclusiveEnd=TRUE, ignoreEpoch)) +# quality="missing", default to NULL +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", + channel="character", starttime="POSIXct", endtime="POSIXct", + quality="missing", inclusiveEnd="logical", ignoreEpoch="logical"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) + getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, + quality=NULL, inclusiveEnd, ignoreEpoch)) + +# inclusiveEnd="missing", default to TRUE, ignoreEpoch="missing", default to FALSE +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", + channel="character", starttime="POSIXct", endtime="POSIXct", + quality="character", inclusiveEnd="missing", ignoreEpoch="missing"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) + getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, + quality, inclusiveEnd=TRUE, ignoreEpoch=FALSE)) + +# quality="missing", default to NULL, inclusiveEnd="missing", default to TRUE +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", + channel="character", starttime="POSIXct", endtime="POSIXct", + quality="missing", inclusiveEnd="missing", ignoreEpoch="logical"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) + getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, + quality=NULL, inclusiveEnd=TRUE, ignoreEpoch)) + +# quality="missing", default to NULL, ignoreEpoch="missing", default to FALSE +setMethod("getDataselect", signature(obj="IrisClient", + network="character", station="character", location="character", + channel="character", starttime="POSIXct", endtime="POSIXct", + quality="missing", inclusiveEnd="logical", ignoreEpoch="missing"), + function(obj, network, station, location, channel, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) + getDataselect.IrisClient(obj, network, station, location, channel, starttime, endtime, + quality=NULL, inclusiveEnd, ignoreEpoch=FALSE)) ################################################################################ # getSNCL method is a convenience wrapper for getDataselect ################################################################################ if (!isGeneric("getSNCL")) { - setGeneric("getSNCL", function(obj, sncl, starttime, endtime, quality) { + setGeneric("getSNCL", function(obj, sncl, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) { standardGeneric("getSNCL") }) } -getSNCL.IrisClient <- function(obj, sncl, starttime, endtime, quality) { +getSNCL.IrisClient <- function(obj, sncl, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) { parts <- unlist(stringr::str_split(sncl,'\\.')) - return( getDataselect.IrisClient(obj, parts[1], parts[2], parts[3], parts[4], starttime, endtime, quality) ) + return( getDataselect.IrisClient(obj, parts[1], parts[2], parts[3], parts[4], starttime, endtime, quality, inclusiveEnd, ignoreEpoch) ) } - # All arguments specified setMethod("getSNCL", signature(obj="IrisClient", sncl="character", starttime="POSIXct", endtime="POSIXct", - quality="missing"), - function(obj, sncl, starttime, endtime, quality) - getSNCL.IrisClient(obj, sncl, starttime, endtime, quality)) + quality="character", inclusiveEnd="logical",ignoreEpoch="logical"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch)) + +# quality="missing", use "B" for "Best", ignoreEpoch="missing", default to FALSE, inclusiveEnd="missing", default to TRUE +setMethod("getSNCL", signature(obj="IrisClient", + sncl="character", starttime="POSIXct", endtime="POSIXct", + quality="missing", inclusiveEnd="missing",ignoreEpoch="missing"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality="B",inclusiveEnd=TRUE,ignoreEpoch=FALSE)) + +# ignoreEpoch="missing", default to FALSE +setMethod("getSNCL", signature(obj="IrisClient", + sncl="character", starttime="POSIXct", endtime="POSIXct", + quality="character", inclusiveEnd="logical",ignoreEpoch="missing"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch=FALSE)) + +# inclusiveEnd="missing", default to TRUE +setMethod("getSNCL", signature(obj="IrisClient", + sncl="character", starttime="POSIXct", endtime="POSIXct", + quality="character", inclusiveEnd="missing",ignoreEpoch="logical"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality,inclusiveEnd=TRUE,ignoreEpoch)) + # quality="missing", use "B" for "Best" setMethod("getSNCL", signature(obj="IrisClient", sncl="character", starttime="POSIXct", endtime="POSIXct", - quality="missing"), - function(obj, sncl, starttime, endtime, quality) - getSNCL.IrisClient(obj, sncl, starttime, endtime, quality="B")) + quality="missing",inclusiveEnd="logical",ignoreEpoch="logical"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality="B",inclusiveEnd,ignoreEpoch)) + +# inclusiveEnd="missing", default to TRUE, ignoreEpoch="missing", default to FALSE +setMethod("getSNCL", signature(obj="IrisClient", + sncl="character", starttime="POSIXct", endtime="POSIXct", + quality="character", inclusiveEnd="missing",ignoreEpoch="missing"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality,inclusiveEnd=TRUE,ignoreEpoch=FALSE)) + +# quality="missing", use "B" for "Best", inclusiveEnd="missing", default to TRUE +setMethod("getSNCL", signature(obj="IrisClient", + sncl="character", starttime="POSIXct", endtime="POSIXct", + quality="missing", inclusiveEnd="missing",ignoreEpoch="logical"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality="B",inclusiveEnd=TRUE,ignoreEpoch)) + +# quality="missing", use "B" for "Best", ignoreEpoch="missing", default to FALSE +setMethod("getSNCL", signature(obj="IrisClient", + sncl="character", starttime="POSIXct", endtime="POSIXct", + quality="missing", inclusiveEnd="logical",ignoreEpoch="missing"), + function(obj, sncl, starttime, endtime, quality,inclusiveEnd,ignoreEpoch) + getSNCL.IrisClient(obj, sncl, starttime, endtime, quality="B",inclusiveEnd,ignoreEpoch=FALSE)) + ################################################################################ @@ -409,7 +501,7 @@ getNetwork.IrisClient <- function(obj, network, station, location, channel, url <- paste(url,"&starttime=",format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&endtime=",format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&includerestricted=",ifelse(includerestricted,"true","false"),sep="") - url <- paste(url,"&output=text",sep="") + url <- paste(url,"&format=text",sep="") # Add optional geographic constraints if they are passed in if (!missing(latitude)) { @@ -440,24 +532,41 @@ getNetwork.IrisClient <- function(obj, network, station, location, channel, # Make webservice request # NOTE: Be sure to set na.strings="" as "NA" is a valid network name - #result <- try( DF <- utils::read.delim(url,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - # silent=TRUE) - # REC -- going with RCurl variation for proper file closure - result <- try( DF <- utils::read.delim(textConnection(RCurl::getURL(url,useragent=obj@useragent)),sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - silent=TRUE) - - # Handle error response + + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error" ) { - err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"204 No Content")) { - stop(paste("getNetwork.IrisClient: No Content:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getNetwork.IrisClient: couldn't connect to host")) + if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getNetwork.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getNetwork.IrisClient: Couldn't connect to host", url)) + } else if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getNetwork.IrisClient: Cannot open connection",url)) } else { - stop(paste("getNetwork.IrisClient:",err_msg)) - } - + stop(paste("getNetwork.IrisClient:",err_msg,url)) + } + } + + header <- gurlc$header + gurlc <- gurlc$body + + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getNetwork.IrisClient: Unexpected http status code", header["status"],url)) + } + if(length(gurlc) == 0) { gurlc <- "" } + + txtcon <- textConnection(gurlc) + on.exit(close(txtcon),add=TRUE) + result <- try( DF <- utils::read.delim(txtcon,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""),silent=TRUE) + + if (class(result) == "try-error" ) { + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getNetwork.IrisClient: Cannot open connection",url)) + } else { + stop(paste("getNetwork.IrisClient:",err_msg, url)) + } } # No errors so proceed @@ -528,8 +637,10 @@ getStation.IrisClient <- function(obj, network, station, location, channel, url <- paste(url,"&cha=",ifelse(channel=="","*",channel),sep="") url <- paste(url,"&starttime=",format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&endtime=",format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") - url <- paste(url,"&includerestricted=",ifelse(includerestricted,"true","false"),sep="") - url <- paste(url,"&output=text",sep="") + if(!is.null(includerestricted)) { + url <- paste(url,"&includerestricted=",includerestricted,sep="") + } + url <- paste(url,"&format=text",sep="") # Add optional geographic constraints if they are passed in if (!missing(latitude)) { @@ -560,22 +671,37 @@ getStation.IrisClient <- function(obj, network, station, location, channel, # Make webservice request # NOTE: Be sure to set na.strings="" as "NA" is a valid network name - #result <- try( DF <- utils::read.delim(url,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - # silent=TRUE) - # REC - going with RCurl variation to ensure proper file closure - result <- try( DF <- utils::read.delim(textConnection(RCurl::getURL(url,useragent=obj@useragent)),sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - silent=TRUE) - # Handle error response + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error" ) { + if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getStation.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getStation.IrisClient: couldn't connect to host", url)) + } else if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getStation.IrisClient: Cannot open connection",url)) + } else { + stop(paste("getStation.IrisClient:",err_msg, url)) + } + } + + header <- gurlc$header + gurlc <- gurlc$body + + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getStation.IrisClient: Unexpected http status code", header["status"],url)) + } + + if(length(gurlc) == 0) { gurlc <- "" } + txtcon <- textConnection(gurlc) + on.exit(close(txtcon),add=TRUE) + result <- try( DF <- utils::read.delim(txtcon,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""),silent=TRUE) if (class(result) == "try-error" ) { - err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"204 No Content")) { - stop(paste("getStation.IrisClient: No Content:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getStation.IrisClient: couldn't connect to host")) + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getStation.IrisClient: Cannot open the connection:",url)) } else { - stop(paste("getStation.IrisClient:",err_msg)) + stop(paste("getStation.IrisClient:",err_msg, url)) } } @@ -602,7 +728,7 @@ setMethod("getStation", signature(obj="IrisClient", getStation.IrisClient(obj, network, station, location, channel, starttime, endtime, includerestricted, latitude, longitude, minradius, maxradius)) -# includerestricted="missing", use FALSE +# includerestricted="missing", use NULL setMethod("getStation", signature(obj="IrisClient", network="character", station="character", location="character", channel="character", starttime="POSIXct", endtime="POSIXct", @@ -611,7 +737,7 @@ setMethod("getStation", signature(obj="IrisClient", function(obj, network, station, location, channel, starttime, endtime, includerestricted, latitude, longitude, minradius, maxradius) getStation.IrisClient(obj, network, station, location, channel, starttime, endtime, - FALSE, latitude, longitude, minradius, maxradius)) + includerestricted=NULL, latitude, longitude, minradius, maxradius)) ################################################################################ @@ -649,8 +775,10 @@ getChannel.IrisClient <- function(obj, network, station, location, channel, url <- paste(url,"&cha=",ifelse(channel=="","*",channel),sep="") url <- paste(url,"&starttime=",format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&endtime=",format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") - url <- paste(url,"&includerestricted=",ifelse(includerestricted,"true","false"),sep="") - url <- paste(url,"&output=text",sep="") + if (!is.null(includerestricted)) { + url <- paste(url,"&includerestricted=",includerestricted,sep="") + } + url <- paste(url,"&format=text",sep="") # Add optional geographic constraints if they are passed in if (!missing(latitude)) { @@ -683,26 +811,43 @@ getChannel.IrisClient <- function(obj, network, station, location, channel, # Make webservice request # NOTE: Be sure to set na.strings="" as "NA" is a valid network name - #result <- try( DF <- utils::read.delim(url,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - # silent=TRUE) - # REC -- going with RCurl variation to ensure proper file closure - result <- try( DF <- utils::read.delim(textConnection(RCurl::getURL(url,useragent=obj@useragent)),sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - silent=TRUE) - # Handle error response - if (class(result) == "try-error" ) { + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error" ) { err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"204 No Content")) { - stop(paste("getChannel.IrisClient: No Content:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getChannel.IrisClient: couldn't connect to host")) + if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getChannel.IrisClient: Couldn't connect to host", url)) + } else if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getChannel.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getChannel.IrisClient: Cannot open connection",url)) } else { - stop(paste("getChannel.IrisClient:",err_msg)) - } - + stop(paste("getChannel.IrisClient:",err_msg, url)) + } } - + + header <- gurlc$header + gurlc <- gurlc$body + + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getChannel.IrisClient: Unexpected http status code", header["status"],url)) + } + + if(length(gurlc) == 0) { gurlc <- "" } + txtcon <- textConnection(gurlc) + on.exit(close(txtcon),add=TRUE) + result <- try( DF <- utils::read.delim(txtcon,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), silent=TRUE) + + if (class(result) == "try-error" ) { + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getChannel.IrisClient: Cannot open connection",url)) + } else { + stop(paste("getChannel.IrisClient:",err_msg, url)) + } + } + # No errors so proceed # Convert " " location codes back into the "" that is used in the miniSEED record @@ -733,7 +878,7 @@ setMethod("getChannel", signature(obj="IrisClient", getChannel.IrisClient(obj, network, station, location, channel, starttime, endtime, includerestricted, latitude, longitude, minradius, maxradius)) -# includerestricted="missing", use FALSE +# includerestricted="missing", use NULL setMethod("getChannel", signature(obj="IrisClient", network="character", station="character", location="character", channel="character", starttime="POSIXct", endtime="POSIXct", @@ -742,7 +887,7 @@ setMethod("getChannel", signature(obj="IrisClient", function(obj, network, station, location, channel, starttime, endtime, includerestricted, latitude, longitude, minradius, maxradius) getChannel.IrisClient(obj, network, station, location, channel, starttime, endtime, - FALSE, latitude, longitude, minradius, maxradius)) + includerestricted=NULL, latitude, longitude, minradius, maxradius)) ################################################################################ @@ -784,7 +929,7 @@ getAvailability.IrisClient <- function(obj, network, station, location, channel, url <- paste(url,"&starttime=",format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&endtime=",format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&includerestricted=",ifelse(includerestricted,"true","false"),sep="") - url <- paste(url,"&output=text",sep="") + url <- paste(url,"&format=text",sep="") # Add optional geographic constraints if they are passed in if (!missing(latitude)) { @@ -819,26 +964,42 @@ getAvailability.IrisClient <- function(obj, network, station, location, channel, # Make webservice request # NOTE: Be sure to set na.strings="" as "NA" is a valid network name - #result <- try( DF <- utils::read.delim(url,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - # silent=TRUE) - # REC -- going with RCurl variation to ensure proper file closure - result <- try( DF <- utils::read.delim(textConnection(RCurl::getURL(url,useragent=obj@useragent)),sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), - silent=TRUE) - - # Handle error response - if (class(result) == "try-error" ) { - + + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error" ) { err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"204 No Content")) { - stop(paste("getAvailability.IrisClient: No Content:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getAvailability.IrisClient: couldn't connect to host")) + if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getAvailability.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getAvailability.IrisClient: couldn't connect to host", url)) + } else if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getAvailability.IrisClient: Cannot open connection",url)) } else { - stop(paste("getAvailability.IrisClient:",err_msg)) - } - + stop(paste("getAvailability.IrisClient:",err_msg, url)) + } } - + + header <- gurlc$header + gurlc <- gurlc$body + + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getAvailability.IrisClient: Unexpected http status code", header["status"],url)) + } + + if(length(gurlc) == 0) { gurlc <- "" } + txtcon <- textConnection(gurlc) + on.exit(close(txtcon), add=TRUE) + result <- try( DF <- utils::read.delim(txtcon,sep="|",col.names=colNames,colClasses=colClasses,na.strings=""), silent=TRUE) + + if (class(result) == "try-error" ) { + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getAvailability.IrisClient: Cannot open connection",url)) + } else { + stop(paste("getAvailability.IrisClient:",err_msg, url)) + } + } + # No errors so proceed # REC Feb 2014 @@ -1006,27 +1167,39 @@ getEvalresp.IrisClient <- function(obj, network, station, location, channel, tim } # Conversion of URL into a data frame is a single line with utils::read.table(). - # REC Jan 2014 -- generate a connection handle for url so we can directly close on error. - # -- use RCurl variation in read.table call so that failed connections close properly - #result <- try ( DF <- utils::read.table(con <- url(url), col.names=colNames), - # silent=TRUE ) - result <- try ( DF <- utils::read.table(textConnection(RCurl::getURL(url,useragent=obj@useragent)), col.names=colNames), - silent=TRUE ) + + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error") { + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getEvalresp.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getEvalresp.IrisClient: Cannot open connection",url)) + } else { + stop(paste("getEvalresp.IrisClient:", err_msg, url)) + } + } + + header <- gurlc$header + gurlc <- gurlc$body - # Handle error response + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getEvalresp.IrisClient: Unexpected http status code", header["status"],url)) + } + + if(length(gurlc) == 0) { gurlc <- "" } + txtcon <- textConnection(gurlc) + on.exit(close(txtcon), add=TRUE) + result <- try ( DF <- utils::read.table(txtcon, col.names=colNames), silent=TRUE ) if (class(result) == "try-error" ) { - # err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"Not Found")) { - stop(paste("getEvalresp.IrisClient: URL Not Found:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getEvalresp.IrisClient: couldn't connect to host")) + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getEvalresp.IrisClient: Cannot open connection",url)) } else { - stop(paste("getEvalresp.IrisClient:",err_msg)) - } - + stop(paste("getEvalresp.IrisClient:", err_msg, url)) + } } - + # No errors so proceed # Return the dataframe, guaranteeing that it is ordered by frequency @@ -1049,7 +1222,7 @@ setMethod("getEvalresp", signature(obj="IrisClient", ################################################################################ # getEvent method returns seismic event data from the event webservice: # -# http://service.iris.edu/fdsnws/event/1/ +# http://earthquake.usgs.gov/fdsnws/event/1/ # # TODO: The getEvent method could be fleshed out with a more complete list # TODO: of arguments to be used as ws-event parameters. @@ -1065,10 +1238,15 @@ if (!isGeneric("getEvent")) { getEvent.IrisClient <- function(obj, starttime, endtime, minmag, maxmag, magtype, mindepth, maxdepth) { - url <- paste(obj@site,"/fdsnws/event/1/query?",sep="") - url <- paste(url,"&starttime=",format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") + if(stringr::str_detect(obj@site,regex("service.*.iris.edu"))) { + url <- paste("http://earthquake.usgs.gov","/fdsnws/event/1/query?",sep="") + } else { + url <- paste(obj@site,"/fdsnws/event/1/query?",sep="") + } + + url <- paste(url,"starttime=",format(starttime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") url <- paste(url,"&endtime=",format(endtime,"%Y-%m-%dT%H:%M:%OS", tz="GMT"),sep="") - url <- paste(url,"&output=text",sep="") + url <- paste(url,"&format=text",sep="") # Add optional arguments if they are non-null if (!missing(minmag)) { @@ -1092,7 +1270,7 @@ getEvent.IrisClient <- function(obj, starttime, endtime, minmag, maxmag, magtype write(paste("URL =",url), stdout()) } - # http://service.iris.edu/fdsnws/event/1/query?starttime=2013-02-01T00:00:00&endtime=2013-02-02T00:00:00&minmag=5&output=text + # http://earthquake.usgs.gov/fdsnws/event/1/query?starttime=2013-02-01T00:00:00&endtime=2013-02-02T00:00:00&minmag=5&format=text # # #EventID | Time | Latitude | Longitude | Depth | Author | Catalog | Contributor | ContributorID | MagType | Magnitude | MagAuthor | EventLocationName # 4075900|2013-02-01T22:18:33|-11.12|165.378|10.0|NEIC|NEIC PDE|NEIC PDE-Q||MW|6.4|GCMT|SANTA CRUZ ISLANDS @@ -1104,33 +1282,42 @@ getEvent.IrisClient <- function(obj, starttime, endtime, minmag, maxmag, magtype "factor","factor","numeric","factor","character") # Conversion of URL into a data frame is a single line with read.table(). - #result <- try ( DF <- utils::read.table(url, sep="|", col.names=colNames, colClasses=colClasses), - # silent=TRUE ) - # REC - going with RCurl variation to ensure proper file closure - result <- try ( DF <- utils::read.table(textConnection(RCurl::getURL(url,useragent=obj@useragent)), sep="|", col.names=colNames, colClasses=colClasses), - silent=TRUE ) + + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error") { + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getEvent.IrisClient: Cannot open connection",url)) + } else if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getEvent.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getEvent.IrisClient: Couldn't connect to host", url)) + } else { + stop(paste("getEvent.IrisClient:",err_msg, url)) + } + } - # Handle error response + header <- gurlc$header + gurlc <- gurlc$body + + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getEvent.IrisClient: Unexpected http status code", header["status"],url)) + } + + if(length(gurlc) == 0) { gurlc <- "" } + txtcon <- textConnection(gurlc) + on.exit(close(txtcon), add=TRUE) + result <- try ( DF <- utils::read.table(txtcon, sep="|", quote="", col.names=colNames, colClasses=colClasses),silent=TRUE ) if (class(result) == "try-error" ) { - err_msg <- geterrmessage() - if (stringr::str_detect(err_msg, "cannot open the connection")) { - stop(paste("getEvent.IrisClient: No Data Found")) - # NOTE: Only "cannot open the connection" has been seen with fdsn web services - - } else if (stringr::str_detect(err_msg, "Not Found")) { - stop(paste("getEvent.IrisClient: URL Not Found:",url)) - } else if (stringr::str_detect(err_msg, "couldn't connect to host")) { - stop(paste("getEvent.IrisClient: couldn't connect to host")) - } else if (stringr::str_detect(err_msg, "No Data Found")) { - stop(paste("getEvent.IrisClient: No Data Found")) + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getEvent.IrisClient: Cannot open connection",url)) } else { - stop(paste("getEvent.IrisClient:",err_msg)) + stop(paste("getEvent.IrisClient:",err_msg,url)) } } - # TODO: improve error trapping so this last check is not needed # Last check to make sure DF is defined if (!exists('DF')) { stop(paste("getEvent.IrisClient: No Data Found")) @@ -1211,26 +1398,39 @@ getTraveltime.IrisClient <- function(obj, latitude, longitude, depth, staLatitud "numeric", "numeric", "numeric", "character", "character") # Conversion of URL into a data frame is a single line with read.table(). - #result <- try ( returnValue <- utils::read.table(url, skip=4, col.names=colNames, colClasses=colClasses), - # silent=TRUE ) - # REC -- going with RCurl variation for better file closure - result <- try ( returnValue <- utils::read.table(textConnection(RCurl::getURL(url,useragent=obj@useragent)), skip=4, col.names=colNames, colClasses=colClasses), - silent=TRUE ) + + result <- try(gurlc <- RCurl::getURLContent(url,useragent=obj@useragent,header=TRUE),silent=TRUE) + if (class(result) == "try-error" ) { + err_msg <- geterrmessage() + if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getTraveltime.IrisClient: Couldn't connect to host", url)) + } else if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { + stop(paste("getTraveltime.IrisClient: URL Not Found",url)) + } else if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getTraveltime.IrisClient: Cannot open connection",url)) + } else { + stop(paste("getTraveltime.IrisClient:",err_msg, url)) + } + } - # Handle error response + header <- gurlc$header + gurlc <- gurlc$body + + if (header["status"] != "200" && header["status"] != "204") { + stop(paste("getTraveltime.IrisClient: Unexpected http status code", header["status"],url)) + } + + if(length(gurlc) == 0) { gurlc <- "" } + txtcon <- textConnection(gurlc) + on.exit(close(txtcon), add=TRUE) + result <- try ( returnValue <- utils::read.table(txtcon, skip=4, col.names=colNames, colClasses=colClasses),silent=TRUE ) if (class(result) == "try-error" ) { - err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"Not Found")) { - stop(paste("getTraveltime.IrisClient: URL Not Found:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getTraveltime.IrisClient: couldn't connect to host")) - } else if (stringr::str_detect(err_msg,"No Data Found")) { - stop(paste("getTraveltime.IrisClient: No Data Found")) + if (stringr::str_detect(err_msg, regex("cannot open the connection",ignore_case=TRUE))) { + stop(paste("getTraveltime.IrisClient: Cannot open connection",url)) } else { - stop(paste("getTraveltime.IrisClient:",err_msg)) + stop(paste("getTraveltime.IrisClient:",err_msg, url)) } - } # No errors so proceed @@ -1288,12 +1488,12 @@ getDistaz.IrisClient <- function(obj, latitude, longitude, staLatitude, staLongi if (class(result) == "try-error" ) { err_msg <- geterrmessage() - if (stringr::str_detect(err_msg,"Not Found")) { + if (stringr::str_detect(err_msg, regex("Not Found",ignore_case=TRUE))) { stop(paste("getDistaz.IrisClient: URL Not Found:",url)) - } else if (stringr::str_detect(err_msg,"couldn't connect to host")) { - stop(paste("getDistaz.IrisClient: couldn't connect to host")) + } else if (stringr::str_detect(err_msg, regex("couldn't connect to host",ignore_case=TRUE))) { + stop(paste("getDistaz.IrisClient: Couldn't connect to host", url)) } else { - stop(paste("getDistaz.IrisClient:",err_msg)) + stop(paste("getDistaz.IrisClient:",err_msg, url)) } } @@ -1307,17 +1507,41 @@ getDistaz.IrisClient <- function(obj, latitude, longitude, staLatitude, staLongi # 0.0 # 14.90407 # - + # See this example for magic with the XML package: # - # http://www.omegahat.org/RSXML/gettingStarted.html + # http://www.omegahat.net/RSXML/gettingStarted.html - doc <- XML::xmlRoot(XML::xmlTreeParse(distazXml)) - namedValueStrings <- XML::xmlSApply(doc,XML::xmlValue) - DF <- as.data.frame(t(as.numeric(namedValueStrings))) - names(DF) <- names(namedValueStrings) + # IRISWS distaz xml output became more compicated on 2016/08/17. Simple leaf node parsing no longer works. + # new output for http://service.iris.edu/irisws/distaz/1/query?stalat=0.0&stalon=0.0&evtlat=15.0&evtlon=0.0 looks like: + # + # + # + # 6378137 + # 1/298.257223563 + # + # 15.0 + # 0.0 + # 0.0 + # 0.0 + # 180.0 + # 0.0 + # 14.90407 + # + + # old code: + # doc <- XML::xmlRoot(XML::xmlTreeParse(distazXml)) + # namedValueStrings <- XML::xmlSApply(doc,XML::xmlValue) + # DF <- as.data.frame(t(as.numeric(namedValueStrings))) + # names(DF) <- names(namedValueStrings) - # Return dataframe with rows ordered by travelTime + # new code: + xmlList <- XML::xmlToList(distazXml) + xmlNames <- names(xmlList) + xmlList <- c(xmlList["ellipsoid"],sapply(xmlList[xmlNames[! xmlNames %in% c("ellipsoid")]], as.numeric)) + DF <- as.data.frame(xmlList) + colnames(DF)[colnames(DF)=="ellipsoid..attrs"] <- "ellipsoid.name" + return(DF) } diff --git a/R/Class-Stream.R b/R/Class-Stream.R index 8cc664a..a238211 100644 --- a/R/Class-Stream.R +++ b/R/Class-Stream.R @@ -340,62 +340,66 @@ getGaps.Stream <- function(x, min_gap) { headers <- lapply(x@traces, slot, "stats") num_headers <- length(headers) - # Sanity check -- single sampling_rate sampling_rates <- sapply(headers, slot, "sampling_rate") - num_rates <- length(unique(round(sampling_rates,digits=4))) - if (num_rates > 1) { - stop(paste("getGaps.Stream:",num_rates,"unique sampling rates encountered in Stream.")) - } - sampling_rate <- sampling_rates[1] - + # Set up arrays for information about gaps/overlaps gaps <- numeric(num_headers+1) - nsamples <- integer(num_headers+1) - - # Set min_gap and make sure it is at least 1/sampling_rate - min_gap <- ifelse(is.null(min_gap), 1/sampling_rate, min_gap) - min_gap <- max(min_gap, 1/sampling_rate) - - # NOTE: The delta we calculate here has a valid datapoint at it's start and end. - # NOTE: For the purposes of calculating the correct number of missing points we - # NOTE: will calculate delta as the total time between points minus 1/sampling rate. - # NOTE: That is how many extra points could be shoved into this gap. - - # Initial gap (no overlap possible) - delta <- as.numeric(difftime(headers[[1]]@starttime, x@requestedStarttime, units="secs")) - 1/sampling_rate - if (abs(delta) > min_gap) { - gaps[1] <- delta - nsamples[1] <- as.integer(round(abs(delta) * sampling_rate)) - } else { - gaps[1] <- 0 - nsamples[1] <- 0 - } - - # Inter-trace gaps and overlaps - if (num_headers > 1) { - for ( i in seq(from=2, to=num_headers) ) { - h1 <- headers[[i-1]] - h2 <- headers[[i]] - delta <- difftime(h2@starttime, h1@endtime, units="secs") - 1/sampling_rate - if (abs(delta) > min_gap) { - gaps[i] <- delta - nsamples[i] <- as.integer(round(abs(delta) * sampling_rate)) - } else { - gaps[i] <- 0 - nsamples[i] <- 0 - } - } - } - - # Final gap (no overlap possible) - delta <- as.numeric(difftime(x@requestedEndtime, headers[[num_headers]]@endtime, units="secs")) - 1/sampling_rate - if (abs(delta) > min_gap) { - gaps[num_headers+1] <- delta - nsamples[num_headers+1] <- as.integer(round(abs(delta) * sampling_rate)) - } else { - gaps[num_headers+1] <- 0 - nsamples[num_headers+1] <- 0 - } + nsamples <- integer(num_headers+1) + + for ( i in seq(from=1, to=num_headers)) { + + # NOTE: The delta we calculate here has a valid datapoint at it's start and end. + # NOTE: For the purposes of calculating the correct number of missing points we + # NOTE: will calculate delta as the total time between points minus 1/sampling rate. + # NOTE: That is how many extra points could be shoved into this gap. + + if (i == 1) { + # Initial gap (no overlap possible) + sampling_rate <- sampling_rates[1] + # Set min_gap and make sure it is at least 1/sampling_rate + min_gap_new <- ifelse(is.null(min_gap), 1/sampling_rate, min_gap) + min_gap_new <- max(min_gap_new, 1/sampling_rate) + + delta <- as.numeric(difftime(headers[[1]]@starttime, x@requestedStarttime, units="secs")) - 1/sampling_rate + if (delta > min_gap_new - 0.5/sampling_rate) { + gaps[1] <- delta + 1/sampling_rate + nsamples[1] <- as.integer(round(gaps[1] * sampling_rate)) + } else { + gaps[1] <- 0 + nsamples[1] <- 0 + } + } else { + # Inter-trace gaps and overlaps + sampling_rate <- sampling_rates[[i-1]] + min_gap_new <- ifelse(is.null(min_gap), 1/sampling_rate, min_gap) + min_gap_new <- max(min_gap_new, 1/sampling_rate) + h1 <- headers[[i-1]] + h2 <- headers[[i]] + delta <- difftime(h2@starttime, h1@endtime, units="secs") - 1/sampling_rate + if (abs(delta) > min_gap_new - 0.5/sampling_rate) { + gaps[i] <- delta + nsamples[i] <- as.integer(round(abs(delta) * sampling_rate)) + } else { + gaps[i] <- 0 + nsamples[i] <- 0 + } + } + if (i == num_headers) { + # Final gap (no overlap possible) + sampling_rate <- sampling_rates[[i]] + min_gap_new <- ifelse(is.null(min_gap), 1/sampling_rate, min_gap) + min_gap_new <- max(min_gap_new, 1/sampling_rate) + delta <- as.numeric(difftime(x@requestedEndtime, headers[[num_headers]]@endtime, units="secs")) - 1/sampling_rate + if (delta > min_gap_new - 0.5/sampling_rate) { + gaps[num_headers+1] <- delta + nsamples[num_headers+1] <- as.integer(round(delta * sampling_rate)) + } else { + gaps[num_headers+1] <- 0 + nsamples[num_headers+1] <- 0 + } + } + + } gap_list <- list(gaps=gaps, nsamples=nsamples) @@ -579,15 +583,21 @@ mergeTraces.Stream <- function(x, fillMethod) { gapInfo <- getGaps(x) num_gaps <- length(gapInfo$nsamples) - # Sanity check + # Sanity checks if (num_gaps != num_traces+1) { stop(paste("mergeTraces.Stream: num_gaps (",num_gaps,") should be one more than num_traces (",num_traces,")", sep="")) } + + headers <- lapply(x@traces, slot, "stats") + sampling_rates <- sapply(headers, slot, "sampling_rate") + num_rates <- length(unique(round(sampling_rates,digits=4))) + #if (num_rates > 1) { + if ( ! all(stats::dist(unique(sampling_rates)) < 0.0002 )) { # slightly more forgiving criteria for acceptable sample rate jitter than round(sampling_rates,digits=4) + stop(paste("mergeTraces.Stream:",num_rates,"unique sampling rates encountered in Stream.")) + } - # NOTE: Increment totalPoints by one to account for point at the very end. - # NOTE: For example, the sequence *--*--*--* has 4 points but only three time intervals. totalSecs <- as.numeric(difftime(x@requestedEndtime, x@requestedStarttime, units="secs")) - totalPoints <- as.integer(round(totalSecs) * x@traces[[1]]@stats@sampling_rate) + 1 + totalPoints <- as.integer(round(totalSecs) * x@traces[[1]]@stats@sampling_rate) # NOTE: Setting up a list of vectors that we will concatenate in one fell swoop at the end. # NOTE: This avoids the incremental growth and copying that is incredibly slow. @@ -643,9 +653,10 @@ mergeTraces.Stream <- function(x, fillMethod) { id <- x@traces[[1]]@id Sensor <- x@traces[[1]]@Sensor InstrumentSensitivity <- x@traces[[1]]@InstrumentSensitivity + SensitivityFrequency <- x@traces[[1]]@SensitivityFrequency InputUnits <- x@traces[[1]]@InputUnits - traces <- list( new("Trace", id, stats, Sensor, InstrumentSensitivity, InputUnits, data=data[1:totalPoints]) ) + traces <- list( new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data=data[1:totalPoints]) ) return( new("Stream", url=x@url, requestedStarttime=x@requestedStarttime, requestedEndtime=x@requestedEndtime, act_flags=x@act_flags, io_flags=x@io_flags, dq_flags=x@dq_flags, timing_qual=x@timing_qual, diff --git a/R/Class-Trace.R b/R/Class-Trace.R index d1b82fe..27a9273 100644 --- a/R/Class-Trace.R +++ b/R/Class-Trace.R @@ -30,10 +30,10 @@ ################################################################################ # Class TraceHeader # -# A container for additional header informaion for each Trace object. +# A container for additional header information for each Trace object. # # A TraceHeader object contains all header information (aka metadata) -# associated with a Trace object. +# associated with a Trace object. # # Documentation for attributes copied verbatim from obspy.core.trace.py: # @@ -61,6 +61,15 @@ # Date and time of the last data sample given in UTC # (default value is "1970-01-01T00:00:00.0Z"). # +# It can also contain some optional metadata, not included in miniseed headers: +# +# ``latitude`` : numeric, optional +# ``longitude`` : numeric, optional +# ``elevation`` : numeric, optional +# ``depth`` : numeric, optional +# ``azimuth`` : numeric, optional +# ``dip`` : numeric, optional +# ################################################################################ setClass("TraceHeader", @@ -76,6 +85,12 @@ setClass("TraceHeader", quality = "character", starttime = "POSIXct", endtime = "POSIXct", + latitude = "numeric", + longitude = "numeric", + elevation = "numeric", + depth = "numeric", + azimuth = "numeric", + dip= "numeric", processing = "character"), # default values for slots prototype(sampling_rate = 1.0, @@ -88,11 +103,16 @@ setClass("TraceHeader", channel = "", quality = "X", starttime = as.POSIXct("1900-01-01T00:00:00",format="%Y-%m-%dT%H:%M:%OS", tz="GMT"), - endtime = as.POSIXct("1900-01-01T00:00:00",format="%Y-%m-%dT%H:%M:%OS", tz="GMT")) + endtime = as.POSIXct("1900-01-01T00:00:00",format="%Y-%m-%dT%H:%M:%OS", tz="GMT"), + latitude = as.numeric(NA), + longitude = as.numeric(NA), + elevation = as.numeric(NA), + depth = as.numeric(NA), + azimuth = as.numeric(NA), + dip= as.numeric(NA)) ) ################################################################################ -# TODO: Improve notes on TraceHeader initialization # # A TraceHeader can be initialized with a headerLine as returned by the IRIS DMC timeseries webservice: # @@ -105,7 +125,8 @@ setClass("TraceHeader", # TODO: I am currently using "Trace" to mean what libmseed calls a "Trace Segment". # TODO: Each libmseed "Trace" has both "numsamples" and "samplecnt". # TODO: Need to figure out what the difference is and when to bail with zero samples in the SNCL. -# TODO: Zero samples in a "Trace Segment" is not necessariliy a deal breaker. +# TODO: Zero samples in a "Trace Segment" is not necessarily a deal breaker. + setMethod("initialize", "TraceHeader", function(.Object, headerList=list(), @@ -124,13 +145,19 @@ setMethod("initialize", "TraceHeader", .Object@npts <- as.integer(strsplit(headerValues[2]," ")[[1]][1]) .Object@sampling_rate <- as.numeric(strsplit(headerValues[3]," ")[[1]][1]) .Object@starttime <- as.POSIXct(headerValues[4], format="%Y-%m-%dT%H:%M:%OS", tz="GMT") # header uses 'T' format + .Object@latitude <- as.numeric(NA) + .Object@longitude <- as.numeric(NA) + .Object@elevation <- as.numeric(NA) + .Object@depth <- as.numeric(NA) + .Object@azimuth <- as.numeric(NA) + .Object@dip <- as.numeric(NA) } # override with anything found in the headerList for (key in names(headerList)) { slot(.Object,key) = headerList[[key]] } - + # guarantee npts and sampling_rate are non-zero numbers if (is.na(.Object@npts)) { err_msg <- paste("initialize.TraceHeader: npts=", .Object@npts, sep="") @@ -161,7 +188,9 @@ setMethod("initialize", "TraceHeader", if (!is.finite(delta)) delta = 0 } - .Object@endtime = .Object@starttime + delta + if (is.na(.Object@endtime)) { + .Object@endtime = .Object@starttime + delta + } # return the newly created .Object return(.Object) @@ -184,6 +213,12 @@ show.TraceHeader <- function(object) { cat (" delta: " , object@delta, "\n") cat (" starttime: " , format(object@starttime), "\n") cat (" endtime: " , format(object@endtime), "\n") + cat (" latitude: " , format(object@latitude), "\n") + cat (" longitude: " , format(object@longitude), "\n") + cat (" elevation: " , format(object@elevation), "\n") + cat (" depth: " , format(object@depth), "\n") + cat (" azimuth: " , format(object@azimuth), "\n") + cat (" dip: " , format(object@dip), "\n") cat (" processing: " , paste(object@processing,collapse="; "), "\n") } # NOTE: method signature must match generic signature for 'show' with argument: 'object' @@ -224,6 +259,10 @@ setMethod("as.headerLine", signature(obj="TraceHeader"), function(obj) as.header # :var id: A SEED compatible identifier of the trace. # :var stats: A container :class:`~obspy.core.trace.Stats` for additional # header information of the trace. +# :var Sensor: +# :var InstrumentSensitivity: +# :var SensitivityFrequency: +# :var InputUnits: # :var data: Data samples in a :class:`~numpy.ndarray` or # :class:`~numpy.ma.MaskedArray` # @@ -235,6 +274,7 @@ setClass("Trace", stats = "TraceHeader", Sensor = "character", InstrumentSensitivity = "numeric", + SensitivityFrequency = "numeric", InputUnits = "character", data = "numeric"), # default values for slots @@ -242,17 +282,19 @@ setClass("Trace", stats = new("TraceHeader"), Sensor = "", InstrumentSensitivity = 1.0, + SensitivityFrequency = 1.0, InputUnits = "", data = c(0)) ) -# initialze method +# initialize method setMethod("initialize", "Trace", function(.Object, id="", stats=new("TraceHeader"), Sensor = "", InstrumentSensitivity = 1.0, + SensitivityFrequency = 1.0, InputUnits = "", data=c(0), ...) { @@ -261,7 +303,8 @@ setMethod("initialize", "Trace", .Object@stats <- stats .Object@Sensor <- Sensor .Object@InstrumentSensitivity <- InstrumentSensitivity - .Object@InputUnits <- InputUnits + .Object@SensitivityFrequency <- SensitivityFrequency + .Object@InputUnits <- InputUnits .Object@data <- data if (.Object@id == "") { @@ -409,11 +452,12 @@ multiplyBy.Trace <- function(x, y) { stats <- x@stats Sensor <- x@Sensor InstrumentSensitivity <- x@InstrumentSensitivity + SensitivityFrequency <- x@SensitivityFrequency InputUnits <- x@InputUnits data <- x@data * y stats@processing <- append(stats@processing,paste("multiply by",y)) - return( new("Trace", id, stats, Sensor, InstrumentSensitivity, InputUnits, data) ) + return( new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data) ) } setMethod("multiplyBy", signature("Trace", y="numeric"), function(x, y) multiplyBy.Trace(x, y=y)) @@ -436,6 +480,7 @@ DDT.Trace <- function(x, demean, detrend, taper) { stats <- x@stats Sensor <- x@Sensor InstrumentSensitivity <- x@InstrumentSensitivity + SensitivityFrequency <- x@SensitivityFrequency InputUnits <- x@InputUnits # NOTE: Use the pracma::detrend() function. @@ -462,7 +507,7 @@ DDT.Trace <- function(x, demean, detrend, taper) { } - return( new("Trace", id, stats, Sensor, InstrumentSensitivity, InputUnits, data) ) + return( new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data) ) } # All parameters specified @@ -491,6 +536,7 @@ butterworth.Trace <- function(x, n, low, high, type) { stats <- x@stats Sensor <- x@Sensor InstrumentSensitivity <- x@InstrumentSensitivity + SensitivityFrequency <- x@SensitivityFrequency InputUnits <- x@InputUnits # If neither demean nor detrend has ever been applied, demean the data @@ -527,7 +573,7 @@ butterworth.Trace <- function(x, n, low, high, type) { data <- as.numeric(ts_filtered) stats@processing <- append(stats@processing,paste0("Butterworth filtered with (",n,",",low,",",high,",",type,")")) - return( new("Trace", id, stats, Sensor, InstrumentSensitivity, InputUnits, data) ) + return( new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data) ) } # All parameters specified @@ -562,6 +608,7 @@ slice.Trace <- function(x, starttime, endtime) { stats <- x@stats Sensor <- x@Sensor InstrumentSensitivity <- x@InstrumentSensitivity + SensitivityFrequency <- x@SensitivityFrequency InputUnits <- x@InputUnits data <- x@data @@ -602,7 +649,7 @@ slice.Trace <- function(x, starttime, endtime) { if (sliced) { stats@npts <- as.integer(end_index - start_index + 1) stats@processing <- append(stats@processing,"slice") - return( new("Trace", id, stats, Sensor, InstrumentSensitivity, InputUnits, data=data[start_index:end_index]) ) + return( new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data=data[start_index:end_index]) ) } else { return(x) } @@ -870,7 +917,7 @@ setMethod("eventWindow", signature(x="Trace", picker="numeric", threshold="missi # NOTE: style functions rather than S4 methods. plot.Trace <- function(x, starttime=x@stats@starttime, endtime=x@stats@endtime, - subsampling=NULL, add=FALSE, ...) { + subsampling=NULL, add=FALSE, ylab="raw", xlab="", ...) { # Create subsampled indices if necessary max_size <- 10000 @@ -893,13 +940,11 @@ plot.Trace <- function(x, starttime=x@stats@starttime, endtime=x@stats@endtime, times <- seq(from=x@stats@starttime, to=x@stats@endtime, length.out=length(x@data)) if (length(x) == length(indices)) { - xlab <- paste("(",length(x), "points )") + xlab <- paste(xlab,"\n"," (",length(x), "points )") } else { - xlab <- paste("( ",length(x), " points, subsampling=", subsampling, " for this plot. )", sep="") + xlab <- paste(xlab, "\n"," ( ",length(x), " points, subsampling=", subsampling, " for this plot. )", sep="") } -### ylab <- paste(x@InputUnits) - ylab <- "raw" - + # # Set up the time range to plot # if (missing(starttime)) { # starttime <- x@stats@starttime diff --git a/R/Utils.R b/R/Utils.R index 3a353f4..6a0fb6f 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -30,7 +30,14 @@ miniseed2Stream <- function(miniseed, requestedEndtime=NULL, sensor="synthetic trace", scale=as.integer(NA), - scaleunits="") { + scalefreq=as.numeric(NA), + scaleunits="", + latitude=as.numeric(NA), + longitude=as.numeric(NA), + elevation=as.numeric(NA), + depth=as.numeric(NA), + azimuth=as.numeric(NA), + dip=as.numeric(NA)) { # Use C code to parse the bytes into a list of lists result <- try( segList <- parseMiniSEED(miniseed), @@ -41,12 +48,12 @@ miniseed2Stream <- function(miniseed, err_msg <- geterrmessage() if (stringr::str_detect(err_msg,"libmseed__zero traces in miniSEED record")) { - stop(paste("miniseed2Stream: No data found.")) + stop(paste("miniseed2Stream: No data found.", url)) } else if ((stringr::str_detect(err_msg,"Data integrity check for Steim"))) { # REC - May 2014 - allowing libmseed Steim coefficient errors to pass unabated - print(paste("miniseed2Stream:",err_msg,"- Using data anyway")) + print(paste("miniseed2Stream:",err_msg,"- Using data anyway",":",url)) } else { - stop(paste("miniseed2Stream:",err_msg)) + stop(paste("miniseed2Stream:",err_msg,":",url)) } } @@ -73,18 +80,27 @@ miniseed2Stream <- function(miniseed, channel=segList[[i]]$channel, quality=segList[[i]]$quality, starttime=as.POSIXct(segList[[i]]$starttime, origin=origin, tz="GMT"), + endtime=as.POSIXct(segList[[i]]$endtime, origin=origin, tz="GMT"), npts=segList[[i]]$npts, - sampling_rate=segList[[i]]$sampling_rate) + sampling_rate=segList[[i]]$sampling_rate, + latitude=latitude, + longitude=longitude, + elevation=elevation, + depth=depth, + azimuth=azimuth, + dip=dip + ) stats <- new("TraceHeader", headerList=headerList) traces[[i]] <- new("Trace", stats=stats, Sensor=sensor, InstrumentSensitivity=scale, + SensitivityFrequency=scalefreq, InputUnits=scaleunits, data=segList[[i]]$data) } - # Each minSEED record has one set of quality flags (currently attached to each element in segList) + # Each miniSEED record has one set of quality flags (currently attached to each element in segList) act_flags <- segList[[1]]$act_flags io_flags <- segList[[1]]$io_flags dq_flags <- segList[[1]]$dq_flags @@ -96,7 +112,7 @@ miniseed2Stream <- function(miniseed, } if (is.null(requestedEndtime)) { last <- length(traces) - requestedEndtime <- traces[[last]]@stats@starttime + traces[[last]]@stats@npts / traces[[last]]@stats@sampling_rate + requestedEndtime <- traces[[last]]@stats@endtime } # Create a new Stream object @@ -113,8 +129,18 @@ miniseed2Stream <- function(miniseed, # # Reads miniSEED bytes from a file and converts them into a Stream object. # -readMiniseedFile <- function(file, sensor, scale, scaleunits) { - +readMiniseedFile <- function(file, + sensor="synthetic trace", + scale=as.integer(NA), + scalefreq=as.numeric(NA), + scaleunits="", + latitude=as.numeric(NA), + longitude=as.numeric(NA), + elevation=as.numeric(NA), + depth=as.numeric(NA), + azimuth=as.numeric(NA), + dip=as.numeric(NA)) { + # Read in the binary data bytes <- file.info(file)$size miniseed <- readBin(file, "raw", n=bytes) @@ -125,7 +151,8 @@ readMiniseedFile <- function(file, sensor, scale, scaleunits) { requestedEndtime <- NULL stream <- miniseed2Stream(miniseed,url,requestedStarttime,requestedEndtime, - sensor,scale,scaleunits) + sensor,scale,scalefreq,scaleunits,latitude,longitude, + elevation,depth,azimuth,dip) return(stream) } @@ -134,34 +161,61 @@ readMiniseedFile <- function(file, sensor, scale, scaleunits) { ############################################################ # rotate2D # -# Rotates non-vertical compoentns of a seismic waveform into a new coordinate system. +# Rotates non-vertical components of a seismic waveform into a new coordinate system. # -# Functionas similarly to the rotation service with "&azimuth=angle&components=ZRT" +# Functions similarly to the rotation service with "&azimuth=angle&components=ZRT" +# expects orthogonal traces # # http://service.iris.edu/irisws/rotation/1/ # http://service.iris.edu/irisws/rotation/docs/1/help/ rotate2D <- function(st1,st2,angle) { - # Make sure that all stream data is in a single trace -# st1 <- mergeTraces(st1) -# st2 <- mergeTraces(st2) - # TODO: mergeTraces can introduce NA's - # TODO: Figure out how to properly merge. - if (length(st1@traces) > 1) { - stop(paste("rotate2D: Stream st1 has more than one trace.")) + stop(paste("Stream st1 has more than one trace.")) } if (length(st2@traces) > 1) { - stop(paste("rotate2D: Stream st2 has more than one trace.")) + stop(paste("Stream st2 has more than one trace.")) } - tr1 <- st1@traces[[1]] - tr2 <- st2@traces[[1]] - # Sanity check - if (length(tr1) != length(tr2)) { - stop(paste("rotate2D: Incoming streams have different data lengths.")) + if (length(st1) != length(st2)) { + stop(paste("Incoming streams have different data lengths.")) + } + + # need to determine which is 1 (lag, usually oriented North) and 2 (lead, usually oriented E) using metadata information + + if (!is.na(st1@traces[[1]]@stats@azimuth) && !is.na(st2@traces[[1]]@stats@azimuth)){ + az1 <- st1@traces[[1]]@stats@azimuth + az2 <- st2@traces[[1]]@stats@azimuth + + azdiff <- az1-az2 + + if ( !( (abs(azdiff) > 87 & abs(azdiff) < 93) || (abs(azdiff) > 267 & abs(azdiff) < 273) ) ) { + stop(paste("Incoming streams are not orthogonal (+/- 3 degrees):", st1@traces[[1]]@id,": azimuth=",st1@traces[[1]]@stats@azimuth,";",st2@traces[[1]]@id,": azimuth=",st2@traces[[1]]@stats@azimuth)) + } + + + if(azdiff < 0 ) { + if ( azdiff <= 180 ) { + tr1 <- st1@traces[[1]] + tr2 <- st2@traces[[1]] + } else { + tr1 <- st2@traces[[1]] + tr2 <- st1@traces[[1]] + } + }else if( azdiff > 0) { + if ( azdiff > 180 ) { + tr1 <- st1@traces[[1]] + tr2 <- st2@traces[[1]] + } else { + tr1 <- st2@traces[[1]] + tr2 <- st1@traces[[1]] + } + } + } else { #if azimuth information is not available, assume that the first trace is lagging and the second is leading. Assume orthogonality. + tr1 <- st1@traces[[1]] + tr2 <- st2@traces[[1]] } # NOTE: The azimuth circle is different from standard geometry!!! @@ -194,10 +248,16 @@ rotate2D <- function(st1,st2,angle) { parts[4] <- stringr::str_replace(parts[4],'.$','R') stR@traces[[1]]@id <- paste(parts,collapse='.') stR@traces[[1]]@stats@channel <- parts[4] - stR@traces[[1]]@Sensor <- paste("synthetic rotation by",formatC(angle,digits=1,format="f"),"degrees") - # TODO: Should we retain InstrumentSensitivity and InputUnits from N/E or 1/2 traces information - stR@traces[[1]]@InstrumentSensitivity <- as.integer(NA) - stR@traces[[1]]@InputUnits <- "" + if ( (tr1@stats@azimuth+angle) < 360 && (tr1@stats@azimuth+angle) >= 0) { + stR@traces[[1]]@stats@azimuth <- tr1@stats@azimuth+angle + } else if ((tr1@stats@azimuth+angle) > 360) { + stR@traces[[1]]@stats@azimuth <- (tr1@stats@azimuth+angle)-360 + } else if ((tr1@stats@azimuth+angle) < 0) { + stR@traces[[1]]@stats@azimuth <- (tr1@stats@azimuth+angle)+360 + } + stR@traces[[1]]@Sensor <- paste("synthetic rotation by",formatC(angle,digits=1,format="f"),"degrees",tr1@Sensor) + stR@traces[[1]]@InstrumentSensitivity <- tr1@InstrumentSensitivity + stR@traces[[1]]@InputUnits <- tr1@InputUnits # Create a new "transverse" Stream object @@ -208,10 +268,16 @@ rotate2D <- function(st1,st2,angle) { parts[4] <- stringr::str_replace(parts[4],'.$','T') stT@traces[[1]]@id <- paste(parts,collapse='.') stT@traces[[1]]@stats@channel <- parts[4] - stT@traces[[1]]@Sensor <- paste("synthetic rotation by",formatC(angle,digits=1,format="f"),"degrees") - # TODO: Should we retain InstrumentSensitivity and InputUnits from N/E or 1/2 traces information - stT@traces[[1]]@InstrumentSensitivity <- as.integer(NA) - stT@traces[[1]]@InputUnits <- "" + if ( (tr2@stats@azimuth+angle) < 360 && (tr2@stats@azimuth+angle) >= 0) { + stT@traces[[1]]@stats@azimuth <- tr2@stats@azimuth+angle + } else if ((tr2@stats@azimuth+angle) > 360) { + stT@traces[[1]]@stats@azimuth <- (tr2@stats@azimuth+angle)-360 + } else if ((tr2@stats@azimuth+angle) < 0) { + stT@traces[[1]]@stats@azimuth <- (tr2@stats@azimuth+angle)+360 + } + stT@traces[[1]]@Sensor <- paste("synthetic rotation by",formatC(angle,digits=1,format="f"),"degrees",tr2@Sensor) + stT@traces[[1]]@InstrumentSensitivity <- tr2@InstrumentSensitivity + stT@traces[[1]]@InputUnits <- tr2@InputUnits # Return the rotated data return(list(stR=stR, stT=stT)) diff --git a/R/spectralUtils.R b/R/spectralUtils.R index d214933..986df22 100644 --- a/R/spectralUtils.R +++ b/R/spectralUtils.R @@ -350,7 +350,7 @@ McNamaraPSD <- function(tr, loFreq=.005, hiFreq=10, alignFreq=0.1, binned=TRUE) # Average sp$spec <- specSum / 13 - # Bin specturm if requested ---------- + # Bin spectrum if requested ---------- if (binned) { df <- data.frame(freq=sp$freq,spec=sp$spec) @@ -399,26 +399,30 @@ psdList <- function(st) { tr_merged <- st_merged@traces[[1]] # Choose chunk size based on the chanel 'band code'(=sampling rate) - # See: http://www.iris.edu/manuals/SEED_appA.htm + # See: http://www.fdsn.org/seed_manual/SEEDManual_V2.4.pdf , Appendix A # NOTE: This choice was recommended by Mary Templeton + # + # get the channel name channel <- st_merged@traces[[1]]@stats@channel - if (stringr::str_detect(channel,"^L")) { + # the high frequency boundary is set to nyquist + hiFreq <- 0.5 * tr_merged@stats@sampling_rate + # Set an alignment frequency from which octaves will be generated + alignFreq <- 0.1 + if (stringr::str_detect(channel,"^V")) { ## added by REC - Apr 2016 + Z <- 24 * 3600 + loFreq <- 0.0001 + alignFreq <- 0.025 # special for VH + } else if (stringr::str_detect(channel,"^L")) { Z <- 3 * 3600 loFreq <- 0.001 - hiFreq <- 0.5 * tr_merged@stats@sampling_rate } else if (stringr::str_detect(channel,"^M")) { Z <- 2 * 3600 loFreq <- 0.0025 - hiFreq <- 0.5 * tr_merged@stats@sampling_rate } else { Z <- 3600 loFreq <- 0.005 - hiFreq <- 0.5 * tr_merged@stats@sampling_rate } - # Set an alignment frequency from which octaves will be generated - alignFreq <- 0.1 - # Initialize psdList <- list() index <- 1 @@ -432,14 +436,14 @@ psdList <- function(st) { while (secsLeft >= 0.99*Z) { # Slice a chunk out of the trace - tr <- slice(tr_merged,start,end) + tr <- IRISSeismic::slice(tr_merged,start,end) if (!isDC(tr)) { # NOTE: Each psd has elements: freq, spec, snclq, starttime, endtime psds <- McNamaraPSD(tr, loFreq, hiFreq, alignFreq) if (! "-Inf" %in% psds$spec) { psdCount <- psdCount + 1 - psdList[[psdCount]] <- McNamaraPSD(tr, loFreq, hiFreq, alignFreq) + psdList[[psdCount]] <- psds } } @@ -461,8 +465,10 @@ psdList <- function(st) { # Create instrument corrected noiseMatrix from psdList # NOTE: The incoming psdList contains raw, uncorrected PSDs. +# NOTE: Optional input evalresp is a dataframe of freq, amp, phase columns +# matching output of getEvalresp function # NOTE: The returned matrix contains corrected PSDs, one per row. -psdList2NoiseMatrix <- function(psdList) { +psdList2NoiseMatrix <- function(psdList, evalresp=NULL) { # Need to ensure that the IrisClient object exists as R has a default # "iris" dataset. See help(iris, package="datasets") @@ -493,11 +499,15 @@ psdList2NoiseMatrix <- function(psdList) { nfreq <- length(psdList[[1]]$freq) units <- "acc" - # Get instrument response - evalresp <- getEvalresp(iris,network,station,location,channel,time=starttime, + # Get instrument response from IRIS webservices if not provided by input evalresp + if (is.null(evalresp)) { + evalresp <- getEvalresp(iris,network,station,location,channel,time=starttime+1, minfreq=minfreq,maxfreq=maxfreq,nfreq=nfreq,units=units) + } - # NOTE: No try block needed as getEvalresp will generate appropriate errors + if (!("amp" %in% colnames(evalresp) & "freq" %in% colnames(evalresp))) { + stop(paste("error evalresp dataframe does not have columns named 'amp' and 'freq'")) + } # NOTE: Because we're operating in dB space we need to think in terms of logarithms. @@ -523,7 +533,7 @@ psdList2NoiseMatrix <- function(psdList) { } - # Create a correciton matrix that matches the dimensions of noiseMatrix + # Create a correction matrix that matches the dimensions of noiseMatrix correctionMatrix <- matrix(rep(correction,times=psdCount), nrow=psdCount, byrow=TRUE) # Dividing by the correction, in dB space, is just subtracting @@ -550,7 +560,7 @@ psdList2NoiseMatrix <- function(psdList) { # ... ... through all frequencies # ... repeated for each PSD, typically associated with an hour long chunk of signal -psdDF2NoiseMatrix <- function(DF) { +psdDF2NoiseMatrix <- function(DF, evalresp=NULL) { # Need to ensure that the IrisClient object exists as R has a default # "iris" dataset. See help(iris, package="datasets") @@ -586,10 +596,15 @@ psdDF2NoiseMatrix <- function(DF) { nfreq <- length(unique(DF$frequency)) units <- "acc" - evalresp <- getEvalresp(iris,network,station,location,channel,time=starttime, + if (is.null(evalresp)) { + evalresp <- getEvalresp(iris,network,station,location,channel,time=starttime+1, minfreq=minfreq,maxfreq=maxfreq,nfreq=nfreq,units=units) - - # NOTE: No try block needed as getEvalresp will generate appropriate errors + } + + if (!("amp" %in% colnames(evalresp) & "freq" %in% colnames(evalresp))) { + stop(paste("error evalresp dataframe does not have columns named 'amp' and 'freq'")) + } + # NOTE: Because we're operating in dB space we need to think in terms of logarithms. @@ -621,7 +636,7 @@ psdDF2NoiseMatrix <- function(DF) { # INPUT: matrix where columns are frequencies and rows are individual, corrected PSDs # OUTPUT: matrix where columns are frequencies and rows are counts of how many input PSDs have that power level -noiseMatrix2PdfMatrix <- function(noiseMatrix, lo=-200, hi=-50, binSize=1) { +noiseMatrix2PdfMatrix <- function(noiseMatrix, lo=-310, hi=55, binSize=1) { # NOTE: Define a function to convert one column of noiseMatrix into a # NOTE: column of histogram counts. However many rows exist in noiseMatrix, @@ -629,8 +644,8 @@ noiseMatrix2PdfMatrix <- function(noiseMatrix, lo=-200, hi=-50, binSize=1) { histCounts <- function(x, lo, hi, binSize) { breaks <- seq(lo,hi,binSize) - nbins <- length(breaks) - 1 - discretizedValue <- .bincode(x, breaks) + nbins <- length(breaks) - 1 + discretizedValue <- .bincode(x, breaks,include.lowest=TRUE) return(tabulate(discretizedValue, nbins)) } @@ -654,7 +669,7 @@ noiseModels <- function(freq) { period <- 1/freq # NOTE: Original New High/Low Noise Models in Peterson paper: - # NOTE: http://www.mttmllr.com/ADS/DATA/peterson_usgs_seismic_noise_ofr93-322.pdf + # NOTE: https://pubs.usgs.gov/of/1993/0322/report.pdf # NOTE: # NOTE: IRIS DMC 2005 paper # NOTE: http://www.earth.northwestern.edu/people/seth/327/HV/McNamaraetal_AmbientNoise.3.0.pdf @@ -664,6 +679,7 @@ noiseModels <- function(freq) { # NOTE: # NOTE: Source code to compare: # NOTE: https://github.com/g2e/seizmo/blob/master/noise/nlnm.m + # NOTE: https://github.com/g2e/seizmo/blob/master/noise/nhnm.m # New Low Noise Model in acceleration: NLNMacc = A + B*log10(T) referred to 1 (m/s^2)^2/Hz # T ( minimum period) @@ -717,7 +733,8 @@ minPeriod A B 7.90 -93.37 -22.42 15.40 73.54 -162.98 20.00 -151.52 10.01 - 100000 -206.66 31.63 + 354.80 -206.66 31.63 + 100000 -206.66 31.63 ") # We create breaks based on the first column of each table to figure out @@ -737,19 +754,19 @@ minPeriod A B # Calculate basic statistics on all PSDs in list or dataframe. -psdStatistics <- function(PSDs) { +psdStatistics <- function(PSDs, evalresp=NULL) { # Get list of frequencies and a noiseMatrix if (class(PSDs) == "list") { freq <- PSDs[[1]]$freq - noiseMatrix <- psdList2NoiseMatrix(PSDs) + noiseMatrix <- psdList2NoiseMatrix(PSDs, evalresp) } else if (class(PSDs) == "data.frame") { freq <- sort(unique(PSDs$freq)) - noiseMatrix <- psdDF2NoiseMatrix(PSDs) + noiseMatrix <- psdDF2NoiseMatrix(PSDs, evalresp) } @@ -775,15 +792,21 @@ psdStatistics <- function(PSDs) { # Calculate pctAboveNHNM and pctBelowNLNM ------------ noiseModels <- noiseModels(freq) - - aboveCount <- rep(0,ncol) - belowCount <- rep(0,ncol) - - for (j in seq(ncol)) { - aboveCount[j] <- length(which(noiseMatrix[,j] > noiseModels$nhnm[j])) + + # calculate percent metrics based on frequencies that have a noise model value and a noiseMatrix value + notNA <- which(!is.na(noiseMatrix[1,]) & !is.na(noiseModels[,1])) # high and low noise models have same upper and lower frequency limits, so only check one model + + aboveCount <- rep(NA,ncol) + aboveCount[notNA] <- 0 + + belowCount <- rep(NA,ncol) + belowCount[notNA] <- 0 + + for (j in notNA) { + aboveCount[j] <- length(which(noiseMatrix[,j] > noiseModels$nhnm[j])) belowCount[j] <- length(which(noiseMatrix[,j] < noiseModels$nlnm[j])) } - + pct_above <- 100 * aboveCount/nrow pct_below <- 100 * belowCount/nrow @@ -793,10 +816,12 @@ psdStatistics <- function(PSDs) { # For mode, we need to convert to the discretized pdfMatrix # that contains in each bin, the count of PSDs that have # that value. The mode is just the bin with the highest count. - - # Default dbBins = seq(-200,-50,1) - pdfBins <- seq(-200,-50,1) - pdfMatrix <- noiseMatrix2PdfMatrix(noiseMatrix) + + lo <- floor(min(noiseMatrix[,which(!is.na(noiseMatrix[1,]))])) + hi <- ceiling(max(noiseMatrix[,which(!is.na(noiseMatrix[1,]))])) + pdfBins <- seq(lo,hi,1) + + pdfMatrix <- noiseMatrix2PdfMatrix(noiseMatrix,lo-0.5,hi+0.5) colMode <- rep(NA,ncol) for (j in seq(ncol)) { @@ -807,7 +832,7 @@ psdStatistics <- function(PSDs) { return(list(noiseMatrix=noiseMatrix, pdfMatrix=pdfMatrix, freq=freq, - pdfBins=pdfBins, + pdfBins=pdfBins, # midpoint of bins max=colMax, min=colMin, mean=colMean, @@ -824,7 +849,8 @@ psdStatistics <- function(PSDs) { # Plot instrument corrected noise values of PSDs in list or dataframe -psdPlot <- function(PSDs, style='psd', showNoiseModel=TRUE, showMaxMin=TRUE, showMode=TRUE, showMean=FALSE, ...) { +psdPlot <- function(PSDs, style='psd', evalresp=NULL, ylo=-200,yhi=-50, showNoiseModel=TRUE, + showMaxMin=TRUE, showMode=TRUE, showMean=FALSE, showMedian=FALSE, ...) { if (class(PSDs) == "list") { @@ -866,11 +892,16 @@ psdPlot <- function(PSDs, style='psd', showNoiseModel=TRUE, showMaxMin=TRUE, sho } # Generate basic statics as well as noiseMatrix and pdfMatrix - - stats <- psdStatistics(PSDs) + + if (is.null(evalresp)) { + stats <- psdStatistics(PSDs) + } else { + stats <- psdStatistics(PSDs, evalresp=evalresp) + } noiseMatrix <- stats$noiseMatrix pdfMatrix <- stats$pdfMatrix + pdfBins <- stats$pdfBins # Choose appropriate limits for period period <- 1/freq @@ -886,9 +917,7 @@ psdPlot <- function(PSDs, style='psd', showNoiseModel=TRUE, showMaxMin=TRUE, sho seq(10,100,length.out=10)) } - # Choose appropriate limits for dB - ylo <- -200 - yhi <- -50 + # Choose appropriate limits for dB, now as input variables ylim <- c(ylo,yhi) horizontalLines <- seq(ylo,yhi,10) @@ -910,7 +939,8 @@ psdPlot <- function(PSDs, style='psd', showNoiseModel=TRUE, showMaxMin=TRUE, sho } else if (style == 'pdf') { # Set up colors and breaks - cols <- c('grey90', rev(grDevices::rainbow(30))[4:30]) + # cols <- c('grey90', rev(grDevices::rainbow(30,alpha=seq(1,0.1,-1/30))[4:30]) + cols <- c('white', rev(grDevices::rainbow(30))[4:30]) breaks <- c(0,seq(0.001,max(pdfMatrix),length.out=length(cols))) # NOTE: To get image() to plot with the same axes as the data displayed as a table you @@ -920,11 +950,13 @@ psdPlot <- function(PSDs, style='psd', showNoiseModel=TRUE, showMaxMin=TRUE, sho # Initial plot graphics::image(t(pdfMatrix)[ncol(pdfMatrix):1,], x=rev(period), - y=seq(ylo,yhi), + y=pdfBins, + xlim=xlim, + ylim=c(ylo,yhi), breaks=breaks, col=cols, las=1, log="x", - xlab="Period (Sec)", ylab="Power dB", + xlab="Period (Sec)", ylab="Power (dB)", main=paste("PDF plot of",psdCount,"corrected, hourly PSDs for",snclq), ...) @@ -954,29 +986,38 @@ psdPlot <- function(PSDs, style='psd', showNoiseModel=TRUE, showMaxMin=TRUE, sho col <- c(col,'blue','red') } if (showMode) { - graphics::points(period, stats$mode, type='l', col='yellow', lwd=3) + graphics::points(period, stats$mode, type='l', col='black', lwd=2) legend <- c(legend,"mode") - lwd <- c(lwd,3) - col <- c(col,'yellow') + lwd <- c(lwd,2) + col <- c(col,'black') } if (showMean) { - graphics::points(period, stats$mean, type='l', col='orange', lwd=4) + graphics::points(period, stats$mean, type='l', col='orange', lwd=2) legend <- c(legend,"mean") - lwd <- c(lwd,4) + lwd <- c(lwd,2) col <- c(col,'orange') } + if (showMedian) { + graphics::points(period, stats$median, type='l', col='yellow4', lwd=2) + legend <- c(legend,"median") + lwd <- c(lwd,2) + col <- c(col,'yellow4') + } + # Add grid lines graphics::abline(h=horizontalLines, v=verticalLines, col='gray50', lty='dotted') # Add a legend - legend("topleft", bg='white', + # legend("topleft", bg='white', + legend("topleft", bty='n', legend=legend, lwd=lwd, col=col) # Add the starttime and endtime - legend("topright", bg='white', + # legend("topright", bg='white', + legend("topright", bty='n', inset=c(0.09,0), legend=c(paste(starttime," start "), paste(endtime," end "))) @@ -1011,6 +1052,23 @@ hilbertFFT <- function(x) { return(hfft) } +#' @author REC +#' if vec represents a set of binned counts of incrementing values (ascending) +#' return a vector of associated bin values with the proper count of each value. +#' @param vec a histogram vector or ordered set of binned counts +#' @param startVal the initial value of the first bin element +#' @param incr the increment rate of each subsequent bin value +#' @return a vector of bin values with appropriate counts of each +unHistogram <- function(vec, startVal=1, incr=1) { + j <- 0 + k <- list() + for (i in vec) { + if (i > 0) k <- c(k,rep(startVal+(j*incr),i)) + j <- j+1 + } + return(unlist(k)) +} + ################################################################################ # END ################################################################################ diff --git a/build/vignette.rds b/build/vignette.rds index f521682..4ba740d 100644 Binary files a/build/vignette.rds and b/build/vignette.rds differ diff --git a/inst/doc/IRISSeismic-intro.R b/inst/doc/IRISSeismic-intro.R index 34d1cf5..697dc2f 100644 --- a/inst/doc/IRISSeismic-intro.R +++ b/inst/doc/IRISSeismic-intro.R @@ -65,7 +65,7 @@ slotNames(st@traces[[1]]) ################################################### -### code chunk number 10: IRISSeismic-intro.Rnw:246-253 +### code chunk number 10: IRISSeismic-intro.Rnw:247-254 ################################################### as.POSIXct("2010-02-27", tz="GMT") # good as.POSIXct("2010-02-27 04:00:00", tz="GMT") # good @@ -77,13 +77,13 @@ as.POSIXct("2010-02-27T04:00:00", tz="GMT") # BAD -- no formatting ################################################### -### code chunk number 11: IRISSeismic-intro.Rnw:268-269 +### code chunk number 11: IRISSeismic-intro.Rnw:269-270 ################################################### help("IRISSeismic",package="IRISSeismic") ################################################### -### code chunk number 12: IRISSeismic-intro.Rnw:297-311 +### code chunk number 12: IRISSeismic-intro.Rnw:298-312 ################################################### starttime <- as.POSIXct("2010-02-27", tz="GMT") endtime <- as.POSIXct("2010-02-28", tz="GMT") @@ -102,7 +102,7 @@ layout(1) # restore original layout ################################################### -### code chunk number 13: IRISSeismic-intro.Rnw:343-350 +### code chunk number 13: IRISSeismic-intro.Rnw:344-351 ################################################### starttime <- as.POSIXct("2002-04-20", tz="GMT") endtime <- as.POSIXct("2002-04-21", tz="GMT") @@ -114,7 +114,7 @@ to <- triggerOnset(tr,picker,threshold) ################################################### -### code chunk number 14: IRISSeismic-intro.Rnw:372-382 +### code chunk number 14: IRISSeismic-intro.Rnw:373-383 ################################################### layout(matrix(seq(3))) # layout a 3x1 matrix closeup1 <- eventWindow(tr,picker,threshold,3600) @@ -129,7 +129,7 @@ layout(1) # restore original layout ################################################### -### code chunk number 15: IRISSeismic-intro.Rnw:400-404 +### code chunk number 15: IRISSeismic-intro.Rnw:401-405 ################################################### starttime <- as.POSIXct("2010-02-27", tz="GMT") endtime <- as.POSIXct("2010-02-28", tz="GMT") @@ -138,7 +138,7 @@ availability ################################################### -### code chunk number 16: IRISSeismic-intro.Rnw:453-501 +### code chunk number 16: IRISSeismic-intro.Rnw:456-504 ################################################### # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -150,7 +150,7 @@ endtime <- starttime + 3600 * 24 * 2 # Find biggest seismic event over these two days -- it's the "Nisqually" events <- getEvent(iris, starttime, endtime, minmag=5.0) bigOneIndex <- which(events$magnitude == max(events$magnitude)) -bigOne <- events[bigOneIndex,] +bigOne <- events[bigOneIndex[1],] # Find US stations that are available within 10 degrees of arc of the # event location during the hour after the event diff --git a/inst/doc/IRISSeismic-intro.Rnw b/inst/doc/IRISSeismic-intro.Rnw index 24b3df7..c1d6715 100644 --- a/inst/doc/IRISSeismic-intro.Rnw +++ b/inst/doc/IRISSeismic-intro.Rnw @@ -216,7 +216,8 @@ The \code{Trace} object has the following slots: \item \code{stats} -- \code{TraceHeader} object (metadata from the trace) \item \code{Sensor} -- instrument equipment name \item \code{InstrumentSensitivity} -- instrument total sensitivity (stage 0 gain) - \item \code{InputUnits} -- instrument data qcquisition input units + \item \code{SensitivityFrequency} -- the frequency (Hz) at which the \code{InstrumentSensitivity} is correct + \item \code{InputUnits} -- instrument data acquisition input units \item \code{data} -- vector of \code{numeric} data (the actual signal) \end{itemize} @@ -414,7 +415,7 @@ The \code{getAvailability()} method accepts the following arguments: \item \code{channel} -- channel code \item \code{starttime} -- POSIXct starttime (GMT) \item \code{endtime} -- POSIXct endtime (GMT) - \item \code{restricted} -- optional flag whether to report on restricted data (default=\code{FALSE}) + \item \code{includerestricted} -- optional flag whether to report on restricted data (default=\code{FALSE}) \item \code{latitude} -- optional latitude when specifying location and radius \item \code{longitude} -- optional longitude when specifying location and radius \item \code{minradius} -- optional minimum radius when specifying location and radius @@ -432,9 +433,11 @@ from web services of the same name. The suite of methods returning dataframes i \begin{itemize} \item \code{getAvailability} \item \code{getChannel} + \item \code{getDataselect} \item \code{getEvalresp} \item \code{getEvent} \item \code{getNetwork} + \item \code{getSNCL} \item \code{getStation} \item \code{getTraveltime} \item \code{getUnavailability} @@ -461,7 +464,7 @@ endtime <- starttime + 3600 * 24 * 2 # Find biggest seismic event over these two days -- it's the "Nisqually" events <- getEvent(iris, starttime, endtime, minmag=5.0) bigOneIndex <- which(events$magnitude == max(events$magnitude)) -bigOne <- events[bigOneIndex,] +bigOne <- events[bigOneIndex[1],] # Find US stations that are available within 10 degrees of arc of the # event location during the hour after the event diff --git a/inst/doc/IRISSeismic-intro.pdf b/inst/doc/IRISSeismic-intro.pdf index 4959eb1..618a592 100644 Binary files a/inst/doc/IRISSeismic-intro.pdf and b/inst/doc/IRISSeismic-intro.pdf differ diff --git a/man/DDT.Rd b/man/DDT.Rd index ac52f79..6177af4 100644 --- a/man/DDT.Rd +++ b/man/DDT.Rd @@ -35,7 +35,7 @@ A new \code{Trace} object is returned. %%\references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -43,6 +43,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} %% \seealso{ } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -63,6 +64,7 @@ mtext("Demean - Detrend - Cosine Taper",line=-2,side=3,adj=0.05,col='red') # Restore default layout layout(1) + } } \keyword{methods} diff --git a/man/IRISSeismic-package.Rd b/man/IRISSeismic-package.Rd index f7efaa5..0c19f7d 100644 --- a/man/IRISSeismic-package.Rd +++ b/man/IRISSeismic-package.Rd @@ -17,7 +17,7 @@ are inspired by the functionality available in the python ObsPy package (\url{ht } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \references{ @@ -27,6 +27,130 @@ IRIS DMC web services: \url{http://service.iris.edu/} } \section{History}{ + version 1.3.8 + \itemize{ + \item{\code{getDataselect} does not add a quality indicator to url by default. IRIS webservices itself defaults to quality="M"} + \item{\code{getStation} and \code{getChannel} do not add includerestricted indicator to url by default. IRIS webservices itself defaults to \code{TRUE}} + \item{better handling of textConnections} + } + version 1.3.7 + \itemize{ + \item{users can now supply instrument response information in the form of frequency, amplitude, phase + to the functions \code{psdStatistics}, \code{psdList2NoiseMatrix}, \code{psdPlot}, + in place of the getEvalresp webservice call. Function argument order for \code{psdPlot} is changed.} + \item{added \code{showMedian} option to \code{psdPlot}} + } + + version 1.3.5 + \itemize{ + \item{added ignoreEpoch option to getDataselect} + } + + version 1.3.4 -- webservices and plotting + \itemize{ + \item{\code{getEvent} forwards http://service.iris.edu/fdsnws/event/1/ calls to http://earthquake.usgs.gov/fdsnws/event/1/} + \item{\code{getDistaz} changes output dataframe column name ellipsoid..attrs to ellipsoid.name} + \item{\code{plot.Trace} allows for user supplied ylab and xlab input} + } + + version 1.3.3 -- documentation + \itemize{ + \item{Updated documentation and corrected outdated links} + } + + version 1.3.2 -- bug fix + \itemize{ + \item{noiseModels(), minor correction to the New High Noise Model} + } + + version 1.3.1 -- bug fixes + \itemize{ + \item{psdStatistics() correctly handles NA values when calculating high and low PDF bin limits + and returns pct_above and pct_below vectors of correct length} + } + + version 1.3.0 -- compatibility with IRIS webservices + \itemize{ + \item{getDistaz() returns new variables from output of http://services.iris.edu/irisws/distaz/1/} + } + + version 1.2.2 -- PDF bug fix + \itemize{ + \item{psdList2NoiseMatrix() adds 1 second to start time in getEvalresp call to work around a quirk in http://services.iris.edu/irisws/evalresp/1/ + webservice that will not return a response if the start time is exactly on a metadata epoch boundary.} + } + + version 1.2.1 -- PDF + \itemize{ + \item{psdPlot() now compatible with changes to psdStatistics() in previous version. Adds ylo, yhi arguments to customize y-axis limits in plot.} + } + + version 1.2.0 -- PDF + \itemize{ + \item{psdStatistics() changes method of setting PDF bins from fixed values to bins based on the high and low PSD values + and shifts bin centers by 0.5 dB. The result now matches output from http://services.iris.edu/mustang/noise-pdf.} + } + + verison 1.1.7 -- improved error handling + \itemize{ + \item{getDataselect(), getNetwork(), getStation(), getChannel(), getAvailability(), getEvalresp(), getTraveltime() error handling + now report unexpected http status codes.} + } + + version 1.1.6 -- bug fixes + \itemize{ + \item{getGaps() fixes issues with multiple sample rates and setting minimum gap length.} + \item{mergeTraces.Stream() relaxes criteria for acceptable sample rate jitter.} + } + + version 1.1.5 -- trace rotation + \itemize{ + \item{rotate2D() changes orthogonality test tolerance from 5 degrees to 3 degrees.} + } + + version 1.1.4 -- trace rotation + \itemize{ + \item{rotate2D() exits if traces are not orthogonal.} + } + + version 1.1.3 -- bug fix + \itemize{ + \item{psdStatistics() fixes bug in calculation of pct_above and pct_below.} + } + + version 1.1.1 -- bug fixes + \itemize{ + \item{getGaps() minor bug fix.} + \item{mergeTraces.Stream() minor bug fix.} + } + + version 1.0.10 -- new data request argument and bug fixes + \itemize{ + \item{Imports seismicRoll (>= 1.1.0).} + \item{getGaps() fixes bugs in calculation of initial and final gap of Trace.} + \item{getDataselect(), getSNCL() adds "inclusiveEnd" argument, a logical that determines whether a data point that falls exactly on the requested endtime + is included in the Trace.} + \item{libmseed change, when multiple sample rates exist in miniseed records use the mode of all sample rates instead of using the sample rate in the first record.} + \item{psdList() added rule for octave generation for channel codes that start with "V".} + } + + version 1.0.9 -- Trace class expansion and bug fixes + \itemize{ + \item{Improved error handling for getAvailability(), getChannel(), getDataselect(), getEvalresp(), miniseed2Stream().} + \item{parseMiniSEED.c, unpackdata.c updated. Fixes protection stack overflow issue.} + \item{getGaps() includes a 0.5/sampling_rate tolerance factor.} + \item{miniseed2Stream() uses endtime from parseMiniSEED instead of calculating from the sample rate.} + \item{Trace class now contains slots for optional metadata "latitude", "longitude", "elevation", "depth", "azimuth", "dip", "SensitivityFrequency".} + \item{rotate2D() uses Trace class "azimuth" slot information to identify channel orientation before rotation instead of assuming lead and lag channel + from trace input order.} + } + + version 1.0.8 -- fixes required by ISPAQ + \itemize{ + \item{Removed 'maps' and 'mapdata' from Suggested: packages.} + \item{Changed URL syntax for FDSN web services to use "format=..." instead of "output=...".} + \item{Fixed bug in getSNCL() so that it works when the "quality" argument is missing.} + } version 1.0.6 -- CRAN updates required \itemize{ @@ -172,6 +296,7 @@ IRIS DMC web services: \url{http://service.iris.edu/} } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient", debug=TRUE) @@ -185,6 +310,7 @@ st <- getDataselect(iris,"IU","ANMO","00","BHZ",starttime,endtime) tr1 <- st@traces[[1]] show(tr1@stats) plot(tr1) + } } %% Optional keywords diff --git a/man/IrisClient-class.Rd b/man/IrisClient-class.Rd index 87888b7..4d56d10 100644 --- a/man/IrisClient-class.Rd +++ b/man/IrisClient-class.Rd @@ -56,12 +56,12 @@ A class for making data and metadata requests from IRIS DMC web services. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ -The \code{IrisClient} object is inspired by the \code{iris.client.Client} class found in the -python ObsPy package (\url{http://docs.obspy.org/packages/autogen/obspy.iris.client.Client.html}). +The \code{IrisClient} object is inspired by the \code{clients.fdsn.client.Client} class found in the +python ObsPy package (\url{http://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.html#obspy.clients.fdsn.client.Client}). } %% ~Make other sections like Warning with \section{Warning }{....} ~ @@ -77,7 +77,7 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) mean(st) } } diff --git a/man/McNamaraBins.Rd b/man/McNamaraBins.Rd index 714f004..b0b8af3 100644 --- a/man/McNamaraBins.Rd +++ b/man/McNamaraBins.Rd @@ -41,7 +41,7 @@ A dataframe containing binned values is returned with the same column names as t } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %%\note{} diff --git a/man/McNamaraPSD.Rd b/man/McNamaraPSD.Rd index d07d94c..f109025 100644 --- a/man/McNamaraPSD.Rd +++ b/man/McNamaraPSD.Rd @@ -84,7 +84,7 @@ and \code{endtime} are single values. } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ diff --git a/man/STALTA.Rd b/man/STALTA.Rd index 7713719..6ed4a39 100644 --- a/man/STALTA.Rd +++ b/man/STALTA.Rd @@ -109,7 +109,7 @@ Bulletin of the Seismological Society of America, Vol. 84, No. 2, pp. 366-376, A } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ @@ -130,8 +130,6 @@ For higher resolution channels, picking an increment of \code{2/sampling_rate} c # Open a connection to IRIS DMC webservices iris <- new("IrisClient") -# NOTE: The following trace has 1.728 million points. -# NOTE: Downloading and calculating STALTA may take a while. starttime <- as.POSIXct("2010-02-27",tz="GMT") endtime <- as.POSIXct("2010-02-28",tz="GMT") diff --git a/man/Stream-class.Rd b/man/Stream-class.Rd index 015e927..db314a7 100644 --- a/man/Stream-class.Rd +++ b/man/Stream-class.Rd @@ -42,7 +42,7 @@ Objects are typically created by calls to \link{getDataselect}. \item{\code{dq_flags}:}{Object of class \code{"integer"}: Accumulators for the dq_flags bits in each miniSEED record. } \item{\code{timing_qual}:}{Object of class \code{"numeric"}: - Average timing quality assicuated with miniSEED records.} + Average timing quality associated with miniSEED records.} \item{\code{traces}:}{Object of class \code{"list"}: List of \code{Trace} objects. } } @@ -104,7 +104,7 @@ Objects are typically created by calls to \link{getDataselect}. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ The \code{Stream} object is inspired by the \code{Stream} class found in the @@ -154,6 +154,7 @@ The following code documentation describes how each of the flags is used within %% \seealso{ } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -161,7 +162,7 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) min(st) median(st) mean(st) @@ -169,6 +170,7 @@ max(st) sd(st) rms(st) rmsVariance(st) + } } \keyword{classes} diff --git a/man/Trace-class.Rd b/man/Trace-class.Rd index 0bae36b..f5543aa 100644 --- a/man/Trace-class.Rd +++ b/man/Trace-class.Rd @@ -31,13 +31,15 @@ populated by calls to \link{getDataselect}. \describe{ \item{\code{id}:}{Object of class \code{"character"}: Unique "SNCL" identifier specifying the Network, Station, Location, Channel and Quality factor - associated with this trace: eg. \code{AK.PIN..VEA.B}. The \code{id} is generated automatically + associated with this trace: eg. \code{AK.PIN..VEA.M}. The \code{id} is generated automatically when the trace is first created and is intended for read only.} \item{\code{Sensor}:}{Object of class \code{"character"}: - Equipment name.} + Instrument name.} \item{\code{InstrumentSensitivity}:}{Object of class \code{"numeric"}: The total sensitivity for a channel, representing the complete acquisition system expressed as a scalar. Equivalent to SEED stage 0 gain.} + \item{\code{SensitivityFrequency}:}{Object of class \code{"numeric"}: + The frequency at which the total sensitivity is correct.} \item{\code{InputUnits}:}{Object of class \code{"character"}: The units of the data as input from the perspective of data acquisition. After correcting data for this response, these would be the resulting units.} @@ -52,7 +54,7 @@ populated by calls to \link{getDataselect}. \item{as.vector}{\code{signature(x="Trace")}: returns the data slot; equivalent to \code{x@data} } \item{DDT}{\code{signature(x="Trace", demean="logical", detrend="logical", taper="numeric")}: - returns a new trace that has been 'cleaned up' for further processing (see \link{DDT}) } + returns a new trace that has been 'cleaned up' for further processing by applying demean, detrend, and taper techniques (see \link{DDT}) } \item{envelope}{\code{signature(x="Trace")}: returns the envelope of the seismic signal (see \link{envelope}) } \item{isDC}{\code{signature(x="Trace")}: @@ -78,7 +80,7 @@ populated by calls to \link{getDataselect}. \item{sd}{\code{signature(x="Trace", na.rm="logical")}: returns the standard deviation of the data; equivalent to \code{sd(x@data)} } \item{slice}{\code{signature(x="Trace", starttime="POSIXct", endtime="POSIXct")}: - returns a new \code{Trace} sliced out of an existing \code{Trace} (see \link{slice})} + returns a new \code{Trace} subset of an existing \code{Trace} (see \link{slice})} \item{STALTA}{\code{signature(x="Trace",staSecs="numeric",ltaSecs="numeric",algorithm="character", ...)}: returns the STALTA picker result (see \link{STALTA})} \item{triggerOnset}{\code{signature(x="Trace", picker="numeric", threshold="numeric", ...)}: @@ -89,7 +91,7 @@ populated by calls to \link{getDataselect}. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ @@ -102,6 +104,7 @@ python ObsPy package (\url{http://docs.obspy.org/packages/autogen/obspy.core.tra %% \seealso{ } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -110,7 +113,7 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) # Get the first trace and generate some statistics tr1 <- st@traces[[1]] @@ -121,6 +124,7 @@ max(tr1) sd(tr1) rms(tr1) rmsVariance(tr1) + } } \keyword{classes} diff --git a/man/TraceHeader-class.Rd b/man/TraceHeader-class.Rd index 62cfd7d..d807d54 100644 --- a/man/TraceHeader-class.Rd +++ b/man/TraceHeader-class.Rd @@ -17,7 +17,7 @@ \title{ Class \code{"TraceHeader"} } \description{ -A container for metadata associated with a \code{Trace} object. +A container for metadata associated with a \code{Trace} object. Originally populated by information in the miniseed trace header; it now has the option of including additional station and channel metadata. } \section{Objects from the Class}{ Objects can be created by calls of the form \code{new("TraceHeader", headerList, headerLine, ...)}. @@ -28,16 +28,22 @@ typically populated by a webservice request. (see \link{IrisClient-class}) \section{Slots}{ \describe{ \item{\code{sampling_rate}:}{Object of class \code{"numeric"}: Sampling rate in hertz.} - \item{\code{delta}:}{Object of class \code{"numeric"}: Sample distance in seconds.} + \item{\code{delta}:}{Object of class \code{"numeric"}: Sample interval in seconds.} \item{\code{calib}:}{Object of class \code{"numeric"}: Calibration factor.} \item{\code{npts}:}{Object of class \code{"integer"}: Number of sample points.} \item{\code{network}:}{Object of class \code{"character"}: Seismic network name.} \item{\code{location}:}{Object of class \code{"character"}: Location code.} \item{\code{station}:}{Object of class \code{"character"}: Station name.} \item{\code{channel}:}{Object of class \code{"character"}: Channel code.} - \item{\code{quality}:}{Object of class \code{"character"} Data quality code} + \item{\code{quality}:}{Object of class \code{"character"}: Data quality code.} \item{\code{starttime}:}{Object of class \code{"POSIXct"}: Start time.} \item{\code{endtime}:}{Object of class \code{"POSIXct"}: End time.} + \item{\code{latitude}:}{Object of class \code{"numeric"}: Latitude.} + \item{\code{longitude}:}{Object of class \code{"numeric"}: Longitude.} + \item{\code{elevation}:}{Object of class \code{"numeric"}: Elevation.} + \item{\code{depth}:}{Object of class \code{"numeric"}: Depth.} + \item{\code{azimuth}:}{Object of class \code{"numeric"}: Azimuth.} + \item{\code{dip}:}{Object of class \code{"numeric"}: Dip.} \item{\code{processing}:}{Object of class \code{"list"}: Information strings describing processing applied to this trace. } } } @@ -45,7 +51,7 @@ typically populated by a webservice request. (see \link{IrisClient-class}) \section{Methods}{ \describe{ \item{as.headerLine}{\code{signature(obj = "TraceHeader")}: - Prints out the information in the \code{TraceHeader} as an ascii header line, e.g. + Prints out the information in the \code{TraceHeader} as an ascii header line, not including any station and channel metadata not found in the miniseed trace header, e.g., \preformatted{ TIMESERIES LD_POTS__HHZ_M, 351 samples, 100.503 sps, \ 2012-01-29T00:00:00.006000, SLIST, INTEGER, COUNTS}} @@ -57,7 +63,7 @@ typically populated by a webservice request. (see \link{IrisClient-class}) %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ @@ -75,6 +81,7 @@ confusion in the context of R. Instead, the name %% \seealso{ } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -82,11 +89,12 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) # Get the first trace and show the associated metadata tr1 <- st@traces[[1]] show(tr1@stats) + } } \keyword{classes} diff --git a/man/basicStats.Rd b/man/basicStats.Rd index 35eb538..540be16 100644 --- a/man/basicStats.Rd +++ b/man/basicStats.Rd @@ -110,7 +110,7 @@ a numeric vector is returned of the same length as \code{Stream@traces}. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ @@ -126,6 +126,7 @@ between \code{Traces} are not counted. %% \seealso{ } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -133,7 +134,7 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) # Get the first trace and generate some statistics tr1 <- st@traces[[1]] @@ -141,6 +142,7 @@ length(tr1) max(tr1) mean(tr1) sd(tr1) + } } % Add one or more standard keywords, see file 'KEYWORDS' in the diff --git a/man/butterworth.Rd b/man/butterworth.Rd index 406b1fc..4926145 100644 --- a/man/butterworth.Rd +++ b/man/butterworth.Rd @@ -42,7 +42,7 @@ A new \code{Trace} object is returned. %%\references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -56,12 +56,14 @@ signal::butter, signal::filter # Open a connection to IRIS DMC webservices iris <- new("IrisClient") -# Trying to match the results in figure 2a of +# Compare to the results in figure 2a of # # "Determination of New Zealand Ocean Bottom Seismometer Orientation # via Rayleigh-Wave Polarization", Stachnik et al. # -# geode.colorado.edu/~sheehan/pdf/StachnikNZorientation-SRL2012.pdf +# http://srl.geoscienceworld.org/content/83/4/704 +# +# (note: since publication, ZU.NZ19..BH1 has been renamed BH2 and ZU.NZ19..BH2 has been renamed BH1) starttime <- as.POSIXct("2009-02-18 22:01:07",tz="GMT") endtime <- starttime + 630 diff --git a/man/crossSpectrum.Rd b/man/crossSpectrum.Rd index 7b8cd3b..c6b80ba 100644 --- a/man/crossSpectrum.Rd +++ b/man/crossSpectrum.Rd @@ -5,7 +5,7 @@ \description{ The crossSpectrum() function is based on \R's spec.pgram() function and attempts to provide -complete results of cross-spectal FFT analysis in a programmer-friendly fashion. +complete results of cross-spectral FFT analysis in a programmer-friendly fashion. } \usage{ @@ -21,7 +21,7 @@ crossSpectrum(x, spans = NULL, kernel = NULL, taper = 0.1, \item{kernel}{alternatively, a kernel smoother of class "tskernel"} \item{taper}{specifies the proportion of data to taper. A split cosine bell taper is applied to this proportion of the data at the beginning and end of the series} \item{pad}{proportion of data to pad. Zeros are added to the end of the series to increase its length by the proportion pad} - \item{fast}{logical; if TRUE, pad the series to a highly composite length} + \item{fast}{logical. if TRUE, pad the series to a highly composite length} \item{demean}{logical. If TRUE, subtract the mean of the series} \item{detrend}{logical. If TRUE, remove a linear trend from the series. This will also remove the mean} \item{na.action}{NA action function} @@ -61,11 +61,11 @@ A dataframe with the following columns: \href{http://octave-signal.sourcearchive.com/documentation/1.0.7/pwelch_8m-source.html}{Octave pwelch() source code} -\href{http://cran.r-project.org/web/packages/psd/vignettes/normalization.pdf}{Normalization of Power Spectral Density estimates} +\href{https://CRAN.R-project.org/package=psd/vignettes/normalization.pdf}{Normalization of Power Spectral Density estimates} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %%\note{} @@ -77,7 +77,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} \examples{ \dontrun{ # Create a new IrisClient -iris <- new("IrisClient", debug=TRUE) +iris <- new("IrisClient") # Get seismic data starttime <- as.POSIXct("2011-05-01", tz="GMT") @@ -88,7 +88,7 @@ st2 <- getDataselect(iris,"CI","PASC","10","BHZ",starttime,endtime) tr1 <- st1@traces[[1]] tr2 <- st2@traces[[1]] -# Both traces have a smpling rate of 40 Hz +# Both traces have a sampling rate of 40 Hz sampling_rate <- tr1@stats@sampling_rate ts1 <- ts(tr1@data,frequency=sampling_rate) diff --git a/man/envelope.Rd b/man/envelope.Rd index 25917bc..fc99bad 100644 --- a/man/envelope.Rd +++ b/man/envelope.Rd @@ -41,7 +41,7 @@ whose data have been replaced with the envelope of the seismic signal. } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ diff --git a/man/eventWindow.Rd b/man/eventWindow.Rd index 77a389f..ec2353b 100644 --- a/man/eventWindow.Rd +++ b/man/eventWindow.Rd @@ -50,7 +50,7 @@ A new \code{Trace} object is returned. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/getAvailability.Rd b/man/getAvailability.Rd index 4bcf363..26d1d87 100644 --- a/man/getAvailability.Rd +++ b/man/getAvailability.Rd @@ -77,14 +77,14 @@ The IRIS DMC station web service: \url{http://service.iris.edu/fdsnws/station/1/} -This implementation was inspired by the functionality in the obspy station() method. +This implementation was inspired by the functionality in the obspy get_stations() method. - \url{http://docs.obspy.org/packages/autogen/obspy.iris.client.Client.station.html} + \url{http://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_stations.html} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -97,6 +97,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} %% max preformatted length --------------------------------------------------- \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -108,8 +109,9 @@ endtime <- starttime + 2*24*3600 events <- getEvent(iris,starttime,endtime,6.0) events -# Event #3 is the Nisqually quake in Washington State -e <- events[3,] +# biggest event is Nisqually +eIndex <- which(events$magnitude == max(events$magnitude)) +e <- events[eIndex[1],] # Find all BHZ channels collecting data at the time of the quake and within # 5 degrees of the quake epicenter @@ -117,6 +119,7 @@ channels <- getAvailability(iris,"*","*","*","BHZ",starttime,endtime, lat=e$latitude,long=e$longitude,maxradius=5) channels + } } \keyword{webservices} diff --git a/man/getChannel.Rd b/man/getChannel.Rd index fc8ca25..494ba44 100644 --- a/man/getChannel.Rd +++ b/man/getChannel.Rd @@ -25,7 +25,7 @@ getChannel(obj, network, station, location, channel, \item{channel}{ character string with the three letter channel code } \item{starttime}{ POSIXct class specifying the starttime (GMT) } \item{endtime}{ POSIXct class specifying the endtime (GMT) } - \item{includerestricted}{ optional logical identifying whether to report on restricted data (default=\code{FALSE}) } + \item{includerestricted}{ optional logical identifying whether to report on restricted data } \item{latitude}{ optional latitude used when specifying a location and radius } \item{longitude}{ optional longitude used when specifying a location and radius } \item{minradius}{ optional minimum radius used when specifying a location and radius } @@ -71,14 +71,14 @@ The IRIS DMC station webservice: \url{http://service.iris.edu/fdsnws/station/1/} -This implementation was inspired by the functionality in the obspy station() method. +This implementation was inspired by the functionality in the obspy get_stations() method. - \url{http://docs.obspy.org/packages/autogen/obspy.iris.client.Client.station.html} + \url{http://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_stations.html} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -90,6 +90,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -101,8 +102,9 @@ endtime <- starttime + 2*24*3600 events <- getEvent(iris,starttime,endtime,6.0) events -# Event #3 is the Nisqually quake in Washington State -e <- events[3,] +# biggest event is Nisqually +eIndex <- which(events$magnitude == max(events$magnitude)) +e <- events[eIndex[1],] # Which stations in the US network are within 5 degrees of the quake epicenter? stations <- getStation(iris,"US","*","*","BHZ",starttime,endtime, @@ -112,6 +114,7 @@ stations # Get some detailed information on any BHZ channels at the "Octopus Mountain" station channels <- getChannel(iris,"US","OCWA","*","BHZ",starttime,endtime) channels + } } \keyword{webservices} diff --git a/man/getDataselect.Rd b/man/getDataselect.Rd index 19b19d0..b69ccf0 100644 --- a/man/getDataselect.Rd +++ b/man/getDataselect.Rd @@ -1,9 +1,13 @@ \name{getDataselect} \alias{getDataselect} -\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,character,logical-method} -\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,character,missing-method} -\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,missing,logical-method} -\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,missing,missing-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,character,logical,logical-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,missing,missing,missing-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,character,logical,missing-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,character,missing,logical-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,missing,logical,logical-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,character,missing,missing-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,missing,missing,logical-method} +\alias{getDataselect,IrisClient,character,character,character,character,POSIXct,POSIXct,missing,logical,missing-method} \title{Retrieve seismic data from IRIS DMC} @@ -15,7 +19,7 @@ have been sorted by start time. \usage{ getDataselect(obj, network, station, location, channel, - starttime, endtime, quality, ignoreEpoch) + starttime, endtime, quality, inclusiveEnd, ignoreEpoch) } \arguments{ @@ -26,12 +30,13 @@ getDataselect(obj, network, station, location, channel, \item{channel}{character string with the three letter channel code} \item{starttime}{POSIXct class specifying the starttime (GMT)} \item{endtime}{POSIXct class specifying the endtime (GMT)} - \item{quality}{optional character string identifying the quality (default = "B")} + \item{quality}{optional character string identifying the quality. IRIS webservices defaults to quality="M".} + \item{inclusiveEnd}{optional logical determining whether the endtime is inclusive (default = TRUE)} \item{ignoreEpoch}{optional logical defining behavior when multiple epochs are encountered (default = FALSE)} } \details{ -This is the primary method for retreiving seismic data. Data requests are made through +This is the primary method for retrieving seismic data. Data requests are made through the \code{dataselect} webservice and returned data are parsed using the internal miniseed2Stream() function. @@ -39,6 +44,11 @@ If the location argument contains an empty string to specify a 'blank' location code of \code{"--"} will be used in the dataselect request URL. (See \href{http://service.iris.edu/fdsnws/dataselect/1/}{dataselect documentation}.) +If \code{inclusiveEnd=FALSE}, then \code{getDataselect} will subtract 0.000001 seconds from the endtime before passing the +value to the \code{dataselect} webservice. An endtime of, e.g., \code{as.POSIXct("2016-01-03", tz="GMT")} will be passed +into \code{dataselect} as \code{end=2016-01-02T23:59:59.999999}. A data sample at time 2016-01-03T00:00:00 will not be returned +unless \code{inclusiveEnd=TRUE}. + Error returns from the webservice will stop evaluation and generate an error message. Sometimes, the station webservice will return multiple records for the same SNCL, @@ -63,11 +73,11 @@ The IRIS DMC dataselect webservice: This implementation is similar in functionality to the obspy dataselect function: - \url{http://docs.obspy.org/_modules/obspy/iris/client.html#Client.dataselect} + \url{http://docs.obspy.org/_modules/obspy/clients/fdsn/client.html#Client.dataselect} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -79,6 +89,25 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} \code{\link{IrisClient-class}} } -%% \examples{ } +\examples{ + \dontrun{ +# Open a connection to IRIS DMC webservices +iris <- new("IrisClient") + +# Use getDataselect to request data for II.JTS.00.BHZ +starttime <- as.POSIXct("2001-02-28",tz="GMT") +endtime <- as.POSIXct("2001-03-01",tz="GMT") + +st <- getDataselect(iris, "II","JTS","00","BHZ",starttime, endtime, + inclusiveEnd=FALSE,ignoreEpoch=TRUE) + +# Display structure of trace(s) + +str(st) + +# Plot trace +plot(st) + } +} \keyword{webservices} diff --git a/man/getDistaz.Rd b/man/getDistaz.Rd index 4cf6750..3370e5e 100644 --- a/man/getDistaz.Rd +++ b/man/getDistaz.Rd @@ -24,8 +24,7 @@ getDistaz(obj, latitude, longitude, staLatitude, staLongitude) \details{ The distance-azimuth service will calculate the great-circle angular distance, azimuth, -and backazimuth between two geographic coordinate pairs. All results are reported in -degrees, with azimuth and backazimuth measured clockwise from North. +and back azimuth between two geographic coordinate pairs. Azimuth and back azimuth are measured clockwise from North. } @@ -33,10 +32,16 @@ degrees, with azimuth and backazimuth measured clockwise from North. A dataframe with the following columns: \preformatted{ -azimuth, backAzimuth, distance +ellipsoid.semiMajorAxis, ellipsoid.flattening, ellipsoid.name, fromlat, fromlon, tolat, tolon, + azimuth,backAzimuth, distance, distanceMeters } +Where \code{fromlat} is the event latitude, \code{fromlon} is the event longitude, \code{tolat} is the station latitude, and \code{tolon} +is the station longitude. \code{azimuth}, \code{backAzimuth}, and \code{distance} are measured in degrees. \code{distanceMeters} is distance in meters. +\code{ellipsoid.semiMajorAxis}, \code{ellipsoid.flattening}, and \code{ellipsoid.name} refer to the World Geodetic System standard coordinate system version used to correct for ellipticity when converting to geocentric latitudes. + Only a single row is returned. + } \references{ @@ -46,7 +51,7 @@ The IRIS DMC distaz webservice: } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/getEvalresp.Rd b/man/getEvalresp.Rd index 0faab7c..944c3a6 100644 --- a/man/getEvalresp.Rd +++ b/man/getEvalresp.Rd @@ -66,7 +66,7 @@ The IRIS DMC evalresp webservice: } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/getEvent.Rd b/man/getEvent.Rd index 2024285..f3f45fb 100644 --- a/man/getEvent.Rd +++ b/man/getEvent.Rd @@ -29,6 +29,8 @@ getEvent(obj, starttime, endtime, minmag, maxmag, magtype, \details{ The \code{getEvent} method uses the event web service to obtain data for all events that meet the criteria defined by the arguments and returns that data in a dataframe. Each row of the dataframe represents a unique event. + +\code{getEvent} calls to the IRIS event webservice now go to http://earthquake.usgs.gov/fdsnws/event/1/. } @@ -51,10 +53,14 @@ used for internal consistency -- all other web services return columns named "la The IRIS DMC event webservice: \url{http://service.iris.edu/fdsnws/event/1/} + +The USGS event webservice: + \url{http://earthquake.usgs.gov/fdsnws/event/1/} + } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/getGaps.Rd b/man/getGaps.Rd index d2ddd6d..82570e4 100644 --- a/man/getGaps.Rd +++ b/man/getGaps.Rd @@ -57,7 +57,7 @@ A list is returned with the following elements: %% \references{} \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{} @@ -66,6 +66,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} %% \seealso{} \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -73,7 +74,7 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) # Save the gap analysis in a variable gapInfo <- getGaps(st) @@ -83,8 +84,9 @@ names(gapInfo) # Look at a histogram of data dropouts hist(gapInfo$nsamples, breaks=50, - main="Data Gaps in AK.PIN..VEA during Jan 24, 2012", + main="Data Gaps in AK.PIN..BHZ Jan 24, 2012", xlab="number of missing samples per gap") + } } \keyword{methods} diff --git a/man/getNetwork.Rd b/man/getNetwork.Rd index 17276d0..5631be9 100644 --- a/man/getNetwork.Rd +++ b/man/getNetwork.Rd @@ -61,14 +61,14 @@ The IRIS DMC station web service: \url{http://service.iris.edu/fdsnws/station/1/} -This implementation was inspired by the functionality in the obspy station() method. +This implementation was inspired by the functionality in the obspy get_stations() method. - \url{http://docs.obspy.org/packages/autogen/obspy.iris.client.Client.station.html} + \url{http://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_stations.html} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -80,6 +80,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -91,13 +92,16 @@ endtime <- starttime + 2*24*3600 events <- getEvent(iris,starttime,endtime,6.0) events -# Event #3 is the Nisqually quake in Washington State -e <- events[3,] +# biggest event is Nisqually +eIndex <- which(events$magnitude == max(events$magnitude)) +e <- events[eIndex[1],] + # Which seismic networks have BHZ stations within 5 degrees of the quake epicenter? networks <- getNetwork(iris,"*","*","*","BHZ",starttime,endtime, lat=e$latitude,lon=e$longitude,maxradius=5) networks + } } \keyword{webservices} diff --git a/man/getRotation.Rd b/man/getRotation.Rd index 36eb289..17f21bb 100644 --- a/man/getRotation.Rd +++ b/man/getRotation.Rd @@ -54,7 +54,7 @@ The IRIS DMC rotation web service: } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/getSNCL.Rd b/man/getSNCL.Rd index 54a2a44..1911e35 100644 --- a/man/getSNCL.Rd +++ b/man/getSNCL.Rd @@ -1,7 +1,14 @@ \name{getSNCL} \docType{methods} \alias{getSNCL} -\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,missing-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,character,logical,logical-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,missing,missing,missing-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,character,logical,missing-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,character,missing,logical-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,missing,logical,logical-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,character,missing,missing-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,missing,missing,logical-method} +\alias{getSNCL,IrisClient,character,POSIXct,POSIXct,missing,logical,missing-method} \title{Retrieve seismic data from IRIS DMC} @@ -12,7 +19,7 @@ have been sorted by start time. } \usage{ -getSNCL(obj, sncl, starttime, endtime, quality) +getSNCL(obj, sncl, starttime, endtime, quality, inclusiveEnd, ignoreEpoch) } \arguments{ @@ -21,6 +28,8 @@ getSNCL(obj, sncl, starttime, endtime, quality) \item{starttime}{ POSIXct class specifying the starttime (GMT) } \item{endtime}{ POSIXct class specifying the endtime (GMT) } \item{quality}{ optional character string identifying the quality (default = "B") } + \item{inclusiveEnd}{optional logical determining whether the endtime is inclusive (default = TRUE)} + \item{ignoreEpoch}{optional logical defining behavior when multiple epochs are encountered (default = FALSE)} } \details{ @@ -41,7 +50,7 @@ The IRIS DMC dataselect web service: } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -53,6 +62,24 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} \code{\link{IrisClient-class}} } -%% \examples{ } +\examples{ + \dontrun{ +# Open a connection to IRIS DMC webservices +iris <- new("IrisClient") +# Use getDataselect to request data for II.JTS.00.BHZ +starttime <- as.POSIXct("2001-02-28",tz="GMT") +endtime <- as.POSIXct("2001-03-01",tz="GMT") + +st <- getSNCL(iris, "II.JTS.00.BHZ",starttime, endtime, quality="M", + inclusiveEnd=FALSE,ignoreEpoch=TRUE) + +# Display structure of trace(s) + +str(st) + +# Plot trace +plot(st) + } +} \keyword{webservices} diff --git a/man/getStation.Rd b/man/getStation.Rd index ce4752a..b83d3a5 100644 --- a/man/getStation.Rd +++ b/man/getStation.Rd @@ -25,7 +25,7 @@ getStation(obj, network, station, location, channel, \item{channel}{ character string with the three letter channel code } \item{starttime}{ POSIXct class specifying the starttime (GMT) } \item{endtime}{ POSIXct class specifying the endtime (GMT) } - \item{includerestricted}{ optional logical identifying whether to report on restricted data (default=\code{FALSE}) } + \item{includerestricted}{ optional logical identifying whether to report on restricted data } \item{latitude}{ optional latitude used when specifying a location and radius } \item{longitude}{ optional longitude used when specifying a location and radius } \item{minradius}{ optional minimum radius used when specifying a location and radius } @@ -62,14 +62,14 @@ The IRIS DMC station web service: \url{http://service.iris.edu/fdsnws/station/1/} -This implementation was inspired by the functionality in the obspy station() method. +This implementation was inspired by the functionality in the obspy get_stations() method. - \url{http://docs.obspy.org/packages/autogen/obspy.iris.client.Client.station.html} + \url{http://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_stations.html} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -81,6 +81,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -92,13 +93,15 @@ endtime <- starttime + 2*24*3600 events <- getEvent(iris,starttime,endtime,6.0) events -# Event #3 is the Nisqually quake in Washington State -e <- events[3,] +# biggest event is Nisqually +eIndex <- which(events$magnitude == max(events$magnitude)) +e <- events[eIndex[1],] # Which stations in the US network are within 5 degrees of the quake epicenter? stations <- getStation(iris,"US","*","*","BHZ",starttime,endtime, lat=e$latitude,long=e$longitude,maxradius=5) stations + } } \keyword{webservices} diff --git a/man/getTraveltime.Rd b/man/getTraveltime.Rd index b1f64ce..f90b07f 100644 --- a/man/getTraveltime.Rd +++ b/man/getTraveltime.Rd @@ -47,7 +47,7 @@ The IRIS DMC traveltime web service: } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -68,7 +68,7 @@ endtime <- starttime + 3600 * 24 * 2 # Find biggest seismic event over these two days -- it's the "Nisqually" events <- getEvent(iris, starttime, endtime, minmag=5.0) bigOneIndex <- which(events$magnitude == max(events$magnitude)) -bigOne <- events[bigOneIndex,] +bigOne <- events[bigOneIndex[1],] # Find US stations that are available within an hour of the event start <- bigOne$time diff --git a/man/getUnavailability.Rd b/man/getUnavailability.Rd index 912a053..c26b17d 100644 --- a/man/getUnavailability.Rd +++ b/man/getUnavailability.Rd @@ -7,7 +7,7 @@ \title{Retrieve Channel metadata from IRIS DMC} \description{ -The \code{getUnavailability} method obtains channel metadata for non-available channels from the IRIS DMC station web service +The \code{getUnavailability} method obtains metadata for channels that are not available from the IRIS DMC station web service and returns it in a dataframe. } @@ -62,14 +62,14 @@ The IRIS DMC station webservice: \url{http://service.iris.edu/fdsnws/station/1/} -This implementation was inspired by the functionality in the obspy station() method. +This implementation was inspired by the functionality in the obspy get_stations() method. - \url{http://docs.obspy.org/packages/autogen/obspy.iris.client.Client.station.html} + \url{http://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_stations.html} } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -81,6 +81,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -92,15 +93,16 @@ endtime <- starttime + 2*24*3600 events <- getEvent(iris,starttime,endtime,6.0) events -# Event #3 is the Nisqually quake in Washington State -e <- events[3,] +# biggest event is Nisqually +eIndex <- which(events$magnitude == max(events$magnitude)) +e <- events[eIndex[1],] # Find all BHZ channels that were NOT collecting data at the time of the quake -# and within 5 degrees of the quake epicenter? +# and within 5 degrees of the quake epicenter (or are otherwise unavailable from IRIS). channels <- getUnavailability(iris,"*","*","*","BHZ",starttime,endtime, lat=e$latitude,long=e$longitude,maxradius=5) channels - + } } \keyword{webservices} diff --git a/man/getUpDownTimes.Rd b/man/getUpDownTimes.Rd index 79bed55..c054d44 100644 --- a/man/getUpDownTimes.Rd +++ b/man/getUpDownTimes.Rd @@ -44,7 +44,7 @@ A vector of \code{POSIXct} datetimes associated with on/off transitions. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -55,6 +55,7 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -62,14 +63,14 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) # Determine up/down transitions, ignoring Traces < 3 min and gaps < 5 min upDownTimes <- getUpDownTimes(st, min_signal=180, min_gap=300) # Or just plot them directly plotUpDownTimes(st, min_signal=180, min_gap=300) - + } } \keyword{methods} diff --git a/man/hilbert.Rd b/man/hilbert.Rd index 3fa1599..4b2042d 100644 --- a/man/hilbert.Rd +++ b/man/hilbert.Rd @@ -35,7 +35,7 @@ whose data have been replaced with the Hilbert transform of the seismic signal. } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ @@ -58,9 +58,6 @@ endtime <- as.POSIXct("2010-02-27 09:00:00", tz="GMT") st <- getDataselect(iris,"IU","ANMO","00","BHZ",starttime,endtime) tr <- st@traces[[1]] -# Demean, detrend, cosine taper -tr <- DDT(tr) - # Create Hilbert transform of the trace trh <- hilbert(tr) diff --git a/man/hilbertFFT.Rd b/man/hilbertFFT.Rd index c26ae72..4f0957e 100644 --- a/man/hilbertFFT.Rd +++ b/man/hilbertFFT.Rd @@ -31,7 +31,7 @@ A complex vector containing the Hilbert FFT of \code{x}. } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ diff --git a/man/mergeTraces.Rd b/man/mergeTraces.Rd index e905a1d..0b8f9a1 100644 --- a/man/mergeTraces.Rd +++ b/man/mergeTraces.Rd @@ -36,7 +36,7 @@ A new \code{Stream} object containing a single \code{Trace} is returned. %% \references{} \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/mergeUpDownTimes.Rd b/man/mergeUpDownTimes.Rd index d959e9b..f2fbc50 100644 --- a/man/mergeUpDownTimes.Rd +++ b/man/mergeUpDownTimes.Rd @@ -41,7 +41,7 @@ A vector of \code{POSIXct} datetimes associated with on/off transitions. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ diff --git a/man/miniseed2Stream.Rd b/man/miniseed2Stream.Rd index ba83f5d..9b4ad19 100644 --- a/man/miniseed2Stream.Rd +++ b/man/miniseed2Stream.Rd @@ -11,7 +11,8 @@ object. \usage{ miniseed2Stream(miniseed,url,requestedStarttime,requestedEndtime, - sensor,scale,scaleunits) + sensor,scale,scalefreq,scaleunits,latitude,longitude, + elevation, depth, azimuth,dip) } \arguments{ @@ -19,13 +20,20 @@ miniseed2Stream(miniseed,url,requestedStarttime,requestedEndtime, \item{url}{character source location (see \link{getDataselect})} \item{requestedStarttime}{\code{POSIXct} time associated with the requested starttime (see \link{getDataselect})} \item{requestedEndtime}{\code{POSIXct} time associated with the requested endtime (see \link{getDataselect})} - \item{sensor}{character description of the Sensor associated with this SNCL (see \link{Trace})} + \item{sensor}{character description of the Sensor type associated with this Station-Network-Channel-Location (SNCL) (see \link{Trace})} \item{scale}{character description of the InstrumentSensitivity associated with this SNCL (see \link{Trace})} + \item{scalefreq}{numeric description of frequency at which the InstrumentSensitivity is correct, the SensitivityFrequency (see \link{Trace})} \item{scaleunits}{character description of the InputUnits associated with this SNCL (see \link{Trace})} + \item{latitude}{numeric latitude associated with this SNCL (see \link{Trace})} + \item{longitude}{numeric longitude associated with this SNCL (see \link{Trace})} + \item{elevation}{numeric elevation associated with this SNCL (see \link{Trace})} + \item{depth}{numeric depth associated with this SNCL (see \link{Trace})} + \item{azimuth}{numeric channel azimuth associated with this SNCL (see \link{Trace})} + \item{dip}{numeric channel dip associated with this SNCL (see \link{Trace})} } \details{ -This function takes raw bytes read in from a file or URL and converts them to a \code{Stream} object. +This function takes raw bytes read in from a file or URL and converts them to a \code{Stream} object. Metadata information is optional. This function is primarily for internal use. } @@ -35,7 +43,7 @@ A \code{Stream} object.} %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/multiplyBy.Rd b/man/multiplyBy.Rd index 6bcc746..f03c85e 100644 --- a/man/multiplyBy.Rd +++ b/man/multiplyBy.Rd @@ -33,7 +33,7 @@ A new \code{Trace} or \code{Stream} object is returned. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -41,24 +41,27 @@ Jonathan Callahan \email{jonathan.s.callahan@gmail.com} %% \seealso{ } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") -starttime <- as.POSIXct("2012-01-24", tz="GMT") -endtime <- as.POSIXct("2012-01-25", tz="GMT") +starttime <- as.POSIXct("2011-01-24", tz="GMT") +endtime <- as.POSIXct("2011-01-25", tz="GMT") # Get the waveform -stRaw <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +stRaw <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) -# obtain an instrument sensitivity value with getChannelMetadata) -# ... - -sensitivityValue <- 10.0 +# obtain an instrument sensitivity value with getChannel metadata) +c <- getChannel(iris, "AK","PIN","","BHZ",starttime, endtime) +sensitivityValue <- c$scale # convert raw data -st <- multiplyBy(stRaw, sensitivityValue) - +st <- multiplyBy(stRaw, 1/sensitivityValue) rmsVariance(st) + +# plot trace +plot(st, ylab=c$scaleunits) + } } \keyword{methods} diff --git a/man/noiseMatrix2PdfMatrix.Rd b/man/noiseMatrix2PdfMatrix.Rd index 5b60b8c..cdb88a0 100644 --- a/man/noiseMatrix2PdfMatrix.Rd +++ b/man/noiseMatrix2PdfMatrix.Rd @@ -38,7 +38,7 @@ A \code{matrix} is returned with one row for each power level (-250:-50 dB) and } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/noiseModels.Rd b/man/noiseModels.Rd index 275aea3..ec0658a 100644 --- a/man/noiseModels.Rd +++ b/man/noiseModels.Rd @@ -7,7 +7,7 @@ \description{ The \code{noiseModels} function returns the New High Noise Model and New Low Noise Model -from the Peterson paper. Values are returned for the specific frequencies specified in +from the Peterson paper referenced below. Values are returned for the specific frequencies specified in the \code{freq} argument. } @@ -27,13 +27,13 @@ low noise models, respectively. } \references{ -\href{http://www.mttmllr.com/ADS/DATA/peterson_usgs_seismic_noise_ofr93-322.pdf}{Observations of Modeling and Seismic Background Noise} (Peterson 1993) +\href{http://pubs.usgs.gov/of/1993/0322/report.pdf}{Observations of Modeling and Seismic Background Noise} (Peterson 1993) \href{http://pubs.usgs.gov/of/2005/1438/pdf/OFR-1438.pdf}{Seismic Noise Analysis System Using Power Spectral Density Probability Density Functions} (McNamara and Boaz 2005) } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/psdDF2NoiseMatrix.Rd b/man/psdDF2NoiseMatrix.Rd index 2daccd4..beedd09 100644 --- a/man/psdDF2NoiseMatrix.Rd +++ b/man/psdDF2NoiseMatrix.Rd @@ -8,16 +8,18 @@ \description{ The \code{psdDF2NoiseMatrix} function uses the \code{snclq} identifier associated with the first PSD in the dataframe to obtain instrument correction information at the specified frequencies -from the \code{\link{getEvalresp}} web service. This correction is applied to every PSD in +from the \code{\link{getEvalresp}} web service if instrumentation correction information is +not supplied as an argument. This correction is applied to every PSD in the dataframe and the now corrected PSD values are returned as a matrix. } \usage{ -psdDF2NoiseMatrix(DF) +psdDF2NoiseMatrix(DF, evalresp=NULL) } \arguments{ \item{DF}{a dataframe of PSDs obtained from the \code{getPSDMeasurements} method of \code{IrisClient}.} + \item{evalresp}{dataframe of freq, amp, phase information matching output of \code{getEvalresp}, optional.} } \details{ @@ -34,11 +36,11 @@ A \code{matrix} is returned with one row for each instrument-corrected PSD and o } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ -The incoming dataframe is checked to make sure that it represents only a single SNCL. An error is generated if +The incoming dataframe is checked to make sure that it represents only a single SNCL (Station-Network-Channel-Location). An error is generated if more than one is found. However, the \code{psdDF} is not checked to make sure that no changes to the instrument correction happened during the time period covered by the \code{psdDF}. This occurs at an 'epoch' boundary when an instrument is replaced. diff --git a/man/psdList.Rd b/man/psdList.Rd index 4bbd0ca..c195716 100644 --- a/man/psdList.Rd +++ b/man/psdList.Rd @@ -6,7 +6,7 @@ \title{Apply McNamara PSD algorithm to a seismic signal} \description{ -The \code{psdList} function breaks up a seismic \code{Stream} object into a series of +The \code{psdList} function subsets a seismic \code{Stream} object into a series of shorter segments with 50\% overlap and uses the \code{McNamaraPSD} method to return a smoothed (aka binned) Power Spectral Density (PSD) for each segment. } @@ -20,7 +20,7 @@ psdList(st) } \details{ -A \code{Stream} will be broken up into segments depending upon the snclq identifier associated +A \code{Stream} will be subset into segments depending upon the channel identifier (@stats@channel) associated with this seismic data. The binning frequencies are also channel dependent as exemplified in this code extract where \code{Z} is the segment length in seconds: @@ -31,7 +31,7 @@ if (stringr::str_detect(channel,"^L")) { Z <- 3 * 3600 loFreq <- 0.001 hiFreq <- 0.5 * tr_merged@stats@sampling_rate -} else if (stringr::str_detect(channel,"M")) { +} else if (stringr::str_detect(channel,"^M")) { Z <- 2 * 3600 loFreq <- 0.0025 hiFreq <- 0.5 * tr_merged@stats@sampling_rate @@ -61,7 +61,7 @@ freq, spec, snclq, starttime, endtime } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/psdList2NoiseMatrix.Rd b/man/psdList2NoiseMatrix.Rd index 429eebe..ed34855 100644 --- a/man/psdList2NoiseMatrix.Rd +++ b/man/psdList2NoiseMatrix.Rd @@ -8,16 +8,18 @@ \description{ The \code{psdList2NoiseMatrix} function uses the \code{snclq} identifier associated with the first PSD in the list to obtain instrument correction information at the specified frequencies -from the \code{\link{getEvalresp}} web service. This correction is applied to every PSD in -the list and the now corrected PSD values are returned as a matrix. +from the \code{\link{getEvalresp}} web service if instrumentation correction information is +not supplied as an argument. This correction is applied to every PSD in +the list and the now corrected PSD values are returned as a matrix. } \usage{ -psdList2NoiseMatrix(psdList) +psdList2NoiseMatrix(psdList, evalresp=NULL) } \arguments{ \item{psdList}{a list of PSDs generated by the \code{psdList} function} + \item{evalresp}{dataframe of freq, amp, phase information matching output of \code{getEvalresp}, optional} } %% \details{ } @@ -31,12 +33,12 @@ A \code{matrix} is returned with one row for each instrument-corrected PSD and o } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ The \code{psdList} function generates a \code{psdList} from a single \code{Stream} of data and should thus only -contain data for a single SNCL. However, the \code{psdList} is not checked to make sure that no changes to the +contain data for a single SNCL (Station-Network-Channel-Location). However, the \code{psdList} is not checked to make sure that no changes to the instrument correction happened during the time period covered by the \code{psdList}. This occurs at an 'epoch' boundary when an instrument is replaced. } diff --git a/man/psdPlot.Rd b/man/psdPlot.Rd index 15e0310..15e712c 100644 --- a/man/psdPlot.Rd +++ b/man/psdPlot.Rd @@ -11,16 +11,21 @@ or \code{psdDF} dataframe. } \usage{ -psdPlot(PSDs, style, showNoiseModel, showMaxMin, showMode, showMean, ...) +psdPlot(PSDs, style='psd', evalresp=NULL, ylo=-200, yhi=-50, showNoiseModel=TRUE, + showMaxMin=TRUE, showMode=TRUE, showMean=FALSE, showMedian=FALSE, ...) } \arguments{ \item{PSDs}{either a list as returned by \code{psdList} or a dataframe of PSD values obtained from the BSS} \item{style}{character identifier of plot type: \code{'psd'} plots PSD lines, \code{'pdf'} plots the \code{pdfMatrix}} + \item{evalresp}{dataframe of freq, amp, phase information matching output of \code{getEvalresp}, optional} + \item{ylo}{numeric setting lower limit of plot y-axis (default=\code{-200})} + \item{yhi}{numeric setting upper limit of plot y-axis (default=\code{-50})} \item{showNoiseModel}{logical controlling plotting of noise model lines (default=\code{TRUE})} \item{showMaxMin}{logical controlling plotting of PSD max and min lines (default=\code{TRUE})} \item{showMode}{logical controlling plotting of PDF mode line (default=\code{TRUE})} \item{showMean}{logical controlling plotting of PSD mean line (default=\code{FALSE})} + \item{showMedian}{logical controlling plotting of PSD median line (default=\code{FALSE})} \item{...}{arguments to be passed to plotting methods} } @@ -36,7 +41,7 @@ mimic the plots presented in the McNamara paper. } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/psdStatistics.Rd b/man/psdStatistics.Rd index 0e7d866..1b8738b 100644 --- a/man/psdStatistics.Rd +++ b/man/psdStatistics.Rd @@ -9,11 +9,12 @@ The \code{psdStatistics} function calculates a variety of information associated } \usage{ -psdStatistics(PSDs) +psdStatistics(PSDs, evalresp=NULL) } \arguments{ \item{PSDs}{either a list as returned by \code{psdList} or a dataframe of PSD values obtained from the BSS} + \item{evalresp}{dataframe of freq, amp, phase information matching output of \code{getEvalresp}, optional} } %%\details{ } @@ -46,11 +47,9 @@ A variety of plots can be generated form the information in this list. } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } -%% \note{ } - \seealso{ \code{\link{McNamaraPSD}}, \code{\link{psdList}}, @@ -75,9 +74,9 @@ stats <- psdStatistics(psdList) # Just for fun plot logPeriod <- log10(1/stats$freq) -plot(logPeriod,stats$max,ylim=c(-250,-50), las=1, - xlab="log10(period)", yalb="Power (dB)", - main="Model 'quiet' area and area of seismic signal.") +plot(logPeriod,stats$max,ylim=c(-200,-50), las=1, + xlab="log10(period)", ylab="Power (dB)", + main="Model 'normal background noise' area and area of seismic signal.") points(logPeriod,stats$min) # Overlay a polygon showing the range between the noise models diff --git a/man/readMiniseedFile.Rd b/man/readMiniseedFile.Rd index a78b82d..8dd5fe8 100644 --- a/man/readMiniseedFile.Rd +++ b/man/readMiniseedFile.Rd @@ -10,18 +10,26 @@ object. } \usage{ -readMiniseedFile(file,sensor,scale,scaleunits) +readMiniseedFile(file,sensor,scale,scalefreq,scaleunits, + latitude,longitude,elevation,depth,azimuth,dip) } \arguments{ \item{file}{character path of a miniSEED file} - \item{sensor}{character description of the Sensor associated with this SNCL (see \link{Trace})} + \item{sensor}{character description of the Sensor associated with this Station-Network-Channel-Location (SNCL) (see \link{Trace})} \item{scale}{character description of the InstrumentSensitivity associated with this SNCL (see \link{Trace})} + \item{scalefreq}{numeric description of frequency at which the InstrumentSensitivity is correct, the SensitivityFrequency (see \link{Trace})} \item{scaleunits}{character description of the InputUnits associated with this SNCL (see \link{Trace})} + \item{latitude}{numeric latitude associated with this SNCL (see \link{Trace})} + \item{longitude}{numeric longitude associated with this SNCL (see \link{Trace})} + \item{elevation}{numeric elevation associated with this SNCL (see \link{Trace})} + \item{depth}{numeric depth associated with this SNCL (see \link{Trace})} + \item{azimuth}{numeric channel azimuth associated with this SNCL (see \link{Trace})} + \item{dip}{numeric channel dip associated with this SNCL (see \link{Trace})} } \details{ -This function reads in a raw miniSEED file and converts it to a \code{Stream} object. +This function reads in a raw miniSEED file and converts it to a \code{Stream} object. Metadata information is optional. } \value{ @@ -30,7 +38,7 @@ A \code{Stream} object.} %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/rmsVariance.Rd b/man/rmsVariance.Rd index 714b563..d3d2902 100644 --- a/man/rmsVariance.Rd +++ b/man/rmsVariance.Rd @@ -19,7 +19,7 @@ \description{ The \code{rms} and \code{rmsVariance} methods of \code{Trace} and \code{Stream} objects compute the -Root Mean Square amplitude or RMS variance of the associated data in each object. +Root Mean Square (RMS) amplitude or RMS variance of the associated data in each object. RMS variance removes the DC level from the seismic signal so that the zero line is consistent. } @@ -73,7 +73,7 @@ A numeric vector is returned for \code{parallelRmsVariance}. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } @@ -89,7 +89,7 @@ starttime <- as.POSIXct("2012-01-24", tz="GMT") endtime <- as.POSIXct("2012-01-25", tz="GMT") # Get the waveform -st <- getDataselect(iris,"AK","PIN","","VEA",starttime,endtime) +st <- getDataselect(iris,"AK","PIN","","BHZ",starttime,endtime) # Get the first trace and generate some statistics tr <- st@traces[[1]] diff --git a/man/rotate2D.Rd b/man/rotate2D.Rd index a11f5dd..0a8972c 100644 --- a/man/rotate2D.Rd +++ b/man/rotate2D.Rd @@ -3,11 +3,11 @@ \alias{rotate2D} -\title{Rotate horiztonal components of a seismic signal} +\title{Rotate horizontal components of a seismic signal} \description{ The rotate2D() function rotates the two horizontal components of a -seismic signal into Radial and Transverse compoents returned as a list of 2 \code{Stream} objects. +seismic signal into Radial and Transverse components returned as a list of 2 \code{Stream} objects. } \usage{ @@ -15,8 +15,8 @@ rotate2D(st1, st2, angle) } \arguments{ - \item{st1}{horizontal ('N' or '1') \code{Stream} from a channel set} - \item{st2}{horizontal ('E' or '2') \code{Stream} from a channel set} + \item{st1}{horizontal \code{Stream} from a channel set (channel name usually ending in "N", "E", "1", or "2")} + \item{st2}{horizontal \code{Stream} from a channel set, complementary to st1} \item{angle}{angle (degrees) of the rotation} } @@ -31,16 +31,22 @@ The rotation service uses the following transformation matrix to change the outp \deqn{ \left[ \begin{array}{c} R \cr T \end{array} \right] = M_{2D} \left[ \begin{array}{c} N \cr E \end{array} \right] } where : - - N and E represent data from the original (horizontal) orientations - R and T represent the Radial and Transverse components. - \eqn{\alpha} is the azimuth angle measured clockwise from north + \preformatted{ + N and E represent data from the original (horizontal) orientations. + R and T represent the Radial and Transverse components. + \eqn{\alpha} is the azimuth angle measured clockwise from north. + } } \value{ A list of two \code{Stream} objects \code{stR} and \code{stT} is returned. } +\note{ +N and E are determined by the \code{Stream} @stats@azimuth values. If \code{Stream} @stats@azimuth values are not defined, +st1 is assumed to be N and st2 is assumed to be E. Orthogonality is also assumed to be correct. +} + \references{ %% ~put references to the literature/web site here ~ IRIS DMC rotation web service documentation: @@ -49,7 +55,7 @@ IRIS DMC rotation web service documentation: } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/slice.Rd b/man/slice.Rd index 5af788c..26b5655 100644 --- a/man/slice.Rd +++ b/man/slice.Rd @@ -36,7 +36,7 @@ A new \code{Trace} or \code{Stream} object is returned. %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/surfaceDistance.Rd b/man/surfaceDistance.Rd index 1078687..4cf7f54 100644 --- a/man/surfaceDistance.Rd +++ b/man/surfaceDistance.Rd @@ -32,7 +32,7 @@ Distance in kilometers } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } %% \note{ } diff --git a/man/triggerOnset.Rd b/man/triggerOnset.Rd index be51646..fab609e 100644 --- a/man/triggerOnset.Rd +++ b/man/triggerOnset.Rd @@ -42,7 +42,7 @@ The returned value wil be a \code{POSIXct} time by defult or a numeric index if %% \references{ } \author{ -Jonathan Callahan \email{jonathan.s.callahan@gmail.com} +Jonathan Callahan \email{jonathan@mazamascience.com} } \note{ @@ -59,6 +59,7 @@ signal. } \examples{ + \dontrun{ # Open a connection to IRIS DMC webservices iris <- new("IrisClient") @@ -75,6 +76,7 @@ to <- triggerOnset(tr,picker) plot(tr) abline(v=to, col='red', lwd=2) + } } \keyword{methods} diff --git a/man/unHistogram.Rd b/man/unHistogram.Rd new file mode 100644 index 0000000..f45b81b --- /dev/null +++ b/man/unHistogram.Rd @@ -0,0 +1,41 @@ +\name{unHistogram} + +\alias{unHistogram} + +\title{Histogram to Vector} + +\description{ +If \code{vec} represents a set of binned counts of incrementing values (ascending) +return a vector of associated bin values with the proper count of each value. Intended +for internal use. +} + +\usage{ +unHistogram(vec, startVal, incr) +} + +\arguments{ + \item{vec}{a histogram vector or ordered set of binned counts} + \item{startVal}{the initial value of the first bin element} + \item{incr}{the increment rate of each subsequent bin value} +} + +%% \details{ } + +\value{ +A vector of bin values with appropriate counts of each. +} + +%% \references{ } + +\author{ +Rob Casey \email{rob@iris.washington.edu} +} + +%% \note{ } + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +%% \seealso{ } + +\keyword{methods} diff --git a/src/Makevars.win b/src/Makevars.win index 9b026fb..2e6549a 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -8,7 +8,7 @@ PKG_LIBS=-Llibmseed -lmseed $(SHLIB): libmseed/libmseed.a libmseed/libmseed.a: - @(cd libmseed && $(MAKE) static AR="$(AR)" CC="$(CC)" CFLAGS="$(CFLAGS)") + @(cd libmseed && $(MAKE) static CC="$(CC)" CFLAGS="$(CFLAGS)" AR="$(AR)") clean: - @(cd libmseed && $(MAKE) clean) + @(cd libmseed && $(MAKE) clean ) diff --git a/src/libmseed/Makefile b/src/libmseed/Makefile index 7f416b3..e7d11e3 100644 --- a/src/libmseed/Makefile +++ b/src/libmseed/Makefile @@ -73,8 +73,8 @@ install: # Standard object building .c.o: - $(CC) $(CFLAGS) -fPIC -c $< -o $@ + $(CC) $(CFLAGS) -Wno-pedantic-ms-format -fPIC -c $< -o $@ # Standard object building for dynamic library components using -fPIC .c.lo: - $(CC) $(CFLAGS) -fPIC -c $< -o $@ + $(CC) $(CFLAGS) -Wno-pedantic-ms-format -fPIC -c $< -o $@ diff --git a/src/libmseed/fileutils.c b/src/libmseed/fileutils.c index af1a9cb..50e47bf 100644 --- a/src/libmseed/fileutils.c +++ b/src/libmseed/fileutils.c @@ -15,6 +15,8 @@ #include #include #include +#include +#include #include "libmseed.h" diff --git a/src/libmseed/genutils.c b/src/libmseed/genutils.c index 4258b12..d09ba51 100644 --- a/src/libmseed/genutils.c +++ b/src/libmseed/genutils.c @@ -15,6 +15,8 @@ #include #include #include +#include +#include #include "lmplatform.h" #include "libmseed.h" diff --git a/src/libmseed/libmseed.h b/src/libmseed/libmseed.h index 0ea3023..c2b426d 100644 --- a/src/libmseed/libmseed.h +++ b/src/libmseed/libmseed.h @@ -450,6 +450,8 @@ typedef struct MSTraceSeg_s { void *datasamples; /* Data samples, 'numsamples' of type 'sampletype'*/ int64_t numsamples; /* Number of data samples in datasamples */ char sampletype; /* Sample type code: a, i, f, d */ + /* REC */ + struct MSSampRateList_s* samprate_list; /* REC - pointer to linked list of sample rate counts */ void *prvtptr; /* Private pointer for general use, unused by libmseed */ struct MSTraceSeg_s *prev; /* Pointer to previous segment */ struct MSTraceSeg_s *next; /* Pointer to next segment */ @@ -483,6 +485,14 @@ typedef struct MSTraceList_s { } MSTraceList; +/* REC -- adding a linked list struct for MSTraceSeg */ +typedef struct MSSampRateList_s { + double samprate; /* sample rate key value for this node */ + int count; /* count value for this sample rate */ + struct MSSampRateList_s* next; /* next node in the linked list */ +} +MSSampRateList; + /* Data selection structure time window definition containers */ typedef struct SelectTime_s { hptime_t starttime; /* Earliest data for matching channels */ @@ -604,6 +614,10 @@ extern void mstl_printtracelist ( MSTraceList *mstl, flag timeformat, extern void mstl_printsynclist ( MSTraceList *mstl, char *dccid, flag subsecond ); extern void mstl_printgaplist (MSTraceList *mstl, flag timeformat, double *mingap, double *maxgap); +/* REC - in tracelist.c */ +MSSampRateList* add_segsamprate(MSSampRateList* srl, double samprate); +double get_segsamprate_mode(MSSampRateList* srl); +/* */ /* Reading Mini-SEED records from files */ typedef struct MSFileParam_s diff --git a/src/libmseed/tracelist.c b/src/libmseed/tracelist.c index 930981c..3fc3188 100644 --- a/src/libmseed/tracelist.c +++ b/src/libmseed/tracelist.c @@ -5,7 +5,7 @@ * * Written by Chad Trabant, IRIS Data Management Center * - * modified: 2015.108 + * modified: 2016/04/28 - Robert Casey, IRIS DMC - (REC) ***************************************************************************/ #include @@ -33,22 +33,22 @@ MSTraceSeg *mstl_addsegtoseg (MSTraceSeg *seg1, MSTraceSeg *seg2); MSTraceList * mstl_init ( MSTraceList *mstl ) { - if ( mstl ) - { - mstl_free (&mstl, 1); - } - - mstl = (MSTraceList *) malloc (sizeof(MSTraceList)); - - if ( mstl == NULL ) - { - ms_log (2, "mstl_init(): Cannot allocate memory\n"); - return NULL; - } - - memset (mstl, 0, sizeof (MSTraceList)); - - return mstl; + if ( mstl ) + { + mstl_free (&mstl, 1); + } + + mstl = (MSTraceList *) malloc (sizeof(MSTraceList)); + + if ( mstl == NULL ) + { + ms_log (2, "mstl_init(): Cannot allocate memory\n"); + return NULL; + } + + memset (mstl, 0, sizeof (MSTraceList)); + + return mstl; } /* End of mstl_init() */ @@ -64,54 +64,63 @@ mstl_init ( MSTraceList *mstl ) void mstl_free ( MSTraceList **ppmstl, flag freeprvtptr ) { - MSTraceID *id = 0; - MSTraceID *nextid = 0; - MSTraceSeg *seg = 0; - MSTraceSeg *nextseg = 0; - - if ( ! ppmstl ) - return; - - if ( *ppmstl ) - { - /* Free any associated traces */ - id = (*ppmstl)->traces; - while ( id ) + MSTraceID *id = 0; + MSTraceID *nextid = 0; + MSTraceSeg *seg = 0; + MSTraceSeg *nextseg = 0; + + if ( ! ppmstl ) + return; + + if ( *ppmstl ) { - nextid = id->next; - - /* Free any associated trace segments */ - seg = id->first; - while ( seg ) - { - nextseg = seg->next; - - /* Free private pointer data if present and requested*/ - if ( freeprvtptr && seg->prvtptr ) - free (seg->prvtptr); - - /* Free data array if allocated */ - if ( seg->datasamples ) - free (seg->datasamples); - - free (seg); - seg = nextseg; - } - - /* Free private pointer data if present and requested*/ - if ( freeprvtptr && id->prvtptr ) - free (id->prvtptr); - - free (id); - id = nextid; + /* Free any associated traces */ + id = (*ppmstl)->traces; + while ( id ) + { + nextid = id->next; + + /* Free any associated trace segments */ + seg = id->first; + while ( seg ) + { + nextseg = seg->next; + + /* Free private pointer data if present and requested*/ + if ( freeprvtptr && seg->prvtptr ) + free (seg->prvtptr); + + /* Free data array if allocated */ + if ( seg->datasamples ) + free (seg->datasamples); + + /* REC - free samprate list if allocated */ + MSSampRateList *sr_ptr, *sr_free; + sr_ptr = sr_free = seg->samprate_list; + while (sr_ptr != NULL) { + sr_ptr = sr_ptr->next; + free (sr_free); + sr_free = sr_ptr; + } + + free (seg); + seg = nextseg; + } + + /* Free private pointer data if present and requested*/ + if ( freeprvtptr && id->prvtptr ) + free (id->prvtptr); + + free (id); + id = nextid; + } + + free (*ppmstl); + + *ppmstl = NULL; } - - free (*ppmstl); - - *ppmstl = NULL; - } - - return; + + return; } /* End of mstl_free() */ @@ -138,468 +147,468 @@ mstl_free ( MSTraceList **ppmstl, flag freeprvtptr ) ***************************************************************************/ MSTraceSeg * mstl_addmsr ( MSTraceList *mstl, MSRecord *msr, flag dataquality, - flag autoheal, double timetol, double sampratetol ) + flag autoheal, double timetol, double sampratetol ) { - MSTraceID *id = 0; - MSTraceID *searchid = 0; - MSTraceID *ltid = 0; - - MSTraceSeg *seg = 0; - MSTraceSeg *searchseg = 0; - MSTraceSeg *segbefore = 0; - MSTraceSeg *segafter = 0; - MSTraceSeg *followseg = 0; - - hptime_t endtime; - hptime_t pregap; - hptime_t postgap; - hptime_t lastgap; - hptime_t firstgap; - hptime_t hpdelta; - hptime_t hptimetol = 0; - hptime_t nhptimetol = 0; - - char srcname[45]; - char *s1, *s2; - flag whence; - flag lastratecheck; - flag firstratecheck; - int mag; - int cmp; - int ltmag; - int ltcmp; - - if ( ! mstl || ! msr ) - return 0; - - /* Calculate end time for MSRecord */ - if ( (endtime = msr_endtime (msr)) == HPTERROR ) - { - ms_log (2, "mstl_addmsr(): Error calculating record end time\n"); - return 0; - } - - /* Generate source name string */ - if ( ! msr_srcname (msr, srcname, dataquality) ) - { - ms_log (2, "mstl_addmsr(): Error generating srcname for MSRecord\n"); - return 0; - } - - /* Search for matching trace ID starting with last accessed ID and - then looping through the trace ID list. */ - if ( mstl->last ) - { - s1 = mstl->last->srcname; - s2 = srcname; - while ( *s1 == *s2++ ) + MSTraceID *id = 0; + MSTraceID *searchid = 0; + MSTraceID *ltid = 0; + + MSTraceSeg *seg = 0; + MSTraceSeg *searchseg = 0; + MSTraceSeg *segbefore = 0; + MSTraceSeg *segafter = 0; + MSTraceSeg *followseg = 0; + + hptime_t endtime; + hptime_t pregap; + hptime_t postgap; + hptime_t lastgap; + hptime_t firstgap; + hptime_t hpdelta; + hptime_t hptimetol = 0; + hptime_t nhptimetol = 0; + + char srcname[45]; + char *s1, *s2; + flag whence; + flag lastratecheck; + flag firstratecheck; + int mag; + int cmp; + int ltmag; + int ltcmp; + + if ( ! mstl || ! msr ) + return 0; + + /* Calculate end time for MSRecord */ + if ( (endtime = msr_endtime (msr)) == HPTERROR ) { - if ( *s1++ == '\0' ) - break; + ms_log (2, "mstl_addmsr(): Error calculating record end time\n"); + return 0; } - cmp = (*s1 - *--s2); - - if ( ! cmp ) + + /* Generate source name string */ + if ( ! msr_srcname (msr, srcname, dataquality) ) { - id = mstl->last; + ms_log (2, "mstl_addmsr(): Error generating srcname for MSRecord\n"); + return 0; } - else + + /* Search for matching trace ID starting with last accessed ID and + then looping through the trace ID list. */ + if ( mstl->last ) { - /* Loop through trace ID list searching for a match, simultaneously - * track the source name which is closest but less than the MSRecord - * to allow for later insertion with sort order. */ - searchid = mstl->traces; - ltcmp = 0; - ltmag = 0; - while ( searchid ) - { - /* Compare source names */ - s1 = searchid->srcname; - s2 = srcname; - mag = 0; - while ( *s1 == *s2++ ) + s1 = mstl->last->srcname; + s2 = srcname; + while ( *s1 == *s2++ ) { - mag++; - if ( *s1++ == '\0' ) - break; + if ( *s1++ == '\0' ) + break; } - cmp = (*s1 - *--s2); - - /* If source names did not match track closest "less than" value - * and continue searching. */ - if ( cmp != 0 ) + cmp = (*s1 - *--s2); + + if ( ! cmp ) { - if ( cmp < 0 ) - { - if ( (ltcmp == 0 || cmp >= ltcmp) && mag >= ltmag ) - { - ltcmp = cmp; - ltmag = mag; - ltid = searchid; - } - else if ( mag > ltmag ) + id = mstl->last; + } + else + { + /* Loop through trace ID list searching for a match, simultaneously + * track the source name which is closest but less than the MSRecord + * to allow for later insertion with sort order. */ + searchid = mstl->traces; + ltcmp = 0; + ltmag = 0; + while ( searchid ) { - ltcmp = cmp; - ltmag = mag; - ltid = searchid; + /* Compare source names */ + s1 = searchid->srcname; + s2 = srcname; + mag = 0; + while ( *s1 == *s2++ ) + { + mag++; + if ( *s1++ == '\0' ) + break; + } + cmp = (*s1 - *--s2); + + /* If source names did not match track closest "less than" value + * and continue searching. */ + if ( cmp != 0 ) + { + if ( cmp < 0 ) + { + if ( (ltcmp == 0 || cmp >= ltcmp) && mag >= ltmag ) + { + ltcmp = cmp; + ltmag = mag; + ltid = searchid; + } + else if ( mag > ltmag ) + { + ltcmp = cmp; + ltmag = mag; + ltid = searchid; + } + } + + searchid = searchid->next; + continue; + } + + /* If we made it this far we found a match */ + id = searchid; + break; } - } - - searchid = searchid->next; - continue; } - - /* If we made it this far we found a match */ - id = searchid; - break; - } - } - } /* Done searching for match in trace ID list */ - - /* If no matching ID was found create new MSTraceID and MSTraceSeg entries */ - if ( ! id ) - { - if ( ! (id = (MSTraceID *) calloc (1, sizeof(MSTraceID))) ) - { - ms_log (2, "mstl_addmsr(): Error allocating memory\n"); - return 0; - } - - /* Populate MSTraceID */ - strcpy (id->network, msr->network); - strcpy (id->station, msr->station); - strcpy (id->location, msr->location); - strcpy (id->channel, msr->channel); - id->dataquality = msr->dataquality; - strcpy (id->srcname, srcname); - - id->earliest = msr->starttime; - id->latest = endtime; - id->numsegments = 1; - - if ( ! (seg = mstl_msr2seg (msr, endtime)) ) - { - return 0; - } - id->first = id->last = seg; - - /* Add new MSTraceID to MSTraceList */ - if ( ! mstl->traces || ! ltid ) - { - id->next = mstl->traces; - mstl->traces = id; - } - else - { - id->next = ltid->next; - ltid->next = id; - } - - mstl->numtraces++; - } - /* Add data coverage to the matching MSTraceID */ - else - { - /* Calculate high-precision sample period */ - hpdelta = (hptime_t) (( msr->samprate ) ? (HPTMODULUS / msr->samprate) : 0.0); - - /* Calculate high-precision time tolerance */ - if ( timetol == -1.0 ) - hptimetol = (hptime_t) (0.5 * hpdelta); /* Default time tolerance is 1/2 sample period */ - else if ( timetol >= 0.0 ) - hptimetol = (hptime_t) (timetol * HPTMODULUS); - - nhptimetol = ( hptimetol ) ? -hptimetol : 0; - - /* last/firstgap are negative when the record overlaps the trace - * segment and positive when there is a time gap. */ - - /* Gap relative to the last segment */ - lastgap = msr->starttime - id->last->endtime - hpdelta; - - /* Gap relative to the first segment */ - firstgap = id->first->starttime - endtime - hpdelta; - - /* Sample rate tolerance checks for first and last segments */ - if ( sampratetol == -1.0 ) - { - lastratecheck = MS_ISRATETOLERABLE (msr->samprate, id->last->samprate); - firstratecheck = MS_ISRATETOLERABLE (msr->samprate, id->first->samprate); - } - else - { - lastratecheck = (ms_dabs (msr->samprate - id->last->samprate) > sampratetol) ? 0 : 1; - firstratecheck = (ms_dabs (msr->samprate - id->first->samprate) > sampratetol) ? 0 : 1; - } - - /* Search first for the simple scenarios in order of likelihood: - * - Record fits at end of last segment - * - Record fits after all coverage - * - Record fits before all coverage - * - Record fits at beginning of first segment - * - * If none of those scenarios are true search the complete segment list. - */ - - /* Record coverage fits at end of last segment */ - if ( lastgap <= hptimetol && lastgap >= nhptimetol && lastratecheck ) - { - if ( ! mstl_addmsrtoseg (id->last, msr, endtime, 1) ) - return 0; - - seg = id->last; - - if ( endtime > id->latest ) - id->latest = endtime; - } - /* Record coverage is after all other coverage */ - else if ( (msr->starttime - hpdelta - hptimetol) > id->latest ) - { - if ( ! (seg = mstl_msr2seg (msr, endtime)) ) - return 0; - - /* Add to end of list */ - id->last->next = seg; - seg->prev = id->last; - id->last = seg; - id->numsegments++; - - if ( endtime > id->latest ) - id->latest = endtime; - } - /* Record coverage is before all other coverage */ - else if ( (endtime + hpdelta + hptimetol) < id->earliest ) - { - if ( ! (seg = mstl_msr2seg (msr, endtime)) ) - return 0; - - /* Add to beginning of list */ - id->first->prev = seg; - seg->next = id->first; - id->first = seg; - id->numsegments++; - - if ( msr->starttime < id->earliest ) - id->earliest = msr->starttime; - } - /* Record coverage fits at beginning of first segment */ - else if ( firstgap <= hptimetol && firstgap >= nhptimetol && firstratecheck ) - { - if ( ! mstl_addmsrtoseg (id->first, msr, endtime, 2) ) - return 0; - - seg = id->first; - - if ( msr->starttime < id->earliest ) - id->earliest = msr->starttime; - } - /* Search complete segment list for matches */ - else + } /* Done searching for match in trace ID list */ + + /* If no matching ID was found create new MSTraceID and MSTraceSeg entries */ + if ( ! id ) { - searchseg = id->first; - segbefore = 0; /* Find segment that record fits before */ - segafter = 0; /* Find segment that record fits after */ - followseg = 0; /* Track segment that record follows in time order */ - while ( searchseg ) - { - if ( msr->starttime > searchseg->starttime ) - followseg = searchseg; - - whence = 0; - - postgap = msr->starttime - searchseg->endtime - hpdelta; - if ( ! segbefore && postgap <= hptimetol && postgap >= nhptimetol ) - whence = 1; - - pregap = searchseg->starttime - endtime - hpdelta; - if ( ! segafter && pregap <= hptimetol && pregap >= nhptimetol ) - whence = 2; - - if ( ! whence ) + if ( ! (id = (MSTraceID *) calloc (1, sizeof(MSTraceID))) ) { - searchseg = searchseg->next; - continue; + ms_log (2, "mstl_addmsr(): Error allocating memory\n"); + return 0; } - - if ( sampratetol == -1.0 ) + + /* Populate MSTraceID */ + strcpy (id->network, msr->network); + strcpy (id->station, msr->station); + strcpy (id->location, msr->location); + strcpy (id->channel, msr->channel); + id->dataquality = msr->dataquality; + strcpy (id->srcname, srcname); + + id->earliest = msr->starttime; + id->latest = endtime; + id->numsegments = 1; + + if ( ! (seg = mstl_msr2seg (msr, endtime)) ) { - if ( ! MS_ISRATETOLERABLE (msr->samprate, searchseg->samprate) ) - { - searchseg = searchseg->next; - continue; - } + return 0; } - else + id->first = id->last = seg; + + /* Add new MSTraceID to MSTraceList */ + if ( ! mstl->traces || ! ltid ) { - if ( ms_dabs (msr->samprate - searchseg->samprate) > sampratetol ) - { - searchseg = searchseg->next; - continue; - } + id->next = mstl->traces; + mstl->traces = id; } - - if ( whence == 1 ) - segbefore = searchseg; - else - segafter = searchseg; - - /* Done searching if not autohealing */ - if ( ! autoheal ) - break; - - /* Done searching if both before and after segments are found */ - if ( segbefore && segafter ) - break; - - searchseg = searchseg->next; - } /* Done looping through segments */ - - /* Add MSRecord coverage to end of segment before */ - if ( segbefore ) - { - if ( ! mstl_addmsrtoseg (segbefore, msr, endtime, 1) ) + else { - return 0; + id->next = ltid->next; + ltid->next = id; } - - /* Merge two segments that now fit if autohealing */ - if ( autoheal && segafter && segbefore != segafter ) + + mstl->numtraces++; + } + /* Add data coverage to the matching MSTraceID */ + else + { + /* Calculate high-precision sample period */ + hpdelta = (hptime_t) (( msr->samprate ) ? (HPTMODULUS / msr->samprate) : 0.0); + + /* Calculate high-precision time tolerance */ + if ( timetol == -1.0 ) + hptimetol = (hptime_t) (0.5 * hpdelta); /* Default time tolerance is 1/2 sample period */ + else if ( timetol >= 0.0 ) + hptimetol = (hptime_t) (timetol * HPTMODULUS); + + nhptimetol = ( hptimetol ) ? -hptimetol : 0; + + /* last/firstgap are negative when the record overlaps the trace + * segment and positive when there is a time gap. */ + + /* Gap relative to the last segment */ + lastgap = msr->starttime - id->last->endtime - hpdelta; + + /* Gap relative to the first segment */ + firstgap = id->first->starttime - endtime - hpdelta; + + /* Sample rate tolerance checks for first and last segments */ + if ( sampratetol == -1.0 ) { - /* Add segafter coverage to segbefore */ - if ( ! mstl_addsegtoseg (segbefore, segafter) ) - { - return 0; - } - - /* Shift last segment pointer if it's going to be removed */ - if ( segafter == id->last ) - id->last = id->last->prev; - - /* Remove segafter from list */ - if ( segafter->prev ) - segafter->prev->next = segafter->next; - if ( segafter->next ) - segafter->next->prev = segafter->prev; - - /* Free data samples, private data and segment structure */ - if (segafter->datasamples) - free (segafter->datasamples); - - if (segafter->prvtptr) - free (segafter->prvtptr); - - free (segafter); + lastratecheck = MS_ISRATETOLERABLE (msr->samprate, id->last->samprate); + firstratecheck = MS_ISRATETOLERABLE (msr->samprate, id->first->samprate); } - - seg = segbefore; - } - /* Add MSRecord coverage to beginning of segment after */ - else if ( segafter ) - { - if ( ! mstl_addmsrtoseg (segafter, msr, endtime, 2) ) + else { - return 0; + lastratecheck = (ms_dabs (msr->samprate - id->last->samprate) > sampratetol) ? 0 : 1; + firstratecheck = (ms_dabs (msr->samprate - id->first->samprate) > sampratetol) ? 0 : 1; } - - seg = segafter; - } - /* Add MSRecord coverage to new segment */ - else - { - /* Create new segment */ - if ( ! (seg = mstl_msr2seg (msr, endtime)) ) + + /* Search first for the simple scenarios in order of likelihood: + * - Record fits at end of last segment + * - Record fits after all coverage + * - Record fits before all coverage + * - Record fits at beginning of first segment + * + * If none of those scenarios are true search the complete segment list. + */ + + /* Record coverage fits at end of last segment */ + if ( lastgap <= hptimetol && lastgap >= nhptimetol && lastratecheck ) { - return 0; + if ( ! mstl_addmsrtoseg (id->last, msr, endtime, 1) ) + return 0; + + seg = id->last; + + if ( endtime > id->latest ) + id->latest = endtime; } - - /* Add new segment as first in list */ - if ( ! followseg ) + /* Record coverage is after all other coverage */ + else if ( (msr->starttime - hpdelta - hptimetol) > id->latest ) { - seg->next = id->first; - if ( id->first ) - id->first->prev = seg; - - id->first = seg; + if ( ! (seg = mstl_msr2seg (msr, endtime)) ) + return 0; + + /* Add to end of list */ + id->last->next = seg; + seg->prev = id->last; + id->last = seg; + id->numsegments++; + + if ( endtime > id->latest ) + id->latest = endtime; + } + /* Record coverage is before all other coverage */ + else if ( (endtime + hpdelta + hptimetol) < id->earliest ) + { + if ( ! (seg = mstl_msr2seg (msr, endtime)) ) + return 0; + + /* Add to beginning of list */ + id->first->prev = seg; + seg->next = id->first; + id->first = seg; + id->numsegments++; + + if ( msr->starttime < id->earliest ) + id->earliest = msr->starttime; } - /* Add new segment after the followseg segment */ - else + /* Record coverage fits at beginning of first segment */ + else if ( firstgap <= hptimetol && firstgap >= nhptimetol && firstratecheck ) { - seg->next = followseg->next; - seg->prev = followseg; - if ( followseg->next ) - followseg->next->prev = seg; - followseg->next = seg; - - if ( followseg == id->last ) - id->last = seg; + if ( ! mstl_addmsrtoseg (id->first, msr, endtime, 2) ) + return 0; + + seg = id->first; + + if ( msr->starttime < id->earliest ) + id->earliest = msr->starttime; } - - id->numsegments++; - } - - /* Track earliest and latest times */ - if ( msr->starttime < id->earliest ) - id->earliest = msr->starttime; - - if ( endtime > id->latest ) - id->latest = endtime; - } /* End of searching segment list */ - } /* End of adding coverage to matching ID */ - - /* Sort modified segment into place, logic above should limit these to few shifts if any */ - while ( seg->next && ( seg->starttime > seg->next->starttime || - (seg->starttime == seg->next->starttime && seg->endtime < seg->next->endtime) ) ) - { - /* Move segment down list, swap seg and seg->next */ - segafter = seg->next; - - if ( seg->prev ) - seg->prev->next = segafter; - - if ( segafter->next ) - segafter->next->prev = seg; - - segafter->prev = seg->prev; - seg->prev = segafter; - seg->next = segafter->next; - segafter->next = seg; - - /* Reset first and last segment pointers if replaced */ - if ( id->first == seg ) - id->first = segafter; - - if ( id->last == segafter ) - id->last = seg; - } - while ( seg->prev && ( seg->starttime < seg->prev->starttime || - (seg->starttime == seg->prev->starttime && seg->endtime > seg->prev->endtime) ) ) - { - /* Move segment up list, swap seg and seg->prev */ - segbefore = seg->prev; - - if ( seg->next ) - seg->next->prev = segbefore; - - if ( segbefore->prev ) - segbefore->prev->next = seg; - - segbefore->next = seg->next; - seg->next = segbefore; - seg->prev = segbefore->prev; - segbefore->prev = seg; - - /* Reset first and last segment pointers if replaced */ - if ( id->first == segbefore ) - id->first = seg; - - if ( id->last == seg ) - id->last = segbefore; - } - - /* Set MSTraceID as last accessed */ - mstl->last = id; - - return seg; + /* Search complete segment list for matches */ + else + { + searchseg = id->first; + segbefore = 0; /* Find segment that record fits before */ + segafter = 0; /* Find segment that record fits after */ + followseg = 0; /* Track segment that record follows in time order */ + while ( searchseg ) + { + if ( msr->starttime > searchseg->starttime ) + followseg = searchseg; + + whence = 0; + + postgap = msr->starttime - searchseg->endtime - hpdelta; + if ( ! segbefore && postgap <= hptimetol && postgap >= nhptimetol ) + whence = 1; + + pregap = searchseg->starttime - endtime - hpdelta; + if ( ! segafter && pregap <= hptimetol && pregap >= nhptimetol ) + whence = 2; + + if ( ! whence ) + { + searchseg = searchseg->next; + continue; + } + + if ( sampratetol == -1.0 ) + { + if ( ! MS_ISRATETOLERABLE (msr->samprate, searchseg->samprate) ) + { + searchseg = searchseg->next; + continue; + } + } + else + { + if ( ms_dabs (msr->samprate - searchseg->samprate) > sampratetol ) + { + searchseg = searchseg->next; + continue; + } + } + + if ( whence == 1 ) + segbefore = searchseg; + else + segafter = searchseg; + + /* Done searching if not autohealing */ + if ( ! autoheal ) + break; + + /* Done searching if both before and after segments are found */ + if ( segbefore && segafter ) + break; + + searchseg = searchseg->next; + } /* Done looping through segments */ + + /* Add MSRecord coverage to end of segment before */ + if ( segbefore ) + { + if ( ! mstl_addmsrtoseg (segbefore, msr, endtime, 1) ) + { + return 0; + } + + /* Merge two segments that now fit if autohealing */ + if ( autoheal && segafter && segbefore != segafter ) + { + /* Add segafter coverage to segbefore */ + if ( ! mstl_addsegtoseg (segbefore, segafter) ) + { + return 0; + } + + /* Shift last segment pointer if it's going to be removed */ + if ( segafter == id->last ) + id->last = id->last->prev; + + /* Remove segafter from list */ + if ( segafter->prev ) + segafter->prev->next = segafter->next; + if ( segafter->next ) + segafter->next->prev = segafter->prev; + + /* Free data samples, private data and segment structure */ + if (segafter->datasamples) + free (segafter->datasamples); + + if (segafter->prvtptr) + free (segafter->prvtptr); + + free (segafter); + } + + seg = segbefore; + } + /* Add MSRecord coverage to beginning of segment after */ + else if ( segafter ) + { + if ( ! mstl_addmsrtoseg (segafter, msr, endtime, 2) ) + { + return 0; + } + + seg = segafter; + } + /* Add MSRecord coverage to new segment */ + else + { + /* Create new segment */ + if ( ! (seg = mstl_msr2seg (msr, endtime)) ) + { + return 0; + } + + /* Add new segment as first in list */ + if ( ! followseg ) + { + seg->next = id->first; + if ( id->first ) + id->first->prev = seg; + + id->first = seg; + } + /* Add new segment after the followseg segment */ + else + { + seg->next = followseg->next; + seg->prev = followseg; + if ( followseg->next ) + followseg->next->prev = seg; + followseg->next = seg; + + if ( followseg == id->last ) + id->last = seg; + } + + id->numsegments++; + } + + /* Track earliest and latest times */ + if ( msr->starttime < id->earliest ) + id->earliest = msr->starttime; + + if ( endtime > id->latest ) + id->latest = endtime; + } /* End of searching segment list */ + } /* End of adding coverage to matching ID */ + + /* Sort modified segment into place, logic above should limit these to few shifts if any */ + while ( seg->next && ( seg->starttime > seg->next->starttime || + (seg->starttime == seg->next->starttime && seg->endtime < seg->next->endtime) ) ) + { + /* Move segment down list, swap seg and seg->next */ + segafter = seg->next; + + if ( seg->prev ) + seg->prev->next = segafter; + + if ( segafter->next ) + segafter->next->prev = seg; + + segafter->prev = seg->prev; + seg->prev = segafter; + seg->next = segafter->next; + segafter->next = seg; + + /* Reset first and last segment pointers if replaced */ + if ( id->first == seg ) + id->first = segafter; + + if ( id->last == segafter ) + id->last = seg; + } + while ( seg->prev && ( seg->starttime < seg->prev->starttime || + (seg->starttime == seg->prev->starttime && seg->endtime > seg->prev->endtime) ) ) + { + /* Move segment up list, swap seg and seg->prev */ + segbefore = seg->prev; + + if ( seg->next ) + seg->next->prev = segbefore; + + if ( segbefore->prev ) + segbefore->prev->next = seg; + + segbefore->next = seg->next; + seg->next = segbefore; + seg->prev = segbefore->prev; + segbefore->prev = seg; + + /* Reset first and last segment pointers if replaced */ + if ( id->first == segbefore ) + id->first = seg; + + if ( id->last == seg ) + id->last = segbefore; + } + + /* Set MSTraceID as last accessed */ + mstl->last = id; + + return seg; } /* End of mstl_addmsr() */ @@ -613,39 +622,45 @@ mstl_addmsr ( MSTraceList *mstl, MSRecord *msr, flag dataquality, MSTraceSeg * mstl_msr2seg (MSRecord *msr, hptime_t endtime) { - MSTraceSeg *seg = 0; - int samplesize; - - if ( ! (seg = (MSTraceSeg *) calloc (1, sizeof(MSTraceSeg))) ) - { - ms_log (2, "mstl_addmsr(): Error allocating memory\n"); - return 0; - } - - /* Populate MSTraceSeg */ - seg->starttime = msr->starttime; - seg->endtime = endtime; - seg->samprate = msr->samprate; - seg->samplecnt = msr->samplecnt; - seg->sampletype = msr->sampletype; - seg->numsamples = msr->numsamples; - - /* Allocate space for and copy datasamples */ - if ( msr->datasamples && msr->numsamples ) - { - samplesize = ms_samplesize (msr->sampletype); - - if ( ! (seg->datasamples = malloc ((size_t) (samplesize * msr->numsamples))) ) + MSTraceSeg *seg = 0; + int samplesize; + + if ( ! (seg = (MSTraceSeg *) calloc (1, sizeof(MSTraceSeg))) ) + { + ms_log (2, "mstl_addmsr(): Error allocating memory\n"); + return 0; + } + + /* Populate MSTraceSeg */ + seg->starttime = msr->starttime; + seg->endtime = endtime; + seg->samprate = msr->samprate; + seg->samplecnt = msr->samplecnt; + seg->sampletype = msr->sampletype; + seg->numsamples = msr->numsamples; + + /* Allocate space for and copy datasamples */ + if ( msr->datasamples && msr->numsamples ) { - ms_log (2, "mstl_msr2seg(): Error allocating memory\n"); - return 0; + samplesize = ms_samplesize (msr->sampletype); + + if ( ! (seg->datasamples = malloc ((size_t) (samplesize * msr->numsamples))) ) + { + ms_log (2, "mstl_msr2seg(): Error allocating memory\n"); + return 0; + } + + /* Copy data samples from MSRecord to MSTraceSeg */ + memcpy (seg->datasamples, msr->datasamples, (size_t) (samplesize * msr->numsamples)); } - - /* Copy data samples from MSRecord to MSTraceSeg */ - memcpy (seg->datasamples, msr->datasamples, (size_t) (samplesize * msr->numsamples)); - } - - return seg; + + /* REC -- adding on a samprate tracking structure so we can get the + * mode (or most common) sample rate to represent as the definitive + * segment sample rate + */ + seg->samprate_list = add_segsamprate(seg->samprate_list,seg->samprate); + + return seg; } /* End of mstl_msr2seg() */ @@ -664,78 +679,85 @@ mstl_msr2seg (MSRecord *msr, hptime_t endtime) MSTraceSeg * mstl_addmsrtoseg (MSTraceSeg *seg, MSRecord *msr, hptime_t endtime, flag whence) { - int samplesize = 0; - void *newdatasamples; - - if ( ! seg || ! msr ) - return 0; - - /* Allocate more memory for data samples if included */ - if ( msr->datasamples && msr->numsamples > 0 ) - { - if ( msr->sampletype != seg->sampletype ) - { - ms_log (2, "mstl_addmsrtoseg(): MSRecord sample type (%c) does not match segment sample type (%c)\n", - msr->sampletype, seg->sampletype); - return 0; - } - - if ( ! (samplesize = ms_samplesize (msr->sampletype)) ) + int samplesize = 0; + void *newdatasamples; + + if ( ! seg || ! msr ) + return 0; + + /* Allocate more memory for data samples if included */ + if ( msr->datasamples && msr->numsamples > 0 ) { - ms_log (2, "mstl_addmsrtoseg(): Unknown sample size for sample type: %c\n", msr->sampletype); - return 0; + if ( msr->sampletype != seg->sampletype ) + { + ms_log (2, "mstl_addmsrtoseg(): MSRecord sample type (%c) does not match segment sample type (%c)\n", + msr->sampletype, seg->sampletype); + return 0; + } + + if ( ! (samplesize = ms_samplesize (msr->sampletype)) ) + { + ms_log (2, "mstl_addmsrtoseg(): Unknown sample size for sample type: %c\n", msr->sampletype); + return 0; + } + + if ( ! (newdatasamples = realloc (seg->datasamples, (size_t) ((seg->numsamples + msr->numsamples) * samplesize))) ) + { + ms_log (2, "mstl_addmsrtoseg(): Error allocating memory\n"); + return 0; + } + + seg->datasamples = newdatasamples; } - - if ( ! (newdatasamples = realloc (seg->datasamples, (size_t) ((seg->numsamples + msr->numsamples) * samplesize))) ) + + /* Add coverage to end of segment */ + if ( whence == 1 ) { - ms_log (2, "mstl_addmsrtoseg(): Error allocating memory\n"); - return 0; + seg->endtime = endtime; + seg->samplecnt += msr->samplecnt; + + if ( msr->datasamples && msr->numsamples > 0 ) + { + memcpy ((char *)seg->datasamples + (seg->numsamples * samplesize), + msr->datasamples, + (size_t) (msr->numsamples * samplesize)); + + seg->numsamples += msr->numsamples; + } } - - seg->datasamples = newdatasamples; - } - - /* Add coverage to end of segment */ - if ( whence == 1 ) - { - seg->endtime = endtime; - seg->samplecnt += msr->samplecnt; - - if ( msr->datasamples && msr->numsamples > 0 ) + /* Add coverage to beginning of segment */ + else if ( whence == 2 ) { - memcpy ((char *)seg->datasamples + (seg->numsamples * samplesize), - msr->datasamples, - (size_t) (msr->numsamples * samplesize)); - - seg->numsamples += msr->numsamples; + seg->starttime = msr->starttime; + seg->samplecnt += msr->samplecnt; + + if ( msr->datasamples && msr->numsamples > 0 ) + { + memmove ((char *)seg->datasamples + (msr->numsamples * samplesize), + seg->datasamples, + (size_t) (seg->numsamples * samplesize)); + + memcpy (seg->datasamples, + msr->datasamples, + (size_t) (msr->numsamples * samplesize)); + + seg->numsamples += msr->numsamples; + } } - } - /* Add coverage to beginning of segment */ - else if ( whence == 2 ) - { - seg->starttime = msr->starttime; - seg->samplecnt += msr->samplecnt; - - if ( msr->datasamples && msr->numsamples > 0 ) + else { - memmove ((char *)seg->datasamples + (msr->numsamples * samplesize), - seg->datasamples, - (size_t) (seg->numsamples * samplesize)); - - memcpy (seg->datasamples, - msr->datasamples, - (size_t) (msr->numsamples * samplesize)); - - seg->numsamples += msr->numsamples; + ms_log (2, "mstl_addmsrtoseg(): unrecognized whence value: %d\n", whence); + return 0; } - } - else - { - ms_log (2, "mstl_addmsrtoseg(): unrecognized whence value: %d\n", whence); - return 0; - } - - return seg; + + /* REC -- adding to a samprate tracking structure so we can get the + * mode (or most common) sample rate to represent as the definitive + * segment sample rate. Use the msr sample rate as seg only recorded + * the first initial sample rate it encountered ( see mstl_msr2seg() ) + */ + seg->samprate_list = add_segsamprate(seg->samprate_list,msr->samprate); + + return seg; } /* End of mstl_addmsrtoseg() */ @@ -749,51 +771,51 @@ mstl_addmsrtoseg (MSTraceSeg *seg, MSRecord *msr, hptime_t endtime, flag whence) MSTraceSeg * mstl_addsegtoseg (MSTraceSeg *seg1, MSTraceSeg *seg2) { - int samplesize = 0; - void *newdatasamples; - - if ( ! seg1 || ! seg2 ) - return 0; - - /* Allocate more memory for data samples if included */ - if ( seg2->datasamples && seg2->numsamples > 0 ) - { - if ( seg2->sampletype != seg1->sampletype ) - { - ms_log (2, "mstl_addsegtoseg(): MSTraceSeg sample types do not match (%c and %c)\n", - seg1->sampletype, seg2->sampletype); - return 0; - } - - if ( ! (samplesize = ms_samplesize (seg1->sampletype)) ) + int samplesize = 0; + void *newdatasamples; + + if ( ! seg1 || ! seg2 ) + return 0; + + /* Allocate more memory for data samples if included */ + if ( seg2->datasamples && seg2->numsamples > 0 ) { - ms_log (2, "mstl_addsegtoseg(): Unknown sample size for sample type: %c\n", seg1->sampletype); - return 0; + if ( seg2->sampletype != seg1->sampletype ) + { + ms_log (2, "mstl_addsegtoseg(): MSTraceSeg sample types do not match (%c and %c)\n", + seg1->sampletype, seg2->sampletype); + return 0; + } + + if ( ! (samplesize = ms_samplesize (seg1->sampletype)) ) + { + ms_log (2, "mstl_addsegtoseg(): Unknown sample size for sample type: %c\n", seg1->sampletype); + return 0; + } + + if ( ! (newdatasamples = realloc (seg1->datasamples, (size_t) ((seg1->numsamples + seg2->numsamples) * samplesize))) ) + { + ms_log (2, "mstl_addsegtoseg(): Error allocating memory\n"); + return 0; + } + + seg1->datasamples = newdatasamples; } - - if ( ! (newdatasamples = realloc (seg1->datasamples, (size_t) ((seg1->numsamples + seg2->numsamples) * samplesize))) ) + + /* Add seg2 coverage to end of seg1 */ + seg1->endtime = seg2->endtime; + seg1->samplecnt += seg2->samplecnt; + + if ( seg2->datasamples && seg2->numsamples > 0 ) { - ms_log (2, "mstl_addsegtoseg(): Error allocating memory\n"); - return 0; + memcpy ((char *)seg1->datasamples + (seg1->numsamples * samplesize), + seg2->datasamples, + (size_t) (seg2->numsamples * samplesize)); + + seg1->numsamples += seg2->numsamples; } - - seg1->datasamples = newdatasamples; - } - - /* Add seg2 coverage to end of seg1 */ - seg1->endtime = seg2->endtime; - seg1->samplecnt += seg2->samplecnt; - - if ( seg2->datasamples && seg2->numsamples > 0 ) - { - memcpy ((char *)seg1->datasamples + (seg1->numsamples * samplesize), - seg2->datasamples, - (size_t) (seg2->numsamples * samplesize)); - - seg1->numsamples += seg2->numsamples; - } - - return seg1; + + return seg1; } /* End of mstl_addsegtoseg() */ @@ -818,125 +840,125 @@ mstl_addsegtoseg (MSTraceSeg *seg1, MSTraceSeg *seg2) int mstl_convertsamples ( MSTraceSeg *seg, char type, flag truncate ) { - int32_t *idata; - float *fdata; - double *ddata; - int64_t idx; - - if ( ! seg ) - return -1; - - /* No conversion necessary, report success */ - if ( seg->sampletype == type ) - return 0; - - if ( seg->sampletype == 'a' || type == 'a' ) - { - ms_log (2, "mstl_convertsamples: cannot convert ASCII samples to/from numeric type\n"); - return -1; - } - - idata = (int32_t *) seg->datasamples; - fdata = (float *) seg->datasamples; - ddata = (double *) seg->datasamples; - - /* Convert to 32-bit integers */ - if ( type == 'i' ) - { - if ( seg->sampletype == 'f' ) /* Convert floats to integers with simple rounding */ + int32_t *idata; + float *fdata; + double *ddata; + int64_t idx; + + if ( ! seg ) + return -1; + + /* No conversion necessary, report success */ + if ( seg->sampletype == type ) + return 0; + + if ( seg->sampletype == 'a' || type == 'a' ) { - for (idx = 0; idx < seg->numsamples; idx++) - { - /* Check for loss of sub-integer */ - if ( ! truncate && (fdata[idx] - (int32_t)fdata[idx]) > 0.000001 ) - { - ms_log (1, "mstl_convertsamples: Warning, loss of precision when converting floats to integers, loss: %g\n", - (fdata[idx] - (int32_t)fdata[idx])); - return -1; - } - - idata[idx] = (int32_t) (fdata[idx] + 0.5); - } + ms_log (2, "mstl_convertsamples: cannot convert ASCII samples to/from numeric type\n"); + return -1; } - else if ( seg->sampletype == 'd' ) /* Convert doubles to integers with simple rounding */ + + idata = (int32_t *) seg->datasamples; + fdata = (float *) seg->datasamples; + ddata = (double *) seg->datasamples; + + /* Convert to 32-bit integers */ + if ( type == 'i' ) { - for (idx = 0; idx < seg->numsamples; idx++) - { - /* Check for loss of sub-integer */ - if ( ! truncate && (ddata[idx] - (int32_t)ddata[idx]) > 0.000001 ) + if ( seg->sampletype == 'f' ) /* Convert floats to integers with simple rounding */ { - ms_log (1, "mstl_convertsamples: Warning, loss of precision when converting doubles to integers, loss: %g\n", - (ddata[idx] - (int32_t)ddata[idx])); - return -1; + for (idx = 0; idx < seg->numsamples; idx++) + { + /* Check for loss of sub-integer */ + if ( ! truncate && (fdata[idx] - (int32_t)fdata[idx]) > 0.000001 ) + { + ms_log (1, "mstl_convertsamples: Warning, loss of precision when converting floats to integers, loss: %g\n", + (fdata[idx] - (int32_t)fdata[idx])); + return -1; + } + + idata[idx] = (int32_t) (fdata[idx] + 0.5); + } } - - idata[idx] = (int32_t) (ddata[idx] + 0.5); - } - - /* Reallocate buffer for reduced size needed */ - if ( ! (seg->datasamples = realloc (seg->datasamples, (size_t) (seg->numsamples * sizeof(int32_t)))) ) - { - ms_log (2, "mstl_convertsamples: cannot re-allocate buffer for sample conversion\n"); - return -1; - } - } - - seg->sampletype = 'i'; - } /* Done converting to 32-bit integers */ - - /* Convert to 32-bit floats */ - else if ( type == 'f' ) - { - if ( seg->sampletype == 'i' ) /* Convert integers to floats */ - { - for (idx = 0; idx < seg->numsamples; idx++) - fdata[idx] = (float) idata[idx]; - } - else if ( seg->sampletype == 'd' ) /* Convert doubles to floats */ - { - for (idx = 0; idx < seg->numsamples; idx++) - fdata[idx] = (float) ddata[idx]; - - /* Reallocate buffer for reduced size needed */ - if ( ! (seg->datasamples = realloc (seg->datasamples, (size_t) (seg->numsamples * sizeof(float)))) ) - { - ms_log (2, "mstl_convertsamples: cannot re-allocate buffer after sample conversion\n"); - return -1; - } - } - - seg->sampletype = 'f'; - } /* Done converting to 32-bit floats */ - - /* Convert to 64-bit doubles */ - else if ( type == 'd' ) - { - if ( ! (ddata = (double *) malloc ((size_t) (seg->numsamples * sizeof(double)))) ) - { - ms_log (2, "mstl_convertsamples: cannot allocate buffer for sample conversion to doubles\n"); - return -1; - } - - if ( seg->sampletype == 'i' ) /* Convert integers to doubles */ + else if ( seg->sampletype == 'd' ) /* Convert doubles to integers with simple rounding */ + { + for (idx = 0; idx < seg->numsamples; idx++) + { + /* Check for loss of sub-integer */ + if ( ! truncate && (ddata[idx] - (int32_t)ddata[idx]) > 0.000001 ) + { + ms_log (1, "mstl_convertsamples: Warning, loss of precision when converting doubles to integers, loss: %g\n", + (ddata[idx] - (int32_t)ddata[idx])); + return -1; + } + + idata[idx] = (int32_t) (ddata[idx] + 0.5); + } + + /* Reallocate buffer for reduced size needed */ + if ( ! (seg->datasamples = realloc (seg->datasamples, (size_t) (seg->numsamples * sizeof(int32_t)))) ) + { + ms_log (2, "mstl_convertsamples: cannot re-allocate buffer for sample conversion\n"); + return -1; + } + } + + seg->sampletype = 'i'; + } /* Done converting to 32-bit integers */ + + /* Convert to 32-bit floats */ + else if ( type == 'f' ) { - for (idx = 0; idx < seg->numsamples; idx++) - ddata[idx] = (double) idata[idx]; - - free (idata); - } - else if ( seg->sampletype == 'f' ) /* Convert floats to doubles */ + if ( seg->sampletype == 'i' ) /* Convert integers to floats */ + { + for (idx = 0; idx < seg->numsamples; idx++) + fdata[idx] = (float) idata[idx]; + } + else if ( seg->sampletype == 'd' ) /* Convert doubles to floats */ + { + for (idx = 0; idx < seg->numsamples; idx++) + fdata[idx] = (float) ddata[idx]; + + /* Reallocate buffer for reduced size needed */ + if ( ! (seg->datasamples = realloc (seg->datasamples, (size_t) (seg->numsamples * sizeof(float)))) ) + { + ms_log (2, "mstl_convertsamples: cannot re-allocate buffer after sample conversion\n"); + return -1; + } + } + + seg->sampletype = 'f'; + } /* Done converting to 32-bit floats */ + + /* Convert to 64-bit doubles */ + else if ( type == 'd' ) { - for (idx = 0; idx < seg->numsamples; idx++) - ddata[idx] = (double) fdata[idx]; - - free (fdata); - } - - seg->datasamples = ddata; - seg->sampletype = 'd'; - } /* Done converting to 64-bit doubles */ - - return 0; + if ( ! (ddata = (double *) malloc ((size_t) (seg->numsamples * sizeof(double)))) ) + { + ms_log (2, "mstl_convertsamples: cannot allocate buffer for sample conversion to doubles\n"); + return -1; + } + + if ( seg->sampletype == 'i' ) /* Convert integers to doubles */ + { + for (idx = 0; idx < seg->numsamples; idx++) + ddata[idx] = (double) idata[idx]; + + free (idata); + } + else if ( seg->sampletype == 'f' ) /* Convert floats to doubles */ + { + for (idx = 0; idx < seg->numsamples; idx++) + ddata[idx] = (double) fdata[idx]; + + free (fdata); + } + + seg->datasamples = ddata; + seg->sampletype = 'd'; + } /* Done converting to 64-bit doubles */ + + return 0; } /* End of mstl_convertsamples() */ @@ -959,125 +981,128 @@ mstl_convertsamples ( MSTraceSeg *seg, char type, flag truncate ) ***************************************************************************/ void mstl_printtracelist ( MSTraceList *mstl, flag timeformat, - flag details, flag gaps ) + flag details, flag gaps ) { - MSTraceID *id = 0; - MSTraceSeg *seg = 0; - char stime[30]; - char etime[30]; - char gapstr[20]; - flag nogap; - double gap; - double delta; - int tracecnt = 0; - int segcnt = 0; - - if ( ! mstl ) - { - return; - } - - /* Print out the appropriate header */ - if ( details > 0 && gaps > 0 ) - ms_log (0, " Source Start sample End sample Gap Hz Samples\n"); - else if ( details <= 0 && gaps > 0 ) - ms_log (0, " Source Start sample End sample Gap\n"); - else if ( details > 0 && gaps <= 0 ) - ms_log (0, " Source Start sample End sample Hz Samples\n"); - else - ms_log (0, " Source Start sample End sample\n"); - - /* Loop through trace list */ - id = mstl->traces; - while ( id ) - { - /* Loop through segment list */ - seg = id->first; - while ( seg ) + MSTraceID *id = 0; + MSTraceSeg *seg = 0; + char stime[30]; + char etime[30]; + char gapstr[20]; + flag nogap; + double gap; + double delta; + int tracecnt = 0; + int segcnt = 0; + + if ( ! mstl ) + { + return; + } + + /* Print out the appropriate header */ + if ( details > 0 && gaps > 0 ) + ms_log (0, " Source Start sample End sample Gap Hz Samples\n"); + else if ( details <= 0 && gaps > 0 ) + ms_log (0, " Source Start sample End sample Gap\n"); + else if ( details > 0 && gaps <= 0 ) + ms_log (0, " Source Start sample End sample Hz Samples\n"); + else + ms_log (0, " Source Start sample End sample\n"); + + /* Loop through trace list */ + id = mstl->traces; + while ( id ) { - /* Create formatted time strings */ - if ( timeformat == 2 ) - { - snprintf (stime, sizeof(stime), "%.6f", (double) MS_HPTIME2EPOCH(seg->starttime) ); - snprintf (etime, sizeof(etime), "%.6f", (double) MS_HPTIME2EPOCH(seg->endtime) ); - } - else if ( timeformat == 1 ) - { - if ( ms_hptime2isotimestr (seg->starttime, stime, 1) == NULL ) - ms_log (2, "Cannot convert trace start time for %s\n", id->srcname); - - if ( ms_hptime2isotimestr (seg->endtime, etime, 1) == NULL ) - ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); - } - else - { - if ( ms_hptime2seedtimestr (seg->starttime, stime, 1) == NULL ) - ms_log (2, "Cannot convert trace start time for %s\n", id->srcname); - - if ( ms_hptime2seedtimestr (seg->endtime, etime, 1) == NULL ) - ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); - } - - /* Print segment info at varying levels */ - if ( gaps > 0 ) - { - gap = 0.0; - nogap = 0; - - if ( seg->prev ) - gap = (double) (seg->starttime - seg->prev->endtime) / HPTMODULUS; - else - nogap = 1; - - /* Check that any overlap is not larger than the trace coverage */ - if ( gap < 0.0 ) + /* Loop through segment list */ + seg = id->first; + while ( seg ) { - delta = ( seg->samprate ) ? (1.0 / seg->samprate) : 0.0; - - if ( (gap * -1.0) > (((double)(seg->endtime - seg->starttime)/HPTMODULUS) + delta) ) - gap = -(((double)(seg->endtime - seg->starttime)/HPTMODULUS) + delta); + /* REC -- use the mode of segment sample rate */ + seg->samprate = get_segsamprate_mode(seg->samprate_list); + + /* Create formatted time strings */ + if ( timeformat == 2 ) + { + snprintf (stime, sizeof(stime), "%.6f", (double) MS_HPTIME2EPOCH(seg->starttime) ); + snprintf (etime, sizeof(etime), "%.6f", (double) MS_HPTIME2EPOCH(seg->endtime) ); + } + else if ( timeformat == 1 ) + { + if ( ms_hptime2isotimestr (seg->starttime, stime, 1) == NULL ) + ms_log (2, "Cannot convert trace start time for %s\n", id->srcname); + + if ( ms_hptime2isotimestr (seg->endtime, etime, 1) == NULL ) + ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); + } + else + { + if ( ms_hptime2seedtimestr (seg->starttime, stime, 1) == NULL ) + ms_log (2, "Cannot convert trace start time for %s\n", id->srcname); + + if ( ms_hptime2seedtimestr (seg->endtime, etime, 1) == NULL ) + ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); + } + + /* Print segment info at varying levels */ + if ( gaps > 0 ) + { + gap = 0.0; + nogap = 0; + + if ( seg->prev ) + gap = (double) (seg->starttime - seg->prev->endtime) / HPTMODULUS; + else + nogap = 1; + + /* Check that any overlap is not larger than the trace coverage */ + if ( gap < 0.0 ) + { + delta = ( seg->samprate ) ? (1.0 / seg->samprate) : 0.0; + + if ( (gap * -1.0) > (((double)(seg->endtime - seg->starttime)/HPTMODULUS) + delta) ) + gap = -(((double)(seg->endtime - seg->starttime)/HPTMODULUS) + delta); + } + + /* Fix up gap display */ + if ( nogap ) + snprintf (gapstr, sizeof(gapstr), " == "); + else if ( gap >= 86400.0 || gap <= -86400.0 ) + snprintf (gapstr, sizeof(gapstr), "%-3.1fd", (gap / 86400)); + else if ( gap >= 3600.0 || gap <= -3600.0 ) + snprintf (gapstr, sizeof(gapstr), "%-3.1fh", (gap / 3600)); + else if ( gap == 0.0 ) + snprintf (gapstr, sizeof(gapstr), "-0 "); + else + snprintf (gapstr, sizeof(gapstr), "%-4.4g", gap); + + if ( details <= 0 ) + ms_log (0, "%-17s %-24s %-24s %-4s\n", + id->srcname, stime, etime, gapstr); + else + ms_log (0, "%-17s %-24s %-24s %-s %-3.3g %-"PRId64"\n", + id->srcname, stime, etime, gapstr, seg->samprate, seg->samplecnt); + } + else if ( details > 0 && gaps <= 0 ) + ms_log (0, "%-17s %-24s %-24s %-3.3g %-"PRId64"\n", + id->srcname, stime, etime, seg->samprate, seg->samplecnt); + else + ms_log (0, "%-17s %-24s %-24s\n", id->srcname, stime, etime); + + segcnt++; + seg = seg->next; } - - /* Fix up gap display */ - if ( nogap ) - snprintf (gapstr, sizeof(gapstr), " == "); - else if ( gap >= 86400.0 || gap <= -86400.0 ) - snprintf (gapstr, sizeof(gapstr), "%-3.1fd", (gap / 86400)); - else if ( gap >= 3600.0 || gap <= -3600.0 ) - snprintf (gapstr, sizeof(gapstr), "%-3.1fh", (gap / 3600)); - else if ( gap == 0.0 ) - snprintf (gapstr, sizeof(gapstr), "-0 "); - else - snprintf (gapstr, sizeof(gapstr), "%-4.4g", gap); - - if ( details <= 0 ) - ms_log (0, "%-17s %-24s %-24s %-4s\n", - id->srcname, stime, etime, gapstr); - else - ms_log (0, "%-17s %-24s %-24s %-s %-3.3g %-"PRId64"\n", - id->srcname, stime, etime, gapstr, seg->samprate, seg->samplecnt); - } - else if ( details > 0 && gaps <= 0 ) - ms_log (0, "%-17s %-24s %-24s %-3.3g %-"PRId64"\n", - id->srcname, stime, etime, seg->samprate, seg->samplecnt); - else - ms_log (0, "%-17s %-24s %-24s\n", id->srcname, stime, etime); - - segcnt++; - seg = seg->next; + + tracecnt++; + id = id->next; } - - tracecnt++; - id = id->next; - } - - if ( tracecnt != mstl->numtraces ) - ms_log (2, "mstl_printtracelist(): number of traces in trace list is inconsistent\n"); - - if ( details > 0 ) - ms_log (0, "Total: %d trace(s) with %d segment(s)\n", tracecnt, segcnt); - - return; + + if ( tracecnt != mstl->numtraces ) + ms_log (2, "mstl_printtracelist(): number of traces in trace list is inconsistent\n"); + + if ( details > 0 ) + ms_log (0, "Total: %d trace(s) with %d segment(s)\n", tracecnt, segcnt); + + return; } /* End of mstl_printtracelist() */ @@ -1097,51 +1122,54 @@ mstl_printtracelist ( MSTraceList *mstl, flag timeformat, void mstl_printsynclist ( MSTraceList *mstl, char *dccid, flag subsecond ) { - MSTraceID *id = 0; - MSTraceSeg *seg = 0; - char starttime[30]; - char endtime[30]; - char yearday[10]; - time_t now; - struct tm *nt; - - if ( ! mstl ) - { - return; - } - - /* Generate current time stamp */ - now = time (NULL); - nt = localtime ( &now ); nt->tm_year += 1900; nt->tm_yday += 1; - snprintf ( yearday, sizeof(yearday), "%04d,%03d", nt->tm_year, nt->tm_yday); - - /* Print SYNC header line */ - ms_log (0, "%s|%s\n", (dccid)?dccid:"DCC", yearday); - - /* Loop through trace list */ - id = mstl->traces; - while ( id ) - { - /* Loop through segment list */ - seg = id->first; - while ( seg ) + MSTraceID *id = 0; + MSTraceSeg *seg = 0; + char starttime[30]; + char endtime[30]; + char yearday[10]; + time_t now; + struct tm *nt; + + if ( ! mstl ) + { + return; + } + + /* Generate current time stamp */ + now = time (NULL); + nt = localtime ( &now ); nt->tm_year += 1900; nt->tm_yday += 1; + snprintf ( yearday, sizeof(yearday), "%04d,%03d", nt->tm_year, nt->tm_yday); + + /* Print SYNC header line */ + ms_log (0, "%s|%s\n", (dccid)?dccid:"DCC", yearday); + + /* Loop through trace list */ + id = mstl->traces; + while ( id ) { - ms_hptime2seedtimestr (seg->starttime, starttime, subsecond); - ms_hptime2seedtimestr (seg->endtime, endtime, subsecond); - - /* Print SYNC line */ - ms_log (0, "%s|%s|%s|%s|%s|%s||%.10g|%"PRId64"|||||||%s\n", - id->network, id->station, id->location, id->channel, - starttime, endtime, seg->samprate, seg->samplecnt, - yearday); - - seg = seg->next; + /* Loop through segment list */ + seg = id->first; + while ( seg ) + { + ms_hptime2seedtimestr (seg->starttime, starttime, subsecond); + ms_hptime2seedtimestr (seg->endtime, endtime, subsecond); + + /* REC -- print mode of segment sample rate */ + seg->samprate = get_segsamprate_mode(seg->samprate_list); + + /* Print SYNC line */ + ms_log (0, "%s|%s|%s|%s|%s|%s||%.10g|%"PRId64"|||||||%s\n", + id->network, id->station, id->location, id->channel, + starttime, endtime, seg->samprate, seg->samplecnt, + yearday); + + seg = seg->next; + } + + id = id->next; } - - id = id->next; - } - - return; + + return; } /* End of mstl_printsynclist() */ @@ -1161,117 +1189,174 @@ mstl_printsynclist ( MSTraceList *mstl, char *dccid, flag subsecond ) ***************************************************************************/ void mstl_printgaplist (MSTraceList *mstl, flag timeformat, - double *mingap, double *maxgap) + double *mingap, double *maxgap) { - MSTraceID *id = 0; - MSTraceSeg *seg = 0; - - char time1[30], time2[30]; - char gapstr[30]; - double gap; - double delta; - double nsamples; - flag printflag; - int gapcnt = 0; - - if ( ! mstl ) - return; - - if ( ! mstl->traces ) - return; - - ms_log (0, " Source Last Sample Next Sample Gap Samples\n"); - - id = mstl->traces; - while ( id ) - { - seg = id->first; - while ( seg->next ) + MSTraceID *id = 0; + MSTraceSeg *seg = 0; + + char time1[30], time2[30]; + char gapstr[30]; + double gap; + double delta; + double nsamples; + flag printflag; + int gapcnt = 0; + + if ( ! mstl ) + return; + + if ( ! mstl->traces ) + return; + + ms_log (0, " Source Last Sample Next Sample Gap Samples\n"); + + id = mstl->traces; + while ( id ) { - /* Skip segments with 0 sample rate, usually from SOH records */ - if ( seg->samprate == 0.0 ) - { - seg = seg->next; - continue; - } - - gap = (double) (seg->next->starttime - seg->endtime) / HPTMODULUS; - - /* Check that any overlap is not larger than the trace coverage */ - if ( gap < 0.0 ) - { - delta = ( seg->next->samprate ) ? (1.0 / seg->next->samprate) : 0.0; - - if ( (gap * -1.0) > (((double)(seg->next->endtime - seg->next->starttime)/HPTMODULUS) + delta) ) - gap = -(((double)(seg->next->endtime - seg->next->starttime)/HPTMODULUS) + delta); - } - - printflag = 1; - - /* Check gap/overlap criteria */ - if ( mingap ) - if ( gap < *mingap ) - printflag = 0; - - if ( maxgap ) - if ( gap > *maxgap ) - printflag = 0; - - if ( printflag ) - { - nsamples = ms_dabs(gap) * seg->samprate; - - if ( gap > 0.0 ) - nsamples -= 1.0; - else - nsamples += 1.0; - - /* Fix up gap display */ - if ( gap >= 86400.0 || gap <= -86400.0 ) - snprintf (gapstr, sizeof(gapstr), "%-3.1fd", (gap / 86400)); - else if ( gap >= 3600.0 || gap <= -3600.0 ) - snprintf (gapstr, sizeof(gapstr), "%-3.1fh", (gap / 3600)); - else if ( gap == 0.0 ) - snprintf (gapstr, sizeof(gapstr), "-0 "); - else - snprintf (gapstr, sizeof(gapstr), "%-4.4g", gap); - - /* Create formatted time strings */ - if ( timeformat == 2 ) - { - snprintf (time1, sizeof(time1), "%.6f", (double) MS_HPTIME2EPOCH(seg->endtime) ); - snprintf (time2, sizeof(time2), "%.6f", (double) MS_HPTIME2EPOCH(seg->next->starttime) ); - } - else if ( timeformat == 1 ) + seg = id->first; + while ( seg->next ) { - if ( ms_hptime2isotimestr (seg->endtime, time1, 1) == NULL ) - ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); - - if ( ms_hptime2isotimestr (seg->next->starttime, time2, 1) == NULL ) - ms_log (2, "Cannot convert next trace start time for %s\n", id->srcname); - } - else - { - if ( ms_hptime2seedtimestr (seg->endtime, time1, 1) == NULL ) - ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); - - if ( ms_hptime2seedtimestr (seg->next->starttime, time2, 1) == NULL ) - ms_log (2, "Cannot convert next trace start time for %s\n", id->srcname); + /* Skip segments with 0 sample rate, usually from SOH records */ + if ( seg->samprate == 0.0 ) + { + seg = seg->next; + continue; + } + + gap = (double) (seg->next->starttime - seg->endtime) / HPTMODULUS; + + /* Check that any overlap is not larger than the trace coverage */ + if ( gap < 0.0 ) + { + delta = ( seg->next->samprate ) ? (1.0 / seg->next->samprate) : 0.0; + + if ( (gap * -1.0) > (((double)(seg->next->endtime - seg->next->starttime)/HPTMODULUS) + delta) ) + gap = -(((double)(seg->next->endtime - seg->next->starttime)/HPTMODULUS) + delta); + } + + printflag = 1; + + /* Check gap/overlap criteria */ + if ( mingap ) + if ( gap < *mingap ) + printflag = 0; + + if ( maxgap ) + if ( gap > *maxgap ) + printflag = 0; + + if ( printflag ) + { + nsamples = ms_dabs(gap) * seg->samprate; + + if ( gap > 0.0 ) + nsamples -= 1.0; + else + nsamples += 1.0; + + /* Fix up gap display */ + if ( gap >= 86400.0 || gap <= -86400.0 ) + snprintf (gapstr, sizeof(gapstr), "%-3.1fd", (gap / 86400)); + else if ( gap >= 3600.0 || gap <= -3600.0 ) + snprintf (gapstr, sizeof(gapstr), "%-3.1fh", (gap / 3600)); + else if ( gap == 0.0 ) + snprintf (gapstr, sizeof(gapstr), "-0 "); + else + snprintf (gapstr, sizeof(gapstr), "%-4.4g", gap); + + /* Create formatted time strings */ + if ( timeformat == 2 ) + { + snprintf (time1, sizeof(time1), "%.6f", (double) MS_HPTIME2EPOCH(seg->endtime) ); + snprintf (time2, sizeof(time2), "%.6f", (double) MS_HPTIME2EPOCH(seg->next->starttime) ); + } + else if ( timeformat == 1 ) + { + if ( ms_hptime2isotimestr (seg->endtime, time1, 1) == NULL ) + ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); + + if ( ms_hptime2isotimestr (seg->next->starttime, time2, 1) == NULL ) + ms_log (2, "Cannot convert next trace start time for %s\n", id->srcname); + } + else + { + if ( ms_hptime2seedtimestr (seg->endtime, time1, 1) == NULL ) + ms_log (2, "Cannot convert trace end time for %s\n", id->srcname); + + if ( ms_hptime2seedtimestr (seg->next->starttime, time2, 1) == NULL ) + ms_log (2, "Cannot convert next trace start time for %s\n", id->srcname); + } + + ms_log (0, "%-17s %-24s %-24s %-4s %-.8g\n", + id->srcname, time1, time2, gapstr, nsamples); + + gapcnt++; + } + + seg = seg->next; } - - ms_log (0, "%-17s %-24s %-24s %-4s %-.8g\n", - id->srcname, time1, time2, gapstr, nsamples); - - gapcnt++; - } - - seg = seg->next; + + id = id->next; } - - id = id->next; - } - - ms_log (0, "Total: %d gap(s)\n", gapcnt); - - return; + + ms_log (0, "Total: %d gap(s)\n", gapcnt); + + return; } /* End of mstl_printgaplist() */ + + +/* REC - non-standard function addition. To the indicated pointer to an associative linked list + * add a sample rate value, which will result in either a new node with a count of 1, or an incremented + * count to an existing node with the same value. The result is a histogram of sample rates. + * + * Return - pointer to head of linked list. + */ +MSSampRateList* add_segsamprate(MSSampRateList *srl, double samprate) { + MSSampRateList *head, *current, *prev = NULL; + head = current = srl; + /*printf("DEBUG: add_segsamprate %f\n",samprate);*/ + while (current != NULL) { + /* find a node with our sample rate */ + /*printf("DEBUG: -- add_segsamprate compare: %f to %f\n",current->samprate,samprate);*/ + if (current->samprate == samprate) { + current->count++; + /*printf("DEBUG: ---- add_segsamprate +1: %f,%d\n",current->samprate,current->count);*/ + break; + } else { + prev = current; + current = current->next; + } + } + if (current == NULL) { + /* create new node with a count of 1 */ + current = (MSSampRateList*) malloc(sizeof(MSSampRateList)); + current->samprate = samprate; + current->count = 1; + /*printf("DEBUG: ---- add_segsamprate new: %f,%d\n",current->samprate,current->count);*/ + current->next = NULL; + if (prev != NULL) prev->next = current; /* add to linked list */ + else head = current; + } + return head; + +} /* End of add_segsamprate */ + +/* REC - non-standard function addition. To the indicated pointer to an associative linked list + * of sample rate counts, return the sample rate with the highest count value, or 0.0 if the + * parameter was null. + */ +double get_segsamprate_mode(MSSampRateList *srl) { + MSSampRateList *current,*max; + current = max = srl; + /*printf("DEBUG: get_segsamprate_mode\n");*/ + /* loop over the sample rate list and find the max count */ + while (current != NULL) { + /* find a node with our sample rate */ + if (max == NULL || current->count > max->count) + max = current; + /*printf("DEBUG: ---- get_segsamprate_mode max: %f\n",max->samprate);*/ + current = current->next; + } + if (max == NULL) return (double) 0.0; + else return (double) max->samprate; +} /* End of get_segsamprate_mode */ diff --git a/src/libmseed/unpackdata.c b/src/libmseed/unpackdata.c index b478950..b6fafdf 100644 --- a/src/libmseed/unpackdata.c +++ b/src/libmseed/unpackdata.c @@ -14,7 +14,7 @@ * (previously) ORFEUS/EC-Project MEREDIAN * (currently) IRIS Data Management Center * - * modified: 2012.357 + * modified: 2010.012 ************************************************************************/ /* @@ -144,7 +144,7 @@ int msr_unpack_float_32 if (req_samples < 0) return 0; for (nd=0; nd 0) *data = *px0; + /* Compute all but first values based on previous value */ prev = data - 1; while (--nr > 0 && --nd > 0) @@ -339,12 +340,14 @@ int msr_unpack_steim1 while (--nd > 0) last_data = *++diff + last_data; - /* Verify that the last value is identical to xn = rev. int. constant */ +/* Verify that the last value is identical to xn = rev. int. constant */ +/** REC -- remove integrity check for now -- we have annoying old data that trips this ** if (last_data != *pxn) { ms_log (1, "%s: Warning: Data integrity check for Steim-1 failed, last_data=%d, xn=%d\n", UNPACK_SRCNAME, last_data, *pxn); } +****/ return ((req_samples < num_samples) ? req_samples : num_samples); } /* End of msr_unpack_steim1() */ @@ -533,11 +536,13 @@ int msr_unpack_steim2 last_data = *++diff + last_data; /* Verify that the last value is identical to xn = rev. int. constant */ +/** REC -- remove integrity check for now -- we have annoying old data that trips this ** if (last_data != *pxn) { ms_log (1, "%s: Warning: Data integrity check for Steim-2 failed, last_data=%d, xn=%d\n", UNPACK_SRCNAME, last_data, *pxn); } +*************/ return ((req_samples < num_samples) ? req_samples : num_samples); } /* End of msr_unpack_steim2() */ diff --git a/src/parseMiniSEED.c b/src/parseMiniSEED.c index 51ebca0..c70f246 100644 --- a/src/parseMiniSEED.c +++ b/src/parseMiniSEED.c @@ -31,16 +31,19 @@ #include #include +typedef void (* R_callback)(char *); + SEXP parseMiniSEED (SEXP buffer) { int debug = 0; - int numPROTECT = 0; + int numPROTECT = 0; // number of protected elements outside the loop + int loopPROTECT = 0; // the number of protected elements inside the loop int bufferLength = 0; char *bufferPtr; char msgPrefix[11] = "libmseed__"; // Redirect libmseed logging messages to Matlab functions - ms_loginit( (void *)&Rprintf, msgPrefix, (void *)&Rf_error, msgPrefix); + ms_loginit( (R_callback)&Rprintf, msgPrefix, (R_callback)&Rf_error, msgPrefix); // Allocate space for the buffer PROTECT(buffer = AS_RAW(buffer)); @@ -162,49 +165,53 @@ SEXP parseMiniSEED (SEXP buffer) { PROTECT(returnList = NEW_LIST(id->numsegments)); numPROTECT++; + // Set up the names that go in the returnList + char *names[14] = {"npts","sampling_rate","network","station","location","channel","quality", + "starttime","endtime","data","act_flags","io_flags","dq_flags","timing_qual"}; + SEXP listNames; + PROTECT(listNames = NEW_CHARACTER(14)); + numPROTECT++; + for (int i=0; i<14; i++) { + SET_STRING_ELT(listNames, i, mkChar(names[i])); + } + // Loop through the segments (continuous blocks of data), populating a segmentList // which is then inserted into the returnList. for (int segIndex=0; segIndex < id->numsegments; segIndex++ ) { - // Set up the names that go in the returnList - char *names[14] = {"npts","sampling_rate","network","station","location","channel","quality", - "starttime","endtime","data","act_flags","io_flags","dq_flags","timing_qual"}; - SEXP listNames; - PROTECT(listNames = NEW_CHARACTER(14)); - numPROTECT++; - for (int i=0; i<14; i++) { - SET_STRING_ELT(listNames, i, mkChar(names[i])); - } - // Now set up the variable vectors that go into the segmentList SEXP npts, sampling_rate, network, station, location, channel, quality, starttime, endtime, data; + + loopPROTECT = 0; PROTECT(npts = NEW_INTEGER(1)); - numPROTECT++; + loopPROTECT++; INTEGER(npts)[0] = (int) seg->samplecnt; PROTECT(sampling_rate = NEW_NUMERIC(1)); - numPROTECT++; - REAL(sampling_rate)[0] = (double) seg->samprate; + loopPROTECT++; + // REC -- modified samprate for segment + // -- get the mode value of a set of sample rates recorded for this segment + REAL(sampling_rate)[0] = get_segsamprate_mode(seg->samprate_list); PROTECT(network = NEW_CHARACTER(1)); - numPROTECT++; + loopPROTECT++; SET_STRING_ELT(network, 0, mkChar(id->network)); PROTECT(station = NEW_CHARACTER(1)); - numPROTECT++; + loopPROTECT++; SET_STRING_ELT(station, 0, mkChar(id->station)); PROTECT(location = NEW_CHARACTER(1)); - numPROTECT++; + loopPROTECT++; SET_STRING_ELT(location, 0, mkChar(id->location)); PROTECT(channel = NEW_CHARACTER(1)); - numPROTECT++; + loopPROTECT++; SET_STRING_ELT(channel, 0, mkChar(id->channel)); PROTECT(quality = NEW_CHARACTER(1)); - numPROTECT++; + loopPROTECT++; // Extra step because R only deals with strings, not single characters. // Create the NULL terminated string by hand. char dataquality[2] = { id->dataquality, '\0' }; @@ -214,17 +221,17 @@ SEXP parseMiniSEED (SEXP buffer) { // NOTE: When these values are used in R, starttime and endtime need to be converted with: // s <- as.POSIXct(header$starttime, tz="GMT", origin=as.POSIXct("1970-01-01T00:00:00Z")) PROTECT(starttime = NEW_NUMERIC(1)); - numPROTECT++; + loopPROTECT++; REAL(starttime)[0] = (double) MS_HPTIME2EPOCH(seg->starttime); PROTECT(endtime = NEW_NUMERIC(1)); - numPROTECT++; + loopPROTECT++; REAL(endtime)[0] = (double) MS_HPTIME2EPOCH(seg->endtime); // TODO: Check for non-int data types // // //Rprintf("Data is of type '%c'\n",seg->sampletype); PROTECT(data = NEW_NUMERIC(seg->samplecnt)); - numPROTECT++; + loopPROTECT++; int32_t *idatasamplesPtr = (int32_t *) seg->datasamples; float *fdatasamplesPtr = (float *) seg->datasamples; @@ -248,7 +255,7 @@ SEXP parseMiniSEED (SEXP buffer) { SEXP act_flags, io_flags, dq_flags, timing_qual; /* PROTECT(act_flags = NEW_LOGICAL(8)); - numPROTECT++; + loopPROTECT++; b = msr->fsdh->act_flags; LOGICAL(act_flags)[0] = (int) bit(b,0x01); LOGICAL(act_flags)[1] = (int) bit(b,0x02); @@ -260,14 +267,14 @@ SEXP parseMiniSEED (SEXP buffer) { LOGICAL(act_flags)[7] = (int) bit(b,0x80); */ PROTECT(act_flags = NEW_INTEGER(8)); - numPROTECT++; + loopPROTECT++; for (int i=0; i<8; i++) { INTEGER(act_flags)[i] = (int) total_act_flags[i]; } /* PROTECT(io_flags = NEW_LOGICAL(8)); - numPROTECT++; + loopPROTECT++; b = msr->fsdh->io_flags; LOGICAL(io_flags)[0] = (int) bit(b,0x01); LOGICAL(io_flags)[1] = (int) bit(b,0x02); @@ -279,14 +286,14 @@ SEXP parseMiniSEED (SEXP buffer) { LOGICAL(io_flags)[7] = (int) bit(b,0x80); */ PROTECT(io_flags = NEW_INTEGER(8)); - numPROTECT++; + loopPROTECT++; for (int i=0; i<8; i++) { INTEGER(io_flags)[i] = (int) total_io_flags[i]; } /* PROTECT(dq_flags = NEW_LOGICAL(8)); - numPROTECT++; + loopPROTECT++; b = msr->fsdh->dq_flags; LOGICAL(dq_flags)[0] = (int) bit(b,0x01); LOGICAL(dq_flags)[1] = (int) bit(b,0x02); @@ -298,13 +305,13 @@ SEXP parseMiniSEED (SEXP buffer) { LOGICAL(dq_flags)[7] = (int) bit(b,0x80); */ PROTECT(dq_flags = NEW_INTEGER(8)); - numPROTECT++; + loopPROTECT++; for (int i=0; i<8; i++) { INTEGER(dq_flags)[i] = (int) total_dq_flags[i]; } PROTECT(timing_qual = NEW_NUMERIC(1)); - numPROTECT++; + loopPROTECT++; if (totalBlkt1001 > 0) { REAL(timing_qual)[0] = (double) totalTimingQuality / (double) totalRecords; } else { @@ -314,7 +321,7 @@ SEXP parseMiniSEED (SEXP buffer) { // All the data have been parsed. Now set up and fill the segmentList SEXP segmentList; PROTECT(segmentList = NEW_LIST(14)); - numPROTECT++; + loopPROTECT++; SET_VECTOR_ELT(segmentList, 0, npts); SET_VECTOR_ELT(segmentList, 1, sampling_rate); SET_VECTOR_ELT(segmentList, 2, network); @@ -330,10 +337,14 @@ SEXP parseMiniSEED (SEXP buffer) { SET_VECTOR_ELT(segmentList, 12, dq_flags); SET_VECTOR_ELT(segmentList, 13, timing_qual); setAttrib(segmentList, R_NamesSymbol, listNames); - + // Put the segmentList inside the returnList; SET_VECTOR_ELT(returnList, segIndex, segmentList); - + + // REC -- now UNPROTECT all elements built up in the loop + // data should now protected by returnList + UNPROTECT(loopPROTECT); // pop the stack down to pre-loop count + // Next segment seg = seg->next; } diff --git a/vignettes/IRISSeismic-intro.Rnw b/vignettes/IRISSeismic-intro.Rnw index 24b3df7..c1d6715 100644 --- a/vignettes/IRISSeismic-intro.Rnw +++ b/vignettes/IRISSeismic-intro.Rnw @@ -216,7 +216,8 @@ The \code{Trace} object has the following slots: \item \code{stats} -- \code{TraceHeader} object (metadata from the trace) \item \code{Sensor} -- instrument equipment name \item \code{InstrumentSensitivity} -- instrument total sensitivity (stage 0 gain) - \item \code{InputUnits} -- instrument data qcquisition input units + \item \code{SensitivityFrequency} -- the frequency (Hz) at which the \code{InstrumentSensitivity} is correct + \item \code{InputUnits} -- instrument data acquisition input units \item \code{data} -- vector of \code{numeric} data (the actual signal) \end{itemize} @@ -414,7 +415,7 @@ The \code{getAvailability()} method accepts the following arguments: \item \code{channel} -- channel code \item \code{starttime} -- POSIXct starttime (GMT) \item \code{endtime} -- POSIXct endtime (GMT) - \item \code{restricted} -- optional flag whether to report on restricted data (default=\code{FALSE}) + \item \code{includerestricted} -- optional flag whether to report on restricted data (default=\code{FALSE}) \item \code{latitude} -- optional latitude when specifying location and radius \item \code{longitude} -- optional longitude when specifying location and radius \item \code{minradius} -- optional minimum radius when specifying location and radius @@ -432,9 +433,11 @@ from web services of the same name. The suite of methods returning dataframes i \begin{itemize} \item \code{getAvailability} \item \code{getChannel} + \item \code{getDataselect} \item \code{getEvalresp} \item \code{getEvent} \item \code{getNetwork} + \item \code{getSNCL} \item \code{getStation} \item \code{getTraveltime} \item \code{getUnavailability} @@ -461,7 +464,7 @@ endtime <- starttime + 3600 * 24 * 2 # Find biggest seismic event over these two days -- it's the "Nisqually" events <- getEvent(iris, starttime, endtime, minmag=5.0) bigOneIndex <- which(events$magnitude == max(events$magnitude)) -bigOne <- events[bigOneIndex,] +bigOne <- events[bigOneIndex[1],] # Find US stations that are available within 10 degrees of arc of the # event location during the hour after the event