阿齐姆派

估计海底地震仪的水平方向

版权所有 (c) 2022 Yasunori Sawaki保留所有权利。

AzimPy是一个开源 Python 包，用于估计海底地震仪的水平方向。该软件包执行瑞利波偏振方法（例如，Stachnik+ 2012；Doran & Laske 2017）。OrientOBS继承自的主要类之一obspy.clients.fdsn.Client作为 Web 客户端搜索远震事件的地震目录，并计算每个事件-台站对的瑞利波极化。该软件包还提供了其他类（例如 、OrientSingle）OrientAnalysis和函数，用于圆形数据的统计分析和绘制具有不确定性的估计方位角。

使用条款

• 引用 Zenodo DOI 来获取AzimPy您发表研究成果或进行演示的特定版本。代表特定版本的 DOI 可能通过Zenodo 页面找到最新版本。
• 这个包正在开发中，所以欢迎任何错误报告和建议！

用例

• Sawaki, Y., Yamashita, Y., Ohyanagi, S., Garcia, ESM, Ito, A., Sugioka, H., Takahashi, T., Shinohara, M., & Ito, Y.，在Geophys 修订。J.诠释。

如何安装

[推荐]AzimPy在PyPI新conda环境中安装

$ conda create -n azimpy-test python=3.9 jupyter astropy "matplotlib>=3.5" "scipy>=1.4" pandas numpy tqdm
$ conda activate azimpy-test
(azimpy-test) $ conda install -c conda-forge "obspy>=1.3" 
(azimpy-test) $ python -m pip install AzimPy

• [替代方案]pip install本地环境中

$ conda create -n azimpy-test python=3.9 jupyter
$ conda activate azimpy-test
(azimpy-test) $ git clone -b v0.2.0 https://github.com/yasuit21/AzimPy.git
(azimpy-test) $ cd AzimPy
(azimpy-test) $ python -m pip install .

可选安装：rpy2

安装R

(azimpy-test) $ conda install r-essentials r-base r-circular

请注意，此安装需要时间。

然后，设置环境变量。

export R_HOME=/path/to/envs/azimpy-test/lib/R
export R_USER=/path/to/envs/azimpy-test/lib/python3.9/site-packages/rpy2

安装rpy2在PyPI

(azimpy-test) $ python -m pip install rpy2
(azimpy-test) $ python -c "import azimpy"

如果没有返回警告或错误，则安装已完成。

用法

计算瑞利波极化

import obspy as ob
from azimpy import OrientOBS

## Initialize web client
obs = OrientOBS(base_url='USGS', timezone=9)

## Query earthquake event catalog
obs.query_events(
    starttime=ob.UTCDateTime('20200401000000'),
    endtime=ob.UTCDateTime('20201001000000'),
    minmagnitude=6.0,
    maxdepth=100,
    orderby='time-asc',
)

## Compute Rayleigh-wave polarization for each event
## Raw SAC data should be located in '/path/to/datadir'
obs.find_stream(
    '/path/to/datadir',
    output_path='/path/to/output/stationA1',
    polezero_fpath='/path/to/polezero/hoge.paz',
    fileformat=f'*.*.%y%m%d%H%M.sac',
    freqmin=1./40, freqmax=1./20,
    max_workers=4,
    vel_surface=4.0,
    time_before_arrival=-20.0,
    time_after_arrival=600.0,
    distmin=5., distmax=120.,
)

然后，输出数据框将被腌制为目录stationA1_020_040.pickle下/path/to/output/stationA1。腌制的数据框可以通过pd.read_pickle().

执行循环统计并绘制图表

单站

该示例使用单个工作站stationA1。

import pandas as pd
from azimpy import OrientSingle, read_chtbl

## Init params
min_CC = 0.5
alpha_CI = 0.05  ## 100(1-a)% CI
bootstrap_iteration = 5000

## The output dataframe of orientations
df_orient = pd.read_pickle(
    '/path/to/output/stationA1/stationA1_020_040.pickle'
)

