In Lesson 06 we will roll up our sleeves and begin doing some real work with seismic data using functions availabe in the seismic package. We will begin by creating a new IrisClient object that will be responsible for all subsequent communication with DMC provided web services.

Accessing Seismic Data

In order to get data from one of the IRIS DMC 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 IRIS DMC. 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.

library(IRISSeismic)

# Create a new IrisClient for access to DMC web services
iris <- new("IrisClient")

# Speicfy a SNCL and time range and get seismic data
starttime <- as.POSIXct("2002-04-20", tz="GMT")
endtime <- as.POSIXct("2002-04-21", tz="GMT")
st <- getDataselect(iris,"US","OXF","","BHZ",starttime,endtime)

# How many individual traces?
length(st@traces)

## [1] 5

# Use the plotUpDownTimes function to show gaps
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.

getGaps(st)

## $gaps
## [1]  0.0000 58.7750 57.0749 47.5771 52.1750  0.0000
## 
## $nsamples
## [1]    0 2351 2283 1903 2087    0

Next we can examine various statistics for each individual Trace with the following parallel~ functions.

parallelLength(st)

## [1]   10733  663953 2163653  308769  300267

parallelMax(st)

## [1]  1218  1356 38406  6139  1305

parallelSd(st)

## [1] 101.14186  97.03080 484.53911 135.00670  93.05572

It 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.

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 
##  processing:

Finally, we can look at the seismic signal with the plot() method of the Trace class.

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, the Trace class plot() method subsamples data before 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.

Closer Examination of the Signal

Once seismic data are in memory, performing mathematical analysis on those data can be very fast. All mathematical operations are performed on every data point.

Note: The plot() method of Stream objects deals with gaps by first calling mergeTraces() to fill all gaps with missing values (NA). Then the single, merged Trace is plotted with the plot() method for Trace objects. Any gaps of a significant size will be now visible in the resulting plot. By default, the plot() method of Trace and Stream objects subsamples the data so that approximately 5,000 points are used in the plot. This dramatically speeds up plotting.

One of the first things you may wish to do with a full day’s worth of seismic signal is clip it to a region of interest. One way to do that would be to modify the starttime and endtime parameters to getDataselect() and then make a data request covering a shorter period of time. A simpler technique, if the signal is already in memory, is to use the slice() method.

starttime <- as.POSIXct("2010-02-27", tz="GMT")
endtime <- as.POSIXct("2010-02-28", tz="GMT")
st <- getDataselect(iris,"IU","ANMO","00","BHZ",starttime,endtime)
start2 <- as.POSIXct("2010-02-27 06:40:00", tz="GMT")
end2 <- as.POSIXct("2010-02-27 07:40:00", tz="GMT")

# Extract trace and trace subset
tr1 <- st@traces[[1]]
tr2 <- slice(tr1, start2, end2)

old_par <- par()              # save old graphics parameters
layout(matrix(seq(2)))        # layout a 2x1 matrix
plot(tr1)                     # top plot
plot(tr2)                     # bottom plot

par(old_par)                  # restore original parameters

Detecting Events with STA/LTA

