#!/usr/bin/env python
# coding: utf-8

# <div style="background-color: #28a745 !important; color: white; min-height: 50px; padding: 10px; margin: auto;">
#     <p style="text-align: right; font-style: italic">STMKGxHMGI Long Course</p>
#     <h1>Sesi: Membaca dan Mengolah Data Rekaman Sederhana</h1>
# </div>
# <div style="background-color: whitesmoke; padding: 10px ">
#     <ul>
#         <li>Tujuan: Memperkenalkan cara membaca dan mengolah sederhana data rekaman gempa</li>
#         <li>Keluaran: Peserta dapat memahami cara pengolahan sederhana data rekaman gempa</li>
#         <li>Sesi Materi</li>
#         <li>Waktu/Tempat: Sabtu, 25 September 2021/ Zoom</li>
#     </ul>
# </div>

# ## Membaca data rekaman gempa
# 
# Modul `obspy` memungkinkan kita membaca data rekaman dalam berbagai macam format seperti `miniseed`,`SAC`, dan yang lain. Pada materi ini kita akan berikan contoh membaca data rekaman gempa dari stasiun GSI jaringan GE yang merekam gempa Mamuju M 6,2. Cara membaca data rekaman ini menggunakan fungsi `read`:

# In[74]:


from obspy import read
st=read("data/2021-01-14-mww62-sulawesi-indonesia-4.miniseed")


# Informasi mengenai data dapat kita lihat dengan fungsi `print`:

# In[75]:


print(st)


# Dari informasi tentang data rekaman kita dapat melihat ada 3 `trace` dengan komponen BHE,BHN, dan BH dengan frekuensi sampling 20 Hz
# 
# Hasil pembacaan berupa obyek `Stream` yang kita beri nama `st`, untuk melihat metode apa saja yang bisa dilakukan kita bisa melihat petunjuknya:

# In[76]:


help(st)


# ## Mengeplot data rekaman
# 
# Melihat petunjuk di atas kita ternyata bisa mengeplot data rekaman menggunakan metode `plot`:

# In[77]:


st.plot()


# ## Melakukan filtering
# 
# `Stream` bisa kita filter dengan metode `filter`, filtering dapat dilakukan baik untuk Lowpass, Highpass, maupun Bandpass:

# In[78]:


st_LP_01 = st.copy()
st_LP_01.filter("lowpass",freq=0.1) #low pass 0.1hz
st_LP_01.plot()


# In[79]:


st_HP_1 = st.copy()
st_HP_1.filter("highpass",freq=1) #high pass 1hz
st_HP_1.plot()


# In[80]:


st_BP_23 = st.copy()
st_BP_23.filter("bandpass",freqmin=0.01, freqmax=0.05) #bandpass 0.01-0.05
st_BP_23.plot()


# ### Memotong stream
# 
# Kita dapat memotong stream menggunakan metode `slice`:

# In[81]:


from obspy import UTCDateTime
start=UTCDateTime("2021-01-14T18:31:00")
stop=start + 2*60 #2 menit
st_slice=st.copy()
st_slice=st_slice.slice(start,stop)
st_slice.plot()


# ### Mengeplot spektrogram
# 
# Spektrogram juga dapat diplot menggunakan metode `spectrogram`:

# In[82]:


st.spectrogram(wlen=50,log=True)


# ## Mengambil `trace` dari `stream`
# 
# Ada beberapa cara mengambil `trace` dari stream, paling mudah adalah menggunakan indeksing:

# In[83]:


tr_BHZ = st[0]
print(tr_BHZ)


# In[84]:


tr_BHZ.plot()


# metode pada `trace` sebagian besar sama dengan metode pada `stream`:

# In[85]:


tr_BHZ.spectrogram(wlen=50, log=True)


# In[86]:


start=UTCDateTime("2021-01-14T18:32:00")
stop=start + 2*60 #2 menit
tr_BHZ_slice = tr_BHZ.slice(start,stop)
tr_BHZ_slice.plot()


# ### Mengambil informasi dari `trace`
# 
# Informasi tentang rekaman dapat dilihat dengan atribut `stats`"

# In[87]:


tr_BHZ.stats


# Jika ingin mengambil informasi spesifik:

