AzimPy 0.2.0

Description:

AzimPy 0.2.0

AzimPy
Estimate horizontal orientation of ocean-bottom seismograph








Copyright (c) 2022 Yasunori Sawaki All rights reserved.
AzimPy is an open-source Python package for estimating the horizontal orientation of ocean-bottom seismographs.
This package performs the Rayleigh-wave polarization method (e.g., Stachnik+ 2012; Doran & Laske 2017).
One of the main classes OrientOBS, inherited from obspy.clients.fdsn.Client, searches an earthquake catalog for teleseismic events as a web client and computes Rayleigh-wave polarizations for each event–station pair.
This package also provides other classes (e.g., OrientSingle, OrientAnalysis) and functions for statistical analysis of circular data and plotting the estimated azimuths with uncertainties.
Terms of use

Cite an article below (Sawaki et al., in press) and a Zenodo DOI for the specific version of AzimPy when you publish your reseach or make a presentation. The DOI representing the specific version is probably found through the Zenodo page for the latest version.
This package is under development, so any bug reports and suggestions are welcome!

Use cases

Sawaki, Y., Yamashita, Y., Ohyanagi, S., Garcia, E.S.M., Ito, A., Sugioka, H., Takahashi, T., Shinohara, M., & Ito, Y., Seafloor Depth Controls Seismograph Orientation Uncertainty, Geophys. J. Int., in press, https://doi.org/10.1093/gji/ggac397

How to install
[Recommended] Install AzimPy from PyPI in a new conda environment
$ 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


[Alternative] pip install locally in the environment

$ 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 .

Optional installation : rpy2
Installation of R
(azimpy-test) $ conda install r-essentials r-base r-circular

Note that this installation will take time.
Then, set environent variables.
export R_HOME=/path/to/envs/azimpy-test/lib/R
export R_USER=/path/to/envs/azimpy-test/lib/python3.9/site-packages/rpy2

Installation of rpy2 in PyPI
(azimpy-test) $ python -m pip install rpy2
(azimpy-test) $ python -c "import azimpy"

If no warning or error is returned, the installation has been completed.
Usage
Compute Rayleigh-wave polarization
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.,
)

Then, the output dataframe will be pickled as stationA1_020_040.pickle under /path/to/output/stationA1 directory. The pickled dataframe can be loaded by pd.read_pickle().
Perform circular statistics and make a plot
Single station
The example uses a single station 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()


Multiple stations
The example uses multiple stations whose names are 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
)


Plot the results using matplotlib.pyplot

Original results w/o bootstrap resampling

fig = oa_raw.plot()


Results of bootstrap analysis

fig = oa_boot.plot()



Note

SAC format is only supported, but you may use some other formats.
The observed data files must be located in one directory, where OrientOBS.find_stream() will try to search for necessary input files. No waveform data in websites and repository are available in this package at this moment.
The author has tested this package in Linux environments (CentOS 7 and WSL Ubuntu 20.04), so it might be incompatible when installed in Windows.
rpy2 is an optional wrapper to run circular in R language, which performs the Kuiper test.

References

Sawaki, Y., Yamashita, Y., Ohyanagi, S., Garcia, E.S.M., Ito, A., Sugioka, H., Takahashi, T., Shinohara, M., & Ito, Y., Seafloor Depth Controls Seismograph Orientation Uncertainty, Geophys. J. Int., in press, https://doi.org/10.1093/gji/ggac397
Stachnik, J.C., Sheehan, A.F., Zietlow, D.W., et al., 2012, Determination of New Zealand ocean bottom seismometer orientation via Rayleigh-wave polarization. Seismol. Res. Lett., 83, 704–713. https://doi.org/10.1785/0220110128
Doran, A.K. & Laske, G., 2017, Ocean‐bottom deismometer instrument orientations via automated Rayleigh‐wave arrival‐angle measurements. Bull. Seismol. Soc. Am., 107, 691–708. https://doi.org/10.1785/0120160165
Takagi, R., Uchida, N., Nakayama, T., et al., 2019, Estimation of the orientations of the S‐net cabled ocean‐bottom sensors. Seismol. Res. Lett., 90, 2175–2187. https://doi.org/10.1785/0220190093
Concept DOI for the latest AzimPy: 10.5281/zenodo.6972713

Acknowledgments
This package makes use of ObsPy v1.3.0 for FDSN web client services and processing seismograms.
License
This project is licensed under the MIT License, see the LICENSE for details.