obspy中文教程（五）

 

Interfacing R from Python（从python对接到R）

rpy2包允许python对接到R。下例展示如何转换numpy.ndarray数据为R矩阵，并对其执行R的summary命令。

>>>from obspy.core import read

>>>import rpy2.robjects as RO

>>>import rpy2.robjects.numpy2ri

>>>r = RO.r

>>>st = read("test/BW.BGLD..EHE.D.2008.001")

>>>M = RO.RMatrix(st[0].data)

>>>print(r.summary(M))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

-1056.0  -409.0  -393.0  -393.7  -378.0   233.0

Coordinate Conversions（坐标转换）

使用pyproj可以方便的进行坐标转换。在查找完源和目标坐标系的EPSG代码后，只需几行代码即可完成坐标转换。下例将两个德国台站的坐标信息转换为区域使用的Gauß-Krüger坐标。

>>> import pyproj

>>> lat = [49.6919, 48.1629]

>>> lon = [11.2217, 11.2752]

>>> proj_wgs84 = pyproj.Proj(init="epsg:4326")

>>> proj_gk4 = pyproj.Proj(init="epsg:31468")

>>> x, y = pyproj.transform(proj_wgs84, proj_gk4, lon, lat)

>>> print(x)

[4443947.179347951, 4446185.667319892]

>>> print(y)

[5506428.401023342, 5336354.054996853]

另一种常见用法是将纬度和经度的位置信息转换到Universal Transverse Mercator coordinate system (UTM)坐标系。这对于小区域中的大密集阵列特别有用。 使用utm包可以轻松完成这种转换。

>>> import utm

>>> utm.from_latlon(51.2, 7.5)

(395201.3103811303, 5673135.241182375, 32, ’U’)

>>> utm.to_latlon(340000, 5710000, 32, ’U’)

(51.51852098408468, 6.693872395145327)

 

 

Hierarchical Clustering（分级聚类）

使用SciPy包中提供的方法可以执行分级聚类。它允许从相似性矩阵构建聚类并制作树状图。以下示例显示了如何对已计算的相似性矩阵执行此操作。相似性数据是根据具有诱发地震活动的区域中的事件计算的（使用obspy.signal中的相关例程），并且可以从我们的示例网络服务器获取：

首先，导入必要的模块并通过我们的网页载入数据：

>>> import io, urllib

>>> import numpy as np

>>> import matplotlib.pyplot as plt

>>> from scipy.cluster import hierarchy

>>> from scipy.spatial import distance



>>> url = "https://examples.obspy.org/dissimilarities.npz"

>>> with io.BytesIO(urllib.urlopen(url).read()) as fh, np.load(fh) as data:

... dissimilarity = data[’dissimilarity’]

现在我们绘制相异矩阵：

>>> plt.subplot(121)

>>> plt.imshow(1 - dissimilarity, interpolation="nearest")

然后我们使用SciPy构建并绘制树形图到图像右侧的子图中：

>>> dissimilarity = distance.squareform(dissimilarity)



>>> threshold = 0.3



>>> linkage = hierarchy.linkage(dissimilarity, method="single")



>>> clusters = hierarchy.fcluster(linkage, threshold, criterion="distance")



>>> plt.subplot(122)



>>> hierarchy.dendrogram(linkage, color_threshold=0.3)



>>> plt.xlabel("Event number")



>>> plt.ylabel("Dissimilarity")



>>> plt.show()

Visualizing Probabilistic Power Spectral Densities（可视化概率功率谱密度）

下列代码展示使用obspy.signal中定义的PPSD类，该例程对于解释现场质量控制检查的噪声测量非常有用。更多信息参考[McNamara2004]。

>>> from obspy import read

>>> from obspy.io.xseed import Parser

>>> from obspy.signal import PPSD

读取数据并选择一条所需的台站/信道轨迹：

>>> st = read("https://examples.obspy.org/BW.KW1..EHZ.D.2011.037")

