#! /usr/bin/env python2

# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>

# <markdowncell>

# ##UTCDateTime
# Documentation: [UTCDateTime](http://docs.obspy.org/packages/autogen/obspy.core.utcdatetime.UTCDateTime.html#obspy.core.utcdatetime.UTCDateTime)

# <codecell>

from obspy.core import UTCDateTime
t = UTCDateTime()
t1 = UTCDateTime("2014-03-14T10:15:00+01:00")
t2= UTCDateTime(2012, 9, 7, 12, 15, 0)
t3= UTCDateTime(1347020100.0)
print t.year, t.julday, t.timestamp, t.weekday
print t+86400

# <markdowncell>

# ##Načítání seismogramů

# <codecell>

from obspy.core import read
st = read('http://examples.obspy.org/RJOB_061005_072159.ehz.new')
print st
tr = st[0]  # assign first and only trace to new variable
print tr

# <codecell>

print tr.stats
tr.stats.station
tr.stats.gse2.datatype

# <codecell>

print tr.data
print tr.data[0:3]
print len(tr)

# <markdowncell>

# ## Arclink

# <codecell>

from obspy.core import UTCDateTime
from obspy.arclink import Client

t = UTCDateTime()
try:
    client = Client(host="arclink.ethz.ch", user="jiri.vackar@sed.ethz.ch")
    st = client.getWaveform('CH', 'ZUR', '', 'HH*', t-300, t)
except:
    print 'stahovani selhalo'
else:
    print st
    st.plot()

# <markdowncell>

# ## Filtrování
# The following script shows how to filter a seismogram. The example uses a zero-phase-shift low-pass filter with a corner frequency of 1 Hz using 2 corners. This is done in two runs forward and backward, so we end up with 4 corners de facto.
# 
# The available filters are:
# 
# * bandpass
# * bandstop
# * lowpass
# * highpass
# 
# Dokumentace: [obspy.core.stream.Stream.filter](http://http://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.filter.html)

# <codecell>

t = UTCDateTime("2013-12-12 0:59:18")
st = client.getWaveform('CH', 'LIENZ', '', 'HH*', t-10, t+60)
st.plot()
st_filt = st.copy()
st_filt.filter('lowpass', freq=1.0, corners=2, zerophase=True)
st_filt.plot()

# <markdowncell>

# ## Obálky seismogramů a vykreslování pomocí matplotlib
# Dokumentace: [Seismogram envelopes](http://docs.obspy.org/packages/autogen/obspy.signal.filter.envelope.html)
# [Matplotlib: plot](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot)

# <codecell>

import numpy as np
import matplotlib.pyplot as plt
from obspy.core import read
import obspy.signal

st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
npts = st[0].stats.npts
samprate = st[0].stats.sampling_rate

st.filter('bandpass', freqmin=1, freqmax=3, corners=2, zerophase=True) # Filtrujeme puvodni zaznam

# Envelope of filtered data
data_envelope = obspy.signal.filter.envelope(st_filt[0].data)

# The plotting, plain matplotlib
t = np.arange(0, npts / samprate, 1 / samprate)
plt.plot(t, st[0].data, 'k')
plt.plot(t, data_envelope, 'r.')
plt.title(st[0].stats.starttime)
plt.ylabel('Filtered Data w/ Envelope')
plt.xlabel('Time [s]')
plt.xlim(80, 90)
plt.show()

# <markdowncell>

# ## Spektrogram

# <codecell>

from obspy.core import read

st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
st.spectrogram(log=True, title='BW.RJOB ' + str(st[0].stats.starttime))

# <markdowncell>

# ## Nuly a póly, přenosová funkce přístroje

# <codecell>

from obspy.core import UTCDateTime
from obspy.arclink import Client
from obspy.signal import pazToFreqResp
import matplotlib.pyplot as plt

t = UTCDateTime("2013-12-12 10:34:21")
client = Client(host="arclink.ethz.ch", user="jiri.vackar@sed.ethz.ch")
st = client.getWaveform('CH', 'VANNI', '', 'HHZ', t, t+10)
tr = st[0]
tr2 = tr.copy()

# stahneme nuly a poly pres ArcLink
paz = client.getPAZ('CH', 'VANNI', '', 'HHZ', t)

# nakreslime frekvencni odezvu
h,f = pazToFreqResp(paz['poles'], paz['zeros'], paz['gain'], 0.002, 2**16, freq=True)
plt.figure()
plt.loglog(f, abs(h))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.show()

# provedeme instrumentalni korekci
tr2.simulate(paz_remove=paz)

# vykreslime
t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, 1 / tr.stats.sampling_rate)
plt.plot(t, tr.data, 'k')
plt2 = plt.twinx()
plt2.plot(t, tr2.data, 'r')
plt.show()