From 745eaf1f55fb07a8b5f61599a85afc15244a4d5d Mon Sep 17 00:00:00 2001
From: nult <leitd@qq.com>
Date: Wed, 18 May 2022 18:14:31 -0400
Subject: [PATCH] update-plot-script

---
 main_hkrec.py | 21 +++++++++++++++++++--
 plt_rf.py     | 44 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 63 insertions(+), 2 deletions(-)
 create mode 100644 plt_rf.py

diff --git a/main_hkrec.py b/main_hkrec.py
index 3c815df..0e5aa5c 100755
--- a/main_hkrec.py
+++ b/main_hkrec.py
@@ -6,9 +6,12 @@
 from hkrec import recfun
 # to find the station active time range: https://ds.iris.edu/SeismiQuery/station.htm (select time range)
 # first test the example from HK package example directory KUL.BH?.01   2003/6/16 (167) 22:08:02.140
+staname='KUL'
+netname='XH'
 #rec=recfun(sta='KUL',net='XH',chan='BH',loc_code='',data_starttime='2002-12-22',data_endtime='2003-12-31')
-rec=recfun(sta='KUL',net='XH',chan='BH',loc_code='',data_starttime='2003-06-15',data_endtime='2003-07-16')
+#rec=recfun(sta='KUL',net='XH',chan='BH',loc_code='',data_starttime='2003-06-15',data_endtime='2003-07-16')
 #rec=recfun(sta='WINE',net='7A',chan='BH',loc_code='',data_starttime='2013-10-01',data_endtime='2016-10-01')
+rec=recfun(sta=staname,net=netname,chan='BH',loc_code='',data_starttime='2003-06-15',data_endtime='2003-07-16')
 rec.event_par.magnitude_range.min_mag=6.5
 rec.get_station_info()
 rec.get_events()
@@ -16,4 +19,18 @@
 # more on data processing: http://eqseis.geosc.psu.edu/cammon/HTML/RftnDocs/prep01.html
 # TODO: more to be added here
 print('selected_event ids:', rec.accepted_event_list)
-pickle.dump(rec,open('rec.pkl','wb'))
\ No newline at end of file
+pickle.dump(rec,open('rec.pkl','wb'))
+
+# output bash script to run hk
+with open(staname+'.'+netname+'.sh', 'w') as fpp:
+    fpp.write('#!/bin/bash\n')
+    fpp.write('f1=3\n')
+    fpp.write('f2=5\n')
+    for ievn in rec.accepted_event_list:
+        fpp.write('./hk/iter_decon -C-2/-10/80 -F1/$f1/-$f2 -N100 -T0.1 %s %s\n'%('./data/'+staname+'.'+netname+\
+                '..BHZ.'+str(ievn)+'.p.sac','./data/'+staname+'.'+netname+'..BHR.'+str(ievn)+'.p.sac'))
+    fpp.write('\n')
+    fpp.write('./hk/k_stack -R30/60/1.6/2.0 -I0.5/0.01 -G%s %s\n'%('./'+staname+'.'+netname+\
+            '.grd','./data/'+staname+'.'+netname+'..BHR.*.p.saci'))
+    fpp.write('./hk/grdmin -D %s\n'%('./'+staname+'.'+netname+'.grd'))
+    fpp.write('./hk/grdmin -D %s | ./hk/hk_plt.pl > %s'%('./'+staname+'.'+netname+'.grd','./'+staname+'.'+netname+'.ps'))
diff --git a/plt_rf.py b/plt_rf.py
new file mode 100644
index 0000000..3b116fa
--- /dev/null
+++ b/plt_rf.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+import os
+import obspy
+import numpy as np
+from matplotlib import pyplot as plt
+
+def get_file_list(basis_dir="./", begin="", end=""):
+  path_list = os.listdir(basis_dir)
+  list_final = []
+  for partial in path_list:
+    if begin and end:
+      if partial[:len(begin)] == begin and partial[-len(end):] == end:
+        list_final.append(partial)
+    elif end:
+      if partial[-len(end):] == end:
+        list_final.append(partial)
+    elif begin:
+      if partial[:len(begin)] == begin:
+        list_final.append(partial)
+    else:
+      list_final.append(partial)
+  return list_final
+
+staname='LBDL'
+netname='7A'
+scale_num=0.005
+file_list = get_file_list("./data",begin=staname+'.'+netname,end=".saci")
+plt.figure()
+for ifile in file_list:
+  print('File: %s'%(ifile))
+  st=obspy.read("./data/"+ifile)
+  print(st[0].stats.sac.get('user0'))
+  plt.plot(np.arange(st[0].stats.sac.get('b'),st[0].stats.sac.get('e'),\
+          st[0].stats.sac.get('delta')),st[0].stats.sac.get('user0')+\
+          scale_num*st[0].data/max(st[0].data),linewidth='1',color='black')
+  plt.fill_between(np.arange(st[0].stats.sac.get('b'),st[0].stats.sac.get('e'),\
+          st[0].stats.sac.get('delta')),st[0].stats.sac.get('user0'),\
+          st[0].stats.sac.get('user0')+scale_num*st[0].data/max(st[0].data),\
+          where=scale_num*st[0].data>0,facecolor='black')
+plt.xlim(-3,20)
+plt.xlabel('Time (s)')
+plt.ylabel('p (s/km)')
+plt.show()