Access to triggering algorithms for detecting events is provided by the STALTA() method of Trace objects. (see A Comparison of Select Trigger Algorithms for Automated Global Seismic Phase and Event Detection The STALTA() method has the following arguments and defaults:

x – Trace being analyzed
staSecs – size of the short window in secs
ltaSecs – size of the long window in secs
algorithm – named algorithm (default=”classic LR”)
demean – should the signal have the mean removed (default=TRUE)
detrend – should the signal have the trend removed (default=TRUE)
taper – portion of the seismic signal to be tapered at each end (default=0.0)
increment – integer increment to use when sliding the averaging windows to the next location (default=1)

The STALTA() method returns a picker, a vector of numeric values, one for every value in the Trace@data slot. Note that this is a fairly compute-intensive operation. This picker can then be used with the triggerOnset() function to return the approximate start of the seismic signal.

Lets test this with our original seismic signal:

starttime <- as.POSIXct("2002-04-20", tz="GMT")
endtime <- as.POSIXct("2002-04-21", tz="GMT")
st <- getDataselect(iris,"US","OXF","","BHZ",starttime,endtime)
tr <- st@traces[[3]]
picker <- STALTA(tr,3,30)
threshold <- quantile(picker,0.99999,na.rm=TRUE)
to <- triggerOnset(tr,picker,threshold)

NOTE: The STALTA() method is intended to be used for crude, automatic event detection, not precise determination of signal arrival. Optimal values for the arguments to the STALTA() method will depend on the details of the seismic signal.

Event Windowing

The eventWindow() method allows you to focus on the region identified by the picker by automatically finding the trigger onset time and then slicing out the region of the trace centered on that time. This method has the following arguments and defaults:

x – Trace being analyzed
picker – picker returned by STALTA()
threshold – threshold value as used in triggerOnset()
windowSecs – total window size (secs)

old_par <- par()              # save old graphics parameters
layout(matrix(seq(3)))        # layout a 3x1 matrix
closeup1 <- eventWindow(tr,picker,threshold,3600)
closeup2 <- eventWindow(tr,picker,threshold,600)
plot(tr)
abline(v=to, col='red', lwd=2)
plot(closeup1)
abline(v=to, col='red', lwd=2)
plot(closeup2)
abline(v=to, col='red', lwd=2)

par(old_par)                  # restore original parameters

Data Availability

The IrisClient also provides functionality for interacting with other web services at the DMC. The getAvailability() method allows users to query what SNCLs are available, obtaining that information from the station webservice.

Information is returned as a dataframe containing all the information returned by the wsavailability web service. Standard DMC webservice wildcards can be used as in the example below which tells us what other ’B’ channels are available at our station of interest during the time of the big quake above:

starttime <- as.POSIXct("2010-02-27", tz="GMT")
endtime <- as.POSIXct("2010-02-28", tz="GMT")
availability <- getAvailability(iris,"IU","ANMO","*","B??",starttime,endtime)
availability

##   network station location channel latitude longitude elevation depth
## 1      IU    ANMO       00     BH1 34.94598 -106.4571    1671.0 145.0
## 2      IU    ANMO       00     BH2 34.94598 -106.4571    1671.0 145.0
## 3      IU    ANMO       00     BHZ 34.94598 -106.4571    1671.0 145.0
## 4      IU    ANMO       10     BH1 34.94591 -106.4571    1767.2  48.8
## 5      IU    ANMO       10     BH2 34.94591 -106.4571    1767.2  48.8
## 6      IU    ANMO       10     BHZ 34.94591 -106.4571    1767.2  48.8
##   azimuth dip                            instrument       scale scalefreq
## 1     328   0 Geotech KS-54000 Borehole Seismometer  3456610000      0.02
## 2      58   0 Geotech KS-54000 Borehole Seismometer  3344370000      0.02
## 3       0 -90 Geotech KS-54000 Borehole Seismometer  3275080000      0.02
## 4      64   0  Guralp CMG3-T Seismometer (borehole) 32805600000      0.02
## 5     154   0  Guralp CMG3-T Seismometer (borehole) 32655000000      0.02
## 6       0 -90  Guralp CMG3-T Seismometer (borehole) 33067200000      0.02
##   scaleunits samplerate           starttime             endtime
## 1        M/S         20 2008-06-30 20:00:00 2011-02-18 19:11:00
## 2        M/S         20 2008-06-30 20:00:00 2011-02-18 19:11:00
## 3        M/S         20 2008-06-30 20:00:00 2011-02-18 19:11:00
## 4        M/S         40 2008-06-30 20:00:00 2011-02-19 06:53:00
## 5        M/S         40 2008-06-30 20:00:00 2011-02-19 06:53:00
## 6        M/S         40 2008-06-30 20:00:00 2011-02-19 06:53:00
##           snclId
## 1 IU.ANMO.00.BH1
## 2 IU.ANMO.00.BH2
## 3 IU.ANMO.00.BHZ
## 4 IU.ANMO.10.BH1
## 5 IU.ANMO.10.BH2
## 6 IU.ANMO.10.BHZ

The getAvailability() method accepts the following arguments:

obj – an IrisClient object
network – network code
station – station code
location – location code
channel – channel code
starttime – POSIXct starttime (GMT)
endtime – POSIXct endtime (GMT)
restricted – optional flag whether to report on restricted data (de- fault=FALSE)
latitude – optional latitude when specifying location and radius
longitude – optional longitude when specifying location and radius
minradius – optional minimum radius when specifying location and radius
maxradius – optional maximum radius when specifying location and ra- dius

Other IRIS DMC web services

Several methods of the IrisClient class work very similarly to the getAvailability() method in that they return dataframes of information obtained from web services of the same name. The suite of methods returning dataframes includes:

getAvailability()
getChannel()
getEvalresp()
getEvent()
getNetwork()
getStation()
getTraveltime()
getUnavailability()

The following example demonstrates the use of several of these services together to do the following:

find seismic events on a particular day
find available US network BHZ channels in the hour after the biggest event that day
determine the easternmost of those channels
get the P and S travel times to that station
plot the seismic signal detected at that station with markers for P and S arrival times

# Open a connection to IRIS DMC webservices
iris <- new("IrisClient")

# Two days around the "Nisqually Quake"
starttime <- as.POSIXct("2001-02-27", tz="GMT")
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,]

