<jupyter_start><jupyter_text>Rotational Seismology            Estimate Event Backazimuth             ![](images/obspy_logo_full_524x179px.png) Seismo-Live: http://seismo-live.org Authors:* <NAME> ([@salve-](https://github.com/salve-))---<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 15, 8<jupyter_output><empty_output><jupyter_text>1 Load pre-processed waveforms and event-xmlWaveforms of transverse acceleration and vertical rotation rate are loaded from the repository. They were pre-processed as described in the **Data Download + Pre-processing**-Tutorial. That means they are: instrument correctedresampled (decimation factor: 4)filtered (lowpass: $f_{cutoff}$ = 2 Hz)We also need the locations of station and event epicenter to obtain the theoretical backazimuth.<jupyter_code>from obspy import read, read_events
from obspy.geodetics.base import gps2dist_azimuth
import numpy as np

AC = read('data/acc_Turkey_preproc.mseed')
RLAS = read('data/rot_Turkey_preproc.mseed')
cat = read_events('data/xml_Turkey.xml')
event = cat[0]

# event location from event info
source_latitude = event.origins[0].latitude
source_longitude = event.origins[0].longitude

# station location (Wettzell)
station_latitude = 49.144001
station_longitude = 12.8782

print(cat)
print(event.event_descriptions[0]['type'], ': ',event.event_descriptions[0]['text'])<jupyter_output>1 Event(s) in Catalog:
2011-10-23T10:41:27.500000Z | +38.640,  +43.400 | 7.1 MW
Flinn-Engdahl region :  TURKEY<jupyter_text>** Obtain the theoretical backazimuth via gps2distazimuth**<jupyter_code># theoretical backazimuth and distance
baz = gps2dist_azimuth(source_latitude, source_longitude, station_latitude, station_longitude)

print('Epicentral distance [m]: ',np.round(baz[0],1))
print('Theoretical azimuth [deg]: ', np.round(baz[1],1))
print('Theoretical backazimuth [deg]: ', np.round(baz[2],1))<jupyter_output>Epicentral distance [m]:  2691458.3
Theoretical azimuth [deg]:  305.6
Theoretical backazimuth [deg]:  104.1<jupyter_text>2 Estimate correlation coefficients and backazimuth for each time windowThe estimation of the backazimuth (BAz = source direction angle) is based on *Igel et al.(2007)* and illustrated as a simple example here:We define the length of the time windows, we want to use and calculate the number of time windows for our downloaded and pre-processed waveform data. For each time-window, we rotate the horizontal acceleration components through all possible source directions (BAz, 0-360°) in **10°-steps**. In each of the steps, we check the similarity (phase-match) between vertical rotation rate and the rotated acceleration by cross-correlation (X-corr).The rotation angle/ BAz-value yielding the highest correlation coefficient for the specific time window rotates the  horizontal acceleration components to radial/ transverse (= into the angular direction of the strongest signal source).That's exactly what *McLeod et al. (1998) and Pancha et al. (2000)* showed for Love waves:**Vertical rotation rate and transverse acceleration are in phase!**<jupyter_code>from obspy.signal.cross_correlation import xcorr
from obspy.signal.rotate import rotate_ne_rt

sampling_rate = int(RLAS[0].stats.sampling_rate)
sec = 60  # window length for correlation
num_windows = len(RLAS[0]) // (int(sampling_rate * sec))

# estimate the Backazimuth for each time window
step = 10
backas = np.linspace(0, 360 - step, 360 / step)
corrbaz = []
ind=None
for i_deg in range(0, len(backas)):
    for i_win in range(0, num_windows):
        corrbazz = xcorr(RLAS[0][sampling_rate * sec * i_win : sampling_rate * sec * (i_win + 1)],
                             rotate_ne_rt(AC.select(component='N')[0].data,
                                          AC.select(component='E')[0].data, backas[i_deg])
                             [1][sampling_rate * sec * i_win : sampling_rate * sec * (i_win + 1)],0)
        corrbaz.append(corrbazz[1])
corrbaz = np.asarray(corrbaz)
corrbaz = corrbaz.reshape(len(backas), num_windows)

maxcorr = []
for l1 in range(0, num_windows):
    maxcor_r = backas[corrbaz[:, l1].argmax()]
    maxcorr.append(maxcor_r)
maxcorr = np.asarray(maxcorr)
X, Y = np.meshgrid(np.arange(0, sec * num_windows, sec), backas)<jupyter_output>/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/ipykernel_launcher.py:10: DeprecationWarning: object of type <class 'float'> cannot be safely interpreted as an integer.
  # Remove the CWD from sys.path while we load stuff.
