You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
233 lines
7.9 KiB
233 lines
7.9 KiB
from __future__ import print_function
|
|
from pyrocko.gui.snuffling import Param, Snuffling, Choice
|
|
|
|
import numpy as num
|
|
|
|
|
|
def p2o_trace(ptrace, station):
|
|
'''Convert Pyrocko trace to ObsPy trace.'''
|
|
from obspy.core import UTCDateTime
|
|
from obspy.core import Trace as oTrace
|
|
|
|
otr = oTrace(
|
|
data=ptrace.get_ydata(),
|
|
header=dict(
|
|
network=ptrace.network,
|
|
station=ptrace.station,
|
|
location=ptrace.location,
|
|
channel=ptrace.channel,
|
|
delta=ptrace.deltat,
|
|
starttime=UTCDateTime(ptrace.tmin),
|
|
coordinates=dict(
|
|
latitude=station.lat,
|
|
longitude=station.lon,
|
|
elevation=station.elevation/1000.)))
|
|
|
|
return otr
|
|
|
|
|
|
class FK(Snuffling):
|
|
'''
|
|
<html>
|
|
<head>
|
|
<style type="text/css">
|
|
body { margin-left:10px };
|
|
</style>
|
|
</head>
|
|
<body>
|
|
<h1 align='center'>FK ANALYSIS</h1>
|
|
<h2 align='center'>based in ObsPy</h2>
|
|
<p>
|
|
<u>Prerequisites</u><br>
|
|
This snuffling requires the ObsPy package which can be found at
|
|
<a href='https://github.com/obspy/obspy/wiki'> ObsPy's github page </a>
|
|
</p>
|
|
<p>
|
|
On way to install this package is to do:
|
|
<PRE>
|
|
git clone git://github.com/obspy/obspy.git
|
|
cd obspy
|
|
sudo python setup.py install
|
|
</PRE>
|
|
<p>
|
|
<u>Usage</u><br>
|
|
|
|
- Load station information at startup <br>
|
|
|
|
- Zoom into the data until you see only data you desire to analyse or
|
|
use extended markers to selected time regions for analysis<br>
|
|
|
|
- Press the 'Run' button <br>
|
|
|
|
</p>
|
|
<p>
|
|
The slowness is given in s/km.
|
|
The assumed location of the geometrical center is printed to the terminal.
|
|
</p>
|
|
<p>
|
|
Further information can be gathered from
|
|
<a href="http://docs.obspy.org/master/tutorial/code_snippets/beamforming_fk_analysis.html"> # noqa
|
|
ObsPy's FK tutorial</a>.
|
|
</p>
|
|
</body>
|
|
</html>
|
|
'''
|
|
|
|
def setup(self):
|
|
self.set_name('FK Analysis')
|
|
self.add_parameter(Param('Slowness range[+-]', 'smax', 0.2, 0., 1.))
|
|
self.add_parameter(Param(
|
|
'Number of slowness divisions', 'divisor', 20, 10, 50))
|
|
self.add_parameter(Param(
|
|
'Number of radial sections', 'numberOfFraction', 32, 4, 50))
|
|
self.add_parameter(Param(
|
|
'Length of Sliding Window [s]', 'window_lenth', 1., 0.5, 10.))
|
|
self.add_parameter(Param(
|
|
'Step fraction of Sliding Window [s]', 'win_frac', 0.05, 0., 10.))
|
|
self.add_parameter(Choice(
|
|
'If sampling rates differ', 'downresample', 'resample',
|
|
['resample', 'downsample', 'downsample to "target dt"']))
|
|
self.add_parameter(Param('target dt', 'target_dt', 0.2, 0., 10))
|
|
|
|
# self.add_parameter(Choice('Units: ','unit','[s/km]',('[s/km]','[s/deg]'))) # noqa
|
|
self.set_live_update(False)
|
|
|
|
def call(self):
|
|
try:
|
|
from obspy.core import UTCDateTime, stream
|
|
from obspy.signal import array_analysis
|
|
from obspy.imaging.cm import obspy_sequential as cmap
|
|
except ImportError as _import_error:
|
|
self.fail('ImportError:\n%s' % _import_error)
|
|
|
|
from matplotlib.colorbar import ColorbarBase
|
|
from matplotlib.colors import Normalize
|
|
import matplotlib.dates as mdates
|
|
self.cleanup()
|
|
viewer = self.get_viewer()
|
|
|
|
if viewer.lowpass is None or viewer.highpass is None:
|
|
self.fail('highpass and lowpass in viewer must be set!')
|
|
|
|
traces = []
|
|
for trs in self.chopper_selected_traces(fallback=True):
|
|
for tr in trs:
|
|
tr.lowpass(2, viewer.lowpass)
|
|
tr.highpass(2, viewer.highpass)
|
|
|
|
traces.extend(trs)
|
|
|
|
if not traces:
|
|
self.fail('no traces selected')
|
|
|
|
if self.downresample == 'resample':
|
|
dt_want = min([t.deltat for t in traces])
|
|
for t in traces:
|
|
t.resample(dt_want)
|
|
|
|
elif self.downresample == 'downsample':
|
|
dt_want = max([t.deltat for t in traces])
|
|
for t in traces:
|
|
t.downsample_to(dt_want)
|
|
|
|
elif self.downresample == 'downsample to "target dt"':
|
|
for t in traces:
|
|
t.downsample_to(float(self.target_dt))
|
|
|
|
tmin = max([t.tmin for t in traces])
|
|
tmax = min([t.tmax for t in traces])
|
|
try:
|
|
obspy_traces = [p2o_trace(
|
|
tr, viewer.get_station(viewer.station_key(tr)))
|
|
for tr in traces]
|
|
|
|
except KeyError:
|
|
self.fail('station information missing')
|
|
|
|
st = stream.Stream(traces=obspy_traces)
|
|
center = array_analysis.get_geometry(st, return_center=True)
|
|
center_lon, center_lat, center_ele = center[len(center)-1]
|
|
|
|
# Execute sonic
|
|
kwargs = dict(
|
|
sll_x=-self.smax, slm_x=self.smax, sll_y=-self.smax,
|
|
slm_y=self.smax, sl_s=self.smax/self.divisor,
|
|
win_len=self.window_lenth, win_frac=self.win_frac,
|
|
frqlow=viewer.highpass, frqhigh=viewer.lowpass, prewhiten=0,
|
|
semb_thres=-1.0e9, vel_thres=-1.0e9, verbose=True,
|
|
timestamp='mlabday', stime=UTCDateTime(tmin),
|
|
etime=UTCDateTime(tmax)
|
|
)
|
|
|
|
try:
|
|
out = array_analysis.array_processing(st, **kwargs)
|
|
except AttributeError:
|
|
from obspy.signal.array_analysis import sonic
|
|
out = sonic(st, **kwargs)
|
|
|
|
pi = num.pi
|
|
|
|
# make output human readable, adjust backazimuth to values between 0
|
|
# and 360
|
|
t, rel_power, abs_power, baz, slow = out.T
|
|
baz[baz < 0.0] += 360.
|
|
|
|
# choose number of fractions in plot (desirably 360 degree/N is an
|
|
# integer!)
|
|
N = int(self.numberOfFraction)
|
|
abins = num.arange(N + 1) * 360. / N
|
|
sbins = num.linspace(0., self.smax, N + 1)
|
|
|
|
# sum rel power in bins given by abins and sbins
|
|
hist, baz_edges, sl_edges = num.histogram2d(
|
|
baz, slow, bins=[abins, sbins], weights=rel_power)
|
|
|
|
# transform to gradient
|
|
baz_edges = baz_edges / 180. * pi
|
|
|
|
fig = self.pylab(get='figure')
|
|
cax = fig.add_axes([0.85, 0.2, 0.05, 0.5])
|
|
ax = fig.add_axes([0.10, 0.1, 0.70, 0.7], polar=True)
|
|
ax.grid(False)
|
|
|
|
dh = abs(sl_edges[1] - sl_edges[0])
|
|
dw = abs(baz_edges[1] - baz_edges[0])
|
|
|
|
# circle through backazimuth
|
|
for i, row in enumerate(hist):
|
|
ax.bar(left=(pi / 2 - (i + 1) * dw) * num.ones(N),
|
|
height=dh * num.ones(N), width=dw,
|
|
bottom=dh * num.arange(N), color=cmap(row / hist.max()))
|
|
|
|
ax.set_xticks([pi / 2, 0, 3. / 2 * pi, pi])
|
|
ax.set_xticklabels(['N', 'E', 'S', 'W'])
|
|
ax.set_ylim(0., self.smax)
|
|
ColorbarBase(cax, cmap=cmap,
|
|
norm=Normalize(vmin=hist.min(), vmax=hist.max()))
|
|
|
|
fig2 = self.pylab(get='figure')
|
|
labels = ['rel.power', 'abs.power', 'baz', 'slow']
|
|
xlocator = mdates.AutoDateLocator()
|
|
ax = None
|
|
for i, lab in enumerate(labels):
|
|
ax = fig2.add_subplot(4, 1, i + 1, sharex=ax)
|
|
ax.scatter(out[:, 0], out[:, i + 1], c=out[:, 1], alpha=0.6,
|
|
edgecolors='none', cmap=cmap)
|
|
ax.set_ylabel(lab)
|
|
ax.set_xlim(out[0, 0], out[-1, 0])
|
|
ax.set_ylim(out[:, i + 1].min(), out[:, i + 1].max())
|
|
ax.xaxis.set_tick_params(which='both', direction='in')
|
|
ax.xaxis.set_major_locator(xlocator)
|
|
ax.xaxis.set_major_formatter(mdates.AutoDateFormatter(xlocator))
|
|
if i != 3:
|
|
ax.set_xticklabels([])
|
|
fig2.subplots_adjust(hspace=0.)
|
|
fig2.canvas.draw()
|
|
fig.canvas.draw()
|
|
|
|
print('Center of Array at latitude %s and longitude %s' %
|
|
(center_lat, center_lon))
|
|
|
|
|
|
def __snufflings__():
|
|
return [FK()]
|