# Find US stations that are available within 10 degrees of arc of the 
# event location during the hour after the event
start <- bigOne$time
end <- start + 3600
av <- getAvailability(iris, "US", "", "", "BHZ", start, end,
                      latitude=bigOne$latitude, longitude=bigOne$longitude,
                      minradius=0, maxradius=10)

# Get the station the furthest East
minLonIndex <- which(av$longitude == max(av$longitude))
snclE <- av[minLonIndex,]

# Get travel times to this station
traveltimes <- getTraveltime(iris, bigOne$latitude, bigOne$longitude, bigOne$depth,
                             snclE$latitude, snclE$longitude)

# Look at the list                             
traveltimes

##   distance depth phaseName travelTime rayParam takeoff incident
## 1     8.89    56         P     125.70   13.680   86.58    45.53
## 2     8.89    56         S     225.12   24.518   84.88    47.81
## 3     8.89    56       PcP     506.70    0.849    3.55     2.54
## 4     8.89    56       ScS     927.91    1.565    3.65     2.71
## 5     8.89    56     PKiKP     987.09    0.199    0.83     0.59
## 6     8.89    56     SKiKS    1405.21    0.222    0.52     0.38
##   puristDistance puristName
## 1           8.89          P
## 2           8.89          S
## 3           8.89        PcP
## 4           8.89        ScS
## 5           8.89      PKiKP
## 6           8.89      SKiKS

# Find the P and S arrival times
pArrival <- start + traveltimes$travelTime[traveltimes$phaseName=="P"]
sArrival <- start + traveltimes$travelTime[traveltimes$phaseName=="S"] 

# Get the BHZ signal for this station
st <- getDataselect(iris,snclE$network,snclE$station,
                    snclE$location,snclE$channel,
                    start,end)

# Check that there is only a single trace
length(st@traces)

## [1] 1

# Plot the seismic trace and mark the "P" and "S" arrival times
tr <- st@traces[[1]]
plot(tr, subsampling=1) # need subsmpling=1 to add vertical lines with abline()
abline(v=pArrival, col='red')
abline(v=sArrival, col='blue')

There are many more features for seismic signal processing available in the IRISSeismic package. You should have enough skills at this point to begin your own investigations of seismic signals.

Task 1: Explore your favorite earthquake

pick a quake you are interested in
find nearby stations
explore this quake’s signal at some of these stations.

‘seismic’ Packages < prev | next > MUSTANG Metrics

Lesson 06 – Working with Seismic Traces

Mazama Science