### obspy.imaging

## obspy.imaging.beachball: Beachball plot

from obspy.imaging.beachball import Beachball

fm2=[318, 88, -169] #strike,dip,rake
Beachball(fm2)
fm=[-0.33, -2.23, 2.55, 0.83, -0.32, 0.32] # m11,m22,m33,m12,m13,m23
Beachball(fm,facecolor='r')



### obspy.taup

## obspy.taup.taup: Travel time calculation tool 

from obspy.taup.taup import travelTimePlot, getTravelTimes
travelTimePlot(max_degree=50,phases=['P','S','PP','PS','SS'])
tt=getTravelTimes(delta=40.,depth=100.,model='iasp91')
print tt[0]
print tt[1]



### obspy.signal

from obspy.core import read
st = read("./event.gse2")
tr = st[0]

print tr.stats



tr.plot()
tr2=tr.copy()
tr2.plot(color='blue',tick_format='%I:%M:%S %p', starttime=st[0].stats.starttime+20, endtime=st[0].stats.endtime-30,outfile='/home/mess/trplot.png')



### Trigger/Picker: obspy.signal.trigger

# Metody z obspy.signal.trigger (mnoho): skusime klasicky STA/LTA - classicSTALTA(array,nsta,nlta)

from obspy.signal.trigger import classicSTALTA
t_sta=1.
t_lta=10.
samp=tr.stats.sampling_rate
cft=classicSTALTA(tr.data,int(t_sta*samp),int(t_lta*samp))
print cft



from obspy.signal.trigger import plotTrigger
thr_on=2. #threshold for switching trigger on
thr_off=.5 #threshold for switching trigger off
plotTrigger(tr,cft,thr_on,thr_off)



from obspy.signal.trigger import recSTALTA, plotTrigger
cft2=recSTALTA(tr.data,int(t_sta*samp),int(t_lta*samp))
plotTrigger(tr,cft2,thr_on,thr_off)



from obspy.signal.trigger import triggerOnset
trg=triggerOnset(cft,thr_on,thr_off)
trg2=triggerOnset(cft2,thr_on,thr_off)
print trg/samp
print trg2/samp



## Coincidence trigger - ziska zoznam prekryvajucich sa triggerov na sieti stanic

from obspy.signal import coincidenceTrigger
st2=st.copy() #kopia
st2=st.copy()+st2 #vyrobime dalsi stream
st2[1].stats.station="RNO2" #premenujeme stanicu na druhom streame, ma to detekciu

coinc_sum=2 #minimalny pocet prekryvov
trig=coincidenceTrigger('classicSTALTA',thr_on,thr_off,st2,coinc_sum,sta=t_sta,lta=t_lta)

print trig



## Baer Picker

from obspy.signal.trigger import pkBaer
tdownmax=25
min_nr_samples=100
sigma=12.
preset=100
p_pick, phase_info = pkBaer(tr.data,samp,tdownmax,min_nr_samples,thr_on,sigma,preset,6)
print p_pick/samp



## obspy.signal.cross_correlation

from obspy.core import read
st1=read("/home/mess/Downloads/BW.UH1..EHZ.D.2010.147.a.slist.gz")
st2=read("/home/mess/Downloads/BW.UH1..EHZ.D.2010.147.b.slist.gz")
#st1 = read("http://examples.obspy.org/BW.UH1..EHZ.D.2010.147.a.slist.gz")
#st2 = read("http://examples.obspy.org/BW.UH1..EHZ.D.2010.147.b.slist.gz")
tr1=st1[0]
tr2=st2[0]
print tr1,tr2
tr1.plot(),tr2.plot()



from obspy.signal.cross_correlation import xcorr, xcorr_max
tr1cp=tr1.copy()
tr2cp=tr2.copy()
#shift_len has to be selected carefully, make it a bit bigger than the highest signal shifts that can ever occur
index, maxim, xcfct = xcorr(tr1.data,tr2.data,shift_len=500,full_xcorr=True)

print index/samp,maxim
print xcfct



from obspy.core import UTCDateTime
from obspy.signal.cross_correlation import xcorrPickCorrection
from obspy.signal.trigger import recSTALTA
from obspy.signal.trigger import triggerOnset,plotTrigger

#estimate pick1:
t_sta=.1
t_lta=1.
samp=tr1cp.stats.sampling_rate
thr_on=4.
thr_off=0.5
cft1=recSTALTA(tr1cp.data,int(t_sta*samp),int(t_lta*samp))
trg1=triggerOnset(cft1,thr_on,thr_off)
pick1=trg1[0][0]/samp
#plotTrigger(tr1cp,cft1,thr_on,thr_off)

t1=tr1.stats.starttime+pick1
t2=tr2.stats.starttime+pick1

dt,coeff = xcorrPickCorrection(t1,tr1,t2,tr2,t_before=0.05,t_after=0.2,cc_maxlag=0.1,plot=True)

print dt,coeff



## PPSD = Probabilistic Power Spectral Densities

