# -*- coding: utf-8 -*-
# 2009-03-12 Moritz

import obspy
from obspy.imaging.spectrogram import spectroGram
import matplotlib.pyplot as plt
import numpy as np
plt.close('all')

# Read 2003 Hokkaido earthquake directly into one trace
# and zoom into a certain part of the 2003 Hokkaido earthquake
tr = obspy.read("data/GR.FUR..LHN.D.2003.268")[0]
t = tr.stats.starttime
tr.trim(t+4600,t+7500) #zoom window

# Uncomment to read in noise
#tr = obspy.read("data/GR.FUR..LHZ.D.2008.131")[0]
#t = tr.stats.endtime
#tr.trim(t-(3600*6),t)

#
# The plotting part
#
plt.figure(1)
npts, df = tr.stats.npts, tr.stats.sampling_rate
# Seismogram
plt.subplot(211)
plt.plot(np.arange(0.0,npts)/df, tr.data)
# Spectrogram
plt.subplot(212)
spectroGram(tr.data, df, log=True)
plt.show()

# # This is what it actually calls
# # more control on the spectogram
# from matplotlib import mlab
# spectogram, freq, time = mlab.specgram(data, Fs=samp_rate,
# # db scale and remove offset
# spectogram = 10 * np.log10(spectogram[1:, :])
# freq = freq[1:]
# spectogram = np.flipud(spectogram)
# extent = 0, np.amax(time), freq[0], freq[-1]
# plt.imshow(spectogram, None, extent=extent)