if (FALSE) {
# 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[1],]
# Find US stations that are available within an hour of the event
start <- bigOne$time
end <- start + 3600
availability <- getAvailability(iris, "US", "", "", "BHZ",
starttime=start, endtime=end,
latitude=bigOne$latitude, longitude=bigOne$longitude,
minradius=0, maxradius=10)
# Get the station the furthest East
minLonIndex <- which(availability$longitude == max(availability$longitude))
snclE <- availability[minLonIndex,]
# Plot the BHZ signal from this station
st <- getDataselect(iris,snclE$network,snclE$station,snclE$location,snclE$channel,
start,end)
# Check that there is only a single trace and then plot it
length(st@traces)
tr <- st@traces[[1]]
plot(tr, subsampling=1) # need subsmpling=1 to add vertical lines with abline()
# Find travel times to this station
traveltimes <- getTraveltime(iris, bigOne$latitude, bigOne$longitude, bigOne$depth,
snclE$latitude, snclE$longitude)
# Look at the list
traveltimes
# mark the P and S arrival times
pArrival <- start + traveltimes$travelTime[traveltimes$phaseName=="P"]
sArrival <- start + traveltimes$travelTime[traveltimes$phaseName=="S"]
abline(v=pArrival, col='red')
abline(v=sArrival, col='blue')
}
Run the code above in your browser using DataLab