Getting Started
For those who are not used to working with R, the Using R series of blog posts offers tips on how to get started and includes links to other introductory documentation.
RStudio
Users of the IRISSeismic package are encouraged to first download and install the RStudio integrated development environment for R. Newcomers to R will find RStudio a much friendlier environment in which to work.
First Example
Once you have an R environment up and running, the first step is to load the IRISSeismic package. Then you can create a new IrisClient object that will be responsible for all subsequent communication with DMC provided web services.
library(IRISSeismic)
iris <- new("IrisClient")In order to get data from one of the EarthScope web services we must specify all the information needed to create a webservice request: network, station, location, channel, starttime, endtime. Each unique combination of these elements is known as a SNCL. These elements are passed to the getDataselect() method of the IrisClient as a series of character strings except for the times which are of type POSIXct. The user is responsible for creating datetime objects of class POSIXct.
The first three commands in the following code chunk use the IrisClient object to communicate with web services and return a Stream object full of data from the EarthScope. The fourth line checks to see how many distinct chunks of seismic data exist. The last line passes this Stream object to a function that will plot the times at which this channel was collecting data.
starttime <- as.POSIXct("2002-04-20", tz="GMT")
endtime <- as.POSIXct("2002-04-21", tz="GMT")
result <- try(st <- getDataselect(iris,"US","OXF","","BHZ",starttime,endtime))
if (inherits(result,"try-error")) {
message(geterrmessage())
} else {
length(st@traces)
plotUpDownTimes(st, min_signal=1, min_gap=1)
}
This station had a few minor data dropouts, causing the data to be broken up into several separate signals that the IRISSeismic package stores in Trace objects.
We can get more information on the gaps between traces with the getGaps() function. The duration (secs) of gaps between traces is displayed along with the number of samples that were missed during the gap.
if (exists("st")){
getGaps(st)
}## $gaps
## [1] 0.0000 58.7750 57.0749 47.5750 52.1750 0.0000
##
## $nsamples
## [1] 0 2351 2283 1903 2087 0Next we can examine various statistics for each individual trace with the following parallel- functions.
if (exists("st")){
parallelLength(st)
parallelMax(st)
parallelSd(st)
}## [1] 101.14186 97.03080 484.53911 135.00670 93.05572It looks like the third trace, with a larger maximum and standard deviation, might have a signal. Metadata for this trace is stored in the stats slot of the Trace object.
if (exists("st")){
tr <- st@traces[[3]]
tr@stats
}## Seismic Trace TraceHeader
## Network: US
## Station: OXF
## Location:
## Channel: BHZ
## Quality: M
## calib: 1
## npts: 2163653
## sampling rate: 40
## delta: 0.025
## starttime: 2002-04-20 04:43:03
## endtime: 2002-04-20 19:44:34
## latitude: 34.5118
## longitude: -89.4092
## elevation: 101
## depth: 0
## azimuth: 0
## dip: -90
## processing:Finally, we can look at the seismic signal with the plot method.
if (exists("tr")){
plot(tr)
}
This small seismic signal was recorded in Oxford, Mississippi and is from a quake that occurred in New York state
Note: By default, data are subsampled before plotting to greatly! improve plotting speed. You can sometimes improve the appearance of a plot by reducing the amount of subsampling used. The plot method accepts a subsampling parameter to specify this.
Understanding Stream and Trace objects
In order to work effectively with the IRISSeismic package you must first understand the structure of the new S4 objects it defines. The package documentation gives a full description of each object but we can also interrogate them using the slotNames() function.
if (exists("st")){
slotNames(st)
}## [1] "url" "requestedStarttime" "requestedEndtime"
## [4] "act_flags" "io_flags" "dq_flags"
## [7] "timing_qual" "traces"The Stream object has the following slots (aka properties or attributes):
url-- full web services URL used to obtain datarequestedStarttime-- POSIXct datetime of the requested startrequestedEndtime-- POSIXct datetime of the requested endact_flags-- activity flags from the miniSEED recordio_flags-- I/O flags from the miniSEED recorddq_flags-- data quality flags from the miniSEED recordtiming_qual-- timing quality from the miniSEED recordtraces-- list ofTraceobjects
When in doubt about what a particular slot contains, it is always a good idea to ask what type of object it is.
if (exists("st")){
class(st@url)
class(st@requestedStarttime)
class(st@traces)
}## [1] "list"The next code chunk examines the first Trace in our Stream.
Note: R uses double square brackets, [[...]] to access list items.
if (exists("st")){
slotNames(st@traces[[1]])
}## [1] "id" "stats" "Sensor"
## [4] "InstrumentSensitivity" "SensitivityFrequency" "InputUnits"
## [7] "data"The Trace object has the following slots:
id-- SNCLQ identifier of the form "US.OXF..BHZ.M"stats--TraceHeaderobject (metadata from the trace)Sensor-- instrument equipment nameInstrumentSensitivity-- instrument total sensitivity (stage 0 gain)SensitivityFrequency-- the frequency (Hz) at which theInstrumentSensitivityis correctInputUnits-- instrument data acquisition input unitsdata-- vector ofnumericdata (the actual signal)
The TraceHeader metadata and the actual signal come from the dataselect webservice. The instrument metadata are obtained from the station webservice.
Be careful with times
Time stamps associated with seismic data should be given as "Universal" or "GMT" times. When specifying times to be used with methods of the IRISSeismic package you must be careful to specify the timezone as R assumes the local timezone by default.
Also, R assumes that datetime strings are formatted with a space separating date and time as opposed to the ISO 8601 'T' separator. If an ISO 8601 character string is provided without specific formatting instructions, the time portion of the string will be lost without any warning! So it is very important to be careful and consistent if you write code that converts ASCII strings into times.
A few examples will demonstrate the issues:
as.POSIXct("2010-02-27", tz="GMT") # good## [1] "2010-02-27 GMT"as.POSIXct("2010-02-27 04:00:00", tz="GMT") # good## [1] "2010-02-27 04:00:00 GMT"as.POSIXct("2010-02-27T04:00:00", tz="GMT",
format="%Y-%m-%dT%H:%M:%OS") # good## [1] "2010-02-27 04:00:00 GMT"as.POSIXct("2010-02-27") # BAD -- no timezone## [1] "2010-02-27 PST"as.POSIXct("2010-02-27T04:00:00", tz="GMT") # BAD -- no formatting## [1] "2010-02-27 GMT"
