3.5. Magnitude calculation

EQcorrscan contains both an automatic amplitude picker and a singular-value decomposition derived magnitude calculation, which is very accurate but requires high levels of event similarity.

3.5.1. Amplitude picker for local magnitudes

Currently this is only implemented for seisan s-files, however we have no plans to extend this to other formats (although it would be simple as seisan files are read in as obspy events not).

This example requires data downloaded from the eqcorrscan github repository.

>>> from eqcorrscan.utils.mag_calc import amp_pick_event
>>> from obspy import read, read_events
>>> from obspy.io.nordic.core import readwavename
>>> import os
>>> testing_path = 'eqcorrscan/tests/test_data'
>>> sfile = os.path.join(testing_path, 'REA', 'TEST_',
...                      '01-0411-15L.S201309')
>>> datapath = os.path.join(testing_path, 'WAV', 'TEST_')
>>> event = read_events(sfile)[0]
>>> st = read(os.path.join(datapath, readwavename(sfile)[0]))
>>> respdir = testing_path
>>> event = amp_pick_event(
...     event=event, st=st, respdir=respdir, chans=['Z'],
...     var_wintype=True, winlen=0.9, pre_pick=0.2, pre_filt=True,
...     lowcut=1.0, highcut=20.0, corners=4) 
Working on ...

3.5.2. Relative moment by singular-value decomposition

This method closely follows the method outlined by Rubinstein & Ellsworth 2010.

This example requires data downloaded from the eqcorrscan github repository.

>>> from eqcorrscan.utils.mag_calc import svd_moments
>>> from obspy import read
>>> import glob
>>> from eqcorrscan.utils.clustering import svd
>>> import numpy as np
>>> # Do the set-up
>>> testing_path = 'eqcorrscan/tests/test_data/similar_events'
>>> stream_files = glob.glob(os.path.join(testing_path, '*'))
>>> stream_list = [read(stream_file) for stream_file in stream_files]
>>> event_list = []
>>> for i, stream in enumerate(stream_list): 
...     st_list = []
...     for tr in stream:
...         # Only use the vertical channels of sites with known high similarity.
...         # You do not need to use this step for your data.
...         if (tr.stats.station, tr.stats.channel) not in\
...                 [('WHAT2', 'SH1'), ('WV04', 'SHZ'), ('GCSZ', 'EHZ')]:
...             stream.remove(tr)
...             continue
...         tr.detrend('simple')
...         tr.filter('bandpass', freqmin=5.0, freqmax=15.0)
...         tr.trim(tr.stats.starttime + 40, tr.stats.endtime - 45)
...         st_list.append(i)
...     event_list.append(st_list)
<obspy.core.stream.Stream object at ...>
>>> event_list = np.asarray(event_list).T.tolist()
>>> SVectors, SValues, Uvectors, stachans = svd(stream_list=stream_list)
>>> M, events_out = svd_moments(u=Uvectors, s=SValues, v=SVectors,
...                             stachans=stachans, event_list=event_list) 
Created Kernel matrix: ...