# In[88]:


tr_BHZ.stats.network


# ### Mengambil data dari `trace`
# 
# Data dapat diambil dengan atribut `.data`:

# In[165]:


import matplotlib.pyplot as plt
data_BHZ=tr_BHZ.data

plt.plot(range(len(data_BHZ)), data_BHZ)


# ### Contoh mengeplot gerak partikel
# 
# Kita akan terlebih dahulu memanggil data untuk komponen BHZ, BHE, dan BHN

# In[90]:


print(st)


# #### Sekitar gelombang P
# 
# Kita memotong di sekitar waktu gelombang P:

# In[112]:


st_pm = st.copy()

P=UTCDateTime("2021-01-14 18:33:08")
start=P-1*60 # 1 menit sebelum
stop=P+1*60 # 1 menit sesudah

st_pm = st_pm.slice(start,stop)
st_pm.filter("lowpass", freq=0.1)
st_pm.plot()


# #### Mengambil data untuk masing-masing komponen

# In[113]:


trdata_BHZ=st_pm[1].data
trdata_BHE=st_pm[0].data
trdata_BHN=st_pm[2].data


# #### Mengeplot untuk Z-N

# In[159]:


plt.plot(trdata_BHZ,trdata_BHN)
plt.xlabel("Z [counts]")
plt.ylabel("N [counts]")
plt.gca().set_aspect('equal')


# #### Mengeplot untuk Z-E

# In[160]:


plt.plot(trdata_BHZ,trdata_BHE)
plt.xlabel("Z [counts]")
plt.ylabel("E [counts]")
plt.gca().set_aspect('equal')


# #### Mengeplot untuk N-E

# In[161]:


plt.plot(trdata_BHN,trdata_BHE)
plt.xlabel("N [counts]")
plt.ylabel("E [counts]")
plt.gca().set_aspect('equal')


# ### Sekitar Gelombang S

# In[143]:


st_sm = st.copy()

S=UTCDateTime("2021-01-14T18:37:08")
start=S-0.5*60 # 0.5 menit sebelum
stop=S+1*60 # 1 menit sesudah

st_sm = st_sm.slice(start,stop)
st_sm.filter("lowpass", freq=0.1)
st_sm.plot()


# In[144]:


Strdata_BHZ=st_sm[1].data
Strdata_BHE=st_sm[0].data
Strdata_BHN=st_sm[2].data


# In[162]:


plt.plot(Strdata_BHZ,Strdata_BHN)
plt.xlabel("Z [counts]")
plt.ylabel("N [counts]")
plt.gca().set_aspect('equal')


# In[163]:


plt.plot(Strdata_BHZ,Strdata_BHE)
plt.xlabel("Z [counts]")
plt.ylabel("E [counts]")
plt.gca().set_aspect('equal')


# In[164]:


plt.plot(Strdata_BHN,Strdata_BHE)
plt.xlabel("N [counts]")
plt.ylabel("E [counts]")
plt.gca().set_aspect('equal')


# ### Sekitar gelombang Love

# In[152]:


st_Lm = st.copy()

L=UTCDateTime("2021-01-14 18:38:30")
start=L-0.5*60 # 1 menit sebelum
stop=L+3*60 # 1 menit sesudah

st_Lm = st_Lm.slice(start,stop)
st_Lm.filter("lowpass", freq=0.1)
st_Lm.plot()


# In[153]:


Ltrdata_BHZ=st_Lm[1].data
Ltrdata_BHE=st_Lm[0].data
Ltrdata_BHN=st_Lm[2].data


# In[156]:


plt.plot(Ltrdata_BHZ,Ltrdata_BHN)
plt.xlabel("Z [counts]")
plt.ylabel("N [counts]")
plt.gca().set_aspect('equal')


# In[157]:


plt.plot(Ltrdata_BHZ,Ltrdata_BHE)
plt.xlabel("Z [counts]")
plt.ylabel("E [counts]")
plt.gca().set_aspect('equal')


# In[158]:


plt.plot(Ltrdata_BHN,Ltrdata_BHE)
plt.xlabel("N [counts]")
plt.ylabel("E [counts]")
plt.gca().set_aspect('equal')


# In[ ]: