obspy中文教程(七)
- Visualize Data Availability of Local Waveform Archive(可视化本地波形存档数据的可用性)
通常,您拥有大量数据并希望知道哪个站点在何时是可用的。对于这种假设的情况,obspy提供了obspy-scan脚本(安装后即可用),它能从文件的数据头检测文件格式(MiniSEED, SAC, SACXY, GSE2, SH-ASC, SH-Q, SEISAN, 等),在间隙处绘制为垂直红线,可用数据的在开始时间绘制十字,数据本身绘为水平线。该脚本可以扫描超过1000个文件(已有被用于扫描30000个文件,耗时约45分钟。),自动绘制年/月范围。它会打开一个可放大的交互式绘图窗口。
从命令提示符执行类似下面的语句,使用通配符匹配文件:
$ obspy-scan /bay_mobil/mobil/20090622/1081019/*_1.*
- Travel Time and Ray Path Plotting(走时和射线路径绘制)
下面的代码展示如何plot_travel_times()函数绘制给定距离和相位的使用iasp91速度模型计算出的走时。
from obspy.taup import plot_travel_times import matplotlib.pyplot as plt fig, ax = plt.subplots() ax = plot_travel_times(source_depth=10, ax=ax, fig=fig, phase_list=['P', 'PP', 'S'], npoints=200)
下面的几行代码展示了如何绘制给定距离和相位的射线路径。射线路径使用iasp91速度模型计算,并使用obspy.taup.tau.Arrivals类的plot_rays()函数的绘制(在笛卡尔坐标系中)。
from obspy.taup import TauPyModel model = TauPyModel(model='iasp91') arrivals = model.get_ray_paths(500, 140, phase_list=['PP', 'SSS']) arrivals.plot_rays(plot_type='cartesian', phase_list=['PP', 'SSS'], plot_all=False, legend=True)
下面的几行代码展示了如何绘制给定距离和相位的射线路径。射线路径使用iasp91速度模型计算,并使用obspy.taup.tau.Arrivals类的plot_rays()函数的绘制(在球形图中)。
下面的几行代码展示了如何绘制有多个震中距和相位的射线路径。射线路径使用iasp91速度模型计算,并使用obspy.taup.tau.Arrivals类的plot_ray_paths()函数的绘制(在球形图中)。
from obspy.taup.tau import plot_ray_paths import matplotlib.pyplot as plt fig, ax = plt.subplots(subplot_kw=dict(polar=True)) ax = plot_ray_paths(source_depth=100, ax=ax, fig=fig, phase_list=['P', 'PKP'], npoints=25)
对于单个震中距离的射线路径示例,请尝试上一节中的plot_rays()方法。以下是一个更高级的示例,其中包含自定义的相位和距离列表:
import numpy as np import matplotlib.pyplot as plt from obspy.taup import TauPyModel PHASES = [ # Phase, distance ('P', 26), ('PP', 60), ('PPP', 94), ('PPS', 155), ('p', 3), ('pPcP', 100), ('PKIKP', 170), ('PKJKP', 194), ('S', 65), ('SP', 85), ('SS', 134.5), ('SSS', 204), ('p', -10), ('pP', -37.5), ('s', -3), ('sP', -49), ('ScS', -44), ('SKS', -82), ('SKKS', -120), ] model = TauPyModel(model='iasp91') fig, ax = plt.subplots(subplot_kw=dict(polar=True)) # Plot all pre-determined phases for phase, distance in PHASES: arrivals = model.get_ray_paths(700, distance, phase_list=[phase]) ax = arrivals.plot_rays(plot_type='spherical', legend=False, label_arrivals=True, plot_all=True, show=False, ax=ax) # Annotate regions ax.text(0, 0, 'Solid\ninner\ncore', horizontalalignment='center', verticalalignment='center', bbox=dict(facecolor='white', edgecolor='none', alpha=0.7)) ocr = (model.model.radius_of_planet - (model.model.s_mod.v_mod.iocb_depth + model.model.s_mod.v_mod.cmb_depth) / 2) ax.text(np.deg2rad(180), ocr, 'Fluid outer core', horizontalalignment='center', bbox=dict(facecolor='white', edgecolor='none', alpha=0.7)) mr = model.model.radius_of_planet - model.model.s_mod.v_mod.cmb_depth / 2 ax.text(np.deg2rad(180), mr, 'Solid mantle', horizontalalignment='center', bbox=dict(facecolor='white', edgecolor='none', alpha=0.7)) plt.show()
- Cross Correlation Pick Correction(交叉相关拾取校正)
该示例展示如何对齐两个地震的起始波形相位,以便纠正在常规分析中无法完全设置一致的原始拾取时间。按照[Deichmann1992]的方法,互相关函数的凹陷部分最大值附近可用抛物线拟合。
为调整参数并验证检查结果,可以选择展示图形或者将其存为图像文件。参见xcorr_pick_correction()。
该示例将打印拾取序列2的时间校正和相应的相关系数,并打开原始和预处理数据相关性的绘图窗口:
No preprocessing: Time correction for pick 2: -0.014459 Correlation coefficient: 0.92 Bandpass prefiltering: Time correction for pick 2: -0.013025 Correlation coefficient: 0.98
from __future__ import print_function import obspy from obspy.signal.cross_correlation import xcorr_pick_correction # read example data of two small earthquakes path = "https://examples.obspy.org/BW.UH1..EHZ.D.2010.147.%s.slist.gz" st1 = obspy.read(path % ("a", )) st2 = obspy.read(path % ("b", )) # select the single traces to use in correlation. # to avoid artifacts from preprocessing there should be some data left and # right of the short time window actually used in the correlation. tr1 = st1.select(component="Z")[0] tr2 = st2.select(component="Z")[0] # these are the original pick times set during routine analysis t1 = obspy.UTCDateTime("2010-05-27T16:24:33.315000Z") t2 = obspy.UTCDateTime("2010-05-27T16:27:30.585000Z") # estimate the time correction for pick 2 without any preprocessing and open # a plot window to visually validate the results dt, coeff = xcorr_pick_correction(t1, tr1, t2, tr2, 0.05, 0.2, 0.1, plot=True) print("No preprocessing:") print(" Time correction for pick 2: %.6f" % dt) print(" Correlation coefficient: %.2f" % coeff) # estimate the time correction with bandpass prefiltering dt, coeff = xcorr_pick_correction(t1, tr1, t2, tr2, 0.05, 0.2, 0.1, plot=True, filter="bandpass", filter_options={'freqmin': 1, 'freqmax': 10}) print("Bandpass prefiltering:") print(" Time correction for pick 2: %.6f" % dt) print(" Correlation coefficient: %.2f" % coeff)