<jupyter_start><jupyter_text>Introduction to obspy The obspy package is very useful to download seismic data and to do some signal processing on them. Most signal processing methods are based on the signal processing method in the Python package scipy. First we import useful packages.<jupyter_code>import obspy
import obspy.clients.earthworm.client as earthworm
import obspy.clients.fdsn.client as fdsn
from obspy import read
from obspy import read_inventory
from obspy import UTCDateTime
from obspy.core.stream import Stream
from obspy.signal.cross_correlation import correlate

import matplotlib.pyplot as plt
import numpy as np
import os
import urllib.request

%matplotlib inline<jupyter_output><empty_output><jupyter_text>We are going to download data from an array of seismic stations.<jupyter_code>network = 'XU'
arrayName = 'BS'
staNames = ['BS01', 'BS02', 'BS03', 'BS04', 'BS05', 'BS06', 'BS11', 'BS20', 'BS21', 'BS22', 'BS23', 'BS24', 'BS25', \
    'BS26', 'BS27']
chaNames = ['SHE', 'SHN', 'SHZ']
staCodes = 'BS01,BS02,BS03,BS04,BS05,BS06,BS11,BS20,BS21,BS22,BS23,BS24,BS25,BS26,BS27'
chans = 'SHE,SHN,SHZ'<jupyter_output><empty_output><jupyter_text>We also need to define the time period for which we want to download data.<jupyter_code>myYear = 2010
myMonth = 8
myDay = 17
myHour = 6
TDUR = 2 * 3600.0
Tstart = UTCDateTime(year=myYear, month=myMonth, day=myDay, hour=myHour)
Tend = Tstart + TDUR<jupyter_output><empty_output><jupyter_text>We start by defining the client for downloading the data<jupyter_code>fdsn_client = fdsn.Client('IRIS')<jupyter_output><empty_output><jupyter_text>Download the seismic data for all the stations in the array.<jupyter_code>Dtmp = fdsn_client.get_waveforms(network=network, station=staCodes, location='--', channel=chans, starttime=Tstart, \
    endtime=Tend, attach_response=True)<jupyter_output><empty_output><jupyter_text>Some stations did not record the entire two hours. We delete these and keep only stations with a complte two hour recording.<jupyter_code>ntmp = []
for ksta in range(0, len(Dtmp)):
    ntmp.append(len(Dtmp[ksta]))
ntmp = max(set(ntmp), key=ntmp.count)
D = Dtmp.select(npts=ntmp)<jupyter_output><empty_output><jupyter_text>This is a function for plotting after each operation on the data.<jupyter_code>def plot_2hour(D, channel, offset, title):
    """ Plot seismograms
        D = Stream
        channel = 'E', 'N', or 'Z'
        offset = Offset between two stations
        title = Title of the figure
    """
    fig, ax = plt.subplots(figsize=(15, 10))
    Dplot = D.select(component=channel)
    t = (1.0 / Dplot[0].stats.sampling_rate) * np.arange(0, Dplot[0].stats.npts)
    for ksta in range(0, len(Dplot)):
        plt.plot(t, ksta * offset + Dplot[ksta].data, 'k')
    plt.xlim(np.min(t), np.max(t))
    plt.ylim(- offset, len(Dplot) * offset)
    plt.title(title, fontsize=24)
    plt.xlabel('Time (s)', fontsize=24)
    ax.set_yticklabels([])
    ax.tick_params(labelsize=20)
plot_2hour(D, 'E', 1200.0, 'Downloaded data')<jupyter_output><empty_output><jupyter_text>We start by detrending the data.<jupyter_code>D
D.detrend(type='linear')
plot_2hour(D, 'E', 1200.0, 'Detrended data')<jupyter_output><empty_output><jupyter_text>We then taper the data.<jupyter_code>D.taper(type='hann', max_percentage=None, max_length=5.0)
plot_2hour(D, 'E', 1200.0, 'Tapered data')<jupyter_output><empty_output><jupyter_text>And we remove the instrment response.<jupyter_code>D.remove_response(output='VEL', pre_filt=(0.2, 0.5, 10.0, 15.0), water_level=80.0)
plot_2hour(D, 'E', 1.0e-6, 'Deconvolving the instrument response')<jupyter_output><empty_output><jupyter_text>Then we filter the data.<jupyter_code>D.filter('bandpass', freqmin=2.0, freqmax=8.0, zerophase=True)
plot_2hour(D, 'E', 1.0e-6, 'Filtered data')<jupyter_output><empty_output><jupyter_text>And we resample the data.<jupyter_code>D.interpolate(100.0, method='lanczos', a=10)
D.decimate(5, no_filter=True)
plot_2hour(D, 'E', 1.0e-6, 'Resampled data')<jupyter_output><empty_output><jupyter_text>We can also compute the envelope of the signal.<jupyter_code>for index in range(0, len(D)):
    D[index].data = obspy.signal.filter.envelope(D[index].data)
