AzimPy 0.2.0a2

Estimate OBS sensor orientation

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 is inherited from obspy.clients.fdsn.Client, which allows us to search for earthquake catalog as a web client and compute Rayleigh-wave polarizations. This module also provides other classes and functions for statistical analysis of circular data and plotting the estimated azimuths with uncertainties.

How to install

[Recommended] Install AzimPy from PyPI in a new conda environment

$ conda create -n azimpy-test python=3.9 ipython 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 locally in the environment

$ conda create -n azimpy-test python=3.9 ipython
$ 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

Note that this installation will take time.

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

Then, set environental 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 via pip

(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/station directory.

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

• The supported format is only SAC, 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 Kuiper test.

Use case

• Sawaki, Y., Yamashita, Y., Ohyanagi, S., Garcia, E.S.M., Ito, A., Sugioka, H., Takahashi, T., Shinohara, M., & Ito, Y., in revision at Geophys. J. Int.

References

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