Maybe we would like to see the body waves more clearly and filter out
the surface waves. And instead of the raw velocity data, lets
apply the full instrument response to correct to displacement.
First, lets switch from the LH? channels, recorded
at one sample per second, to the HH? channels, recorded at 100 sps
to give us a wider frequency band to play with. Since we will need
the response information as well, we will need to change from
using
queryChannels()
to
queryResponses()
in the
stationQuery
. This will return stationXML with the full response filled in for each
channel. We can also the time window used in our dataselect query
to use offsets from the predicted arrivals, windowing in on just the
range we need around the P wave.
This is all possible by modifying the functionality
shown in the previous tutorial. However, because getting seismograms
for an earthquake at a channel based on predicted arrival times is
such a common activity, there is a tool,
SeismogramLoader
that does this for us.
The first steps will look a lot like the previous tutorial, where we set up the stationQuery and eventQuery with our search parameters.
let queryTimeWindow = seisplotjs.luxon.Interval.fromDateTimes(
seisplotjs.util.isoToDateTime('2019-07-01'),
seisplotjs.util.isoToDateTime('2019-07-31'));
let eventQuery = new seisplotjs.fdsnevent.EventQuery()
.timeWindow(queryTimeWindow)
.minMag(7)
.latitude(35).longitude(-118)
.maxRadius(3);
let stationQuery = new seisplotjs.fdsnstation.StationQuery()
.networkCode('CO')
.stationCode('HODGE')
.locationCode('00')
.channelCode('HH?')
.timeWindow(queryTimeWindow);
Of course 100 sps means that we will have 100 times as many samples, so lets reduce the time window to just the region where the P arrival is, say 300 seconds instead of 1800.
let loader = new seisplotjs.seismogramloader.SeismogramLoader(stationQuery, eventQuery);
loader.startOffset = -30;
loader.endOffset = 270;
loader.markedPhaseList = "origin,PP,pP";
loader.withResponse = true;
let loaderSDDPromise = loader.loadSeismograms();
We are not displaying maps here, but if we were then we can add a then function on loader.stationList and loader.quakeList, which are both Promises set up by the loadSeismograms method. These are useful when we need direct access to intermediate values as SeismogramLoader works across the various remote call. We can set the text in the span elements here just as an example. Note that the quake and channel are saved within the list of SeismogramDisplayData objects we will get at the end, so they are not lost even if we do not make use of these extra then methods.
let stationsPromise = loader.networkList.then(networkList => {
let staText = "";
for (let s of seisplotjs.stationxml.allStations(networkList)) {
staText += s.codes();
}
seisplotjs.d3.select('span#stationCode').text(staText);
});
let quakePromise = loader.quakeList.then( quakeList => {
let quakeText="";
for (const q of quakeList) {
quakeText+=quakeList[0].description+" ";
}
seisplotjs.d3.select('span#earthquakeDescription').text(quakeText);
});
Now we insert the filtering code after the seismograms have arrived
but before we plot them. We will insert another
then()
call on the actual loader Promise. We create a butterworth
filter using the sampling rate of the seismogram, with a passband
of 0.5 to 10 Hz. Removing the mean is usually a good idea, then
we apply the filter. Tapering is important before a deconvolution
and then we
transfer
the instrument response for the channel.
loaderSDDPromise.then((seismogramDataList) => {
seismogramDataList.forEach(sdd => {
let butterworth = seisplotjs.filter.createButterworth(
2, // poles
seisplotjs.filter.BAND_PASS,
.5, // low corner
10, // high corner
1/sdd.seismogram.sampleRate // delta (period)
);
let rmeanSeis = seisplotjs.filter.rMean(sdd.seismogram);
let filteredSeis = seisplotjs.filter.applyFilter(butterworth, rmeanSeis);
let taperSeis = seisplotjs.taper.taper(filteredSeis);
let correctedSeis = seisplotjs.transfer.transfer(taperSeis,
sdd.channel.response, .001, .02, 250, 500);
sdd.seismogram = correctedSeis;
});
return seismogramDataList ;
}).then( seismogramDataList => {
let div = document.querySelector('div#myseismograph');
let seisConfig = new seisplotjs.seismographconfig.SeismographConfig();
seisConfig.linkedTimeScale = new seisplotjs.scale.LinkedTimeScale();
seisConfig.linkedAmplitudeScale = new seisplotjs.scale.LinkedAmplitudeScale();
seisConfig.wheelZoom = false;
Then, just to make sure we don"t correct the data twice, we disable the gain correction and create the plots.
seisConfig.doGain = false;
seisConfig.centeredAmp = false;
for( let sdd of seismogramDataList) {
let graph = new seisplotjs.seismograph.Seismograph([sdd],
seisConfig);
div.appendChild(graph);
}
return seismogramDataList;
We can also plot the amplitude spectra for the three seismograms. We need to add an additional div to hold them.
<div id="fftplot">
</div>
And an additional style to size it.
div#fftplot {
height: 600px;
}
Then we calculate the fft and plot it.
}).then(seismogramDataList => {
let div = seisplotjs.d3.select('div#fftplot');
let fftList = seismogramDataList.map(sdd => {
if (sdd.seismogram.isContiguous()) {
return seisplotjs.fft.fftForward(sdd.seismogram)
} else {
return null; // can't do fft for non-contiguouus
}
}).filter(x => x); // to remove nulls
let seisConfig = new seisplotjs.seismographconfig.SeismographConfig();
let fftPlot = new seisplotjs.spectraplot.SpectraPlot(fftList, seisConfig);
document.querySelector('div#fftplot').appendChild(fftPlot);
return seismogramDataList;
});
Then just a little bit of cleanup to make sure everything finished correctly.
Promise.all( [ quakePromise, stationsPromise, loaderSDDPromise ] )
.catch( function(error) {
seisplotjs.d3.select("div#myseismograph").append('p').html("Error loading data." +error);
console.assert(false, error);
});
Previous: Predicted phase arrival times
Next: Helicorder