plot_2hour(D, 'E', 1.0e-6, 'Envelope')<jupyter_output><empty_output><jupyter_text>You can also download the instrument response separately:<jupyter_code>network = 'XQ'
station = 'ME12'
channels = 'BHE,BHN,BHZ'
location = '01'<jupyter_output><empty_output><jupyter_text>This is to download the instrument response.<jupyter_code>fdsn_client = fdsn.Client('IRIS')
inventory = fdsn_client.get_stations(network=network, station=station, level='response')
inventory.write('response/' + network + '_' + station + '.xml', format='STATIONXML')<jupyter_output><empty_output><jupyter_text>We then read the data and start precessing the signal as we did above.<jupyter_code>fdsn_client = fdsn.Client('IRIS')
Tstart = UTCDateTime(year=2008, month=4, day=1, hour=4, minute=49)
Tend = UTCDateTime(year=2008, month=4, day=1, hour=4, minute=50)
D = fdsn_client.get_waveforms(network=network, station=station, location=location, channel=channels, starttime=Tstart, endtime=Tend, attach_response=False)
D.detrend(type='linear')
D.taper(type='hann', max_percentage=None, max_length=5.0)<jupyter_output><empty_output><jupyter_text>But we now use the xml file that contains the instrment response to remove it from the signal.<jupyter_code>filename = 'response/' + network + '_' + station + '.xml'
inventory = read_inventory(filename, format='STATIONXML')
D.attach_response(inventory)
D.remove_response(output='VEL', pre_filt=(0.2, 0.5, 10.0, 15.0), water_level=80.0)<jupyter_output><empty_output><jupyter_text>We resume signal processing.<jupyter_code>D.filter('bandpass', freqmin=2.0, freqmax=8.0, zerophase=True)
D.interpolate(100.0, method='lanczos', a=10)
D.decimate(5, no_filter=True)<jupyter_output><empty_output><jupyter_text>And we plot.<jupyter_code>t = (1.0 / D[0].stats.sampling_rate) * np.arange(0, D[0].stats.npts)
plt.plot(t, D[0].data, 'k')
plt.xlim(np.min(t), np.max(t))
plt.title('Single waveform', fontsize=18)
plt.xlabel('Time (s)', fontsize=18)<jupyter_output><empty_output><jupyter_text>Not all seismic data are stored on IRIS. This is an example of how to download data from the Northern California Earthquake Data Center (NCEDC).<jupyter_code>network = 'BK'
station = 'WDC'
channels = 'BHE,BHN,BHZ'
location = '--'<jupyter_output><empty_output><jupyter_text>This is to download the instrument response.<jupyter_code>url = 'http://service.ncedc.org/fdsnws/station/1/query?net=' + network + '&sta=' + station + '&level=response&format=xml&includeavailability=true'
s = urllib.request.urlopen(url)
contents = s.read()
file = open('response/' + network + '_' + station + '.xml', 'wb')
file.write(contents)
file.close()<jupyter_output><empty_output><jupyter_text>And this is to download the data.<jupyter_code>Tstart = UTCDateTime(year=2007, month=2, day=12, hour=1, minute=11, second=54)
Tend = UTCDateTime(year=2007, month=2, day=12, hour=1, minute=12, second=54)
request = 'waveform_' + station + '.request'
file = open(request, 'w')
message = '{} {} {} {} '.format(network, station, location, channels) + \
    '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d} '.format( \
    Tstart.year, Tstart.month, Tstart.day, Tstart.hour, Tstart.minute, Tstart.second) + \
    '{:04d}-{:02d}-{:02d}T{:02d}:{:02d}:{:02d}\n'.format( \
    Tend.year, Tend.month, Tend.day, Tend.hour, Tend.minute, Tend.second)
file.write(message)
file.close()
miniseed = 'station_' + station + '.miniseed'
request = 'curl -s --data-binary @waveform_' + station + '.request -o ' + miniseed + ' http://service.ncedc.org/fdsnws/dataselect/1/query'
os.system(request)
D = read(miniseed)
D.detrend(type='linear')
D.taper(type='hann', max_percentage=None, max_length=5.0)
filename = 'response/' + network + '_' + station + '.xml'
inventory = read_inventory(filename, format='STATIONXML')
D.attach_response(inventory)
D.remove_response(output='VEL', pre_filt=(0.2, 0.5, 10.0, 15.0), water_level=80.0)
D.filter('bandpass', freqmin=2.0, freqmax=8.0, zerophase=True)
D.interpolate(100.0, method='lanczos', a=10)
D.decimate(5, no_filter=True)
t = (1.0 / D[0].stats.sampling_rate) * np.arange(0, D[0].stats.npts)
plt.plot(t, D[0].data, 'k')
plt.xlim(np.min(t), np.max(t))
plt.title('Single waveform', fontsize=18)
plt.xlabel('Time (s)', fontsize=18)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15474209650582363, 'fraction_numerical': 0.035298312336581886, 'mean_word_length': 6.466725820763088, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 89, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '15569', 'license': 'MIT', 'n_tokens_mistral': 3382, 'n_tokens_neox': 3224, 'n_words': 659, 'path': 'tutorials/week3_pandas/obspy.ipynb', 'repo_name': 'geo-smart/ESS490-590-Spr21', 'size': '8414'}