Commit 3a3cf435 authored by Mijian Xu's avatar Mijian Xu 😷
Browse files

ver 1.1.15

parent 15707865
import math
import numpy as np
def sind(deg):
rad = math.radians(deg)
return math.sin(rad)
......
......@@ -92,18 +92,19 @@ class eq(object):
this_st = self.st.copy()
this_st.filter('bandpass', freqmin=0.03, freqmax=0.5)
this_st.trim(this_st[0].stats.starttime+p_arr-time_b, this_st[0].stats.starttime+p_arr+time_e)
bazs = np.arange(bazi-offset, bazi+offset)
bazs = np.mod(np.arange(bazi-offset, bazi+offset), 360)
ampt = np.zeros_like(bazs)
for i, b in enumerate(bazs):
t, _ = seispy.geo.rotateSeisENtoTR(this_st[0].data, this_st[1].data, b)
ampt[i] = seispy.geo.rssq(t)
idx = seispy.geo.extrema(ampt, opt='min')
if len(idx) == 0:
return np.nan
corr_baz = np.nan
elif len(idx) > 1:
return None
corr_baz = None
else:
return bazs[seispy.geo.extrema(ampt, opt='min')[0]] - bazi
corr_baz = bazs[idx[0]] - bazi
return corr_baz, ampt
def fix_channel_name(self):
if self.st.select(channel='??1') and self.st.select(channel='??Z') and hasattr(self.st.select(channel='*1')[0].stats.sac, 'cmpaz'):
......@@ -154,10 +155,18 @@ class eq(object):
st_noise = self.trim(length, 0, phase=phase, isreturn=True)
st_signal = self.trim(0, length, phase=phase, isreturn=True)
try:
snr_R = seispy.geo.snr(st_signal[1].data, st_noise[1].data)
snr_E = seispy.geo.snr(st_signal[0].data, st_noise[0].data)
except IndexError:
snr_R = 0
return snr_R
snr_E = 0
try:
snr_N = seispy.geo.snr(st_signal[1].data, st_noise[1].data)
except IndexError:
snr_N = 0
try:
snr_Z = seispy.geo.snr(st_signal[2].data, st_noise[2].data)
except IndexError:
snr_Z = 0
return snr_E, snr_N, snr_Z
def get_time_offset(self, event_time):
if not isinstance(event_time, obspy.core.utcdatetime.UTCDateTime):
......
......@@ -2,36 +2,43 @@ import numpy as np
from numpy import pi, mod
from seispy import distaz
from scipy import interpolate
import math
def sind(deg):
rad = np.radians(deg)
return np.sin(rad)
def cosd(deg):
rad = np.radians(deg)
return np.cos(rad)
def tand(deg):
rad = np.radians(deg)
return np.tan(rad)
def cotd(deg):
rad = np.radians(deg)
return np.cos(rad) / np.sin(rad)
def asind(x):
rad = np.arcsin(x)
return np.degrees(rad)
def acosd(x):
rad = np.arccos(x)
return np.degrees(rad)
def atand(x):
rad = np.arctan(x)
return np.degrees(rad)
def km2deg(km):
radius = 6371
circum = 2*pi*radius
......@@ -39,6 +46,7 @@ def km2deg(km):
deg = km / conv
return deg
def deg2km(deg):
radius = 6371
circum = 2*pi*radius
......@@ -46,22 +54,27 @@ def deg2km(deg):
km = deg * conv
return km
def rad2deg(rad):
deg = rad*(360/(2*pi))
return deg
def skm2sdeg(skm):
sdeg = skm * deg2km(1)
return sdeg
def sdeg2skm(sdeg):
skm = sdeg / deg2km(1)
return skm
def srad2skm(srad):
sdeg = srad * ((2*pi)/360)
return sdeg / deg2km(1)
def skm2srad(skm):
sdeg = skm * deg2km(1)
return rad2deg(sdeg)
......@@ -115,35 +128,41 @@ def rotateSeisENtoTR(E, N, BAZ):
T = E*cosd(angle) - N*sind(angle)
return T, R
def rssq(x):
return np.sqrt(np.sum(x**2)/len(x))
def snr(x, y):
spow = rssq(x)**2
npow = rssq(y)**2
return 10 * np.log10(spow / npow)
if npow == 0:
return np.nan
else:
return 10 * np.log10(spow / npow)
def latlon_from(lat1, lon1, azimuth, gcarc_dist):
lat2 = asind ((sind (lat1) * cosd (gcarc_dist)) + (cosd (lat1) * sind (gcarc_dist) * cosd (azimuth)))
lat2 = asind((sind(lat1) * cosd(gcarc_dist)) + (cosd(lat1) * sind(gcarc_dist) * cosd(azimuth)))
if isinstance(gcarc_dist, np.ndarray):
lon2 = np.zeros_like(lat2)
for n in range(len(gcarc_dist)):
if cosd (gcarc_dist[n]) >= (cosd (90 - lat1) * cosd (90 - lat2[n])):
lon2[n] = lon1 + asind (sind (gcarc_dist[n]) * sind (azimuth) / cosd (lat2[n]))
if cosd(gcarc_dist[n]) >= (cosd(90 - lat1) * cosd(90 - lat2[n])):
lon2[n] = lon1 + asind(sind(gcarc_dist[n]) * sind(azimuth) / cosd(lat2[n]))
else:
lon2[n] = lon1 + asind (sind (gcarc_dist[n]) * sind (azimuth) / cosd (lat2[n])) + 180
lon2[n] = lon1 + asind(sind(gcarc_dist[n]) * sind(azimuth) / cosd(lat2[n])) + 180
elif isinstance(azimuth, np.ndarray):
lon2 = np.zeros_like(lat2)
for n in range(len(azimuth)):
if cosd (gcarc_dist) >= (cosd (90 - lat1) * cosd (90 - lat2[n])):
lon2[n] = lon1 + asind (sind (gcarc_dist) * sind (azimuth[n]) / cosd (lat2[n]))
if cosd(gcarc_dist) >= (cosd(90 - lat1) * cosd(90 - lat2[n])):
lon2[n] = lon1 + asind(sind(gcarc_dist) * sind(azimuth[n]) / cosd(lat2[n]))
else:
lon2[n] = lon1 + asind (sind (gcarc_dist) * sind (azimuth[n]) / cosd (lat2[n])) + 180
lon2[n] = lon1 + asind(sind(gcarc_dist) * sind(azimuth[n]) / cosd(lat2[n])) + 180
else:
if ( cosd (gcarc_dist) >= (cosd (90 - lat1) * cosd (90 - lat2))):
lon2 = lon1 + asind (sind (gcarc_dist) * sind (azimuth) / cosd (lat2))
if (cosd(gcarc_dist) >= (cosd(90 - lat1) * cosd(90 - lat2))):
lon2 = lon1 + asind(sind(gcarc_dist) * sind(azimuth) / cosd(lat2))
else:
lon2 = lon1 + asind (sind (gcarc_dist) * sind (azimuth) / cosd (lat2)) + 180
lon2 = lon1 + asind(sind(gcarc_dist) * sind(azimuth) / cosd(lat2)) + 180
return lat2, lon2
......@@ -165,6 +184,7 @@ def extrema(x, opt='max'):
raise ImportError('Wrong Options!!!')
return idx
def slantstack(seis, timeaxis, rayp_range, tau_range, ref_dis, dis):
EV_num = seis.shape[0]
tmp = np.zeros([EV_num, tau_range.shape[0]])
......
import numpy as np
def mccc(seis, dt, twin=0):
data = np.zeros([len(seis), seis[0].data.shape[0]])
for i in range(len(seis)):
......
import sys
import os
import random
import matplotlib
import argparse
# matplotlib.use("Qt5Agg")
from PyQt5 import QtCore
from PyQt5.QtCore import Qt
from PyQt5.QtGui import QIcon, QKeySequence
from PyQt5.QtWidgets import QApplication, QMainWindow, QVBoxLayout, QSizePolicy, QWidget, QDesktopWidget, \
QPushButton, QHBoxLayout, QLineEdit, QFileDialog, QAction, QShortcut
from numpy import arange, sin, pi
from PyQt5.QtWidgets import QApplication, QMainWindow, QVBoxLayout, \
QSizePolicy, QWidget, QDesktopWidget, \
QPushButton, QHBoxLayout, QFileDialog, \
QAction, QShortcut
from os.path import exists, dirname, join
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT as NavigationToolbar
from matplotlib.figure import Figure
import matplotlib.pyplot as plt
from seispy.pickfigure import RFFigure
class MyMplCanvas(FigureCanvas):
def __init__(self, parent=None, rfpath='', width=21, height=11, dpi=100):
......@@ -63,7 +58,7 @@ class MatplotlibWidget(QMainWindow):
self._define_global_shortcuts()
self.setWindowTitle('PickRF')
self.setWindowIcon(QIcon(join(dirname(__file__), 'data', 'seispy.png')))
def on_click(self, event):
self.mpl.rffig.onclick(event)
self.mpl.draw()
......@@ -71,7 +66,7 @@ class MatplotlibWidget(QMainWindow):
def previous_connect(self):
self.mpl.rffig.butprevious()
self.mpl.draw()
def next_connect(self):
self.mpl.rffig.butnext()
self.mpl.draw()
......@@ -79,11 +74,11 @@ class MatplotlibWidget(QMainWindow):
def enlarge(self):
self.mpl.rffig.enlarge()
self.mpl.draw()
def reduce(self):
self.mpl.rffig.reduce()
self.mpl.draw()
def finish(self):
self.mpl.rffig.finish()
QApplication.quit()
......@@ -92,8 +87,8 @@ class MatplotlibWidget(QMainWindow):
self.mpl.rffig.plot()
def plot_save(self):
fileName_choose, filetype = QFileDialog.getSaveFileName(self,
"Save the figure",
fileName_choose, filetype = QFileDialog.getSaveFileName(self,
"Save the figure",
os.path.join(os.getcwd(), self.mpl.rffig.staname + 'RT_bazorder'),
"PDF Files (*.pdf);;Images (*.png);;All Files (*)")
......@@ -116,8 +111,8 @@ class MatplotlibWidget(QMainWindow):
self.setGeometry(0, 0, frame_width, frame_height)
self.move((screen_width / 2) - (self.frameSize().width() / 2),
(screen_height / 2) - (self.frameSize().height() / 2))
(screen_height / 2) - (self.frameSize().height() / 2))
def _define_global_shortcuts(self):
self.key_c = QShortcut(QKeySequence('c'), self)
self.key_c.activated.connect(self.next_connect)
......@@ -167,7 +162,8 @@ def main():
app = QApplication(sys.argv)
ui = MatplotlibWidget(rfpath)
ui.show()
sys.exit(app.exec_())
sys.exit(app.exec_())
if __name__ == '__main__':
main()
......@@ -26,8 +26,9 @@ def init_figure():
return h, axr, axt, axb, axr_sum, axt_sum
def read_process_data(lst):
def read_process_data(lst, resamp_dt=0.1):
stadata = SACStation(lst)
stadata.resample(resamp_dt)
idx = np.argsort(stadata.bazi)
stadata.event = stadata.event[idx]
stadata.bazi = stadata.bazi[idx]
......
......@@ -4,19 +4,21 @@ import obspy
from obspy.io.sac import SACTrace
import re
from seispy.io import wsfetch
from os.path import dirname, join, expanduser, exists
from os.path import join
from seispy.para import para
import seispy
from seispy.eq import eq
from seispy.setuplog import setuplog
import glob
import numpy as np
from datetime import timedelta, datetime
from datetime import timedelta
import pandas as pd
from obspy.taup import TauPyModel
import configparser
import argparse
import sys
from tkinter import messagebox
import matplotlib.pyplot as plt
def datestr2regex(datestr):
......@@ -164,6 +166,13 @@ def CfgModify(cfg_file, session, key, value):
cf.write(open(cfg_file, 'w'))
def _plotampt(x, y, ampt):
xx, yy = np.meshgrid(x, y)
f = plt.figure(figsize=(8, 8))
plt.pcolor(xx, yy, ampt)
return f
class RF(object):
def __init__(self, cfg_file=None, log=None):
if cfg_file is None:
......@@ -299,15 +308,32 @@ class RF(object):
def baz_correct(self, time_b=10, time_e=20, offset=90):
self.logger.RFlog.info('correct back-azimuth')
shift_all = []
# ampt_all = []
for i, row in self.eqs.iterrows():
shift_all.append(row['data'].search_baz(row['bazi'], time_b=time_b, time_e=time_e, offset=offset))
curr_baz, _ = row['data'].search_baz(row['bazi'], time_b=time_b, time_e=time_e, offset=offset)
shift_all.append(curr_baz)
# ampt_all.append(ampt)
if None in shift_all:
self.logger.RFlog.error('Range of searching bazi is too small.')
sys.exit(1)
shift_all = np.array(shift_all)
# ampt_all = np.array(ampt_all)
baz_shift = np.mean(shift_all[np.where(np.logical_not(np.isnan(shift_all)))])
self.logger.RFlog.info('Average {:.4f} deg offset in back-azimuth'.format(baz_shift))
self.eqs['bazi'] = self.eqs['bazi'] + baz_shift
# self._baz_confirm(offset, ampt_all)
self.logger.RFlog.info('Average {:.1f} deg offset in back-azimuth'.format(baz_shift))
self.eqs['bazi'] = np.mod(self.eqs['bazi'] + baz_shift, 360)
# def _baz_confirm(self, offset, ampt_all):
# y = np.arange(self.eqs.shape[0])
# x = np.linspace(-offset, offset, ampt_all.shape[1])
# fig = _plotampt(x, y, ampt_all)
# if messagebox.askyesno('Back-azimuth correction',
# 'Would you like to keep this correction?'):
# fig.close()
# else:
# self.logger.RFlog.error('Manual interruption.')
# fig.close()
# sys.exit(1)
def rotate(self, method='NE->RT'):
self.logger.RFlog.info('Rotate {0} phase {1}'.format(self.para.phase, method))
......@@ -316,7 +342,7 @@ class RF(object):
try:
row['data'].rotate(row['bazi'], method=method, phase=self.para.phase)
except Exception as e:
self.logger.RFlog.error('{}'.format(e))
self.logger.RFlog.error('{}: {}'.format(row['data'].datestr, e))
drop_idx.append(i)
self.eqs.drop(drop_idx, inplace=True)
......@@ -326,8 +352,9 @@ class RF(object):
self.logger.RFlog.info('Reject data record with SNR less than {0}'.format(self.para.noisegate))
drop_lst = []
for i, row in self.eqs.iterrows():
snr_R = row['data'].snr(length=length)
if snr_R < self.para.noisegate:
snr_E, snr_N, snr_Z = row['data'].snr(length=length)
if (np.nan in (snr_E, snr_N, snr_Z) or snr_E < self.para.noisegate
or snr_N < self.para.noisegate or snr_Z < self.para.noisegate):
drop_lst.append(i)
self.eqs.drop(drop_lst, inplace=True)
self.logger.RFlog.info('{0} events left after SNR calculation'.format(self.eqs.shape[0]))
......@@ -381,6 +408,7 @@ class RF(object):
evla=row['evla'], evlo=row['evlo'], evdp=row['evdp'], baz=row['bazi'],
mag=row['mag'], gcarc=row['dis'], gauss=self.para.gauss, only_r=self.para.only_r)
good_lst.append(i)
self.logger.RFlog.info('{} PRFs are saved.'.format(len(good_lst)))
self.eqs = self.eqs.loc[good_lst]
......@@ -388,18 +416,22 @@ def prf():
parser = argparse.ArgumentParser(description="Calculating RFs for single station")
parser.add_argument('cfg_file', type=str, help='Path to RF configure file')
parser.add_argument('-l', help="use local catalog. Default is false", dest='islocal', action='store_true')
parser.add_argument('-r', help='Reverse components: EN, E or N', default=None, type=str)
parser.add_argument('-b', help='Correct back-azimuth via minimal energy of T component', default=None, type=str)
parser.add_argument('-r', help='Reverse components: EN, E or N', dest='comp', default=None, type=str)
parser.add_argument('-s', help='Switch East and North components', dest='isswitch', action='store_true')
parser.add_argument('-b', help='Correct back-azimuth with minimal '
'energy of T component. "baz" is specified '
'as half-searching range. Default value is 90 deg',
dest='baz', nargs='?', const=90, type=float)
arg = parser.parse_args()
if arg.r is not None:
arg.r = arg.r.upper()
if arg.r == 'NE' or arg.r == 'EN':
if arg.comp is not None:
arg.comp = arg.comp.upper()
if arg.comp == 'NE' or arg.comp == 'EN':
reverseE = True
reverseN = True
elif arg.r == 'E':
elif arg.comp == 'E':
reverseE = True
reverseN = False
elif arg.r == 'N':
elif arg.comp == 'N':
reverseE = False
reverseN = True
else:
......@@ -407,21 +439,21 @@ def prf():
else:
reverseN = False
reverseE = False
if arg.b is not None:
try:
baz_corr = [float(val) for val in arg.b.split('/')]
except:
raise ValueError('Format of baz correction must be as 10/20/45')
# if arg.baz is not None:
# try:
# baz_corr = [float(val) for val in arg.b.split('/')]
# except:
# raise ValueError('Format of baz correction must be as 10/20/45')
pjt = RF(cfg_file=arg.cfg_file)
pjt.load_stainfo()
pjt.search_eq(local=arg.islocal)
pjt.match_eq(reverseN=reverseN, reverseE=reverseE)
pjt.match_eq(switchEN=arg.isswitch, reverseN=reverseN, reverseE=reverseE)
pjt.detrend()
pjt.filter()
pjt.cal_phase()
pjt.drop_eq_snr()
if arg.b is not None:
pjt.baz_correct(baz_corr[0], baz_corr[1], baz_corr[2])
if arg.baz is not None:
pjt.baz_correct(offset=arg.baz)
pjt.trim()
pjt.rotate()
pjt.deconv()
......@@ -437,6 +469,7 @@ def setpar():
arg = parser.parse_args()
CfgModify(arg.cfg_file, arg.session, arg.key, arg.value)
if __name__ == '__main__':
# get_events_test()
prf()
......
from obspy.io.sac.sactrace import SACTrace
import numpy as np
from scipy.interpolate import interp1d, interpn
from scipy.signal import resample
from os.path import dirname, join, exists
from seispy.geo import skm2srad, sdeg2skm, rad2deg, latlon_from, \
asind, tand, srad2skm, km2deg
......@@ -17,6 +18,7 @@ class SACStation(object):
:param evt_lst: event list in RF data dir. Column as date string, phase, evt lat, evt lon, evt dep,
distance, back-azimuth, ray parameter, magnitude, gauss factor.
"""
self.only_r = only_r
data_path = dirname(evt_lst)
dtype = {'names': ('evt', 'phase', 'evlat', 'evlon', 'evdep', 'dis', 'bazi', 'rayp', 'mag', 'f0'),
'formats': ('U20', 'U20', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4', 'f4')}
......@@ -45,6 +47,14 @@ class SACStation(object):
sac = SACTrace.read(join(data_path, evt + '_' + ph + '_R.sac'))
self.datar[_i] = sac.data
def resample(self, dt):
npts = int(self.RFlength * (self.sampling / dt)) + 1
self.datar = resample(self.datar, npts, axis=1)
if not self.only_r:
self.datat = resample(self.datat, npts, axis=1)
self.sampling = dt
self.RFlength = npts
class DepModel(object):
def __init__(self, YAxisRange, velmod='iasp91'):
......
......@@ -41,21 +41,20 @@ def smooth(x, half_len=5, window='flat'):
if x.size < window_len:
raise ValueError("Input vector needs to be bigger than window size.")
if window_len<3:
if window_len < 3:
return x
if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
raise ValueError("Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'")
if window not in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
raise ValueError("Window is on of 'flat', 'hanning', 'hamming',"
"'bartlett', 'blackman'")
s = np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
#print(len(s))
if window == 'flat': #moving average
w = np.ones(window_len,'d')
if window == 'flat':
w = np.ones(window_len, 'd') # moving average
else:
w = eval('np.'+window+'(window_len)')
y = np.convolve(w/w.sum(), s, mode='valid')
return y[half_len:-half_len]
return y[half_len:-half_len]
def whiten(data, Nfft, delta, f1, f2, f3, f4):
......
......@@ -2,7 +2,7 @@
from setuptools import find_packages, setup
packages = find_packages()
VERSION = "1.1.14"
VERSION = "1.1.15"
setup(name='seispy',
version=VERSION,
author='Mijian Xu',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment