obspy中文教程(二)
下面的脚本展示了如何对地震记录进行下采样。目前支持简单的整数下采样。如果未明确禁用,则在下采样之前应先进行低通滤波以防产生混叠。为了比较,也绘制了过滤过的未下采样数据。应用的处理步骤记录在每个Trace的trace.stats.processing中。请注意处理开始的相位,因为默认情况下应用的滤波器不是零相位类型。这可以通过手动应用零相滤波器并在下下采样期间停用自动滤波来避免(no_filter = True)。
import numpy as np import matplotlib.pyplot as plt import obspy # Read the seismogram st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new") # There is only one trace in the Stream object, let’s work on that trace... tr = st[0] # Decimate the 200 Hz data by a factor of 4 to 50 Hz. Note that this # automatically includes a lowpass filtering with corner frequency 20 Hz. # We work on a copy of the original data just to demonstrate the effects of # downsampling. tr_new = tr.copy() tr_new.decimate(factor=4, strict_length=False) # For comparison also only filter the original data (same filter options as in # automatically applied filtering during downsampling, corner frequency # 0.4 * new sampling rate) tr_filt = tr.copy() tr_filt.filter(’lowpass’, freq=0.4 * tr.stats.sampling_rate / 4.0) # Now let’s plot the raw and filtered data... t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta) t_new = np.arange(0, tr_new.stats.npts / tr_new.stats.sampling_rate,tr_new.stats.delta) plt.plot(t, tr.data, ’k’, label=’Raw’, alpha=0.3) plt.plot(t, tr_filt.data, ’b’, label=’Lowpassed’, alpha=0.7) plt.plot(t_new, tr_new.data, ’r’, label=’Lowpassed/Downsampled’, alpha=0.7) plt.xlabel(’Time [s]’) plt.xlim(82, 83.5) plt.suptitle(tr.stats.starttime) plt.legend() plt.show()
下面的例子显示了如何合并然后绘制三个重叠的地震图,最长的一个被认为是正确的。 请参阅merge()方法的文档。
import numpy as np import matplotlib.pyplot as plt import obspy # Read in all files starting with dis. st = obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE") st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.1") st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.2") # sort st.sort([’starttime’]) # start time in plot equals 0 dt = st[0].stats.starttime.timestamp # Go through the stream object, determine time range in julian seconds # and plot the data with a shared x axis ax = plt.subplot(4, 1, 1) # dummy for tying axis for i in range(3): plt.subplot(4, 1, i + 1, sharex=ax) t = np.linspace(st[i].stats.starttime.timestamp - dt,st[i].stats.endtime.timestamp - dt,st[i].stats.npts) plt.plot(t, st[i].data) # Merge the data together and show plot in a similar way st.merge(method=1) plt.subplot(4, 1, 4, sharex=ax) t = np.linspace(st[0].stats.starttime.timestamp - dt, st[0].stats.endtime.timestamp - dt, st[0].stats.npts) plt.plot(t, st[0].data, ’r’) plt.show()
- Beamforming - FK Analysis (FK分析)
下列代码展示了如何使用obspy进行FK分析。这些数据来自····
我们应用下列设置执行array_processing():
- 设置慢度网络拐点-3.0:3.0s/km,步长sl_s=0.03
- 窗口长1s,步长0.05
- 数据经过1~8Hz的带通滤波,不使用预白化
- semb_thres和vel_thres设置为无穷小的数字,不得更改。
- 时间戳timestamp设置为“mlabday”他可以被我们的绘图程序直接读取
- stime和etime定为UTCDate Time格式
输出被存储在out中。
下半部分展示了如何绘制输出。我们使用array_processing()输出的out,它是包含了时间戳、相对功率、绝对功率、反量和慢度的numpy ndarray。色条代表相对功率。
import matplotlib.pyplot as plt import matplotlib.dates as mdates import obspy from obspy.core.util import AttribDict from obspy.imaging.cm import obspy_sequential from obspy.signal.invsim import corn_freq_2_paz from obspy.signal.array_analysis import array_processing # Load data st = obspy.read("https://examples.obspy.org/agfa.mseed") # Set PAZ and coordinates for all 5 channels st[0].stats.paz = AttribDict({ 'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)], 'zeros': [0j, 0j], 'sensitivity': 205479446.68601453, 'gain': 1.0}) st[0].stats.coordinates = AttribDict({ 'latitude': 48.108589, 'elevation': 0.450000, 'longitude': 11.582967}) st[1].stats.paz = AttribDict({ 'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)], 'zeros': [0j, 0j], 'sensitivity': 205479446.68601453, 'gain': 1.0}) st[1].stats.coordinates = AttribDict({ 'latitude': 48.108192, 'elevation': 0.450000, 'longitude': 11.583120}) st[2].stats.paz = AttribDict({ 'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)], 'zeros': [0j, 0j], 'sensitivity': 250000000.0, 'gain': 1.0}) st[2].stats.coordinates = AttribDict({ 'latitude': 48.108692, 'elevation': 0.450000, 'longitude': 11.583414}) st[3].stats.paz = AttribDict({ 'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j)], 'zeros': [0j, 0j], 'sensitivity': 222222228.10910088, 'gain': 1.0}) st[3].stats.coordinates = AttribDict({ 'latitude': 48.108456, 'elevation': 0.450000, 'longitude': 11.583049}) st[4].stats.paz = AttribDict({ 'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j), (-2.105 + 0j)], 'zeros': [0j, 0j, 0j], 'sensitivity': 222222228.10910088, 'gain': 1.0}) st[4].stats.coordinates = AttribDict({ 'latitude': 48.108730, 'elevation': 0.450000, 'longitude': 11.583157}) # Instrument correction to 1Hz corner frequency paz1hz = corn_freq_2_paz(1.0, damp=0.707) st.simulate(paz_remove='self', paz_simulate=paz1hz) # Execute array_processing stime = obspy.UTCDateTime("20080217110515") etime = obspy.UTCDateTime("20080217110545") kwargs = dict( # slowness grid: X min, X max, Y min, Y max, Slow Step sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03, # sliding window properties win_len=1.0, win_frac=0.05, # frequency properties frqlow=1.0, frqhigh=8.0, prewhiten=0, # restrict output semb_thres=-1e9, vel_thres=-1e9, timestamp='mlabday', stime=stime, etime=etime ) out = array_processing(st, **kwargs) # Plot labels = ['rel.power', 'abs.power', 'baz', 'slow'] xlocator = mdates.AutoDateLocator() fig = plt.figure() for i, lab in enumerate(labels): ax = fig.add_subplot(4, 1, i + 1) ax.scatter(out[:, 0], out[:, i + 1], c=out[:, 1],alpha=0.6, edgecolors='none', cmap=obspy_sequential) 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_major_locator(xlocator) ax.xaxis.set_major_formatter(mdates.AutoDateFormatter(xlocator)) fig.suptitle('AGFA skyscraper blasting in Munich %s' % ( stime.strftime('%Y-%m-%d'), )) fig.autofmt_xdate() fig.subplots_adjust(left=0.15, top=0.95, right=0.95, bottom=0.2, hspace=0) plt.show()
也可以使用极坐标图表示,其将网格箱中的相对功率相加,每个网格箱都由背调角和要分析的部分信号的慢度定义。反向量从北向开始顺时针计数,慢度限制可以手动设定。
import numpy as np import matplotlib.pyplot as plt from matplotlib.colorbar import ColorbarBase from matplotlib.colors import Normalize import obspy from obspy.core.util import AttribDict from obspy.imaging.cm import obspy_sequential from obspy.signal.invsim import corn_freq_2_paz from obspy.signal.array_analysis import array_processing # Load data st = obspy.read("https://examples.obspy.org/agfa.mseed") # Set PAZ and coordinates for all 5 channels st[0].stats.paz = AttribDict({ 'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)], 'zeros': [0j, 0j], 'sensitivity': 205479446.68601453, 'gain': 1.0}) st[0].stats.coordinates = AttribDict({ 'latitude': 48.108589, 'elevation': 0.450000, 'longitude': 11.582967}) st[1].stats.paz = AttribDict({ 'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)], 'zeros': [0j, 0j], 'sensitivity': 205479446.68601453, 'gain': 1.0}) st[1].stats.coordinates = AttribDict({ 'latitude': 48.108192, 'elevation': 0.450000, 'longitude': 11.583120}) st[2].stats.paz = AttribDict({ 'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)], 'zeros': [0j, 0j], 'sensitivity': 250000000.0, 'gain': 1.0}) st[2].stats.coordinates = AttribDict({ 'latitude': 48.108692, 'elevation': 0.450000, 'longitude': 11.583414}) st[3].stats.paz = AttribDict({ 'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j)], 'zeros': [0j, 0j], 'sensitivity': 222222228.10910088, 'gain': 1.0}) st[3].stats.coordinates = AttribDict({ 'latitude': 48.108456, 'elevation': 0.450000, 'longitude': 11.583049}) st[4].stats.paz = AttribDict({ 'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j), (-2.105 + 0j)], 'zeros': [0j, 0j, 0j], 'sensitivity': 222222228.10910088, 'gain': 1.0}) st[4].stats.coordinates = AttribDict({ 'latitude': 48.108730, 'elevation': 0.450000, 'longitude': 11.583157}) # Instrument correction to 1Hz corner frequency paz1hz = corn_freq_2_paz(1.0, damp=0.707) st.simulate(paz_remove='self', paz_simulate=paz1hz) # Execute array_processing kwargs = dict( # slowness grid: X min, X max, Y min, Y max, Slow Step sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03, # sliding window properties win_len=1.0, win_frac=0.05, # frequency properties frqlow=1.0, frqhigh=8.0, prewhiten=0, # restrict output semb_thres=-1e9, vel_thres=-1e9, stime=obspy.UTCDateTime("20080217110515"), etime=obspy.UTCDateTime("20080217110545") ) out = array_processing(st, **kwargs) # Plot cmap = obspy_sequential # 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 = 36 N2 = 30 abins = np.arange(N + 1) * 360. / N sbins = np.linspace(0, 3, N2 + 1) # sum rel power in bins given by abins and sbins hist, baz_edges, sl_edges = \ np.histogram2d(baz, slow, bins=[abins, sbins], weights=rel_power) # transform to radian baz_edges = np.radians(baz_edges) # add polar and colorbar axes fig = plt.figure(figsize=(8, 8)) 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.set_theta_direction(-1) ax.set_theta_zero_location("N") 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): bars = ax.bar(left=(i * dw) * np.ones(N2), height=dh * np.ones(N2), width=dw, bottom=dh * np.arange(N2), color=cmap(row / hist.max())) ax.set_xticks(np.linspace(0, 2 * np.pi, 4, endpoint=False)) ax.set_xticklabels(['N', 'E', 'S', 'W']) # set slowness limits ax.set_ylim(0, 3) [i.set_color('grey') for i in ax.get_yticklabels()] ColorbarBase(cax, cmap=cmap, norm=Normalize(vmin=hist.min(), vmax=hist.max())) plt.show()
- Seismogram Envelopes(信号包络)
以下脚本显示如何过滤地震数据并将其与其包络信号一起绘制图像。本示例使用零相移带通滤波器(1~3Hz)进行滤波。然后我们计算出包络并将其与Trace一起绘制。
import numpy as np import matplotlib.pyplot as plt import obspy import obspy.signal st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new") data = st[0].data npts = st[0].stats.npts samprate = st[0].stats.sampling_rate # Filtering the Stream object st_filt = st.copy() st_filt.filter('bandpass', freqmin=1, freqmax=3, corners=2, zerophase=True) # Envelope of filtered data data_envelope = obspy.signal.filter.envelope(st_filt[0].data) # The plotting, plain matplotlib t = np.arange(0, npts / samprate, 1 / samprate) plt.plot(t, st_filt[0].data, 'k') plt.plot(t, data_envelope, 'k:') plt.title(st[0].stats.starttime) plt.ylabel('Filtered Data w/ Envelope') plt.xlabel('Time [s]') plt.xlim(80, 90) plt.show()
- Plotting Spectrograms(绘制频谱图)
下列代码展示了如何使用obspy Stream对象绘制频谱图。其中的很多选项都可以自定义,参考spectrogram()获取更多细节。比如,可以很轻松的通过调用matplotlib.crm中预定义的色彩图用来调整所需,这里提供几个教程:
- http://www.physics.ox.ac.uk/users/msshin/science/code/matplotlib_cm/
- http://matplotlib.org/examples/color/colormaps_reference.html
- https://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
- import obspy st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new") st.spectrogram(log=True, title='BW.RJOB ' + str(st[0].stats.starttime))