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. We will need the stationXML with the full response filled in for each channel.
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 = sp.util.startEnd('2019-07-01', '2019-07-31');
let eventQuery = new sp.fdsnevent.EventQuery()
.timeRange(queryTimeWindow)
.minMag(7)
.latitude(35).longitude(-118)
.maxRadius(3);
let stationQuery = new sp.fdsnstation.StationQuery()
.networkCode('CO')
.stationCode('HODGE')
.locationCode('00')
.channelCode('HH?')
.timeRange(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. Notice we are passing the station and event queries directly to the loader, it can take either queries or the returned arrays. We will also mark the origin time as if it were a predicted travel time.
let loader = new sp.seismogramloader.SeismogramLoader(stationQuery, eventQuery);
loader.startOffset = -30;
loader.endOffset = 270;
loader.markedPhaseList = "origin,PP,pP";
loader.withResponse = true;
loader.load().then(dataset => {
We can set the text in the span elements here just as an example. Note that the quake and channel are also saved within the list of SeismogramDisplayData objects so they are not lost even if we do not make use of these extra inventory and catalog fields in the dataset.
let staText = "";
for (let s of sp.stationxml.allStations(dataset.inventory)) {
staText += s.codes();
}
document.querySelector('span#stationCode').textContent=staText;
let quakeText="";
for (const q of dataset.catalog) {
quakeText+=q.description+" ";
}
document.querySelector('span#earthquakeDescription').textContent =quakeText;
Now we insert the filtering code after the seismograms have arrived
but before we plot them. 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. These are saved as the
processedWaveforms
of the dataset.
dataset.processedWaveforms = dataset.waveforms.map(sdd => {
let butterworth = sp.filter.createButterworth(
2, // poles
sp.filter.BAND_PASS,
.5, // low corner
10, // high corner
1/sdd.seismogram.sampleRate // delta (period)
);
let rmeanSeis = sp.filter.rMean(sdd.seismogram);
let filteredSeis = sp.filter.applyFilter(butterworth, rmeanSeis);
let taperSeis = sp.taper.taper(filteredSeis);
let correctedSeis = sp.transfer.transfer(taperSeis,
sdd.channel.response, .001, .02, 250, 500);
sdd.seismogram = correctedSeis;
return sdd;
});
let div = document.querySelector('div#myseismograph');
let seisConfig = new sp.seismographconfig.SeismographConfig();
seisConfig.linkedTimeScale = new sp.scale.LinkedTimeScale();
seisConfig.linkedAmplitudeScale = new sp.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;
for( let sdd of dataset.processedWaveforms) {
let graph = new sp.seismograph.Seismograph([sdd], seisConfig);
div.appendChild(graph);
}
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.
let fftList = dataset.processedWaveforms.map(sdd => {
if (sdd.seismogram.isContiguous()) {
return sp.fft.fftForward(sdd.seismogram)
} else {
return null; // can't do fft for non-contiguouus
}
}).filter(x => x); // to remove nulls
let fftSeisConfig = new sp.seismographconfig.SeismographConfig();
let fftPlot = new sp.spectraplot.SpectraPlot(fftList, fftSeisConfig);
document.querySelector('div#fftplot').appendChild(fftPlot);
return dataset;
}).catch( function(error) {
const div = document.querySelector('div#myseismograph');
div.innerHTML = `<p>Error loading data. ${error}</p>`;
console.assert(false, error);
});
Previous: Predicted phase arrival times
Next: Helicorder