User contributed plug-ins for Pyrocko's seismic waveform browser Snuffler.
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

4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
4 years ago
6 years ago
6 years ago
4 years ago
4 years ago
  1. from __future__ import print_function
  2. from pyrocko.gui.snuffling import Param, Snuffling, Choice
  3. import numpy as num
  4. def p2o_trace(ptrace, station):
  5. '''Convert Pyrocko trace to ObsPy trace.'''
  6. from obspy.core import UTCDateTime
  7. from obspy.core import Trace as oTrace
  8. otr = oTrace(
  9. data=ptrace.get_ydata(),
  10. header=dict(
  11. network=ptrace.network,
  12. station=ptrace.station,
  13. location=ptrace.location,
  14. channel=ptrace.channel,
  15. delta=ptrace.deltat,
  16. starttime=UTCDateTime(ptrace.tmin),
  17. coordinates=dict(
  18. latitude=station.lat,
  19. longitude=station.lon,
  20. elevation=station.elevation/1000.)))
  21. return otr
  22. class FK(Snuffling):
  23. '''
  24. <html>
  25. <head>
  26. <style type="text/css">
  27. body { margin-left:10px };
  28. </style>
  29. </head>
  30. <body>
  31. <h1 align='center'>FK ANALYSIS</h1>
  32. <h2 align='center'>based in ObsPy</h2>
  33. <p>
  34. <u>Prerequisites</u><br>
  35. This snuffling requires the ObsPy package which can be found at
  36. <a href='https://github.com/obspy/obspy/wiki'> ObsPy's github page </a>
  37. </p>
  38. <p>
  39. On way to install this package is to do:
  40. <PRE>
  41. git clone git://github.com/obspy/obspy.git
  42. cd obspy
  43. sudo python setup.py install
  44. </PRE>
  45. <p>
  46. <u>Usage</u><br>
  47. - Load station information at startup <br>
  48. - Zoom into the data until you see only data you desire to analyse or
  49. use extended markers to selected time regions for analysis<br>
  50. - Press the 'Run' button <br>
  51. </p>
  52. <p>
  53. The slowness is given in s/km.
  54. The assumed location of the geometrical center is printed to the terminal.
  55. </p>
  56. <p>
  57. Further information can be gathered from
  58. <a href="http://docs.obspy.org/master/tutorial/code_snippets/beamforming_fk_analysis.html"> # noqa
  59. ObsPy's FK tutorial</a>.
  60. </p>
  61. </body>
  62. </html>
  63. '''
  64. def setup(self):
  65. self.set_name('FK Analysis')
  66. self.add_parameter(Param('Slowness range[+-]', 'smax', 0.2, 0., 1.))
  67. self.add_parameter(Param(
  68. 'Number of slowness divisions', 'divisor', 20, 10, 50))
  69. self.add_parameter(Param(
  70. 'Number of radial sections', 'numberOfFraction', 32, 4, 50))
  71. self.add_parameter(Param(
  72. 'Length of Sliding Window [s]', 'window_lenth', 1., 0.5, 10.))
  73. self.add_parameter(Param(
  74. 'Step fraction of Sliding Window [s]', 'win_frac', 0.05, 0., 10.))
  75. self.add_parameter(Choice(
  76. 'If sampling rates differ', 'downresample', 'resample',
  77. ['resample', 'downsample', 'downsample to "target dt"']))
  78. self.add_parameter(Param('target dt', 'target_dt', 0.2, 0., 10))
  79. # self.add_parameter(Choice('Units: ','unit','[s/km]',('[s/km]','[s/deg]'))) # noqa
  80. self.set_live_update(False)
  81. def call(self):
  82. try:
  83. from obspy.core import UTCDateTime, stream
  84. from obspy.signal import array_analysis
  85. from obspy.imaging.cm import obspy_sequential as cmap
  86. except ImportError as _import_error:
  87. self.fail('ImportError:\n%s' % _import_error)
  88. from matplotlib.colorbar import ColorbarBase
  89. from matplotlib.colors import Normalize
  90. import matplotlib.dates as mdates
  91. self.cleanup()
  92. viewer = self.get_viewer()
  93. if viewer.lowpass is None or viewer.highpass is None:
  94. self.fail('highpass and lowpass in viewer must be set!')
  95. traces = []
  96. for trs in self.chopper_selected_traces(fallback=True):
  97. for tr in trs:
  98. tr.lowpass(2, viewer.lowpass)
  99. tr.highpass(2, viewer.highpass)
  100. traces.extend(trs)
  101. if not traces:
  102. self.fail('no traces selected')
  103. if self.downresample == 'resample':
  104. dt_want = min([t.deltat for t in traces])
  105. for t in traces:
  106. t.resample(dt_want)
  107. elif self.downresample == 'downsample':
  108. dt_want = max([t.deltat for t in traces])
  109. for t in traces:
  110. t.downsample_to(dt_want)
  111. elif self.downresample == 'downsample to "target dt"':
  112. for t in traces:
  113. t.downsample_to(float(self.target_dt))
  114. tmin = max([t.tmin for t in traces])
  115. tmax = min([t.tmax for t in traces])
  116. try:
  117. obspy_traces = [p2o_trace(
  118. tr, viewer.get_station(viewer.station_key(tr)))
  119. for tr in traces]
  120. except KeyError:
  121. self.fail('station information missing')
  122. st = stream.Stream(traces=obspy_traces)
  123. center = array_analysis.get_geometry(st, return_center=True)
  124. center_lon, center_lat, center_ele = center[len(center)-1]
  125. # Execute sonic
  126. kwargs = dict(
  127. sll_x=-self.smax, slm_x=self.smax, sll_y=-self.smax,
  128. slm_y=self.smax, sl_s=self.smax/self.divisor,
  129. win_len=self.window_lenth, win_frac=self.win_frac,
  130. frqlow=viewer.highpass, frqhigh=viewer.lowpass, prewhiten=0,
  131. semb_thres=-1.0e9, vel_thres=-1.0e9, verbose=True,
  132. timestamp='mlabday', stime=UTCDateTime(tmin),
  133. etime=UTCDateTime(tmax)
  134. )
  135. try:
  136. out = array_analysis.array_processing(st, **kwargs)
  137. except AttributeError:
  138. from obspy.signal.array_analysis import sonic
  139. out = sonic(st, **kwargs)
  140. pi = num.pi
  141. # make output human readable, adjust backazimuth to values between 0
  142. # and 360
  143. t, rel_power, abs_power, baz, slow = out.T
  144. baz[baz < 0.0] += 360.
  145. # choose number of fractions in plot (desirably 360 degree/N is an
  146. # integer!)
  147. N = int(self.numberOfFraction)
  148. abins = num.arange(N + 1) * 360. / N
  149. sbins = num.linspace(0., self.smax, N + 1)
  150. # sum rel power in bins given by abins and sbins
  151. hist, baz_edges, sl_edges = num.histogram2d(
  152. baz, slow, bins=[abins, sbins], weights=rel_power)
  153. # transform to gradient
  154. baz_edges = baz_edges / 180. * pi
  155. fig = self.pylab(get='figure')
  156. cax = fig.add_axes([0.85, 0.2, 0.05, 0.5])
  157. ax = fig.add_axes([0.10, 0.1, 0.70, 0.7], polar=True)
  158. ax.grid(False)
  159. dh = abs(sl_edges[1] - sl_edges[0])
  160. dw = abs(baz_edges[1] - baz_edges[0])
  161. # circle through backazimuth
  162. for i, row in enumerate(hist):
  163. ax.bar(left=(pi / 2 - (i + 1) * dw) * num.ones(N),
  164. height=dh * num.ones(N), width=dw,
  165. bottom=dh * num.arange(N), color=cmap(row / hist.max()))
  166. ax.set_xticks([pi / 2, 0, 3. / 2 * pi, pi])
  167. ax.set_xticklabels(['N', 'E', 'S', 'W'])
  168. ax.set_ylim(0., self.smax)
  169. ColorbarBase(cax, cmap=cmap,
  170. norm=Normalize(vmin=hist.min(), vmax=hist.max()))
  171. fig2 = self.pylab(get='figure')
  172. labels = ['rel.power', 'abs.power', 'baz', 'slow']
  173. xlocator = mdates.AutoDateLocator()
  174. ax = None
  175. for i, lab in enumerate(labels):
  176. ax = fig2.add_subplot(4, 1, i + 1, sharex=ax)
  177. ax.scatter(out[:, 0], out[:, i + 1], c=out[:, 1], alpha=0.6,
  178. edgecolors='none', cmap=cmap)
  179. ax.set_ylabel(lab)
  180. ax.set_xlim(out[0, 0], out[-1, 0])
  181. ax.set_ylim(out[:, i + 1].min(), out[:, i + 1].max())
  182. ax.xaxis.set_tick_params(which='both', direction='in')
  183. ax.xaxis.set_major_locator(xlocator)
  184. ax.xaxis.set_major_formatter(mdates.AutoDateFormatter(xlocator))
  185. if i != 3:
  186. ax.set_xticklabels([])
  187. fig2.subplots_adjust(hspace=0.)
  188. fig2.canvas.draw()
  189. fig.canvas.draw()
  190. print('Center of Array at latitude %s and longitude %s' %
  191. (center_lat, center_lon))
  192. def __snufflings__():
  193. return [FK()]