obspy中文教程（五）_obspy ppsd

 

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()