# Class to compile probabilistic power spectral densities for one combination of network/station/location/channel/sampling_rate

#nacitat data
from obspy.core import read
from obspy.xseed import Parser
from obspy.signal import PPSD
st = read("/home/mess/Downloads/BW.KW1..EHZ.D.2011.037")
print st



#poles&zeroes
parser = Parser("/home/mess/Downloads/dataless.seed.BW_KW1")
paz=parser.getPAZ("BW.KW1..EHZ")
print paz



#tr2=st[0]
st2=st.copy()
tr2=st2[0]
#st2=st[0:1] #berieme len prvy trace, dlho to trva
#st2.decimate(factor=2) #downsample na polovicnu frekvenciu
#print st2

#inicializuje ppsd
ppsd=PPSD(tr2.stats,paz)
#pridat data do ppsd:
ppsd.add(st2)
print ppsd.times
#rozdeli data do 1hodinovych segmentov, preprocessing(demean,taper..)



print len(ppsd.times)
#a dalsie data 
st3=read("/home/mess/Downloads/BW.KW1..EHZ.D.2011.038")
ppsd.add(st3)

print len(ppsd.times)



#plot
ppsd.plot()



### obspy.signal.tf_misfit

## Continuous wavelet transform

from obspy.signal.tf_misfit import cwt
st3=read()#"./G.SCZ..BHE.sac")
tr3=st3[0]
#help(cwt)
dt=tr3.stats.delta
fmin=0.05
fmax=50
w0=6
tr3.plot()
scalogram=cwt(tr3.data,dt,w0,fmin,fmax)



#vykreslenie pomocou matplotlib...
import numpy as np
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
t0=0.
tmax=dt*tr3.stats.npts
ax.imshow(np.abs(scalogram)[-1::-1],extent=[t0,tmax,fmin,fmax],
          interpolation='nearest')
ax.set_xlabel('time')
ax.set_ylabel('frequency')
ax.set_yscale('log')
plt.show



## TF reprezentacia

#vykreslit TF reprezentaciu
from obspy.signal.tf_misfit import plotTfr
plotTfr(tr3.data,dt,fmin=50,fmax=.1)



tr3_filt=tr3.copy()
tr3_filt.filter('highpass',freq=1.,zerophase=True) #odfiltrujeme dlhovlnny signal
plotTfr(tr3_filt.data,dt,fmin=50.,fmax=0.5)



## TF Misfits

from obspy.signal.tf_misfit import plotTfMisfits,em,pm
plotTfMisfits(tr3.data,tr3_filt.data,fmin=50.,fmax=0.01)
print em(tr3.data,tr3_filt.data,fmin=50.,fmax=0.01)
print pm(tr3.data,tr3_filt.data,fmin=50.,fmax=0.01)



tr3_filt2=tr3.copy()
tr3_filt2.filter('lowpass',freq=1.,zerophase=False) #vnesieme fazovy posun na dlhovlnnejsej zlozke
plotTfMisfits(tr3.data,tr3_filt2.data,fmin=50.,fmax=0.01)
print em(tr3.data,tr3_filt2.data,fmin=50.,fmax=0.01)
print pm(tr3.data,tr3_filt2.data,fmin=50.,fmax=0.01)



## obspy.signal.invsim.estimateMagnitude

from obspy.core import read, UTCDateTime
from obspy.core.util.geodetics import gps2DistAzimuth
from obspy.xseed import Parser
from obspy.signal.invsim import estimateMagnitude

# First read the data. Note that you pass a URL.
st5 = read("http://localhost/data/Advanced%20ObsPy%20Exercise/LKBD.MSEED")


parser = Parser("http://localhost/data/Advanced%20ObsPy%20Exercise/LKBD.dataless")
paz = parser.getPAZ("CH.LKBD..EHZ")

# Get the maximum amplitude of the North component
trn = st5.select(component="N")[0]
amplmax_n = max(trn.data)
amplmin_n = min(trn.data)
# Same for the East component
tre = st5.select(component="E")[0]
amplmax_e = max(tre.data)
amplmin_e = min(tre.data)

trz = st5.select(component="z")[0]
amplmax_z = max(trz.data)
amplmin_z = min(trz.data)

ampl = [amplmax_n-amplmin_n, amplmax_e-amplmin_e,amplmax_z-amplmin_z]
samp=st5[0].stats.sampling_rate

sta_lat = 46.38703
sta_lon = 7.62714
event_lat = 46.218
event_lon = 7.706
event_depth=5.2

epi_dist, az, baz= gps2DistAzimuth(event_lat, event_lon, sta_lat, sta_lon)
epi_dist = epi_dist / 1000

epi_dist= sqrt(epi_dist**2+event_depth**2)

#help(estimateMagnitude)

#st5.plot(starttime=st5[0].stats.starttime+500,endtime=st5[0].stats.starttime+550,number_of_ticks=10,tick_format='%S')
print estimateMagnitude([paz,paz,paz], ampl, h_dist=epi_dist,timespan=0.5)