In [1]:
%matplotlib inline
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import os,sys
sys.path.append('../')
from event_functions import *
from seispy import event
from seispy.trace import Trace
from astropy.table import Table
import os
sta='RRDG'
freq=0.4
/home/meyers/opt/seispy/lib/python2.7/site-packages/matplotlib/__init__.py:1401: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)

Grab data¶

We're taking the 38th event in Gary's database that we have from a long time ago. We're going to look at relative phase between radial and vertical channels in RRDG station at 0.4 Hz (period = 2.5s)

In [2]:
tab = event.EventTable.read_seis_db('../db/gary.db','../db/pat.window')
myevent = event.Event(*tab[9])
tableHHR, envelopesHHR, dataHHR = get_event_table(myevent, [sta], [freq],
                        filter_chan='HHR', write_file='events_38')
tableHHZ, envelopesHHZ, dataHHZ = get_event_table(myevent, [sta], [freq],
                        filter_chan='HHZ', write_file='events_38')
/home/meyers/opt/seispy/lib/python2.7/site-packages/astropy/table/column.py:264: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  return self.data.__eq__(other)
In [3]:
datHHR = dataHHR[sta][freq]['HHR']
datHHZ = dataHHZ[sta][freq]['HHZ']
In [4]:
plt.plot(datHHR.times.value, datHHR.value, label='HHR')
ax = plt.gca()
ax.plot(datHHZ.times.value, datHHZ.value, label='HHZ')
plt.title('Event 38 Time Series, $r$-wave chunk', fontsize=14)
plt.ylabel('Velocity [m/s]')
plt.legend()
plt.show()

Check cross-correlation¶

First, let's check the relative phase between the channels in the first 20s of the data. We do this by maximizing the time-series correlation as a function of time. The time-lag between stations tells us something about their relative phase. Specifically, the phase in radians should be:

$$ \phi = 2 \pi f t_{lag}$$
In [5]:
ARR1 = np.hstack((np.zeros(datHHR[:2000].size), datHHR[:2000].value, np.zeros(datHHR[:2000].size)))
corr = np.correlate(ARR1, datHHZ.value[:2000])[1:-1]
times = np.arange(datHHZ[:2000].size)[1:]
times = np.hstack((-times[::-1],[0],times))
In [6]:
plt.plot(times/100.,corr)
plt.xlabel('Time Lags')
plt.show()
In [7]:
print 'Time lag is: %f' % (times[np.argmax(corr)]/100.)
print 'This corresponds to a phase of: %f' % np.degrees(2*np.pi*(freq*times[np.argmax(corr)]/100.))
Time lag is: -0.650000
This corresponds to a phase of: -93.600000

Check phase during whole waveform on 2T chunks¶

Next we check the phase in each chunk that is twice the filter period. So we check on 5-second time-scales. We see that the phase tends to sit around 90 degrees (which is what we expec)

In [8]:
chunk = 500
idx0=0
phases = []
sts = []
while idx0+chunk < datHHR.size:
    ARR1 = np.hstack((np.zeros(datHHR[idx0:idx0+chunk].size), datHHR[idx0:idx0+chunk].value, np.zeros(datHHR[idx0:idx0+chunk].size)))
    corr = np.correlate(ARR1, datHHZ.value[idx0:idx0+chunk])[1:-1]
    times = np.arange(datHHZ[idx0:idx0+chunk].size)[1:]
    times = np.hstack((-times[::-1],[0],times))/100.
    phases.append(np.degrees(2*np.pi*(freq*times[np.argmax(corr)])))
    sts.append(idx0/100)
    idx0+=chunk
In [9]:
plt.plot(sts, phases)
plt.xlabel('start time [s]')
plt.ylabel('$\phi$ between radial and vertical')
plt.show()
In [ ]: