wip: update magma plot

dsar.py

@@ -1,38 +1,49 @@
 #!/usr/bin/env python
 # coding: utf-8
 
-# In[35]:
+# In[137]:
 
 
 import pandas as pd
 import os
 import matplotlib.pyplot as plt
 from statsmodels.tsa.filters.hp_filter import hpfilter
+from datetime import datetime
 
 
-# In[48]:
+# In[138]:
 
 
 network = "VG"
-station = "PSAG"
+station = "TMKS"
 location = "00"
 channel = "EHZ"
 
 nslc = "{}.{}.{}.{}".format(network, station, location, channel)
 
 
-# In[47]:
+# In[139]:
+
+
+bands: dict[str, list[float]] = {
+    'HF' : [0.1, 8.0, 16.0],
+    'LF' : [0.1, 4.5, 8.0],
+}
+
+
+# In[143]:
 
 
 current_dir: str = os.getcwd()
+input_directory: str = os.path.join(current_dir, "output")
 output_directory: str = os.path.join(current_dir, "output", "dsar")
 os.makedirs(output_directory, exist_ok=True)
 
-combined_HF_csv: str = r'D:\Projects\dsar\output\HF\combined_0.1-8.0-16.0Hz_VG.PSAG.00.EHZ.csv'
-combined_LF_csv: str = r'D:\Projects\dsar\output\LF\combined_0.1-4.5-8.0Hz_VG.PSAG.00.EHZ.csv'
+combined_HF_csv: str = os.path.join(input_directory, "HF", 'combined_{}Hz_{}.csv'.format('-'.join(map(str,bands['HF'])), nslc))
+combined_LF_csv: str = os.path.join(input_directory, "LF", 'combined_{}Hz_{}.csv'.format('-'.join(map(str,bands['LF'])), nslc))
 
 
-# In[37]:
+# In[144]:
 
 
 HF = pd.read_csv(combined_HF_csv, names=["datetime", "values"], 
@@ -41,87 +52,119 @@
                  index_col='datetime', parse_dates=True)
 
 
-# In[38]:
+# In[145]:
 
 
 HF
 
 
-# In[39]:
+# In[146]:
 
 
 LF
 
 
-# In[40]:
+# In[147]:
 
 
 df = pd.DataFrame()
 
 
-# In[41]:
+# In[148]:
 
 
 df['LF'] = LF['values']
 df['HF'] = HF['values']
 df['DSAR'] = (LF['values']/HF['values'])
+df['DSAR_365'] = (LF['values']/HF['values']).rolling(6).median()
 
 
-# In[42]:
+# In[149]:
 
 
 df
 
 
-# In[43]:
+# In[150]:
 
 
 df = df.dropna()
 
 
-# In[44]:
+# In[151]:
 
 
-df = df.apply(lambda col: col.drop_duplicates())
+df.loc[~df.index.duplicated(), :]
 
 
-# In[45]:
+# In[152]:
 
 
-df = df.interpolate('time').interpolate()
+# df = df.apply(lambda col: col.drop_duplicates())
 
 
-# In[46]:
+# In[153]:
 
 
-df
+df = df.interpolate('time').interpolate()
 
 
-# In[49]:
+# In[154]:
 
 
-filename = os.path.join(output_directory, "DSAR_{}.csv".format(nslc))
-df.to_csv(filename)
+df
 
 
-# In[50]:
+# In[155]:
 
 
-plt.subplot(111)
+filename = os.path.join(output_directory, "DSAR_{}.csv".format(nslc))
+df.to_csv(filename)
 
 
-# In[61]:
+# In[157]:
 
 
-plt.scatter(df.index, df.DSAR, c= 'k', alpha=0.3, s=10, label='{} (10 minutes)'.format(nslc))
+import matplotlib.dates as mdates
 
 # HP filter documentation https://www.statsmodels.org/stable/generated/statsmodels.tsa.filters.hp_filter.hpfilter.html
-_,trend = hpfilter(df.DSAR, 1600)
-plt.plot(df.index, trend, c='lightseagreen', label='{} smoothed (HP Filter)'.format(nslc), alpha=0.8)
-
-plt.legend(loc=2)
-plt.ylabel('DSAR')
-plt.ylim(1,5)
+cycle,trend = hpfilter(df.DSAR, 1000000)
+
+fig, axs = plt.subplots(figsize=(12, 5), layout="constrained")
+scatter = axs.scatter(df.index, df.DSAR, c= 'k', alpha=0.3, s=10, label='{} (10 minutes)'.format(nslc))
+smoothed_using_HP_filter = axs.plot(df.index, trend, c='red', label='{} smoothed (HP Filter)'.format(nslc), alpha=1)
+
+axs.legend(loc=2)
+axs.set_ylabel('DSAR')
+axs.set_xlabel('Date')
+axs.xaxis.set_major_locator(mdates.DayLocator(interval=7))
+axs.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))
+axs.set_ylim(0,6.5)
+
+# Agung Eruption
+continous_eruptions = [
+    ['2017-11-21', '2017-11-29'],
+    ['2018-06-27', '2018-07-16'],
+    ['2018-07-24', '2018-07-27'],
+]
+
+for continous in continous_eruptions:
+    axs.axvspan(continous[0], continous[1], alpha=0.4, color='orange')
+
+single_eruptions = [
+    '2018-05-29',
+    '2018-06-10',
+    '2018-06-13',
+    '2018-06-15',
+    '2018-07-21',
+    '2018-07-24',
+]
+
+for date in single_eruptions:
+    axs.axvline(datetime.strptime(date, '%Y-%m-%d'), alpha=0.4, color='orange')
+
+for label in axs.get_xticklabels(which='major'):
+    label.set(rotation=30, horizontalalignment='right')
 
 
 # In[ ]: