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_timesimport matplotlib.pyplot as pltfig, 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 TauPyModelmodel = 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_pathsimport matplotlib.pyplot as pltfig, 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 npimport matplotlib.pyplot as pltfrom obspy.taup import TauPyModelPHASES = [# 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 phasesfor 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 regionsax.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 / 2ax.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.014459Correlation coefficient: 0.92Bandpass prefiltering:Time correction for pick 2: -0.013025Correlation coefficient: 0.98

from __future__ import print_functionimport obspyfrom obspy.signal.cross_correlation import xcorr_pick_correction# read example data of two small earthquakespath = "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 analysist1 = 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 resultsdt, 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 prefilteringdt, 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)