/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/obspy/signal/cross_correlation.py:220: ObsPyDeprecationWarning: Call to deprecated function xcorr(). Please use the correlate and xcorr_max functions.
  warnings.warn(msg, ObsPyDeprecationWarning)<jupyter_text>Get theoretical **P-** and **S-Wave first arrivals** from a 1D velocity model using **TauPy**:<jupyter_code>from obspy.taup import TauPyModel

TauPy_model = TauPyModel('ak135')
arrivals_p = TauPy_model.get_travel_times(distance_in_degree=0.001 * baz[0] / 111.11,
                                        source_depth_in_km=event.origins[0].depth*0.001,
                                       phase_list=["P","p","Pdiff","PP","PKiKP","PKIKP","Pn","Pg"])
arrivals_s = TauPy_model.get_travel_times(distance_in_degree=0.001 * baz[0] / 111.11,
                                        source_depth_in_km=event.origins[0].depth*0.001,
                                       phase_list=["S","s","Sdiff","SS","SKiKS","SKIKS","Sn","Sg"])
tiemp = []
tiems = []
for i in range(0,len(arrivals_p)): tiemp.append(arrivals_p[i].time)
for ii in range(0,len(arrivals_s)): tiems.append(arrivals_s[ii].time)

# first arrivals
arriv_p = min(tiemp)
arriv_s = min(tiems)<jupyter_output>/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/obspy/taup/tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)<jupyter_text>3 Plot backazimuth estimation**Top:** Vertical rotation rate**Middle:** Transverse acceleration (rotated from N-E component acceleration, see above)**Bottom: Backazimuth estimation**Green displays anti-correlation, red shows positive correlation/ phase match.Yellow indicates no correlation.Black dots highlight the BAz values for the maximum correlation coefficient in that time window.The theoretical backazimuth is plotted for comparison.<jupyter_code>import matplotlib as mpl

time = np.linspace(0, len(AC[0].data)/sampling_rate,len(AC[0].data)) ## Time-array
# vertical rotation rate
plt.subplot2grid((3, 30), (0, 0), colspan=29)
plt.plot(time, RLAS[0].data,label='vertical rotation rate')
plt.xlim(0, time[-1])
plt.ylabel('vert. rot. rate \n[nrad/s]')
plt.legend()

# add P- and S-wave arrivals
plt.axvline(arriv_p);plt.annotate('P-arrival', xy=(arriv_p+20,np.max(RLAS[0].data)),xycoords='data');
plt.axvline(arriv_s);plt.annotate('S-arrival', xy=(arriv_s+20,np.max(RLAS[0].data)),xycoords='data');

# transverse acceleration
plt.subplot2grid((3, 30), (1, 0), colspan=29)
plt.plot(time, AC[0].data, 'k',label='transverse acceleration')
plt.xlim(0, time[-1])
plt.ylabel('transv. acc. \n[$nm/s^2$]')
plt.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
plt.legend()
plt.axvline(arriv_p)
plt.axvline(arriv_s)

# backazimuth estimation plot
plt.subplot2grid((3, 30), (2, 0), colspan=29)
im = plt.pcolor(X, Y, corrbaz, cmap=plt.cm.RdYlGn_r, vmin=-1,vmax=1)
plt.plot(np.arange(sec/2., sec * num_windows, sec), maxcorr, '.k')
plt.xlim(0, time[-1])
plt.ylim(0, 360)
plt.ylabel(u'estimated \nbackazimuth [°]')
plt.xlabel('time [s]')

# plot theoretical Backazimuth for comparison
xx = np.arange(0, sec * num_windows + 1, sec)
tba = np.ones(len(xx)) * baz[2]
plt.plot(xx, tba, c='.5', lw=1.5)
plt.text(2600, 106, u'Theor. BAz = '+str(round(baz[2],2))+'°', color='k',fontsize=12,backgroundcolor='.9')


# add colorbar
fig = plt.subplot2grid((3, 30), (2, 29))
norm = mpl.colors.Normalize(vmin=-1, vmax=1)
cb1 = mpl.colorbar.ColorbarBase(fig, cmap=plt.cm.RdYlGn_r, norm=norm, orientation='vertical')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13334961533764805, 'fraction_numerical': 0.03321528880205153, 'mean_word_length': 5.017634092578986, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 3}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '129425', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 3066, 'n_tokens_neox': 2790, 'n_words': 683, 'path': 'notebooks/Rotational Seismology/estimate_backazimuth.ipynb', 'repo_name': 'krischer/seismo_live_build', 'size': '8203'}