Merge branch 'master' of github.com:sstaehler/dailyspec

dailyspec/data/string_bottom.txt

@@ -5,7 +5,7 @@
 <br>
 
 <script>
-var slideIndex = 6;
+var slideIndex = 1;
 showSlides(slideIndex);
 magnify("myimage", 3);
 

dailyspec/data/string_head.txt

@@ -7,6 +7,14 @@
 * {box-sizing: border-box}
 body {font-family: Verdana, sans-serif; margin:0}
 
+h1 {
+   font-family: "Trebuchet MS", Arial, Helvetica, sans-serif;
+   color: navy;
+   margin-left: 14px;
+   text-align: center;
+}
+
+
 .img-magnifier-container {position: relative;}
 
 .img-magnifier-glass {

dailyspec/get_events.py

@@ -3,11 +3,11 @@
 from obspy.clients.fdsn import Client
 
 client = Client("ETH")
-t0 = obspy.UTCDateTime('2022-01-01T00:00:00Z')
+t0 = obspy.UTCDateTime('2021-01-01T00:00:00Z')
 cat = client.get_events(starttime=t0, minmagnitude=2.0)
 
 client = Client("IRIS")
-t0 = obspy.UTCDateTime('2022-01-01T00:00:00Z')
+t0 = obspy.UTCDateTime('2021-01-01T00:00:00Z')
 
 cat += client.get_events(starttime=t0, minmagnitude=4.5)
 

dailyspec/plot_spec.py

@@ -54,6 +54,9 @@ def define_arguments():
                         help='Open a matplotlib window instead of saving to '
                              'disk')
 