## Init OrientSingle for circular statistics
orientsingle_raw = OrientSingle(
    df_orient, 'stationA1', 
    if_selection=False,  # w/o bootstrap analysis
    min_CC=min_CC, weight_CC=True,
)
orientsingle_boot = OrientSingle(
    df_orient, 'stationA1', 
    if_selection=True,  # perform bootstrap analysis
    min_CC=min_CC, weight_CC=True, K=5.0,
    bootstrap_iteration=bootstrap_iteration, alpha_CI=alpha_CI
)

## Init a figure with subfigures
fig = plt.figure(figsize=[8,4])
subfigs = fig.subfigures(nrows=1, ncols=2).flatten()

## Plot for `orientsingle_raw`
orientsingle_raw.plot(
    polar=True, 
    fig=subfigs[0], in_parentheses='BB',
    add_cbar=True
)
subfigs[0].legend(loc=1, bbox_to_anchor=(1,1.15), fontsize='small')

## Plot for `orientsingle_boot`
orientsingle_boot.plot(
    fig=subfigs[1], in_parentheses='BB',
)
subfigs[1].legend(loc=1, bbox_to_anchor=(1,1.15), fontsize='small')

## Show or save the figure
fig.savefig()
plt.show()

多站

该示例使用名称为 的多个电台stationAX。

from azimpy import OrientAnalysis

stationList = ['stationA1','stationA2','stationA3','stationA4']

## Channeltable including above stations' info
df_chtbl = read_chtbl('/path/to/channeltable.txt')
df_chtbl = df_chtbl.query('comp.str.endswith("U")')

## Init OrientAnalysis for circular statistics
oa_raw = OrientAnalysis(
    if_selection=False,  # w/o bootstrap analysis
    df_chtbl=df_chtbl, 
    min_CC=min_CC, 
)
oa_boot = OrientAnalysis(
    if_selection=True,  # perform bootstrap analysis
    df_chtbl=df_chtbl, 
    min_CC=min_CC, alpha_CI=alpha_CI, 
    bootstrap_iteration=bootstrap_iteration, 
)

for stationName in stationList:
    period = df_chtbl.at[stationName,'period']
    df_orient = pd.read_pickle(
        f'/path/to/output/{stationName}/{stationName}_020_040.pickle'
    )
    ## Add the dataframe in `oa_raw`
    ## This is actually passed to `OrientSingle`
    oa_raw.add_station(
        df_orient, stationName, 
        period=period
    )
    ## Add the dataframe in `oa_boot`
    oa_boot.add_station(
        df_orient, stationName, 
        period=period
    )

• 绘制结果使用matplotlib.pyplot
  • 没有自举重采样的原始结果

  fig = oa_raw.plot()
  

  • 自举分析的结果

  fig = oa_boot.plot()
  

笔记

• SAC仅支持格式，但您可以使用其他一些格式。
• 观察到的数据文件必须位于一个目录中，该目录OrientOBS.find_stream()将尝试搜索必要的输入文件。目前此软件包中没有网站和存储库中的波形数据。
• 作者已经在Linux环境（CentOS 7和WSL Ubuntu 20.04）中测试过这个包，所以安装在Windows.
• rpy2circular是在语言中运行的可选包装器R，它执行 Kuiper 测试。

参考

• Stachnik, JC, Sheehan, AF, Zietlow, DW, 等人，2012，通过瑞利波偏振确定新西兰海底地震仪方向。地震。水库。莱特。, 83, 704–713。https://doi.org/10.1785/0220110128
• Doran, AK & Laske, G.，2017，通过自动瑞利波到达角测量的海底防震仪仪器方向。公牛。地震。社会党。是。, 107, 691–708。https://doi.org/10.1785/0120160165
• Takagi, R.、Uchida, N.、Nakayama, T. 等人，2019 年，S-net 电缆海底传感器方向的估计。地震。水库。莱特。, 90, 2175–2187。https://doi.org/10.1785/0220190093
• 最新的概念 DOI AzimPy：10.5281/zenodo.6972713

致谢

该软件包ObsPy v1.3.0用于FDSN Web 客户端服务和处理地震图。

执照

该项目在 MIT License 下获得许可，详情请参阅LICENSE。