>>> tr = st.select(id="BW.KW1..EHZ")[0]

元数据可以是：StationXML文件或FDSN网络服务提供的Inventory；无数据SEED文件的Parser；本地RESP文件的文件名；传统的零极点字典。

初始化一个新的PPSD语句。ppsd对象将确保随后只有适当的数据进入概率psd统计。

>>> parser = Parser("https://examples.obspy.org/dataless.seed.BW_KW1")

>>> ppsd = PPSD(tr.stats, metadata=parser)

然后添加数据（Trace或Stream对象）做PPSD估计。

>>> ppsd.add(st)
True

我们可以检查ppsd估计中表示的时间范围。它包含计算psds的一小时长切片的开始时间的排序列表（此处仅打印显示前两个）。

>>> print(ppsd.times[:2])

[UTCDateTime(2011, 2, 6, 0, 0, 0, 935000), UTCDateTime(2011, 2, 6, 0, 30, 0, 935000)]

>>> print("number of psd segments:", len(ppsd.times))

number of psd segments: 47

再次添加同样的stream无效（返回值False），ppsd对象确保没有重复的数据段进入ppsd估计。

>>> ppsd.add(st)

False

>>> print("number of psd segments:", len(ppsd.times))

number of psd segments: 47

可逐步添加来自其他文件/来源的更多信息。

>>> st = read("https://examples.obspy.org/BW.KW1..EHZ.D.2011.038")

>>> ppsd.add(st)

True

ppsd的图形表示可以显示在matplotlib窗口中。

>>> ppsd.plot()

或者保存为图片文件。

>>> ppsd.plot("/tmp/ppsd.png")  

>>> ppsd.plot("/tmp/ppsd.pdf")

对于每个频率仓，还可以显示累积的直方图：

>>> ppsd.plot(cumulative=True)

 

要使用PQLX /[McNamara2004]使用的色彩映射，可以从obspy.imaging.cm导入并使用该色彩映射：

>>> from obspy.imaging.cm import pqlx

>>> ppsd.plot(cmap=pqlx)

实际PPSD（[McNamara2004]）下方的是可视化PPSD的数据基础。第一行显示输入PPSD的数据，绿条表示可用数据，红条表示添加到PPSD的流中的间隙。蓝色的底行显示进入直方图的单个psd测量值。默认处理方法用零填充间隙，然后这些数据段显示为单个偏离psd行。

Note：从无数据SEED或StationXML文件提供元数据比指定静态零极点信息更安全（参见PPSD）。

psd值的时间序列也可以通过访问psd_values从PPSD中提取并使用plot_temporal()绘制。

>>> ppsd.plot_temporal([0.1, 1, 10])

使用plot_spectrogram()绘制频谱图：

>>> ppsd.plot_spectrogram()

 

Array Response Function（数组响应函数）

下面的代码块展示了如何使用ObsPy的obspy.signal.array_analysis.array_transff_wavenumber()函数绘制波束形成的数组传递函数（波数的函数）。

import numpy as np

import matplotlib.pyplot as plt



from obspy.imaging.cm import obspy_sequential

from obspy.signal.array_analysis import array_transff_wavenumber





# generate array coordinates

coords = np.array([[10., 60., 0.], [200., 50., 0.], [-120., 170., 0.],

                   [-100., -150., 0.], [30., -220., 0.]])



# coordinates in km

coords /= 1000.



# set limits for wavenumber differences to analyze

klim = 40.

kxmin = -klim

kxmax = klim

kymin = -klim

kymax = klim

kstep = klim / 100.



# compute transfer function as a function of wavenumber difference

transff = array_transff_wavenumber(coords, klim, kstep, coordsys='xy')



# plot

plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2.,

           np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2.,

           transff.T, cmap=obspy_sequential)



plt.colorbar()

plt.clim(vmin=0., vmax=1.)

plt.xlim(kxmin, kxmax)

plt.ylim(kymin, kymax)

plt.show()