+    parser.add_argument('--unit', default='VEL',
+                        help='Unit of input data. Options: ACC, VEL, DIS. Plot is in acceleration.')
+
     parser.add_argument('--kind', default='spec',
                         help='Calculate spectrogram (spec) or continuous '
                              'wavelet transfort (cwt, much slower)? Default: '
@@ -81,8 +84,8 @@ def main():
         st.decimate(2)
 
     if args.tstart is not None:
-        t0 = obspy.UTCDateTime(args.tstart) - 120.
-        t1 = obspy.UTCDateTime(args.tend) + 120.
+        t0 = obspy.UTCDateTime(args.tstart) - args.winlen*1.5
+        t1 = obspy.UTCDateTime(args.tend) + args.winlen*1.5
         st.trim(t0, t1)
 
     if args.inventory_file is not None:
@@ -93,21 +96,27 @@ def main():
             tr.stats.longitude = coords['longitude']
             tr.stats.elevation = coords['elevation']
         st.remove_response(inventory=inv, output='ACC')
+    else:
+        if args.unit=='VEL':
+            st.differentiate()
+        if args.unit=='DIS':
+            st.differentiate()
+            st.differentiate()
 
     if args.catalog_file is not None:
         cat = obspy.read_events(args.catalog_file)
     else:
         cat = None
 
-
     # The computation of the LF spectrograms with long time windows or even CWT
     # can be REALLY slow, thus, decimate it to anything larger 2.5 Hz
     st_LF = st.copy()
     while st_LF[0].stats.sampling_rate > 4.:
-        # print('LF samp rate ', st_LF[0].stats.sampling_rate, ' decimating')
         st_LF.decimate(2)
 
+    #if not args.fnam:
     fnam = "spec_{network}.{station}.{location}.{channel}.png".format(**st_LF[0].stats)
+
     calc_specgram_dual(st_LF=st_LF,
                        st_HF=st.copy(),
                        fnam=fnam, show=args.interactive,

dailyspec/process_dir.py

@@ -62,7 +62,7 @@ def main():
         raise ValueError('Either seed_id or list_seed_ids needs to be set')
 
     from .spectrogram import calc_specgram_dual
-    from .produce_cycler import produce_cycler
+    from .produce_cycler import produce_cycler_per_channel, produce_cycler_per_day
     import obspy
     from os.path import join as pjoin
     from os.path import exists as pexists
@@ -76,52 +76,64 @@ def main():
         seed_ids = [args.seed_id]
 
     for seed_id in seed_ids:
-        dir_out = f'fig_{seed_id}'
+        dir_out = pjoin('by_channel', f'fig_{seed_id}')
         mkdir(dir_out, exist_ok=True)
         network, station, location, channel = seed_id.split('.')
 
         iday_start = 366
         iday_end = 0
+        iday_end_all_channels = 0
         for iday in tqdm(range(args.jday_start, args.jday_end)):
             fnam_mseed = pjoin(args.directory, f'{args.year:4d}', network, station,
                                f'{channel}.D',
                                f'{network}.{station}.{location}.{channel}.D.{args.year:4d}.{iday:03d}')
             fnam_out = pjoin(dir_out,
                              f"spec_{network}.{station}.{location}.{channel}_{args.year}.{iday:03d}.png")
 
-            if pexists(fnam_mseed):
-                iday_start = min(iday, iday_start)
-                iday_end = max(iday, iday_end)
-
             if pexists(fnam_mseed) and not pexists(fnam_out):
                 st = obspy.read(fnam_mseed)
                 st.merge(method=1, fill_value='interpolate')
-                samp_rate_original = st[0].stats.sampling_rate
-                if samp_rate_original > args.fmax * 10 and samp_rate_original % 5 == 0.:
-                    st.decimate(5)
-                while st[0].stats.sampling_rate > 4. * args.fmax:
-                    st.decimate(2)
 
                 inv = obspy.read_inventory(args.inventory_file)
                 if args.catalog_file is not None:
                     cat = obspy.read_events(args.catalog_file)
                 else:
                     cat = None
 
-                for tr in st:
-                    coords = inv.get_coordinates(tr.get_id())
-                    tr.stats.latitude = coords['latitude']
-                    tr.stats.longitude = coords['longitude']
-                    tr.stats.elevation = coords['elevation']
-
-                st.remove_response(inventory=inv, output='ACC')
-
+                if args.inventory_file is not None:
+                    for tr in st:
+                        coords = inv.get_coordinates(tr.get_id())
+                        tr.stats.latitude = coords['latitude']
+                        tr.stats.longitude = coords['longitude']
+                        tr.stats.elevation = coords['elevation']
+
+                    try:
+                        st.remove_response(inventory=inv, output='ACC',
+                                pre_filt=(1./args.winlen * 0.5,
+                                          1./args.winlen,
+                                          args.fmax*0.95,
+                                          args.fmax))
+                    except ValueError:
+                        print(f'Problem with response for station {tr.stats.station}')
+                else:
+                    st.differentiate()
+
+                # samp_rate_original = st[0].stats.sampling_rate
+                # if samp_rate_original > args.fmax * 10 and samp_rate_original % 5 == 0.:
+                #     st.decimate(5)
+                # while st[0].stats.sampling_rate > 4. * args.fmax:
+                #     st.decimate(2)
+                st.interpolate(sampling_rate=args.fmax * 2.,
+                               method='nearest')
                 # The computation of the LF spectrograms with long time windows or even CWT
                 # can be REALLY slow, thus, decimate it to anything larger 2.5 Hz
                 st_LF = st.copy()
-                while st_LF[0].stats.sampling_rate > 4.:
-                    # print('LF samp rate ', st_LF[0].stats.sampling_rate, ' decimating')
-                    st_LF.decimate(2)
+                st_LF.filter('lowpass', freq=0.95, corners=16)
+                # while st_LF[0].stats.sampling_rate > 4.:
+                #     # print('LF samp rate ', st_LF[0].stats.sampling_rate, ' decimating')
+                #     st_LF.decimate(2)
+                st_LF.interpolate(sampling_rate=2.,
+                               method='nearest')
 
                 tstart = float(obspy.UTCDateTime(f'{args.year:4d}0101T00:00:00Z')) \
                          + 86400. * (iday - 1) - 20.
@@ -140,11 +152,23 @@ def main():
                                    catalog=cat,
                                    winlen_sec_HF=4,
                                    winlen_sec_LF=args.winlen)
-        produce_cycler(year=args.year,
-                       jday_start=iday_start,
-                       jday_end=iday_end,
-                       seed_id=seed_id,
-                       fnam_out=pjoin(dir_out, 'overview.html'))
+            if pexists(fnam_out):
+                iday_start = min(iday, iday_start)
+                iday_end = max(iday, iday_end)
+
+        produce_cycler_per_channel(year=args.year,
+                                   jday_start=iday_start,
+                                   jday_end=iday_end,
+                                   seed_id=seed_id,
+                                   fnam_out=pjoin(dir_out, 'overview.html'))
+        iday_end_all_channels = max(iday_end_all_channels, iday_end)
+
+    produce_cycler_per_day(year=args.year,
+                           jday_start=iday_start,
+                           jday_end=iday_end_all_channels,
+                           seed_ids=seed_ids,
+                           dir_image='../../by_channel',
+                           dir_out='./by_day')
 
 if __name__ == '__main__':
     main()

dailyspec/produce_cycler.py

@@ -34,7 +34,36 @@ def define_arguments():
 
     return args
 
-def produce_cycler(year, jday_start, jday_end, seed_id, fnam_out):
+def produce_cycler_per_day(year, jday_start, jday_end, seed_ids,
+        dir_image, dir_out):
+    from os.path import join as pjoin
+    import os.path as path
+    from os import makedirs as mkdir
+
+    mydir = path.dirname(path.abspath(__file__))
+    with open(pjoin(mydir, 'data/string_head.txt'), 'r') as f:
+        str_head = f.read()
+    with open(pjoin(mydir, 'data/string_bottom.txt'), 'r') as f:
+        str_bottom = f.read()
+    for iday in range(jday_start, jday_end):
+        str_all = str_head
+        for seed_id in seed_ids:
+            network, station, location, channel = seed_id.split('.')
+            fnam = pjoin(dir_image, f'fig_{seed_id}',
+                    f"spec_{network}.{station}.{location}.{channel}_{year}.{iday:03d}.png")
+            str_fnam = f'<div class="mySlides img-magnifier-container">\n' \
+                       f'<H1>{seed_id} {year}.{iday:03d}</H1>\n' \
+                       f' <img id="myimage" src="{fnam}" style="width:100%">\n' \
+                       f'</div>\n\n'
+            str_all += str_fnam
+        str_all += str_bottom
+        dir_out_sol = pjoin(dir_out, f'{year:04d}_jday_{iday:03d}')
+        mkdir(dir_out_sol, exist_ok=True)
+        with open(pjoin(dir_out_sol, 'index.html'), 'w') as f:
+            f.write(str_all)
+
+
+def produce_cycler_per_channel(year, jday_start, jday_end, seed_id, fnam_out):
     import os.path as path
 
     mydir = path.dirname(path.abspath(__file__))
@@ -46,9 +75,10 @@ def produce_cycler(year, jday_start, jday_end, seed_id, fnam_out):
     for iday in range(jday_start, jday_end):
         network, station, location, channel = seed_id.split('.')
         fnam = f"spec_{network}.{station}.{location}.{channel}_{year}.{iday:03d}.png"
-        str_fnam = f'<div class="mySlides img-magnifier-container">' \
-                   f' <img id="myimage" src="{fnam}" style="width:100%">' \
-                   f' </div>'
+        str_fnam = f'<div class="mySlides img-magnifier-container">\n' \
+                   f'<H1>{seed_id} {year}.{iday:03d}</H1>\n' \
+                   f' <img id="myimage" src="{fnam}" style="width:100%">\n' \
+                   f'</div>\n\n'
         str_all += str_fnam
     str_all += str_bottom