From 068001a6c2018a4bb9f36e4e562e4d6c495de1a7 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 22 Jul 2021 13:41:51 -0700 Subject: [PATCH 01/54] utils.py: Show waveform download exception --- rfpy/utils.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/rfpy/utils.py b/rfpy/utils.py index 72a98b8..df118c8 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -552,9 +552,13 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, print("* - Z12 Data Downloaded") else: st = None - except: + except Exception as e: + print("* Met exception:") + print("* " + e.__repr__()) st = None - except: + except Exception as e: + print("* Met exception:") + print("* " + e.__repr__()) st = None # Break if we successfully obtained 3 components in st From 21cc21496005b1a364be0db51a93daaa96c86d04 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Fri, 6 Aug 2021 07:10:51 -0700 Subject: [PATCH 02/54] parse_local_data: abbreviate search for local data code, catch missing sac user9 kw --- rfpy/utils.py | 190 ++++++++++++++++++-------------------------------- 1 file changed, 67 insertions(+), 123 deletions(-) diff --git a/rfpy/utils.py b/rfpy/utils.py index df118c8..f4d5857 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -84,7 +84,9 @@ def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='SAC', altnet=[]) # Keep paths only for those matching the station for sstring in sstrings: for filename in filter(filenames, sstring): - fpathmatch.append(join(root, filename)) + # No hidden directories + if not '/.' in root: + fpathmatch.append(join(root, filename)) fpathmatch.sort() @@ -135,49 +137,34 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, ("* {0:2s}{1:1s} - Checking Disk".format(sta.channel.upper(), comp.upper()))) + # Possible naming conventions + f1 = '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}' + f2 = '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}' + f3 = '*/{0:4s}.{1:3s}.*.{2:s}.{3:s}.*.{4:2s}{5:1s}*.{6:s}' + # Time Window Spans Single Day if stjd == edjd: - # Format 1 - lclfiles = list(filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'.format( - styr, stjd, sta.network.upper( - ), sta.station.upper(), sta.channel.upper()[0:2], - comp.upper(), dtype))) - # Format 2 - if len(lclfiles) == 0: - lclfiles = list(filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format( - styr, stjd, sta.network.upper(), sta.station.upper(), - comp.upper(), dtype))) - # Alternate Nets (for CN/PO issues) Format 1 - if len(lclfiles) == 0: - lclfiles = [] - for anet in sta.altnet: - lclfiles.extend( - list( - filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.' + - '{4:2s}{5:1s}.{6:s}'.format( - styr, stjd, anet.upper(), sta.station.upper(), - sta.channel.upper()[0:2], comp.upper(), dtype)))) - - # Alternate Nets (for CN/PO issues) Format 2 - if len(lclfiles) == 0: - # Check Alternate Networks - lclfiles = [] - for anet in sta.altnet: - lclfiles.extend( - list( - filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*' + - '{4:1s}.{5:s}'.format( - styr, stjd, sta.network.upper(), - sta.station.upper(), comp.upper(), dtype)))) + lclfiles = [] + nets = [sta.network] + sta.altnet + for net in nets: + # Start day + s1 = f1.format(styr, stjd, net.upper(), sta.station.upper(), + sta.channel.upper()[0:2], comp.upper(), dtype) + s2 = f2.format(styr, stjd, net.upper(), sta.station.upper(), + comp.upper(), dtype) + s3 = f3.format(styr, stjd, net.upper(), sta.station.upper(), + sta.channel.upper()[0:2], comp.upper(), dtype) + + print("* Trying formats:") + print("* " + s1) + print("* " + s2) + print("* " + s3) + print("* ") + + lclfiles.extend(list(filter(stdata, s1))) + lclfiles.extend(list(filter(stdata, s2))) + lclfiles.extend(list(filter(stdata, s3))) # If still no Local files stop if len(lclfiles) == 0: @@ -209,7 +196,10 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, eddt = False # Check for NoData and convert to NaN if a SAC file if dtype.upper() == 'SAC': - stnd = st[0].stats.sac['user9'] + try: + stnd = st[0].stats.sac['user9'] + except KeyError: + stnd = 0.0 if (not stnd == 0.0) and (not stnd == -12345.0): st[0].data[st[0].data == stnd] = ndval eddt = True @@ -244,83 +234,33 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Time Window spans Multiple days else: - # Day 1 Format 1 - lclfiles1 = list( - filter(stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'.format( - styr, stjd, sta.network.upper(), sta.station.upper(), - sta.channel.upper()[0:2], comp.upper(), dtype))) - # Day 1 Format 2 - if len(lclfiles1) == 0: - lclfiles1 = list( - filter(stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format( - styr, stjd, sta.network.upper(), - sta.station.upper(), comp.upper(), dtype))) - # Day 1 Alternate Nets (for CN/PO issues) Format 1 - if len(lclfiles1) == 0: - lclfiles1 = [] - for anet in sta.altnet: - lclfiles1.extend( - list( - filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.' + - '{4:2s}{5:1s}.{6:s}'.format( - styr, stjd, anet.upper(), sta.station.upper( - ), sta.channel.upper()[0:2], - comp.upper(), dtype)))) - # Day 1 Alternate Nets (for CN/PO issues) Format 2 - if len(lclfiles1) == 0: - lclfiles1 = [] - for anet in sta.altnet: - lclfiles1.extend( - list( - filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format( - styr, stjd, anet.upper(), - sta.station.upper(), comp.upper(), dtype)))) - - # Day 2 Format 1 - lclfiles2 = list( - filter(stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.{4:2s}{5:1s}.{6:s}'.format( - edyr, edjd, sta.network.upper( - ), sta.station.upper(), sta.channel.upper()[0:2], - comp.upper(), dtype))) - # Day 2 Format 2 - if len(lclfiles2) == 0: - lclfiles2 = list( - filter(stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*' + - '{4:1s}.{5:s}'.format( - edyr, edjd, sta.network.upper(), - sta.station.upper(), - comp.upper(), dtype))) - # Day 2 Alternate Nets (for CN/PO issues) Format 1 - if len(lclfiles2) == 0: - lclfiles2 = [] - for anet in sta.altnet: - lclfiles2.extend( - list( - filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.' + - '{4:2s}{5:1s}.{6:s}'.format( - edyr, edjd, anet.upper(), sta.station.upper(), - sta.channel.upper()[0:2], comp.upper(), dtype)))) - # Day 2 Alternate Nets (for CN/PO issues) Format 2 - if len(lclfiles2) == 0: - lclfiles2 = [] - for anet in sta.altnet: - lclfiles2.extend( - list( - filter( - stdata, - '*/{0:4s}.{1:3s}.{2:s}.{3:s}.*.*{4:1s}.{5:s}'.format( - edyr, edjd, anet.upper(), sta.station.upper(), - comp.upper(), dtype)))) + lclfiles1 = [] + lclfiles2 = [] + nets = [sta.network] + sta.altnet + for net in nets: + # Start day + s1 = f1.format(styr, stjd, net.upper(), sta.station.upper(), + sta.channel.upper()[0:2], comp.upper(), dtype) + s2 = f2.format(styr, stjd, net.upper(), sta.station.upper(), + comp.upper(), dtype) + s3 = f3.format(styr, stjd, net.upper(), sta.station.upper(), + sta.channel.upper()[0:2], comp.upper(), dtype) + + lclfiles1.extend(list(filter(stdata, s1))) + lclfiles1.extend(list(filter(stdata, s2))) + lclfiles1.extend(list(filter(stdata, s3))) + + # End day + s1 = f1.format(edyr, edjd, net.upper(), sta.station.upper(), + sta.channel.upper()[0:2], comp.upper(), dtype) + s2 = f2.format(edyr, edjd, net.upper(), sta.station.upper(), + comp.upper(), dtype) + s3 = f3.format(edyr, edjd, net.upper(), sta.station.upper(), + sta.channel.upper()[0:2], comp.upper(), dtype) + + lclfiles2.extend(list(filter(stdata, s1))) + lclfiles2.extend(list(filter(stdata, s2))) + lclfiles2.extend(list(filter(stdata, s3))) # If still no Local files stop if len(lclfiles1) == 0 and len(lclfiles2) == 0: @@ -334,8 +274,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, st1 = read(sacf1) if dtype.upper() == 'MSEED': if len(st1) > 1: - st1.merge(method=1, interpolation_samples=- - 1, fill_value=-123456789) + st1.merge(method=1, interpolation_samples=-1, + fill_value=-123456789) # Loop over second day file options for sacf2 in lclfiles2: @@ -343,7 +283,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, if dtype.upper() == 'MSEED': if len(st2) > 1: st2.merge( - method=1, interpolation_samples=-1, fill_value=-123456789) + method=1, interpolation_samples=-1, + fill_value=-123456789) # Check time overlap of the two files. if st1[0].stats.endtime >= \ @@ -383,7 +324,10 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, eddt = False # Check for NoData and convert to NaN if a SAC file if dtype.upper() == 'SAC': - stnd = st[0].stats.sac['user9'] + try: + stnd = st[0].stats.sac['user9'] + except KeyError: + stnd = 0.0 if (not stnd == 0.0) and (not stnd == -12345.0): st[0].data[st[0].data == stnd] = ndval eddt = True From b54c6d5c110785d31341ec918de158df11cb9b6a Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Fri, 6 Aug 2021 15:48:54 -0700 Subject: [PATCH 03/54] More logging funcitonality --- rfpy/rfdata.py | 11 +++--- rfpy/utils.py | 98 +++++++++++++++++++++++++------------------------- 2 files changed, 56 insertions(+), 53 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 4caf39a..651d928 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -22,6 +22,7 @@ # -*- coding: utf-8 -*- from math import ceil +import logging import numpy as np from obspy import Trace, Stream, UTCDateTime from rfpy import utils @@ -114,9 +115,9 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'): source_depth_in_km=self.dep, phase_list=[phase]) if len(arrivals) > 1: - print("arrival has many entries: ", len(arrivals)) + logging.warning("arrival has many entries: ", len(arrivals)) elif len(arrivals) == 0: - print("no arrival found") + logging.warning("no arrival found") self.accept = False return @@ -343,9 +344,9 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, tend = self.meta.time + self.meta.ttime + dts # Get waveforms - print("* Requesting Waveforms: ") - print("* Startime: " + tstart.strftime("%Y-%m-%d %H:%M:%S")) - print("* Endtime: " + tend.strftime("%Y-%m-%d %H:%M:%S")) + logging.info("* Requesting Waveforms: ") + logging.info("* Startime: " + tstart.strftime("%Y-%m-%d %H:%M:%S")) + logging.info("* Endtime: " + tend.strftime("%Y-%m-%d %H:%M:%S")) # Download data err, stream = utils.download_data( diff --git a/rfpy/utils.py b/rfpy/utils.py index f4d5857..d352f7c 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -1,4 +1,5 @@ import math +import logging from obspy import UTCDateTime from numpy import nan, isnan, abs import numpy as np @@ -133,7 +134,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Intialize to default positive error erd = True - print( + logging.info( ("* {0:2s}{1:1s} - Checking Disk".format(sta.channel.upper(), comp.upper()))) @@ -156,11 +157,11 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, s3 = f3.format(styr, stjd, net.upper(), sta.station.upper(), sta.channel.upper()[0:2], comp.upper(), dtype) - print("* Trying formats:") - print("* " + s1) - print("* " + s2) - print("* " + s3) - print("* ") + logging.debug("* Trying formats:") + logging.debug("* " + s1) + logging.debug("* " + s2) + logging.debug("* " + s3) + logging.debug("* ") lclfiles.extend(list(filter(stdata, s1))) lclfiles.extend(list(filter(stdata, s2))) @@ -168,7 +169,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # If still no Local files stop if len(lclfiles) == 0: - print("* - Data Unavailable") + logging.warning("* - Data Unavailable") return erd, None # Process the local Files @@ -206,18 +207,18 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Check for Nan in stream for SAC if True in isnan(st[0].data): - print( + logging.warning( "* !!! Missing Data Present !!! " + "Skipping (NaNs)") # Check for ND Val in stream for MSEED elif -123456789 in st[0].data: - print( + logging.warning( "* !!! Missing Data Present !!! " + "Skipping (MSEED fill)") else: if eddt and (ndval == 0.0): if any(st[0].data == 0.0): - print( + logging.warning( "* !!! Missing Data Present " + "!!! (Set to Zero)") @@ -227,7 +228,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, tloc = "--" # Processed succesfully...Finish - print(("* {1:3s}.{2:2s} - From Disk".format( + logging.info(("* {1:3s}.{2:2s} - From Disk".format( st[0].stats.station, st[0].stats.channel.upper(), tloc))) return False, st @@ -264,7 +265,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # If still no Local files stop if len(lclfiles1) == 0 and len(lclfiles2) == 0: - print("* - Data Unavailable") + logging.warning("* - Data Unavailable") return erd, None # Now try to merge the two separate day files @@ -313,8 +314,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Should only be one component, otherwise keep # reading If more than 1 component, error if len(st) != 1: - print(st) - print("merge failed?") + logging.error(st) + logging.error("merge failed?") else: if (st[0].stats.starttime <= start and @@ -334,18 +335,18 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Check for Nan in stream for SAC if True in isnan(st[0].data): - print( + logging.warning( "* !!! Missing Data " + "Present !!! Skipping (NaNs)") # Check for ND Val in stream for MSEED elif -123456789 in st[0].data: - print( + logging.warning( "* !!! Missing Data Present !!! " + "Skipping (MSEED fill)") else: if (eddt1 or eddt2) and (ndval == 0.0): if any(st[0].data == 0.0): - print( + logging.warning( "* !!! Missing " + "Data Present !!! (Set " + "to Zero)") @@ -356,7 +357,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, tloc = "--" # Processed succesfully...Finish - print(("* {1:3s}.{2:2s} - " + + logging.info(( + "* {1:3s}.{2:2s} - " + "From Disk".format( st[0].stats.station, st[0].stats.channel.upper(), @@ -367,14 +369,14 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, pass else: st2ot = st2[0].stats.endtime-st2[0].stats.delta - print("* - Merge Failed: No " + + logging.warning("* - Merge Failed: No " + "Overlap {0:s} - {1:s}".format( st1[0].stats.endtime.strftime( "%Y-%m-%d %H:%M:%S"), st2ot.strftime("%Y-%m-%d %H:%M:%S"))) # If we got here, we did not get the data. - print("* - Data Unavailable") + logging.warning("* - Data Unavailable") return erd, None @@ -424,7 +426,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, from math import floor # Output - print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, + logging.info(("* {0:s}.{1:2s} - ZNE:".format(sta.station, sta.channel.upper()))) # Set Error Default to True @@ -475,7 +477,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, channel=channelsZNE, starttime=start, endtime=end+1., attach_response=False) if len(st) == 3: - print("* - ZNE Data Downloaded") + logging.info("* - ZNE Data Downloaded") # It's possible if len(st)==1 that data is Z12 else: @@ -485,7 +487,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, sta.channel.upper() + '2' msg = "* {1:2s}[Z12].{2:2s} - Checking Network".format( sta.station, sta.channel.upper(), tloc) - print(msg) + logging.info(msg) try: st = client.get_waveforms( network=sta.network, @@ -493,16 +495,16 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, channel=channelsZ12, starttime=start, endtime=end+1., attach_response=False) if len(st) == 3: - print("* - Z12 Data Downloaded") + logging.info("* - Z12 Data Downloaded") else: st = None except Exception as e: - print("* Met exception:") - print("* " + e.__repr__()) + logging.warning("* Met exception:") + logging.warning("* " + e.__repr__()) st = None except Exception as e: - print("* Met exception:") - print("* " + e.__repr__()) + logging.warning("* Met exception:") + logging.warning("* " + e.__repr__()) st = None # Break if we successfully obtained 3 components in st @@ -512,8 +514,8 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Check the correct 3 components exist if st is None: - print("* Error retrieving waveforms") - print("**************************************************") + logging.warning("* Error retrieving waveforms") + logging.warning("**************************************************") return True, None # Three components successfully retrieved @@ -525,12 +527,12 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Check start times if not np.all([tr.stats.starttime == start for tr in st]): - print("* Start times are not all close to true start: ") - [print("* "+tr.stats.channel+" " + + logging.warning("* Start times are not all close to true start: ") + [logging.warning("* "+tr.stats.channel+" " + str(tr.stats.starttime)+" " + str(tr.stats.endtime)) for tr in st] - print("* True start: "+str(start)) - print("* -> Shifting traces to true start") + logging.warning("* True start: "+str(start)) + logging.warning("* -> Shifting traces to true start") delay = [tr.stats.starttime - start for tr in st] st_shifted = Stream( traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)]) @@ -540,39 +542,39 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, sr = st[0].stats.sampling_rate sr_round = float(floor_decimal(sr, 0)) if not sr == sr_round: - print("* Sampling rate is not an integer value: ", sr) - print("* -> Resampling") + logging.warning("* Sampling rate is not an integer value: ", sr) + logging.warning("* -> Resampling") st.resample(sr_round, no_filter=False) # Try trimming try: st.trim(start, end) except: - print("* Unable to trim") - print("* -> Skipping") - print("**************************************************") + logging.error("* Unable to trim") + logging.error("* -> Skipping") + logging.error("**************************************************") return True, None # Check final lengths - they should all be equal if start times # and sampling rates are all equal and traces have been trimmed if not np.allclose([tr.stats.npts for tr in st[1:]], st[0].stats.npts): - print("* Lengths are incompatible: ") - [print("* "+str(tr.stats.npts)) for tr in st] - print("* -> Skipping") - print("**************************************************") + logging.error("* Lengths are incompatible: ") + [logging.error("* "+str(tr.stats.npts)) for tr in st] + logging.error("* -> Skipping") + logging.error("**************************************************") return True, None elif not np.allclose([st[0].stats.npts], int((end - start)*sr), atol=1): - print("* Length is too short: ") - print("* "+str(st[0].stats.npts) + + logging.error("* Length is too short: ") + logging.error("* "+str(st[0].stats.npts) + " ~= "+str(int((end - start)*sr))) - print("* -> Skipping") - print("**************************************************") + logging.error("* -> Skipping") + logging.error("**************************************************") return True, None else: - print("* Waveforms Retrieved...") + logging.info("* Waveforms Retrieved...") return False, st From 329011128e870f594d4e5aa0e4434afd047b10ed Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 4 Oct 2021 15:56:08 -0700 Subject: [PATCH 04/54] rfdata.py: Add option to write out deconvolution wavelets --- rfpy/rfdata.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 651d928..f957e96 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -23,6 +23,7 @@ # -*- coding: utf-8 -*- from math import ceil import logging +import pickle import numpy as np from obspy import Trace, Stream, UTCDateTime from rfpy import utils @@ -178,7 +179,6 @@ def __init__(self, sta): if sta == 'demo' or sta == 'Demo': print("Uploading demo data - station NY.MMPY") import os - import pickle sta = pickle.load( open(os.path.join( os.path.dirname(__file__), @@ -270,7 +270,6 @@ def add_data(self, stream, returned=False, new_sr=5.): # Load demo data if stream == 'demo' or stream == 'Demo': import os - import pickle file = open(os.path.join( os.path.dirname(__file__), "examples/data", "ZNE_Data.pkl"), "rb") @@ -616,7 +615,7 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.): def deconvolve(self, phase='P', vp=None, vs=None, align=None, method='wiener', pre_filt=None, - gfilt=None, wlevel=0.01): + gfilt=None, wlevel=0.01, writeto=None): """ Deconvolves three-component data using one component as the source wavelet. The source component is always taken as the dominant compressional @@ -641,6 +640,8 @@ def deconvolve(self, phase='P', vp=None, vs=None, Center frequency of Gaussian filter (Hz). wlevel : float Water level used in ``method='water'``. + writeto : str or None + Write wavelets for deconvolution to file. Attributes ---------- @@ -756,6 +757,11 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): gauss = np.ones(npad) gnorm = 1. + # Is this correct? + #parent de-noised = np.fft.ifftshift(np.real(np.fft.ifft(gauss*Sdenom))/gnorm) + #daughter1 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd1p))/gnorm) + #daughter2 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd2p))/gnorm) + # Copy traces rfp = parent.copy() rfd1 = daughter1.copy() @@ -796,6 +802,12 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): trNl = self.data.select(component=cL)[0].copy() trNq = self.data.select(component=cQ)[0].copy() + trL.stats.channel = 'WV' + self.meta.align[0] + trQ.stats.channel = 'WV' + self.meta.align[1] + trT.stats.channel = 'WV' + self.meta.align[2] + trNl.stats.channel = 'WV' + self.meta.align[0] + trNq.stats.channel = 'WV' + self.meta.align[1] + if phase == 'P' or 'PP': # Get signal length (i.e., seismogram to deconvolve) from trace length @@ -847,6 +859,10 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): corners=2, zerophase=True) for tr in [trL, trQ, trT, trNl, trNq]] + if writeto: + with open(writeto, 'wb') as f: + pickle.dump(Stream(traces=[trL, trQ, trT]), f) + # Deconvolve if phase == 'P' or 'PP': rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method) @@ -950,7 +966,6 @@ def save(self, file): """ - import pickle output = open(file, 'wb') pickle.dump(self, output) output.close() From d27059e966467a55012c7ac849da06c84a52b972 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Tue, 19 Oct 2021 09:40:53 -0700 Subject: [PATCH 05/54] rfdata.py: add _wavelet method to manipulate de-convolution wavelet --- rfpy/rfdata.py | 117 ++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 106 insertions(+), 11 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index f957e96..c59a571 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -614,7 +614,8 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.): def deconvolve(self, phase='P', vp=None, vs=None, - align=None, method='wiener', pre_filt=None, + align=None, method='wiener', wavelet='complete', + envelope_threshold=0.05, time=5, pre_filt=None, gfilt=None, wlevel=0.01, writeto=None): """ Deconvolves three-component data using one component as the source wavelet. @@ -633,6 +634,13 @@ def deconvolve(self, phase='P', vp=None, vs=None, method : str Method for deconvolution. Options are 'wiener', 'water' or 'multitaper' + wavelet : str + Type of wavelet for deconvolution. Options are 'complete', 'time' or + 'envelope' + envelope_threshold : float + Threshold [0-1] used in ``wavelet='envelope'``. + time : float + Window length used in ``wavelet='time'``. pre_filt : list of 2 floats Low and High frequency corners of bandpass filter applied before deconvolution @@ -671,6 +679,84 @@ def _gauss_filt(dt, nft, f0): gauss[nft21:] = np.flip(gauss[1:nft21-1]) return gauss + def _wavelet(parent, method='complete', overhang=5, + envelope_threshold=0.05, time=5): + + """ + Select wavelet from the parent function for deconvolution using method. + parent: obspy.Trace + wavefrom to extract the wavelet from + method: str + 'complete' use complete parent signal after P arrival (current + implementation) + 'envelope' use only the part of the parent signal after the + P arrival where + envelope > envelope_threshold*max(envelope) + fall back to 'complete' if condition not reached + 'time' use only this many seonds after P arrival` + fall back to 'complete' if longer than parent + overhang: float + seconds before start and after end of wavelet to be used for + tapering + envelope_threshold: float + fraction of the envelope that defines wavelet (for + method='envelope') + time: float + window (seconds) that defines the wavelet (for method='time') + + Return: + left, right: (obspy.UTCDateTime) start and end time of wavelet + """ + import obspy.signal.filter as osf + + if method not in ['complete', 'envelope', 'time', 'noise']: + msg = 'Unknow method for wavelet extraction: ' + method + raise NotImplementedError(msg) + + errflag = False + + if method == 'envelope': + split = int((self.meta.time + self.meta.ttime - + parent.stats.starttime) * parent.stats.sampling_rate) + env = osf.envelope(parent.data) + env = env[split:] # only look after P + env /= max(env) # normalize + try: + i = np.nonzero( + np.diff( + np.array( + env > envelope_threshold, dtype=int))==-1)[0][0] + except IndexError: + i = len(parent.data)-1 + dts = i * parent.stats.delta + left = self.meta.time + self.meta.ttime - overhang + right = self.meta.time + self.meta.ttime + dts + overhang + if left < parent.stats.starttime or right > parent.stats.endtime: + logging.info('Envelope wavelet longer than trace.') + logging.info('Falling back to "complete" wavelet') + errflag = True + # TODO: Test me! + + if method == 'time': + left = self.meta.time + self.meta.ttime - overhang + right = self.meta.time + self.meta.ttime + time + overhang + if left < parent.stats.starttime or right > parent.stats.endtime: + logging.info('Time wavelet longer than trace.') + logging.info('Falling back to "complete" wavelet') + errflag = True + + if method == 'complete' or errflag: + dts = len(parent.data)*parent.stats.delta/2. + left = self.meta.time + self.meta.ttime - overhang + right = self.meta.time + self.meta.ttime + dts - 2*overhang + + if method == 'noise': + dts = len(parent.data)*parent.stats.delta/2. + left = self.meta.time + self.meta.ttime - dts + right = self.meta.time + self.meta.ttime - overhang + + return left, right + def _decon(parent, daughter1, daughter2, noise, nn, method): # Get length, zero padding parameters and frequencies @@ -811,24 +897,32 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): if phase == 'P' or 'PP': # Get signal length (i.e., seismogram to deconvolve) from trace length - dts = len(trL.data)*trL.stats.delta/2. - nn = int(round((dts-5.)*trL.stats.sampling_rate)) + 1 + over = 5 + dtsqt = len(trL.data)*trL.stats.delta/2. + + # Traces will be paded to this length (samples) + nn = int(round((dtsqt+over)*trL.stats.sampling_rate)) + 1 + + sig_left, sig_right = _wavelet(trL, method='envelope', + envelope_threshold=envelope_threshold, overhang=over) + + # Trim wavelet + trL.trim(sig_left, sig_right, nearest_sample=False, pad=True, + fill_value=0.) - # Signal window (-5. to dts-10 sec) - sig_left = self.meta.time+self.meta.ttime-5. - sig_right = self.meta.time+self.meta.ttime+dts-10. + # Signal window (-5. to dtsqt-10 sec) + sig_left, sig_right = _wavelet(trQ, method='complete', overhang=over) # Trim signal traces [tr.trim(sig_left, sig_right, nearest_sample=False, - pad=nn, fill_value=0.) for tr in [trL, trQ, trT]] + pad=True, fill_value=0.) for tr in [trQ, trT]] - # Noise window (-dts to -5. sec) - noise_left = self.meta.time+self.meta.ttime-dts - noise_right = self.meta.time+self.meta.ttime-5. + # Noise window (-dtsqt to -5. sec) + noise_left, noise_right = _wavelet(trQ, method='noise', overhang=over) # Trim noise traces [tr.trim(noise_left, noise_right, nearest_sample=False, - pad=nn, fill_value=0.) for tr in [trNl, trNq]] + pad=True, fill_value=0.) for tr in [trNl, trNq]] elif phase == 'S' or 'SKS': @@ -850,6 +944,7 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): self.meta.time+self.meta.ttime-dts/2.) # Taper traces - only necessary processing after trimming + # TODO: What does this to the multitaper method [tr.taper(max_percentage=0.05, max_length=2.) for tr in [trL, trQ, trT, trNl, trNq]] From 4898a6c8081ce125c236be605f6ae668a9c0c219 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 25 Oct 2021 15:47:23 -0700 Subject: [PATCH 06/54] rfdata.py: _Pwavelet function: change envelope to minimum-length-envelope --- rfpy/rfdata.py | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index c59a571..c346482 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -641,6 +641,7 @@ def deconvolve(self, phase='P', vp=None, vs=None, Threshold [0-1] used in ``wavelet='envelope'``. time : float Window length used in ``wavelet='time'``. + Minimum window length for ``wavelet='envelope'``. pre_filt : list of 2 floats Low and High frequency corners of bandpass filter applied before deconvolution @@ -679,7 +680,7 @@ def _gauss_filt(dt, nft, f0): gauss[nft21:] = np.flip(gauss[1:nft21-1]) return gauss - def _wavelet(parent, method='complete', overhang=5, + def _Pwavelet(parent, method='complete', overhang=5, envelope_threshold=0.05, time=5): """ @@ -703,10 +704,12 @@ def _wavelet(parent, method='complete', overhang=5, method='envelope') time: float window (seconds) that defines the wavelet (for method='time') + minimum time (seconds) of the wavelet (for method='envelope') Return: left, right: (obspy.UTCDateTime) start and end time of wavelet """ + import obspy.signal.filter as osf if method not in ['complete', 'envelope', 'time', 'noise']: @@ -716,11 +719,12 @@ def _wavelet(parent, method='complete', overhang=5, errflag = False if method == 'envelope': - split = int((self.meta.time + self.meta.ttime - - parent.stats.starttime) * parent.stats.sampling_rate) + split = int((self.meta.time + self.meta.ttime + + time - parent.stats.starttime ) * + parent.stats.sampling_rate) env = osf.envelope(parent.data) - env = env[split:] # only look after P env /= max(env) # normalize + env = env[split:] # only look after P + time try: i = np.nonzero( np.diff( @@ -728,21 +732,20 @@ def _wavelet(parent, method='complete', overhang=5, env > envelope_threshold, dtype=int))==-1)[0][0] except IndexError: i = len(parent.data)-1 - dts = i * parent.stats.delta + dts = i * parent.stats.delta + time left = self.meta.time + self.meta.ttime - overhang right = self.meta.time + self.meta.ttime + dts + overhang if left < parent.stats.starttime or right > parent.stats.endtime: - logging.info('Envelope wavelet longer than trace.') - logging.info('Falling back to "complete" wavelet') + print('Envelope wavelet longer than trace.') + print('Falling back to "complete" wavelet') errflag = True - # TODO: Test me! if method == 'time': left = self.meta.time + self.meta.ttime - overhang right = self.meta.time + self.meta.ttime + time + overhang if left < parent.stats.starttime or right > parent.stats.endtime: - logging.info('Time wavelet longer than trace.') - logging.info('Falling back to "complete" wavelet') + print('Time wavelet longer than trace.') + print('Falling back to "complete" wavelet') errflag = True if method == 'complete' or errflag: @@ -900,25 +903,25 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): over = 5 dtsqt = len(trL.data)*trL.stats.delta/2. - # Traces will be paded to this length (samples) + # Traces will be zero-paded to this length (samples) nn = int(round((dtsqt+over)*trL.stats.sampling_rate)) + 1 - sig_left, sig_right = _wavelet(trL, method='envelope', - envelope_threshold=envelope_threshold, overhang=over) + sig_left, sig_right = _Pwavelet(trL, method=wavelet, + envelope_threshold=envelope_threshold, time=time, overhang=over) # Trim wavelet trL.trim(sig_left, sig_right, nearest_sample=False, pad=True, fill_value=0.) # Signal window (-5. to dtsqt-10 sec) - sig_left, sig_right = _wavelet(trQ, method='complete', overhang=over) + sig_left, sig_right = _Pwavelet(trQ, method='complete', overhang=over) # Trim signal traces [tr.trim(sig_left, sig_right, nearest_sample=False, pad=True, fill_value=0.) for tr in [trQ, trT]] # Noise window (-dtsqt to -5. sec) - noise_left, noise_right = _wavelet(trQ, method='noise', overhang=over) + noise_left, noise_right = _Pwavelet(trQ, method='noise', overhang=over) # Trim noise traces [tr.trim(noise_left, noise_right, nearest_sample=False, From 31e83ea1b3f8a09c85c31a51f4480e1a1b064ebc Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 25 Oct 2021 15:51:02 -0700 Subject: [PATCH 07/54] Revert "More logging funcitonality" This reverts commit b54c6d5c110785d31341ec918de158df11cb9b6a. Remove logging module --- rfpy/rfdata.py | 11 +++--- rfpy/utils.py | 98 +++++++++++++++++++++++++------------------------- 2 files changed, 53 insertions(+), 56 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index c346482..946227d 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -22,7 +22,6 @@ # -*- coding: utf-8 -*- from math import ceil -import logging import pickle import numpy as np from obspy import Trace, Stream, UTCDateTime @@ -116,9 +115,9 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'): source_depth_in_km=self.dep, phase_list=[phase]) if len(arrivals) > 1: - logging.warning("arrival has many entries: ", len(arrivals)) + print("arrival has many entries: ", len(arrivals)) elif len(arrivals) == 0: - logging.warning("no arrival found") + print("no arrival found") self.accept = False return @@ -343,9 +342,9 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, tend = self.meta.time + self.meta.ttime + dts # Get waveforms - logging.info("* Requesting Waveforms: ") - logging.info("* Startime: " + tstart.strftime("%Y-%m-%d %H:%M:%S")) - logging.info("* Endtime: " + tend.strftime("%Y-%m-%d %H:%M:%S")) + print("* Requesting Waveforms: ") + print("* Startime: " + tstart.strftime("%Y-%m-%d %H:%M:%S")) + print("* Endtime: " + tend.strftime("%Y-%m-%d %H:%M:%S")) # Download data err, stream = utils.download_data( diff --git a/rfpy/utils.py b/rfpy/utils.py index d352f7c..f4d5857 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -1,5 +1,4 @@ import math -import logging from obspy import UTCDateTime from numpy import nan, isnan, abs import numpy as np @@ -134,7 +133,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Intialize to default positive error erd = True - logging.info( + print( ("* {0:2s}{1:1s} - Checking Disk".format(sta.channel.upper(), comp.upper()))) @@ -157,11 +156,11 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, s3 = f3.format(styr, stjd, net.upper(), sta.station.upper(), sta.channel.upper()[0:2], comp.upper(), dtype) - logging.debug("* Trying formats:") - logging.debug("* " + s1) - logging.debug("* " + s2) - logging.debug("* " + s3) - logging.debug("* ") + print("* Trying formats:") + print("* " + s1) + print("* " + s2) + print("* " + s3) + print("* ") lclfiles.extend(list(filter(stdata, s1))) lclfiles.extend(list(filter(stdata, s2))) @@ -169,7 +168,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # If still no Local files stop if len(lclfiles) == 0: - logging.warning("* - Data Unavailable") + print("* - Data Unavailable") return erd, None # Process the local Files @@ -207,18 +206,18 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Check for Nan in stream for SAC if True in isnan(st[0].data): - logging.warning( + print( "* !!! Missing Data Present !!! " + "Skipping (NaNs)") # Check for ND Val in stream for MSEED elif -123456789 in st[0].data: - logging.warning( + print( "* !!! Missing Data Present !!! " + "Skipping (MSEED fill)") else: if eddt and (ndval == 0.0): if any(st[0].data == 0.0): - logging.warning( + print( "* !!! Missing Data Present " + "!!! (Set to Zero)") @@ -228,7 +227,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, tloc = "--" # Processed succesfully...Finish - logging.info(("* {1:3s}.{2:2s} - From Disk".format( + print(("* {1:3s}.{2:2s} - From Disk".format( st[0].stats.station, st[0].stats.channel.upper(), tloc))) return False, st @@ -265,7 +264,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # If still no Local files stop if len(lclfiles1) == 0 and len(lclfiles2) == 0: - logging.warning("* - Data Unavailable") + print("* - Data Unavailable") return erd, None # Now try to merge the two separate day files @@ -314,8 +313,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Should only be one component, otherwise keep # reading If more than 1 component, error if len(st) != 1: - logging.error(st) - logging.error("merge failed?") + print(st) + print("merge failed?") else: if (st[0].stats.starttime <= start and @@ -335,18 +334,18 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, # Check for Nan in stream for SAC if True in isnan(st[0].data): - logging.warning( + print( "* !!! Missing Data " + "Present !!! Skipping (NaNs)") # Check for ND Val in stream for MSEED elif -123456789 in st[0].data: - logging.warning( + print( "* !!! Missing Data Present !!! " + "Skipping (MSEED fill)") else: if (eddt1 or eddt2) and (ndval == 0.0): if any(st[0].data == 0.0): - logging.warning( + print( "* !!! Missing " + "Data Present !!! (Set " + "to Zero)") @@ -357,8 +356,7 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, tloc = "--" # Processed succesfully...Finish - logging.info(( - "* {1:3s}.{2:2s} - " + + print(("* {1:3s}.{2:2s} - " + "From Disk".format( st[0].stats.station, st[0].stats.channel.upper(), @@ -369,14 +367,14 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, pass else: st2ot = st2[0].stats.endtime-st2[0].stats.delta - logging.warning("* - Merge Failed: No " + + print("* - Merge Failed: No " + "Overlap {0:s} - {1:s}".format( st1[0].stats.endtime.strftime( "%Y-%m-%d %H:%M:%S"), st2ot.strftime("%Y-%m-%d %H:%M:%S"))) # If we got here, we did not get the data. - logging.warning("* - Data Unavailable") + print("* - Data Unavailable") return erd, None @@ -426,7 +424,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, from math import floor # Output - logging.info(("* {0:s}.{1:2s} - ZNE:".format(sta.station, + print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, sta.channel.upper()))) # Set Error Default to True @@ -477,7 +475,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, channel=channelsZNE, starttime=start, endtime=end+1., attach_response=False) if len(st) == 3: - logging.info("* - ZNE Data Downloaded") + print("* - ZNE Data Downloaded") # It's possible if len(st)==1 that data is Z12 else: @@ -487,7 +485,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, sta.channel.upper() + '2' msg = "* {1:2s}[Z12].{2:2s} - Checking Network".format( sta.station, sta.channel.upper(), tloc) - logging.info(msg) + print(msg) try: st = client.get_waveforms( network=sta.network, @@ -495,16 +493,16 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, channel=channelsZ12, starttime=start, endtime=end+1., attach_response=False) if len(st) == 3: - logging.info("* - Z12 Data Downloaded") + print("* - Z12 Data Downloaded") else: st = None except Exception as e: - logging.warning("* Met exception:") - logging.warning("* " + e.__repr__()) + print("* Met exception:") + print("* " + e.__repr__()) st = None except Exception as e: - logging.warning("* Met exception:") - logging.warning("* " + e.__repr__()) + print("* Met exception:") + print("* " + e.__repr__()) st = None # Break if we successfully obtained 3 components in st @@ -514,8 +512,8 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Check the correct 3 components exist if st is None: - logging.warning("* Error retrieving waveforms") - logging.warning("**************************************************") + print("* Error retrieving waveforms") + print("**************************************************") return True, None # Three components successfully retrieved @@ -527,12 +525,12 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Check start times if not np.all([tr.stats.starttime == start for tr in st]): - logging.warning("* Start times are not all close to true start: ") - [logging.warning("* "+tr.stats.channel+" " + + print("* Start times are not all close to true start: ") + [print("* "+tr.stats.channel+" " + str(tr.stats.starttime)+" " + str(tr.stats.endtime)) for tr in st] - logging.warning("* True start: "+str(start)) - logging.warning("* -> Shifting traces to true start") + print("* True start: "+str(start)) + print("* -> Shifting traces to true start") delay = [tr.stats.starttime - start for tr in st] st_shifted = Stream( traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)]) @@ -542,39 +540,39 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, sr = st[0].stats.sampling_rate sr_round = float(floor_decimal(sr, 0)) if not sr == sr_round: - logging.warning("* Sampling rate is not an integer value: ", sr) - logging.warning("* -> Resampling") + print("* Sampling rate is not an integer value: ", sr) + print("* -> Resampling") st.resample(sr_round, no_filter=False) # Try trimming try: st.trim(start, end) except: - logging.error("* Unable to trim") - logging.error("* -> Skipping") - logging.error("**************************************************") + print("* Unable to trim") + print("* -> Skipping") + print("**************************************************") return True, None # Check final lengths - they should all be equal if start times # and sampling rates are all equal and traces have been trimmed if not np.allclose([tr.stats.npts for tr in st[1:]], st[0].stats.npts): - logging.error("* Lengths are incompatible: ") - [logging.error("* "+str(tr.stats.npts)) for tr in st] - logging.error("* -> Skipping") - logging.error("**************************************************") + print("* Lengths are incompatible: ") + [print("* "+str(tr.stats.npts)) for tr in st] + print("* -> Skipping") + print("**************************************************") return True, None elif not np.allclose([st[0].stats.npts], int((end - start)*sr), atol=1): - logging.error("* Length is too short: ") - logging.error("* "+str(st[0].stats.npts) + + print("* Length is too short: ") + print("* "+str(st[0].stats.npts) + " ~= "+str(int((end - start)*sr))) - logging.error("* -> Skipping") - logging.error("**************************************************") + print("* -> Skipping") + print("**************************************************") return True, None else: - logging.info("* Waveforms Retrieved...") + print("* Waveforms Retrieved...") return False, st From ad8f5b42ee6e56e76ae8a2dff369c1d00da9e404 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Tue, 26 Oct 2021 10:02:09 -0700 Subject: [PATCH 08/54] binning: updated documentation --- rfpy/binning.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 4ed72b0..9899e67 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -45,10 +45,13 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False): a single trace. stream2 : :class:`~obspy.core.Stream` Optionally stack a second stream in the same operation. - dbaz : int - Number of bazk-azimuth samples in bins - dslow : int - Number of slowness samples in bins + typ: str + Attribute to bin + 'baz': backazimuth (degree) + 'slow': Horizontal slowness (s/km) + 'dist': epicentral distance (degree) + nbin : int + Number of equally sized bins within data range pws : bool Whether or not to perform phase-weighted stacking @@ -157,9 +160,9 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): a single trace. stream2 : :class:`~obspy.core.Stream` Optionally stack a second stream in the same operation. - dbaz : int + nbaz : int Number of bazk-azimuth samples in bins - dslow : int + nslow : int Number of slowness samples in bins pws : bool Whether or not to perform phase-weighted stacking @@ -236,6 +239,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): return final_stream + def bin_all(stream1, stream2=None, pws=False): """ Function to bin all streams into a single trace. From da613882fbb571c1a942da54d5ef105dc51bad25 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 15 Oct 2021 13:28:27 -0400 Subject: [PATCH 09/54] fixed setup for scripts package --- setup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 5fdb266..c9f78dc 100644 --- a/setup.py +++ b/setup.py @@ -1,3 +1,4 @@ +import setuptools import os.path from os import listdir import re @@ -33,7 +34,7 @@ def find_version(*paths): 'Programming Language :: Python :: 3.9'], install_requires=['numpy', 'obspy', 'stdb>=0.2.0'], python_requires='>=3.6', - packages=['rfpy'], + packages=setuptools.find_packages(), entry_points={'console_scripts':[ 'rfpy_calc=rfpy.scripts.rfpy_calc:main', 'rfpy_recalc=rfpy.scripts.rfpy_recalc:main', From 448e51018c50bae3eadc851725ffd8c237b41f1a Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Tue, 26 Oct 2021 10:02:09 -0700 Subject: [PATCH 10/54] binning: updated documentation --- rfpy/binning.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 4ed72b0..9899e67 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -45,10 +45,13 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False): a single trace. stream2 : :class:`~obspy.core.Stream` Optionally stack a second stream in the same operation. - dbaz : int - Number of bazk-azimuth samples in bins - dslow : int - Number of slowness samples in bins + typ: str + Attribute to bin + 'baz': backazimuth (degree) + 'slow': Horizontal slowness (s/km) + 'dist': epicentral distance (degree) + nbin : int + Number of equally sized bins within data range pws : bool Whether or not to perform phase-weighted stacking @@ -157,9 +160,9 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): a single trace. stream2 : :class:`~obspy.core.Stream` Optionally stack a second stream in the same operation. - dbaz : int + nbaz : int Number of bazk-azimuth samples in bins - dslow : int + nslow : int Number of slowness samples in bins pws : bool Whether or not to perform phase-weighted stacking @@ -236,6 +239,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): return final_stream + def bin_all(stream1, stream2=None, pws=False): """ Function to bin all streams into a single trace. From ef558b0aa516b2549a9c2376ac31485bc43eac19 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Tue, 26 Oct 2021 10:02:09 -0700 Subject: [PATCH 11/54] binning: updated documentation --- rfpy/binning.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 9899e67..2a3eda7 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -161,7 +161,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): stream2 : :class:`~obspy.core.Stream` Optionally stack a second stream in the same operation. nbaz : int - Number of bazk-azimuth samples in bins + Number of back-azimuth samples in bins nslow : int Number of slowness samples in bins pws : bool From 8129104ee440e49daa317cc5684c4c584d292382 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 1 Nov 2021 16:03:50 -0700 Subject: [PATCH 12/54] Meta: initialize floats with np.nan, instead of None --- rfpy/rfdata.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 946227d..3f3f12f 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -130,9 +130,9 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'): self.phase = phase self.accept = True else: - self.ttime = None - self.slow = None - self.inc = None + self.ttime = np.nan + self.slow = np.nan + self.inc = np.nan self.phase = None self.accept = False @@ -143,9 +143,9 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'): # Attributes that get updated as analysis progresses self.rotated = False - self.snr = None - self.snrh = None - self.cc = None + self.snr = np.nan + self.snrh = np.nan + self.cc = np.nan class RFData(object): @@ -551,9 +551,10 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.): if not self.meta.accept: return - if self.meta.snr: - print("SNR already calculated - continuing") - return + if self.meta.snr: # False if None (for backward compatibility) + if np.isfinite(self.meta.snr): # False if nan (None throws error) + print("SNR already calculated - continuing") + return t1 = self.meta.time + self.meta.ttime @@ -869,7 +870,8 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): print("Warning: Data have not been rotated yet - rotating now") self.rotate(vp=vp, vs=vs, align=align) - if not self.meta.snr: + # v--True if None v--True if nan, error if None + if not self.meta.snr or not np.isfinite(self.meta.snr): print("Warning: SNR has not been calculated - " + "calculating now using default") self.calc_snr() From c90d31049bc2551eab27ea042e2f51349915335d Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 1 Nov 2021 16:58:29 -0700 Subject: [PATCH 13/54] binning: Added option 'include_empty' to return streams that have no data. Useful for plotting data gaps. --- rfpy/binning.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 2a3eda7..dba288e 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -32,7 +32,8 @@ from scipy.signal import hilbert -def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False): +def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, + include_empty=False): """ Function to stack receiver functions into (baz or slow) bins This can be done using a linear stack (i.e., simple @@ -54,6 +55,8 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False): Number of equally sized bins within data range pws : bool Whether or not to perform phase-weighted stacking + include_empty : bool + Return empty bins as null traces (default omits them) Returns ------- @@ -113,7 +116,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False): continue - if nb > 0: + if nb > 0 or include_empty: # Average and update stats array /= nb @@ -147,7 +150,8 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False): return final_stream -def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): +def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, + include_empty=False): """ Function to stack receiver functions into back-azimuth and slowness bins. This can be done using a linear stack (i.e., simple @@ -166,6 +170,8 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): Number of slowness samples in bins pws : bool Whether or not to perform phase-weighted stacking + include_empty : bool + Return empty bins as null traces (default omits them) Returns ------- @@ -216,7 +222,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False): continue - if nbin > 0: + if nbin > 0 or include_empty: # Average and update stats array /= nbin From 052f96ecd09bcccf0ea0b49041e2d9d708b416c7 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Fri, 5 Nov 2021 08:48:42 -0700 Subject: [PATCH 14/54] utils: Print download exception only if verbose --- rfpy/utils.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/rfpy/utils.py b/rfpy/utils.py index f4d5857..2e9f615 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -497,12 +497,14 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, else: st = None except Exception as e: - print("* Met exception:") - print("* " + e.__repr__()) + if verbose: + print("* Met exception:") + print("* " + e.__repr__()) st = None except Exception as e: - print("* Met exception:") - print("* " + e.__repr__()) + if verbose: + print("* Met exception:") + print("* " + e.__repr__()) st = None # Break if we successfully obtained 3 components in st From 171fb0179d8d6360301a1a66732b103fa46b1cb4 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Fri, 5 Nov 2021 13:03:41 -0700 Subject: [PATCH 15/54] rfdata, utils: added option to remove instrument response from donwloaded data --- rfpy/rfdata.py | 10 ++- rfpy/utils.py | 203 ++++++++++++++++++++++++------------------------- 2 files changed, 108 insertions(+), 105 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 3f3f12f..49cb98b 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -299,7 +299,8 @@ def add_data(self, stream, returned=False, new_sr=5.): return self.meta.accept def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, - new_sr=5.,dts=120., returned=False, verbose=False): + new_sr=5.,dts=120., remove_response=False, + returned=False, verbose=False): """ Downloads seismograms based on event origin time and P phase arrival. @@ -316,8 +317,13 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, Time duration (sec) stdata : List Station list + remove_response : bool + Remove instrument response from seismogram and resitute to true ground + velocity (m/s) using obspy.core.trace.Trace.remove_response() returned : bool Whether or not to return the ``accept`` attribute + verbose : bool + Output diagnostics to screen Returns ------- @@ -350,7 +356,7 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, err, stream = utils.download_data( client=client, sta=self.sta, start=tstart, end=tend, stdata=stdata, dtype=dtype, ndval=ndval, new_sr=new_sr, - verbose=verbose) + remove_response=remove_response, verbose=verbose) # Store as attributes with traces in dictionary try: diff --git a/rfpy/utils.py b/rfpy/utils.py index 2e9f615..1bc4928 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -379,7 +379,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, - stdata=[], dtype='SAC', ndval=nan, new_sr=0., verbose=False): + stdata=[], dtype='SAC', ndval=nan, new_sr=0., verbose=False, + remove_response=False): """ Function to build a stream object for a seismogram in a given time window either by downloading data from the client object or alternatively first checking if the @@ -403,6 +404,9 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, Station list ndval : float or nan Default value for missing data + remove_response : bool + Remove instrument response from seismogram and resitute to true ground + velocity (m/s) using obspy.core.trace.Trace.remove_response() Returns ------- @@ -460,56 +464,45 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Construct location name if len(tloc) == 0: tloc = "--" + # Construct Channel List - channelsZNE = sta.channel.upper() + 'Z,' + sta.channel.upper() + \ - 'N,' + sta.channel.upper() + 'E' - print(("* {1:2s}[ZNE].{2:2s} - Checking Network".format( - sta.station, sta.channel.upper(), tloc))) - - # Get waveforms, with extra 1 second to avoid - # traces cropped too short - traces are trimmed later - try: - st = client.get_waveforms( - network=sta.network, - station=sta.station, location=loc, - channel=channelsZNE, starttime=start, - endtime=end+1., attach_response=False) - if len(st) == 3: - print("* - ZNE Data Downloaded") - - # It's possible if len(st)==1 that data is Z12 + chaZNE = (sta.channel.upper() + 'Z,' + + sta.channel.upper() + 'N,' + + sta.channel.upper() + 'E') + msgZNE = "* {1:2s}[ZNE].{2:2s} - Checking Network".format( + sta.station, sta.channel.upper(), tloc) + + chaZ12 = (sta.channel.upper() + 'Z,' + + sta.channel.upper() + '1,' + + sta.channel.upper() + '2') + msgZ12 = "* {1:2s}[Z12].{2:2s} - Checking Network".format( + sta.station, sta.channel.upper(), tloc) + + # Loop over possible channels + for channel, msg in zip([chaZNE, chaZ12], [msgZNE, msgZ12]): + print(msg) + + # Get waveforms, with extra 1 second to avoid + # traces cropped too short - traces are trimmed later + try: + st = client.get_waveforms( + network=sta.network, + station=sta.station, location=loc, + channel=channel, starttime=start, + endtime=end+1., attach_response=remove_response) + except Exception as e: + if verbose: + print("* Met exception:") + print("* " + e.__repr__()) + st = None else: - # Construct Channel List - channelsZ12 = sta.channel.upper() + 'Z,' + \ - sta.channel.upper() + '1,' + \ - sta.channel.upper() + '2' - msg = "* {1:2s}[Z12].{2:2s} - Checking Network".format( - sta.station, sta.channel.upper(), tloc) - print(msg) - try: - st = client.get_waveforms( - network=sta.network, - station=sta.station, location=loc, - channel=channelsZ12, starttime=start, - endtime=end+1., attach_response=False) - if len(st) == 3: - print("* - Z12 Data Downloaded") - else: - st = None - except Exception as e: - if verbose: - print("* Met exception:") - print("* " + e.__repr__()) - st = None - except Exception as e: - if verbose: - print("* Met exception:") - print("* " + e.__repr__()) - st = None + if len(st) == 3: + # It's possible if len(st)==1 that data is Z12 + print("* - ZNE Data Downloaded") + break # Break if we successfully obtained 3 components in st if not erd: - break # Check the correct 3 components exist @@ -519,62 +512,66 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, return True, None # Three components successfully retrieved - else: + if remove_response: + st.remove_response() + print("*") + print("* Restituted stream to true ground velocity.") + print("*") + + # Detrend and apply taper + st.detrend('demean').detrend('linear').taper( + max_percentage=0.05, max_length=5.) + + # Check start times + if not np.all([tr.stats.starttime == start for tr in st]): + print("* Start times are not all close to true start: ") + [print("* "+tr.stats.channel+" " + + str(tr.stats.starttime)+" " + + str(tr.stats.endtime)) for tr in st] + print("* True start: "+str(start)) + print("* -> Shifting traces to true start") + delay = [tr.stats.starttime - start for tr in st] + st_shifted = Stream( + traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)]) + st = st_shifted.copy() + + # Check sampling rate + sr = st[0].stats.sampling_rate + sr_round = float(floor_decimal(sr, 0)) + if not sr == sr_round: + print("* Sampling rate is not an integer value: ", sr) + print("* -> Resampling") + st.resample(sr_round, no_filter=False) + + # Try trimming + try: + st.trim(start, end) + except: + print("* Unable to trim") + print("* -> Skipping") + print("**************************************************") + return True, None - # Detrend and apply taper - st.detrend('demean').detrend('linear').taper( - max_percentage=0.05, max_length=5.) - - # Check start times - if not np.all([tr.stats.starttime == start for tr in st]): - print("* Start times are not all close to true start: ") - [print("* "+tr.stats.channel+" " + - str(tr.stats.starttime)+" " + - str(tr.stats.endtime)) for tr in st] - print("* True start: "+str(start)) - print("* -> Shifting traces to true start") - delay = [tr.stats.starttime - start for tr in st] - st_shifted = Stream( - traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)]) - st = st_shifted.copy() - - # Check sampling rate - sr = st[0].stats.sampling_rate - sr_round = float(floor_decimal(sr, 0)) - if not sr == sr_round: - print("* Sampling rate is not an integer value: ", sr) - print("* -> Resampling") - st.resample(sr_round, no_filter=False) - - # Try trimming - try: - st.trim(start, end) - except: - print("* Unable to trim") - print("* -> Skipping") - print("**************************************************") - return True, None - - # Check final lengths - they should all be equal if start times - # and sampling rates are all equal and traces have been trimmed - if not np.allclose([tr.stats.npts for tr in st[1:]], st[0].stats.npts): - print("* Lengths are incompatible: ") - [print("* "+str(tr.stats.npts)) for tr in st] - print("* -> Skipping") - print("**************************************************") - - return True, None - - elif not np.allclose([st[0].stats.npts], int((end - start)*sr), - atol=1): - print("* Length is too short: ") - print("* "+str(st[0].stats.npts) + - " ~= "+str(int((end - start)*sr))) - print("* -> Skipping") - print("**************************************************") - - return True, None + # Check final lengths - they should all be equal if start times + # and sampling rates are all equal and traces have been trimmed + if not np.allclose([tr.stats.npts for tr in st[1:]], st[0].stats.npts): + print("* Lengths are incompatible: ") + [print("* "+str(tr.stats.npts)) for tr in st] + print("* -> Skipping") + print("**************************************************") - else: - print("* Waveforms Retrieved...") - return False, st + return True, None + + elif not np.allclose([st[0].stats.npts], int((end - start)*sr), + atol=1): + print("* Length is too short: ") + print("* "+str(st[0].stats.npts) + + " ~= "+str(int((end - start)*sr))) + print("* -> Skipping") + print("**************************************************") + + return True, None + + else: + print("* Waveforms Retrieved...") + return False, st From 86e5e84258eb6fabb5ab639d061a84af00b87284 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Wed, 10 Nov 2021 10:34:56 -0800 Subject: [PATCH 16/54] Option to use locally saved response information --- rfpy/rfdata.py | 5 ++++- rfpy/utils.py | 40 +++++++++++++++++++++++++++++++++++++--- 2 files changed, 41 insertions(+), 4 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 49cb98b..20c06cf 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -300,6 +300,7 @@ def add_data(self, stream, returned=False, new_sr=5.): def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, new_sr=5.,dts=120., remove_response=False, + local_response_dir='', returned=False, verbose=False): """ Downloads seismograms based on event origin time and @@ -356,7 +357,9 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, err, stream = utils.download_data( client=client, sta=self.sta, start=tstart, end=tend, stdata=stdata, dtype=dtype, ndval=ndval, new_sr=new_sr, - remove_response=remove_response, verbose=verbose) + remove_response=remove_response, + local_response_dir=local_response_dir, + verbose=verbose) # Store as attributes with traces in dictionary try: diff --git a/rfpy/utils.py b/rfpy/utils.py index 1bc4928..b77f709 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -2,7 +2,8 @@ from obspy import UTCDateTime from numpy import nan, isnan, abs import numpy as np -from obspy.core import Stream, read +from obspy import Stream, Inventory +from obspy import read, read_inventory def floor_decimal(n, decimals=0): @@ -378,9 +379,29 @@ def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, return erd, None +def attach_local_response_data(stream, local_response_dir): + """ + Function to restrieve response data from locally stored dataless seed + + Parameters + ---------- + stream : obspy.Stream + Event seismogram + local_response_dir: str + Directory holding response information. All files containing the station + name are read. + """ + stations = [t.stats.station for t in stream.traces] + inventory = Inventory() + for station in stations: + inventory += read_inventory( + '{:}/*{:}*'.format(local_response_dir, station)) + stream.attach_response(inventories=inventory) + + def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, stdata=[], dtype='SAC', ndval=nan, new_sr=0., verbose=False, - remove_response=False): + remove_response=False, local_response_dir=''): """ Function to build a stream object for a seismogram in a given time window either by downloading data from the client object or alternatively first checking if the @@ -401,12 +422,16 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, end : :class:`~obspy.core.utcdatetime.UTCDateTime` End time for request stdata : List - Station list + List of directories holding local waveform data ndval : float or nan Default value for missing data remove_response : bool Remove instrument response from seismogram and resitute to true ground velocity (m/s) using obspy.core.trace.Trace.remove_response() + local_response_dir: str + Directory holding response files to be read by obspy.read_inventory(). + Required when ``remove_response`` and using locally stored waveform data + via ``stdata``. Returns ------- @@ -454,6 +479,15 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, if not erd: # Combine Data st = stZ + stN + stE + if remove_response: + if local_response_dir: + attach_local_response_data(st, local_response_dir) + else: + print('*') + print('* Warning: No local_response_dir given.') + print('* Warning: Continuing without.') + print('*') + # No local data? Request using client if erd: From 8c442de064ee8542fadf7109d1a6213e145fb307 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 15 Nov 2021 12:51:28 -0800 Subject: [PATCH 17/54] download: Only prefilter traces when downsampling. Add some verbosity. --- rfpy/rfdata.py | 45 ++++++++++++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 20c06cf..cd420c1 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -338,6 +338,16 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, """ + def _resample(tr, new_sr): + + if tr.stats.sampling_rate > new_sr: + # only filter when downsampling + tr.filter('lowpass', freq=0.5*new_sr, corners=2, zerophase=True) + + tr.resample(new_sr, no_filter=True) + return tr + + if self.meta is None: raise(Exception("Requires event data as attribute - aborting")) @@ -366,26 +376,28 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, trE = stream.select(component='E')[0] trN = stream.select(component='N')[0] trZ = stream.select(component='Z')[0] - self.data = Stream(traces=[trZ, trN, trE]) - # Filter Traces and resample - self.data.filter('lowpass', freq=0.5*new_sr, - corners=2, zerophase=True) - self.data.resample(new_sr, no_filter=True) + for tr in [trE, trN, trZ]: + tr = _resample(tr, new_sr) + + self.data = Stream(traces=[trZ, trN, trE]) # If there is no ZNE, perhaps there is Z12? - except: + except Exception as e: + if verbose: + print('* Met exception: ') + print('* ' + str(e)) + print('* Trying to access Z12 data. ') try: tr1 = stream.select(component='1')[0] tr2 = stream.select(component='2')[0] trZ = stream.select(component='Z')[0] - self.data = Stream(traces=[trZ, tr1, tr2]) - # Filter Traces and resample - self.data.filter('lowpass', freq=0.5*new_sr, - corners=2, zerophase=True) - self.data.resample(new_sr, no_filter=True) + for tr in [tr1, tr2, trZ]: + tr = _resample(tr, new_sr) + + self.data = Stream(traces=[trZ, tr1, tr2]) # Save Z12 components in case it's necessary for later self.dataZ12 = self.data.copy() @@ -393,7 +405,11 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, # Rotate from Z12 to ZNE using StDb azcorr attribute self.rotate(align='ZNE') - except: + except Exception as e: + if verbose: + print('* Met exception: ') + print('* ' + str(e)) + print('* Cannot process ') self.meta.accept = False if returned: @@ -855,11 +871,6 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): gauss = np.ones(npad) gnorm = 1. - # Is this correct? - #parent de-noised = np.fft.ifftshift(np.real(np.fft.ifft(gauss*Sdenom))/gnorm) - #daughter1 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd1p))/gnorm) - #daughter2 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd2p))/gnorm) - # Copy traces rfp = parent.copy() rfd1 = daughter1.copy() From 5257cadabdcbb337d145cf54874bc99c6ac9b472 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 17 Nov 2021 10:47:55 -0500 Subject: [PATCH 18/54] First big commit for simdec --- rfpy/binning.py | 83 +-- rfpy/rfdata.py | 282 +++++--- rfpy/rfstreams.py | 560 ++++++++++++++++ rfpy/scripts/rfpy_calc.py | 15 +- rfpy/scripts/rfpy_calc_simdec.py | 481 +++++++++++++ .../{rfpy_recalc.py => rfpy_calc_sse.py} | 0 rfpy/scripts/rfpy_get_data.py | 631 ++++++++++++++++++ setup.py | 6 +- 8 files changed, 1907 insertions(+), 151 deletions(-) create mode 100644 rfpy/rfstreams.py create mode 100644 rfpy/scripts/rfpy_calc_simdec.py rename rfpy/scripts/{rfpy_recalc.py => rfpy_calc_sse.py} (100%) create mode 100644 rfpy/scripts/rfpy_get_data.py diff --git a/rfpy/binning.py b/rfpy/binning.py index dba288e..0aa6530 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -194,54 +194,59 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, islow = np.digitize(slow, slow_bins) final_stream = [] + if stream2: + stream1 = [stream1, stream2] + else: + stream1 = [stream1] + + for stream in stream1: + # try: + # Define empty streams + binned_stream = Stream() + + # Loop through baz_bins + for i in range(nbaz): + for j in range(nslow): + + nbin = 0 + array = np.zeros(len(stream[0].data), dtype=type(stream[0].data)) + weight = np.zeros(len(stream[0].data), dtype=complex) - for stream in [stream1, stream2]: - try: - # Define empty streams - binned_stream = Stream() - - # Loop through baz_bins - for i in range(nbaz): - for j in range(nslow): - - nbin = 0 - array = np.zeros(len(stream[0].data)) - weight = np.zeros(len(stream[0].data), dtype=complex) - - # Loop all traces - for k, tr in enumerate(stream): + # Loop all traces + for k, tr in enumerate(stream): - # If index of baz_bins is equal to ibaz - if i == ibaz[k] and j == islow[k]: + # If index of baz_bins is equal to ibaz + if i == ibaz[k] and j == islow[k]: - nbin += 1 - array += tr.data - hilb = hilbert(tr.data) - phase = np.arctan2(hilb.imag, hilb.real) - weight += np.exp(1j*phase) - - continue + nbin += 1 + array += tr.data + hilb = hilbert(np.real(tr.data)) + phase = np.arctan2(hilb.imag, hilb.real) + weight += np.exp(1j*phase) + + continue - if nbin > 0 or include_empty: + if nbin > 0 or include_empty: - # Average and update stats - array /= nbin - weight = np.real(abs(weight/nbin)) + # Average and update stats + array /= nbin + weight = np.real(abs(weight/nbin)) - trace = Trace(header=stream[0].stats) - trace.stats.baz = baz_bins[i] - trace.stats.slow = slow_bins[j] - trace.stats.nbin = nbin + trace = Trace(header=stream[0].stats) + trace.stats.baz = baz_bins[i] + trace.stats.slow = slow_bins[j] + trace.stats.nbin = nbin - if not pws: - weight = np.ones(len(stream[0].data)) - trace.data = weight*array - binned_stream.append(trace) + if not pws: + weight = np.ones(len(stream[0].data)) + trace.data = weight*array + binned_stream.append(trace) - final_stream.append(binned_stream) + final_stream.append(binned_stream) - except: - continue + # except: + # final_stream = [Stream()] + # continue return final_stream diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index cd420c1..9c393ec 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -152,7 +152,7 @@ class RFData(object): """ A RFData object contains Class attributes that associate station information with a single event (i.e., earthquake) - metadata, corresponding raw and rotated seismograms and + metadata, corresponding raw and rotated seismograms, spectra and receiver functions. Note @@ -338,16 +338,6 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, """ - def _resample(tr, new_sr): - - if tr.stats.sampling_rate > new_sr: - # only filter when downsampling - tr.filter('lowpass', freq=0.5*new_sr, corners=2, zerophase=True) - - tr.resample(new_sr, no_filter=True) - return tr - - if self.meta is None: raise(Exception("Requires event data as attribute - aborting")) @@ -376,40 +366,34 @@ def _resample(tr, new_sr): trE = stream.select(component='E')[0] trN = stream.select(component='N')[0] trZ = stream.select(component='Z')[0] - - for tr in [trE, trN, trZ]: - tr = _resample(tr, new_sr) - self.data = Stream(traces=[trZ, trN, trE]) + # Filter Traces and resample + self.data.filter('lowpass', freq=0.5*new_sr, + corners=2, zerophase=True) + self.data.resample(new_sr, no_filter=True) + # If there is no ZNE, perhaps there is Z12? - except Exception as e: - if verbose: - print('* Met exception: ') - print('* ' + str(e)) - print('* Trying to access Z12 data. ') + except: try: tr1 = stream.select(component='1')[0] tr2 = stream.select(component='2')[0] trZ = stream.select(component='Z')[0] - - for tr in [tr1, tr2, trZ]: - tr = _resample(tr, new_sr) - self.data = Stream(traces=[trZ, tr1, tr2]) + # Filter Traces and resample + self.data.filter('lowpass', freq=0.5*new_sr, + corners=2, zerophase=True) + self.data.resample(new_sr, no_filter=True) + # Save Z12 components in case it's necessary for later self.dataZ12 = self.data.copy() # Rotate from Z12 to ZNE using StDb azcorr attribute self.rotate(align='ZNE') - except Exception as e: - if verbose: - print('* Met exception: ') - print('* ' + str(e)) - print('* Cannot process ') + except: self.meta.accept = False if returned: @@ -637,13 +621,13 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.): # Calculate signal/noise ratio in dB self.meta.snrh = 10*np.log10(srms*srms/nrms/nrms) - - def deconvolve(self, phase='P', vp=None, vs=None, + def calc_spectra(self, phase='P', vp=None, vs=None, align=None, method='wiener', wavelet='complete', - envelope_threshold=0.05, time=5, pre_filt=None, - gfilt=None, wlevel=0.01, writeto=None): + envelope_threshold=0.05, time=5., pre_filt=None, + writeto=False): """ - Deconvolves three-component data using one component as the source wavelet. + Calculate the numerators and denominator of the spectral division for receiver + function calculation using one component as the source wavelet. The source component is always taken as the dominant compressional component, which can be either 'Z', 'L', or 'P'. @@ -679,8 +663,8 @@ def deconvolve(self, phase='P', vp=None, vs=None, Attributes ---------- - rf : :class:`~obspy.core.Stream` - Stream containing the receiver function traces + specs : :class:`~obspy.core.Stream` + Stream containing the cross spectral quantities """ @@ -705,8 +689,8 @@ def _gauss_filt(dt, nft, f0): gauss[nft21:] = np.flip(gauss[1:nft21-1]) return gauss - def _Pwavelet(parent, method='complete', overhang=5, - envelope_threshold=0.05, time=5): + def _Pwavelet(parent, method='complete', overhang=5., + envelope_threshold=0.05, time=5.): """ Select wavelet from the parent function for deconvolution using method. @@ -785,7 +769,7 @@ def _Pwavelet(parent, method='complete', overhang=5, return left, right - def _decon(parent, daughter1, daughter2, noise, nn, method): + def _calc_specs(parent, daughter1, daughter2, noise, nn, method): # Get length, zero padding parameters and frequencies dt = parent.stats.delta @@ -809,14 +793,6 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): Sd2p = Fd2*np.conjugate(Fp) Snn = np.real(Fn*np.conjugate(Fn)) - # Final processing depends on method - if method == 'wiener': - Sdenom = Spp + Snn - elif method == 'water': - phi = np.amax(Spp)*wlevel - Sdenom = Spp - Sdenom[Sdenom < phi] = phi - # Multitaper deconvolution elif method == 'multitaper': @@ -856,35 +832,12 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): Sd2p = np.sum(Fd2*np.conjugate(Fp), axis=0) Snn = np.sum(np.real(Fn*np.conjugate(Fn)), axis=0) - # Denominator - Sdenom = Spp + Snn - else: print("Method not implemented") pass - # Apply Gaussian filter? - if gfilt: - gauss = _gauss_filt(dt, npad, gfilt) - gnorm = np.sum(gauss)*(freqs[1]-freqs[0])*dt - else: - gauss = np.ones(npad) - gnorm = 1. - - # Copy traces - rfp = parent.copy() - rfd1 = daughter1.copy() - rfd2 = daughter2.copy() - - # Spectral division and inverse transform - rfp.data = np.fft.ifftshift(np.real(np.fft.ifft( - gauss*Spp/Sdenom))/gnorm) - rfd1.data = np.fft.ifftshift(np.real(np.fft.ifft( - gauss*Sd1p/Sdenom))/gnorm) - rfd2.data = np.fft.ifftshift(np.real(np.fft.ifft( - gauss*Sd2p/Sdenom))/gnorm) + return Spp, Sd1p, Sd2p, Snn - return rfp, rfd1, rfd2 if not self.meta.rotated: print("Warning: Data have not been rotated yet - rotating now") @@ -948,31 +901,32 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): [tr.trim(noise_left, noise_right, nearest_sample=False, pad=True, fill_value=0.) for tr in [trNl, trNq]] - elif phase == 'S' or 'SKS': + # Can't deal with S phases right now #### + # elif phase == 'S' or 'SKS': - # Get signal length (i.e., seismogram to deconvolve) from trace length - dts = len(trL.data)*trL.stats.delta/2. - - # Trim signal traces (-5. to dts-10 sec) - trL.trim(self.meta.time+self.meta.ttime+25.-dts/2., - self.meta.time+self.meta.ttime+25.) - trQ.trim(self.meta.time+self.meta.ttime+25.-dts/2., - self.meta.time+self.meta.ttime+25.) - trT.trim(self.meta.time+self.meta.ttime+25.-dts/2., - self.meta.time+self.meta.ttime+25.) + # # Get signal length (i.e., seismogram to deconvolve) from trace length + # dts = len(trL.data)*trL.stats.delta/2. + + # # Trim signal traces (-5. to dts-10 sec) + # trL.trim(self.meta.time+self.meta.ttime+25.-dts/2., + # self.meta.time+self.meta.ttime+25.) + # trQ.trim(self.meta.time+self.meta.ttime+25.-dts/2., + # self.meta.time+self.meta.ttime+25.) + # trT.trim(self.meta.time+self.meta.ttime+25.-dts/2., + # self.meta.time+self.meta.ttime+25.) - # Trim noise traces (-dts to -5 sec) - trNl.trim(self.meta.time+self.meta.ttime-dts, - self.meta.time+self.meta.ttime-dts/2.) - trNq.trim(self.meta.time+self.meta.ttime-dts, - self.meta.time+self.meta.ttime-dts/2.) + # # Trim noise traces (-dts to -5 sec) + # trNl.trim(self.meta.time+self.meta.ttime-dts, + # self.meta.time+self.meta.ttime-dts/2.) + # trNq.trim(self.meta.time+self.meta.ttime-dts, + # self.meta.time+self.meta.ttime-dts/2.) # Taper traces - only necessary processing after trimming # TODO: What does this to the multitaper method [tr.taper(max_percentage=0.05, max_length=2.) for tr in [trL, trQ, trT, trNl, trNq]] - # Pre-filter waveforms before deconvolution + # Pre-filter waveforms before calculating spectra if pre_filt: [tr.filter('bandpass', freqmin=pre_filt[0], freqmax=pre_filt[1], corners=2, zerophase=True) @@ -982,17 +936,114 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): with open(writeto, 'wb') as f: pickle.dump(Stream(traces=[trL, trQ, trT]), f) - # Deconvolve + # Calculate spectra if phase == 'P' or 'PP': - rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method) + spp, sd1p, sd2p, snn = _calc_specs(trL, trQ, trT, trNl, nn, method) + + # Can't deal with S phases right now #### + # elif phase == 'S' or 'SKS': + # rfQ, rfL, rfT = _calc_specs(trQ, trL, trT, trNq, nn, method) + + # Copy traces + sLL = trL.copy() + sQL = trQ.copy() + sTL = trT.copy() + sNN = trNl.copy() + sLL.data = spp + sQL.data = sd1p + sTL.data = sd2p + sNN.data = snn + + # Update stats of streams + sLL.stats.channel = 'S' + self.meta.align[0] + self.meta.align[0] + sQL.stats.channel = 'S' + self.meta.align[1] + self.meta.align[0] + sTL.stats.channel = 'S' + self.meta.align[2] + self.meta.align[0] + sNN.stats.channel = 'SNN' + + self.specs = Stream(traces=[sLL, sQL, sTL, sNN]) - elif phase == 'S' or 'SKS': - rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method) + + def deconvolve(self, align=None, method='wiener', + gfilt=None, wlevel=0.01, writeto=None): + """ + Deconvolves three-component data using one component as the source wavelet. + The source component is always taken as the dominant compressional + component, which can be either 'Z', 'L', or 'P'. + + Parameters + ---------- + align : str + Alignment of coordinate system for rotation + ('ZRT', 'LQT', or 'PVH') + method : str + Method for deconvolution. Options are 'wiener', 'water' or + 'multitaper' + gfilt : float + Center frequency of Gaussian filter (Hz). + wlevel : float + Water level used in ``method='water'``. + writeto : str or None + Write wavelets for deconvolution to file. + + Attributes + ---------- + rf : :class:`~obspy.core.Stream` + Stream containing the receiver function traces + + """ + + try: + len(self.specs) > 0 + except: + Exception("spectra have not been calculated") + + + # Make everything explicit + SLL = self.specs[0].copy() + SQL = self.specs[1].copy() + STL = self.specs[2].copy() + SNN = self.specs[3].copy() + + dt = SLL.stats.delta + npad = SLL.stats.npts + freqs = np.fft.fftfreq(npad, d=dt) + + # Wiener or Water level deconvolution + if method in ['wiener', 'multitaper']: + + # Denominator (Spp + Snn) + Sdenom = SLL.data + SNN.data + + elif method == 'water': + phi = np.amax(SLL)*wlevel + Sdenom = SLL.data + Sdenom[Sdenom < phi] = phi + + # Apply Gaussian filter? + if gfilt: + gauss = _gauss_filt(dt, npad, gfilt) + gnorm = np.sum(gauss)*(freqs[1]-freqs[0])*dt + else: + gauss = np.ones(npad) + gnorm = 1. + + # Copy traces + rfL = SLL.copy() + rfQ = SQL.copy() + rfT = STL.copy() + + # Spectral division and inverse transform + rfL.data = np.fft.ifftshift(np.real(np.fft.ifft( + gauss*SLL.data/Sdenom))/gnorm) + rfQ.data = np.fft.ifftshift(np.real(np.fft.ifft( + gauss*SQL.data/Sdenom))/gnorm) + rfT.data = np.fft.ifftshift(np.real(np.fft.ifft( + gauss*STL.data/Sdenom))/gnorm) # Update stats of streams - rfL.stats.channel = 'RF' + self.meta.align[0] - rfQ.stats.channel = 'RF' + self.meta.align[1] - rfT.stats.channel = 'RF' + self.meta.align[2] + rfL.stats.channel = 'RF' + align[0] + rfQ.stats.channel = 'RF' + align[1] + rfT.stats.channel = 'RF' + align[2] self.rf = Stream(traces=[rfL, rfQ, rfT]) @@ -1034,7 +1085,7 @@ def calc_cc(self): self.meta.cc = np.corrcoef(obs_Q.data, pred_Q.data)[0][1] - def to_stream(self): + def to_stream(self, store=None): """ Method to switch from RFData object to Stream object. This allows easier manipulation of the receiver functions @@ -1045,7 +1096,7 @@ def to_stream(self): if not self.meta.accept: return - def _add_rfstats(trace): + def _add_stats(trace, store): trace.stats.snr = self.meta.snr trace.stats.snrh = self.meta.snrh trace.stats.cc = self.meta.cc @@ -1059,18 +1110,38 @@ def _add_rfstats(trace): trace.stats.vp = self.meta.vp trace.stats.vs = self.meta.vs trace.stats.phase = self.meta.phase - trace.stats.is_rf = True - nn = self.rf[0].stats.npts - sr = self.rf[0].stats.sampling_rate + trace.stats.align = self.meta.align + nn = self.data[0].stats.npts + sr = self.data[0].stats.sampling_rate trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn) + + if store == 'rf': + trace.stats.is_rf = True + elif store == 'specs': + trace.stats.is_specs = True + return trace - if not hasattr(self, 'rf'): - raise(Exception("Warning: Receiver functions are not available")) + if store is None: + stream = self.data + for tr in stream: + tr = _add_stats(tr, store) - stream = self.rf - for tr in stream: - tr = _add_rfstats(tr) + if store == 'rf': + if not hasattr(self, 'rf'): + raise(Exception("Warning: Receiver functions are not available")) + else: + stream = self.rf + for tr in stream: + tr = _add_stats(tr, store) + + elif store == 'specs': + if not hasattr(self, 'specs'): + raise(Exception("Warning: Spectra are not available")) + else: + stream = self.specs + for tr in stream: + tr = _add_stats(tr, store) return stream @@ -1088,3 +1159,4 @@ def save(self, file): output = open(file, 'wb') pickle.dump(self, output) output.close() + diff --git a/rfpy/rfstreams.py b/rfpy/rfstreams.py new file mode 100644 index 0000000..22249ba --- /dev/null +++ b/rfpy/rfstreams.py @@ -0,0 +1,560 @@ +def RFStreams(object): + + def _init__(self, rfdata=None): + + self.rfdata = [] + + if isinstance(rfdata, RFData): + rfdata = [rfdata] + if rfdata: + self.rfdata.extend(rfdata) + + def __add__(self, other): + """ + Add two `:class:`~plateflex.classes.Grid` objects or a + :class:`~plateflex.classes.Project` object with a single grid. + + """ + if isinstance(other, RFData): + other = RFStreams([other]) + if not isinstance(other, RFStreams): + raise TypeError + rfdata = self.rfdata + other.rfdata + return self.__class__(rfdata=rfdata) + + def __iter__(self): + """ + Return a robust iterator for :class:`~plateflex.classes.Grid` + objects + + """ + return list(self.rfdata).__iter__() + + def append(self, rfdata): + """ + Append a single :class:`~plateflex.classes.Grid` object to the + current `:class:`~plateflex.classes.Project` object. + + :type grid: :class:`~plateflex.classes.Grid` + :param grid: object to append to project + + .. rubric:: Example + + >>> import numpy as np + >>> from plateflex import Grid, Project + >>> nn = 200; dd = 10. + >>> grid = Grid(np.random.randn(nn, nn), dd, dd) + >>> project = Project() + >>> project.append(grid) + """ + + if isinstance(rfdata, RFData): + self.rfdata.append(rfdata) + else: + msg = 'Append only supports a single RFData object as an argument.' + raise TypeError(msg) + + return self + + def extend(self, rfdata_list): + """ + Extend the current Project object with a list of Grid objects. + + :param trace_list: list of :class:`~plateflex.classes.Grid` objects or + :class:`~plateflex.classes.Project`. + + .. rubric:: Example + + >>> import numpy as np + >>> from plateflex import Grid, Project + >>> nn = 200; dd = 10. + >>> grid1 = Grid(np.random.randn(nn, nn), dd, dd) + >>> grid2 = Grid(np.random.randn(nn, nn), dd, dd) + >>> project = Project() + >>> project.extend(grids=[grid1, grid2]) + + """ + if isinstance(rfdata_list, list): + for _i in rfdata_list: + # Make sure each item in the list is a Grid object. + if not isinstance(_i, RFData): + msg = 'Extend only accepts a list of RFData objects.' + raise TypeError(msg) + self.rfdata.extend(rfdata_list) + elif isinstance(rfdata_list, RFStreams): + self.rfdata.extend(rfdata_list.rfdata) + else: + msg = 'Extend only supports a list of RFData objects as argument.' + raise TypeError(msg) + return self + + + def calc_spectra(self, wavelet='complete',): + + def _npow2(x): + return 1 if x == 0 else 2**(x-1).bit_length() + + def _pad(array, n): + tmp = np.zeros(n) + tmp[:array.shape[0]] = array + return tmp + + def _Pwavelet(parent, method='complete', overhang=5, + envelope_threshold=0.05, time=5): + + """ + Select wavelet from the parent function for deconvolution using method. + parent: obspy.Trace + wavefrom to extract the wavelet from + method: str + 'complete' use complete parent signal after P arrival (current + implementation) + 'envelope' use only the part of the parent signal after the + P arrival where + envelope > envelope_threshold*max(envelope) + fall back to 'complete' if condition not reached + 'time' use only this many seonds after P arrival` + fall back to 'complete' if longer than parent + overhang: float + seconds before start and after end of wavelet to be used for + tapering + envelope_threshold: float + fraction of the envelope that defines wavelet (for + method='envelope') + time: float + window (seconds) that defines the wavelet (for method='time') + minimum time (seconds) of the wavelet (for method='envelope') + + Return: + left, right: (obspy.UTCDateTime) start and end time of wavelet + """ + + import obspy.signal.filter as osf + + if method not in ['complete', 'envelope', 'time', 'noise']: + msg = 'Unknow method for wavelet extraction: ' + method + raise NotImplementedError(msg) + + errflag = False + + if method == 'envelope': + split = int((self.meta.time + self.meta.ttime + + time - parent.stats.starttime ) * + parent.stats.sampling_rate) + env = osf.envelope(parent.data) + env /= max(env) # normalize + env = env[split:] # only look after P + time + try: + i = np.nonzero( + np.diff( + np.array( + env > envelope_threshold, dtype=int))==-1)[0][0] + except IndexError: + i = len(parent.data)-1 + dts = i * parent.stats.delta + time + left = self.meta.time + self.meta.ttime - overhang + right = self.meta.time + self.meta.ttime + dts + overhang + if left < parent.stats.starttime or right > parent.stats.endtime: + print('Envelope wavelet longer than trace.') + print('Falling back to "complete" wavelet') + errflag = True + + if method == 'time': + left = self.meta.time + self.meta.ttime - overhang + right = self.meta.time + self.meta.ttime + time + overhang + if left < parent.stats.starttime or right > parent.stats.endtime: + print('Time wavelet longer than trace.') + print('Falling back to "complete" wavelet') + errflag = True + + if method == 'complete' or errflag: + dts = len(parent.data)*parent.stats.delta/2. + left = self.meta.time + self.meta.ttime - overhang + right = self.meta.time + self.meta.ttime + dts - 2*overhang + + if method == 'noise': + dts = len(parent.data)*parent.stats.delta/2. + left = self.meta.time + self.meta.ttime - dts + right = self.meta.time + self.meta.ttime - overhang + + return left, right + + def _specs(trL, trQ, trT, trN, nn, method): + + # Get length, zero padding parameters and frequencies + dt = parent.stats.delta + + # npad = _npow2(nn*2) + npad = nn + freqs = np.fft.fftfreq(npad, d=dt) + + # Fourier transform + Fl = np.fft.fft(trL.data, n=npad) + Fq = np.fft.fft(trQ.data, n=npad) + Ft = np.fft.fft(trT.data, n=npad) + Fn = np.fft.fft(trN.data, n=npad) + + # Copy traces + Sll = trL.copy() + Sql = trQ.copy() + Stl = trT.copy() + Snn = trN.copy() + + # Auto and cross spectra + Sll.data = np.real(Fl*np.conjugate(Fl)) + Sql.data = Fq*np.conjugate(Fl) + Stl.data = Ft*np.conjugate(Fl) + Snn.data = np.real(Fn*np.conjugate(Fn)) + + return Sll, Sql, Stl, Snn + + + for rfdata in self.rfdata: + + # Get the name of components (order is critical here) + cL = stream.meta.align[0] + cQ = stream.meta.align[1] + cT = stream.meta.align[2] + + # Define signal and noise + trL = stream.data.select(component=cL)[0].copy() + trQ = stream.data.select(component=cQ)[0].copy() + trT = stream.data.select(component=cT)[0].copy() + trN = stream.data.select(component=cL)[0].copy() + + # Standardize all trace data + meanN = np.mean(trN.data) + stdN = np.std(trN.data) + [tr.data = (tr.data - meanN)/stdN for tr in [trL, trQ, trT, trN]] + + trL.stats.channel = 'WV' + stream.meta.align[0] + trQ.stats.channel = 'WV' + stream.meta.align[1] + trT.stats.channel = 'WV' + stream.meta.align[2] + trN.stats.channel = 'WV' + stream.meta.align[0] + + if phase == 'P' or 'PP': + + # Get signal length (i.e., seismogram to deconvolve) from trace length + over = 5 + dtsqt = len(trL.data)*trL.stats.delta/2. + + # Traces will be zero-paded to this length (samples) + nn = int(round((dtsqt+over)*trL.stats.sampling_rate)) + 1 + + sig_left, sig_right = _Pwavelet(trL, method=wavelet, + envelope_threshold=envelope_threshold, time=time, overhang=over) + + # Trim wavelet + trL.trim(sig_left, sig_right, nearest_sample=False, pad=True, + fill_value=0.) + + # Signal window (-5. to dtsqt-10 sec) + sig_left, sig_right = _Pwavelet(trQ, method='complete', overhang=over) + + # Trim signal traces + [tr.trim(sig_left, sig_right, nearest_sample=False, + pad=True, fill_value=0.) for tr in [trQ, trT]] + + # Noise window (-dtsqt to -5. sec) + noise_left, noise_right = _Pwavelet(trQ, method='noise', overhang=over) + + # Trim noise trace + trN.trim(noise_left, noise_right, nearest_sample=False, + pad=True, fill_value=0.) + + elif phase == 'S' or 'SKS': + + # Get signal length (i.e., seismogram to deconvolve) from trace length + dts = len(trL.data)*trL.stats.delta/2. + + # Trim signal traces (-5. to dts-10 sec) + trL.trim(self.meta.time+self.meta.ttime+25.-dts/2., + self.meta.time+self.meta.ttime+25.) + trQ.trim(self.meta.time+self.meta.ttime+25.-dts/2., + self.meta.time+self.meta.ttime+25.) + trT.trim(self.meta.time+self.meta.ttime+25.-dts/2., + self.meta.time+self.meta.ttime+25.) + + # Trim noise traces (-dts to -5 sec) + trN.trim(self.meta.time+self.meta.ttime-dts, + self.meta.time+self.meta.ttime-dts/2.) + + # Taper traces - only necessary processing after trimming + # TODO: What does this to the multitaper method + [tr.taper(max_percentage=0.05, max_length=2.) + for tr in [trL, trQ, trT, trN]] + + # Spectra + if phase == 'P' or 'PP': + sll, sql, stl, snn = _specs(trL, trQ, trT, trN, nn) + + rfdata.specs = Stream(traces=[sll, sql, stl, snn]) + + + # elif phase == 'S' or 'SKS': + # rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method) + + + def bin_baz_slow_specs(self, nbaz=72, nslow=1): + + """ + Function to stack auto- and cross-component spectra into back-azimuth + and slowness bins. + + Parameters + ---------- + nbaz : int + Number of back-azimuth samples in bins + nslow : int + Number of slowness samples in bins + include_empty : bool + Return empty bins as null traces (default omits them) + + Returns + ------- + stack : :class:`~obspy.core.Stream` + Stream containing one or two stacked traces, + depending on the number of input streams + + """ + + # Define back-azimuth and slowness bins + baz_bins = np.linspace(0, 360, nbaz) + slow_bins = np.linspace(0.04, 0.08, nslow) + + # Extract baz and slowness + baz = [stream[i].stats.baz for i in range(len(stream))] + slow = [stream[i].stats.slow for i in range(len(stream))] + + # Digitize baz and slowness + ibaz = np.digitize(baz, baz_bins) + islow = np.digitize(slow, slow_bins) + + final_stream = [] + + for stream in self.streams: + try: + # Define empty streams + binned_stream = Stream() + + # Loop through baz_bins + for i in range(nbaz): + for j in range(nslow): + + nbin = 0 + array = np.zeros(len(stream[0].specs.data)) + + # Loop all traces + for k, tr in enumerate(stream): + + # If index of baz_bins is equal to ibaz + if i == ibaz[k] and j == islow[k]: + + nbin += 1 + array += tr.data + + continue + + if nbin > 0 or include_empty: + + # Average and update stats + array /= nbin + weight = np.real(abs(weight/nbin)) + + trace = Trace(header=stream[0].stats) + trace.stats.baz = baz_bins[i] + trace.stats.slow = slow_bins[j] + trace.stats.nbin = nbin + + if not pws: + weight = np.ones(len(stream[0].data)) + trace.data = weight*array + binned_stream.append(trace) + + final_stream.append(binned_stream) + + except: + continue + + return final_stream + + + def deconvolve(self, phase='P', vp=None, vs=None, + align=None, method='wiener', + envelope_threshold=0.05, time=5, pre_filt=None, + gfilt=None, wlevel=0.01, writeto=None): + """ + Deconvolves three-component data using one component as the source wavelet. + The source component is always taken as the dominant compressional + component, which can be either 'Z', 'L', or 'P'. + + Parameters + ---------- + vp : float + P-wave velocity at surface (km/s) + vs : float + S-wave velocity at surface (km/s) + align : str + Alignment of coordinate system for rotation + ('ZRT', 'LQT', or 'PVH') + method : str + Method for deconvolution. Options are 'wiener', 'water' or + 'multitaper' + wavelet : str + Type of wavelet for deconvolution. Options are 'complete', 'time' or + 'envelope' + envelope_threshold : float + Threshold [0-1] used in ``wavelet='envelope'``. + time : float + Window length used in ``wavelet='time'``. + Minimum window length for ``wavelet='envelope'``. + pre_filt : list of 2 floats + Low and High frequency corners of bandpass filter applied + before deconvolution + gfilt : float + Center frequency of Gaussian filter (Hz). + wlevel : float + Water level used in ``method='water'``. + writeto : str or None + Write wavelets for deconvolution to file. + + Attributes + ---------- + rf : :class:`~obspy.core.Stream` + Stream containing the receiver function traces + + """ + + if not self.meta.accept: + return + + + def _gauss_filt(dt, nft, f0): + df = 1./(nft*dt) + nft21 = int(0.5*nft + 1) + f = df*np.arange(nft21) + w = 2.*np.pi*f + gauss = np.zeros(nft) + gauss[:nft21] = np.exp(-0.25*(w/f0)**2.)/dt + gauss[nft21:] = np.flip(gauss[1:nft21-1]) + return gauss + + + def _decon(parent, daughter1, daughter2, noise, nn, method): + + + # Final processing depends on method + if method == 'wiener': + Sdenom = Spp + Snn + elif method == 'water': + phi = np.amax(Spp)*wlevel + Sdenom = Spp + Sdenom[Sdenom < phi] = phi + + # Apply Gaussian filter? + if gfilt: + gauss = _gauss_filt(dt, npad, gfilt) + gnorm = np.sum(gauss)*(freqs[1]-freqs[0])*dt + else: + gauss = np.ones(npad) + gnorm = 1. + + # Is this correct? + #parent de-noised = np.fft.ifftshift(np.real(np.fft.ifft(gauss*Sdenom))/gnorm) + #daughter1 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd1p))/gnorm) + #daughter2 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd2p))/gnorm) + + # Copy traces + rfp = parent.copy() + rfd1 = daughter1.copy() + rfd2 = daughter2.copy() + + # Spectral division and inverse transform + rfp.data = np.fft.ifftshift(np.real(np.fft.ifft( + gauss*Spp/Sdenom))/gnorm) + rfd1.data = np.fft.ifftshift(np.real(np.fft.ifft( + gauss*Sd1p/Sdenom))/gnorm) + rfd2.data = np.fft.ifftshift(np.real(np.fft.ifft( + gauss*Sd2p/Sdenom))/gnorm) + + return rfp, rfd1, rfd2 + + if not self.meta.rotated: + print("Warning: Data have not been rotated yet - rotating now") + self.rotate(vp=vp, vs=vs, align=align) + + # v--True if None v--True if nan, error if None + if not self.meta.snr or not np.isfinite(self.meta.snr): + print("Warning: SNR has not been calculated - " + + "calculating now using default") + self.calc_snr() + + if hasattr(self, 'rf'): + print("Warning: Data have been deconvolved already - passing") + return + + + # Pre-filter waveforms before deconvolution + if pre_filt: + [tr.filter('bandpass', freqmin=pre_filt[0], freqmax=pre_filt[1], + corners=2, zerophase=True) + for tr in [trL, trQ, trT, trNl, trNq]] + + if writeto: + with open(writeto, 'wb') as f: + pickle.dump(Stream(traces=[trL, trQ, trT]), f) + + # Deconvolve + if phase == 'P' or 'PP': + rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method) + + elif phase == 'S' or 'SKS': + rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method) + + # Update stats of streams + rfL.stats.channel = 'RF' + self.meta.align[0] + rfQ.stats.channel = 'RF' + self.meta.align[1] + rfT.stats.channel = 'RF' + self.meta.align[2] + + self.rf = Stream(traces=[rfL, rfQ, rfT]) + + + def to_stream(self): + """ + Method to switch from RFData object to Stream object. + This allows easier manipulation of the receiver functions + for post-processing. + + """ + + if not self.meta.accept: + return + + def _add_specstats(trace): + trace.stats.snr = self.meta.snr + trace.stats.snrh = self.meta.snrh + trace.stats.cc = self.meta.cc + trace.stats.slow = self.meta.slow + trace.stats.baz = self.meta.baz + trace.stats.gac = self.meta.gac + trace.stats.stlo = self.sta.longitude + trace.stats.stla = self.sta.latitude + trace.stats.evlo = self.meta.lon + trace.stats.evla = self.meta.lat + trace.stats.vp = self.meta.vp + trace.stats.vs = self.meta.vs + trace.stats.phase = self.meta.phase + trace.stats.is_rf = True + nn = self.rf[0].stats.npts + sr = self.rf[0].stats.sampling_rate + trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn) + return trace + + if not hasattr(self, 'rf'): + raise(Exception("Warning: Receiver functions are not available")) + + stream = self.rf + for tr in stream: + tr = _add_rfstats(tr) + + return stream diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 32507f8..251fde5 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -740,12 +740,17 @@ def main(): print("*"*50) continue - # Deconvolve data - rfdata.deconvolve( + # Calculate spectra + rfdata.calc_spectra( vp=args.vp, vs=args.vs, align=args.align, method=args.method, - gfilt=args.gfilt, wlevel=args.wlevel, - pre_filt=args.pre_filt) + wavelet='complete', envelope_threshold=0.05, + time=5) + + # Deconvolve data + rfdata.deconvolve( + method=args.method, gfilt=args.gfilt, + wlevel=args.wlevel) # Get cross-correlation QC rfdata.calc_cc() @@ -753,7 +758,7 @@ def main(): print("* CC: {}".format(rfdata.meta.cc)) # Convert to Stream - rfstream = rfdata.to_stream() + rfstream = rfdata.to_stream(store='rf') # Save event meta data pickle.dump(rfdata.meta, open(metafile, "wb")) diff --git a/rfpy/scripts/rfpy_calc_simdec.py b/rfpy/scripts/rfpy_calc_simdec.py new file mode 100644 index 0000000..b9bf65f --- /dev/null +++ b/rfpy/scripts/rfpy_calc_simdec.py @@ -0,0 +1,481 @@ +#!/usr/bin/env python + +# Copyright 2019 Pascal Audet +# +# This file is part of RfPy. +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +# Import modules and functions +from argparse import ArgumentParser +from os.path import exists as exist +import numpy as np +import pickle +import stdb +from obspy import Stream +from rfpy import RFData, binning, plotting +from pathlib import Path + + +def get_simdec_arguments(argv=None): + """ + Get Options from :class:`~optparse.OptionParser` objects. + + This function is used for data processing on-the-fly (requires web connection) + + """ + + parser = ArgumentParser( + usage="%(prog)s [arguments] ", + description="Script used to calculate receiver functions " + + "using simultaneous deconvolution. The stations are processed one " + + "by one and the data are stored to disk. ") + + # General Settings + parser.add_argument( + "indb", + help="Station Database to process from.", + type=str) + parser.add_argument( + "--keys", + action="store", + type=str, + dest="stkeys", + default="", + help="Specify a comma separated list of station keys for " + + "which to perform the analysis. These must be " + + "contained within the station database. Partial keys will " + + "be used to match against those in the dictionary. For " + + "instance, providing IU will match with all stations in " + + "the IU network [Default processes all stations in the database]") + parser.add_argument( + "-v", "-V", "--verbose", + action="store_true", + dest="verb", + default=False, + help="Specify to increase verbosity.") + parser.add_argument( + "-L", "--long-name", + action="store_true", + dest="lkey", + default=False, + help="Force folder names to use long-key form (NET.STN.CHN). " + + "Default behaviour uses short key form (NET.STN) for the folder " + + "names, regardless of the key type of the database." + ) + + # Constants Settings + ConstGroup = parser.add_argument_group( + title='Parameter Settings', + description="Miscellaneous default values and settings") + ConstGroup.add_argument( + "--Z12", + action="store_true", + dest="Z12", + default=False, + help="Use Z12 data if available. [Default uses ZNE data]") + ConstGroup.add_argument( + "--phase", + action="store", + type=str, + dest="phase", + default='allP', + help="Specify the phase name to use. Be careful with the distance. " + + "setting. Options are 'P', 'PP', 'allP', 'S', 'SKS' or 'allS'. " + + "[Default 'allP']") + ConstGroup.add_argument( + "--resample", + action="store", + type=float, + dest="resample", + default=None, + help="Specify the new sampling-rate for the receiver functions. " + + "Note the sampling rate of the original data (ZNE or Z12) stored " + + "on disk is unchanged. [Default None]") + ConstGroup.add_argument( + "--align", + action="store", + type=str, + dest="align", + default=None, + help="Specify component alignment key. Can be either " + + "ZRT, LQT, or PVH. [Default ZRT]") + ConstGroup.add_argument( + "--vp", + action="store", + type=float, + dest="vp", + default=6.0, + help="Specify near-surface Vp to use with --align=PVH (km/s). " + + "[Default 6.0]") + ConstGroup.add_argument( + "--vs", + action="store", + type=float, + dest="vs", + default=3.5, + help="Specify near-surface Vs to use with --align=PVH (km/s). " + + "[Default 3.5]") + ConstGroup.add_argument( + "--dt-snr", + action="store", + type=float, + dest="dt_snr", + default=30., + help="Specify the window length over which to calculate " + + "the SNR in sec. [Default 30.]") + ConstGroup.add_argument( + "--pre-filt", + action="store", + type=str, + dest="pre_filt", + default=None, + help="Specify two floats with low and high frequency corners for " + + "pre-filter (before deconvolution). [Default None]") + ConstGroup.add_argument( + "--fmin", + action="store", + type=float, + dest="fmin", + default=0.05, + help="Specify the minimum frequency corner for SNR " + + "filter (Hz). [Default 0.05]") + ConstGroup.add_argument( + "--fmax", + action="store", + type=float, + dest="fmax", + default=1.0, + help="Specify the maximum frequency corner for SNR " + + "filter (Hz). [Default 1.0]") + + # Deconvolution Settings + DeconGroup = parser.add_argument_group( + title='Deconvolution Settings', + description="Parameters for deconvolution") + DeconGroup.add_argument( + "--method", + action="store", + dest="method", + type=str, + default="wiener", + help="Specify the deconvolution method. Available methods " + + "include 'wiener', 'water' and 'multitaper'. [Default 'wiener']") + DeconGroup.add_argument( + "--gfilt", + action="store", + dest="gfilt", + type=float, + default=None, + help="Specify the Gaussian filter width in Hz. " + + "[Default None]") + DeconGroup.add_argument( + "--wlevel", + action="store", + dest="wlevel", + type=float, + default=0.01, + help="Specify the water level, used in the 'water' method. " + + "[Default 0.01]") + + PreGroup = parser.add_argument_group( + title='Pre-processing Settings', + description="Options for pre-processing " + + "data before stacking") + PreGroup.add_argument( + "--snr", + action="store", + type=float, + dest="snr", + default=-9999., + help="Specify the vertical component SNR threshold for " + + "extracting receiver functions. [Default 5.]") + PreGroup.add_argument( + "--snrh", + action="store", + type=float, + dest="snrh", + default=-9999., + help="Specify the horizontal component SNR threshold for " + + "extracting receiver functions. [Default None]") + PreGroup.add_argument( + "--nbaz", + action="store", + dest="nbaz", + type=int, + default=36, + help="Specify integer number of back-azimuth bins to consider " + + "(typically 36 or 72). If not None, the plot will show receiver " + + "functions sorted by back-azimuth values. [Default 36]") + PreGroup.add_argument( + "--nslow", + action="store", + dest="nslow", + type=int, + default=16, + help="Specify integer number of slowness bins to consider " + + "(range of slowness values is 0.04 to 0.08). " + + "If not None, the plot will show receiver " + + "functions sorted by slowness values. [Default 16]") + + + args = parser.parse_args(argv) + + # Check inputs + if not exist(args.indb): + parser.error("Input file " + args.indb + " does not exist") + + # create station key list + if len(args.stkeys) > 0: + args.stkeys = args.stkeys.split(',') + + if args.phase not in ['P', 'PP', 'allP', 'S', 'SKS', 'allS']: + parser.error( + "Error: choose between 'P', 'PP', 'allP', 'S', 'SKS' and 'allS'.") + if args.phase == 'allP': + args.listphase = ['P', 'PP'] + elif args.phase == 'allS': + args.listphase = ['S', 'SKS'] + else: + args.listphase = [args.phase] + + if args.align is None: + args.align = 'ZRT' + elif args.align not in ['ZRT', 'LQT', 'PVH']: + parser.error( + "Error: Incorrect alignment specifier. Should be " + + "either 'ZRT', 'LQT', or 'PVH'.") + + if args.method not in ['wiener', 'water', 'multitaper']: + parser.error( + "Error: 'method' should be either 'wiener', 'water' or " + + "'multitaper'") + + if args.pre_filt is not None: + args.pre_filt = [float(val) for val in args.pre_filt.split(',')] + args.pre_filt = sorted(args.pre_filt) + if (len(args.pre_filt)) != 2: + parser.error( + "Error: --pre-filt should contain 2 " + + "comma-separated floats") + + return args + + +def main(): + + print() + print("########################################################") + print("# #") + print("# __ _ #") + print("# _ __ / _|_ __ _ _ _ __ ___ ___ __ _| | ___ #") + print("# | '__| |_| '_ \| | | | | '__/ _ \/ __/ _` | |/ __| #") + print("# | | | _| |_) | |_| | | | | __/ (_| (_| | | (__ #") + print("# |_| |_| | .__/ \__, |___|_| \___|\___\__,_|_|\___| #") + print("# |_| |___/_____| #") + print("# #") + print("########################################################") + print() + + # Run Input Parser + args = get_simdec_arguments() + + # Load Database + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + # Track processed folders + procfold = [] + + # Loop over station keys + for stkey in list(stkeys): + + # Extract station information from dictionary + sta = db[stkey] + + # Construct Folder Name + stfld = stkey + if not args.lkey: + stfld = stkey.split('.')[0]+"."+stkey.split('.')[1] + + # Define path to see if it exists + if args.phase in ['P', 'PP', 'allP']: + datapath = Path('P_DATA') / stfld + elif args.phase in ['S', 'SKS', 'allS']: + datapath = Path('S_DATA') / stfld + if not datapath.is_dir(): + print('Path to ' + str(datapath) + ' doesn`t exist - continuing') + continue + + # Temporary print locations + tlocs = sta.location + if len(tlocs) == 0: + tlocs = [''] + for il in range(0, len(tlocs)): + if len(tlocs[il]) == 0: + tlocs[il] = "--" + sta.location = tlocs + + # Update Display + print(" ") + print(" ") + print("|===============================================|") + print("| {0:>8s} |".format( + sta.station)) + print("|===============================================|") + print("| Station: {0:>2s}.{1:5s} |".format( + sta.network, sta.station)) + print("| Channel: {0:2s}; Locations: {1:15s} |".format( + sta.channel, ",".join(tlocs))) + print("| Lon: {0:7.2f}; Lat: {1:6.2f} |".format( + sta.longitude, sta.latitude)) + print("|-----------------------------------------------|") + + # # Check for folder already processed + # if stfld in procfold: + # print(' {0} already processed...skipping '.format(stfld)) + # continue + + LLstreams = Stream() + QLstreams = Stream() + TLstreams = Stream() + NNstreams = Stream() + + datafiles = [x for x in datapath.iterdir() if x.is_dir()] + for folder in datafiles: + + # Skip hidden folders + if folder.name.startswith('.'): + continue + + # Re-initialize RFData object + rfdata = RFData(sta) + + # Load meta data + metafile = folder / "Meta_Data.pkl" + if not metafile.is_file(): + continue + rfdata.meta = pickle.load(open(metafile, 'rb')) + + # Skip data not in list of phases + if rfdata.meta.phase not in args.listphase: + continue + + if args.verb: + print("* Station: {0}; folder: {1}".format(stkey, folder)) + + if args.Z12: + try: + # Try loading Z12_Data + Z12file = folder / "Z12_Data.pkl" + rfdata.data = pickle.load(open(Z12file, 'rb')) + # Remove rotated flag and snr flag + rfdata.meta.rotated = False + rfdata.rotate(align='ZNE') + except: + + print("Z12_Data.pkl not available - using ZNE_Data.pkl") + # Load ZNE data + ZNEfile = folder / "ZNE_Data.pkl" + rfdata.data = pickle.load(open(ZNEfile, 'rb')) + else: + # Load ZNE data + ZNEfile = folder / "ZNE_Data.pkl" + rfdata.data = pickle.load(open(ZNEfile, 'rb')) + + # Resample if requested + if args.resample: + rfdata.data.resample( + args.resample, no_filter=False) + + # Remove rotated flag and snr flag + rfdata.meta.rotated = False + rfdata.meta.snr = None + + # Rotate from ZNE to 'align' ('ZRT', 'LQT', or 'PVH') + rfdata.rotate(vp=args.vp, vs=args.vs, align=args.align) + + # Calculate SNR + rfdata.calc_snr(dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax) + + # QC Thresholding + if rfdata.meta.snrh < args.snrh: + continue + if rfdata.meta.snr < args.snr: + continue + + # Calculate spectra + rfdata.calc_spectra( + vp=args.vp, vs=args.vs, + align=args.align, method=args.method, + wavelet='complete', envelope_threshold=0.05, + time=5) + + # Convert to Stream + rfstream = rfdata.to_stream(store='specs') + + # Add to Stream objects + LLstreams.append(rfstream[0]) + QLstreams.append(rfstream[1]) + TLstreams.append(rfstream[2]) + NNstreams.append(rfstream[3]) + + # Bin into baz and slowness bins + LL_tmp, = binning.bin_baz_slow(LLstreams, nbaz=args.nbaz, nslow=args.nslow, + pws=False, include_empty=False) + QL_tmp, = binning.bin_baz_slow(QLstreams, nbaz=args.nbaz, nslow=args.nslow, + pws=False, include_empty=False) + TL_tmp, = binning.bin_baz_slow(TLstreams, nbaz=args.nbaz, nslow=args.nslow, + pws=False, include_empty=False) + NN_tmp, = binning.bin_baz_slow(NNstreams, nbaz=args.nbaz, nslow=args.nslow, + pws=False, include_empty=False) + print(QL_tmp) + + # Initialize new empty streams + RFL = Stream() + RFQ = Stream() + RFT = Stream() + + # Cycle through binned spectra + for i in range(len(LL_tmp)): + + # Initialize new RFData object with empty Meta header + rfdata = RFData(sta) + + # Extract the specs from the binned streams + rfdata.specs = Stream(traces=[LL_tmp[i], QL_tmp[i], TL_tmp[i], NN_tmp[i]]) + + # Deconvolve + rfdata.deconvolve( + align=args.align, method=args.method, + gfilt=args.gfilt, wlevel=args.wlevel) + + # Append RFs to streams + RFL.append(rfdata.rf[0]) + RFQ.append(rfdata.rf[1]) + RFT.append(rfdata.rf[2]) + + # # Update processed folders + # procfold.append(stfld) + + +if __name__ == "__main__": + + # Run main program + main() diff --git a/rfpy/scripts/rfpy_recalc.py b/rfpy/scripts/rfpy_calc_sse.py similarity index 100% rename from rfpy/scripts/rfpy_recalc.py rename to rfpy/scripts/rfpy_calc_sse.py diff --git a/rfpy/scripts/rfpy_get_data.py b/rfpy/scripts/rfpy_get_data.py new file mode 100644 index 0000000..f1005f5 --- /dev/null +++ b/rfpy/scripts/rfpy_get_data.py @@ -0,0 +1,631 @@ +#!/usr/bin/env python + +# Copyright 2019 Pascal Audet +# +# This file is part of RfPy. +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +# Import modules and functions +import numpy as np +import pickle +import stdb +from obspy.clients.fdsn import Client +from rfpy import utils, RFData +from pathlib import Path +from argparse import ArgumentParser +from os.path import exists as exist +from obspy import UTCDateTime +from numpy import nan + + +def get_data_arguments(argv=None): + """ + Get Options from :class:`~optparse.OptionParser` objects. + + This function is used for data processing on-the-fly (requires web connection) + + """ + + parser = ArgumentParser( + usage="%(prog)s [arguments] ", + description="Script used to download and pre-process " + + "three-component ('Z', 'N', and 'E'), seismograms for individual " + + "events and calculate teleseismic P-wave receiver functions" + + "This version requests data on the fly for a given date " + + "range. Data are requested from the internet using the " + + "client services framework. The stations are processed one " + + "by one and the data are stored to disk.") + + # General Settings + parser.add_argument( + "indb", + help="Station Database to process from.", + type=str) + parser.add_argument( + "--keys", + action="store", + type=str, + dest="stkeys", + default="", + help="Specify a comma separated list of station keys for " + + "which to perform the analysis. These must be " + + "contained within the station database. Partial keys will " + + "be used to match against those in the dictionary. For " + + "instance, providing IU will match with all stations in " + + "the IU network [Default processes all stations in the database]") + parser.add_argument( + "-v", "-V", "--verbose", + action="store_true", + dest="verb", + default=False, + help="Specify to increase verbosity.") + parser.add_argument( + "-O", "--overwrite", + action="store_true", + dest="ovr", + default=False, + help="Force the overwriting of pre-existing data. " + + "[Default False]") + parser.add_argument( + "-L", "--long-name", + action="store_true", + dest="lkey", + default=False, + help="Force folder names to use long-key form (NET.STN.CHN). " + + "Default behaviour uses short key form (NET.STN) for the folder " + + "names, regardless of the key type of the database." + ) + + # Server Settings + ServerGroup = parser.add_argument_group( + title="Server Settings", + description="Settings associated with which " + "datacenter to log into.") + ServerGroup.add_argument( + "-S", "--Server", + action="store", + type=str, + dest="Server", + default="IRIS", + help="Specify the server to connect to. Options include: " + + "BGR, ETH, GEONET, GFZ, INGV, IPGP, IRIS, KOERI, " + + "LMU, NCEDC, NEIP, NERIES, ODC, ORFEUS, RESIF, SCEDC, USGS, USP. " + + "[Default IRIS]") + ServerGroup.add_argument( + "-U", "--User-Auth", + action="store", + type=str, + dest="UserAuth", + default="", + help="Enter your IRIS Authentification Username and Password " + + "(--User-Auth='username:authpassword') to " + + "access and download restricted data. " + + "[Default no user and password]") + + # Database Settings + DataGroup = parser.add_argument_group( + title="Local Data Settings", + description="Settings associated with defining " + + "and using a local data base of pre-downloaded " + + "day-long SAC or MSEED files.") + DataGroup.add_argument( + "--local-data", + action="store", + type=str, + dest="localdata", + default=None, + help="Specify a comma separated list of paths containing " + + "day-long sac or mseed files of data already downloaded. " + + "If data exists for a seismogram is already present on disk, " + + "it is selected preferentially over downloading " + + "the data using the Client interface") + DataGroup.add_argument( + "--dtype", + action="store", + type=str, + dest="dtype", + default='SAC', + help="Specify the data archive file type, either SAC " + + " or MSEED. Note the default behaviour is to search for " + + "SAC files. Local archive files must have extensions of '.SAC' "+ + " or '.MSEED. These are case dependent, so specify the correct case"+ + "here.") + DataGroup.add_argument( + "--no-data-zero", + action="store_true", + dest="ndval", + default=False, + help="Specify to force missing data to be set as zero, rather " + + "than default behaviour which sets to nan. Note this is applied " + + "only to the SAC data") + DataGroup.add_argument( + "--no-local-net", + action="store_false", + dest="useNet", + default=True, + help="Specify to prevent using the Network code in the " + + "search for local data (sometimes for CN stations " + + "the dictionary name for a station may disagree with that " + + "in the filename. [Default Network used]") + DataGroup.add_argument( + "--save-Z12", + action="store_true", + dest="saveZ12", + default=False, + help="Specify to save Z12 (un-rotated) components. [Default " + + "False]") + + # Event Selection Criteria + EventGroup = parser.add_argument_group( + title="Event Settings", + description="Settings associated with refining " + + "the events to include in matching event-station pairs") + EventGroup.add_argument( + "--start", + action="store", + type=str, + dest="startT", + default="", + help="Specify a UTCDateTime compatible string representing " + + "the start time for the event search. This will override any " + + "station start times. [Default start date of station]") + EventGroup.add_argument( + "--end", + action="store", + type=str, + dest="endT", + default="", + help="Specify a UTCDateTime compatible string representing " + + "the end time for the event search. This will override any " + + "station end times [Default end date of station]") + EventGroup.add_argument( + "--reverse", "-R", + action="store_true", + dest="reverse", + default=False, + help="Reverse order of events. Default behaviour starts at " + + "oldest event and works towards most recent. Specify reverse " + + "order and instead the program will start with the most recent " + + "events and work towards older") + EventGroup.add_argument( + "--minmag", + action="store", + type=float, + dest="minmag", + default=6.0, + help="Specify the minimum magnitude of event for which to search. " + + "[Default 6.0]") + EventGroup.add_argument( + "--maxmag", + action="store", + type=float, + dest="maxmag", + default=9.0, + help="Specify the maximum magnitude of event for which to search. " + + "[Default None, i.e. no limit]") + + # Geometry Settings + PhaseGroup = parser.add_argument_group( + title="Geometry Settings", + description="Settings associatd with the " + "event-station geometries for the specified phase") + PhaseGroup.add_argument( + "--phase", + action="store", + type=str, + dest="phase", + default='P', + help="Specify the phase name to use. Be careful with the distance. " + + "setting. Options are 'P' or 'PP'. [Default 'P']") + PhaseGroup.add_argument( + "--mindist", + action="store", + type=float, + dest="mindist", + default=None, + help="Specify the minimum great circle distance (degrees) between " + + "the station and event. [Default depends on phase]") + PhaseGroup.add_argument( + "--maxdist", + action="store", + type=float, + dest="maxdist", + default=None, + help="Specify the maximum great circle distance (degrees) between " + + "the station and event. [Default depends on phase]") + + # Constants Settings + ConstGroup = parser.add_argument_group( + title='Parameter Settings', + description="Miscellaneous default values and settings") + ConstGroup.add_argument( + "--sampling-rate", + action="store", + type=float, + dest="new_sampling_rate", + default=10., + help="Specify new sampling rate in Hz. [Default 10.]") + ConstGroup.add_argument( + "--dts", + action="store", + type=float, + dest="dts", + default=150., + help="Specify the window length in sec (symmetric about arrival " + + "time). [Default 150.]") + + args = parser.parse_args(argv) + + # Check inputs + if not exist(args.indb): + parser.error("Input file " + args.indb + " does not exist") + + # create station key list + if len(args.stkeys) > 0: + args.stkeys = args.stkeys.split(',') + + # construct start time + if len(args.startT) > 0: + try: + args.startT = UTCDateTime(args.startT) + except: + parser.error( + "Cannot construct UTCDateTime from start time: " + + args.startT) + else: + args.startT = None + + # construct end time + if len(args.endT) > 0: + try: + args.endT = UTCDateTime(args.endT) + except: + parser.error( + "Cannot construct UTCDateTime from end time: " + + args.endT) + else: + args.endT = None + + # Parse User Authentification + if not len(args.UserAuth) == 0: + tt = args.UserAuth.split(':') + if not len(tt) == 2: + parser.error( + "Error: Incorrect Username and Password " + + "Strings for User Authentification") + else: + args.UserAuth = tt + else: + args.UserAuth = [] + + # Parse Local Data directories + if args.localdata is not None: + args.localdata = args.localdata.split(',') + else: + args.localdata = [] + + # Check Datatype specification + if (not args.dtype.upper() == 'MSEED') and (not args.dtype.upper() == 'SAC'): + parser.error( + "Error: Local Data Archive must be of types 'SAC'" + + "or MSEED. These must match the file extensions for " + + " the archived data.") + + # Check NoData Value + if args.ndval: + args.ndval = 0.0 + else: + args.ndval = nan + + # Check distances for selected phase + if args.phase not in ['P', 'PP', 'S', 'SKS']: + parser.error( + "Error: choose between 'P', 'PP', 'S' and 'SKS'.") + if args.phase == 'P': + if not args.mindist: + args.mindist = 30. + if not args.maxdist: + args.maxdist = 100. + if args.mindist < 30. or args.maxdist > 100.: + parser.error( + "Distances should be between 30 and 100 deg. for " + + "teleseismic 'P' waves.") + elif args.phase == 'PP': + if not args.mindist: + args.mindist = 100. + if not args.maxdist: + args.maxdist = 180. + if args.mindist < 100. or args.maxdist > 180.: + parser.error( + "Distances should be between 100 and 180 deg. for " + + "teleseismic 'PP' waves.") + elif args.phase == 'S': + if not args.mindist: + args.mindist = 55. + if not args.maxdist: + args.maxdist = 85. + if args.mindist < 55. or args.maxdist > 85.: + parser.error( + "Distances should be between 55 and 85 deg. for " + + "teleseismic 'S' waves.") + elif args.phase == 'SKS': + if not args.mindist: + args.mindist = 85. + if not args.maxdist: + args.maxdist = 115. + if args.mindist < 85. or args.maxdist > 115.: + parser.error( + "Distances should be between 85 and 115 deg. for " + + "teleseismic 'SKS' waves.") + + return args + + +def main(): + + print() + print("##############################################") + print("# __ _ #") + print("# _ __ / _|_ __ _ _ ___ __ _| | ___ #") + print("# | '__| |_| '_ \| | | | / __/ _` | |/ __| #") + print("# | | | _| |_) | |_| | | (_| (_| | | (__ #") + print("# |_| |_| | .__/ \__, |___\___\__,_|_|\___| #") + print("# |_| |___/_____| #") + print("# #") + print("##############################################") + print() + + # Run Input Parser + args = get_data_arguments() + + # Load Database + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + # Loop over station keys + for stkey in list(stkeys): + + # Extract station information from dictionary + sta = db[stkey] + + # Construct Folder Name + stfld = stkey + if not args.lkey: + stfld = stkey.split('.')[0]+"."+stkey.split('.')[1] + + # Define path to see if it exists + if args.phase in ['P', 'PP']: + datapath = Path('P_DATA') / stfld + elif args.phase in ['S', 'SKS']: + datapath = Path('S_DATA') / stfld + if not datapath.exists(): + print('Path to '+str(datapath)+' doesn`t exist - creating it') + datapath.mkdir(parents=True) + + # Establish client for data + if len(args.UserAuth) == 0: + data_client = Client(args.Server) + else: + data_client = Client(args.Server, user=args.UserAuth[0], + password=args.UserAuth[1]) + + # Establish client for events + event_client = Client() + + # Get catalogue search start time + if args.startT is None: + tstart = sta.startdate + else: + tstart = args.startT + + # Get catalogue search end time + if args.endT is None: + tend = sta.enddate + else: + tend = args.endT + if tstart > sta.enddate or tend < sta.startdate: + continue + + # Temporary print locations + tlocs = sta.location + if len(tlocs) == 0: + tlocs = [''] + for il in range(0, len(tlocs)): + if len(tlocs[il]) == 0: + tlocs[il] = "--" + sta.location = tlocs + + # Update Display + print(" ") + print(" ") + print("|===============================================|") + print("|===============================================|") + print("| {0:>8s} |".format( + sta.station)) + print("|===============================================|") + print("|===============================================|") + print("| Station: {0:>2s}.{1:5s} |".format( + sta.network, sta.station)) + print("| Channel: {0:2s}; Locations: {1:15s} |".format( + sta.channel, ",".join(tlocs))) + print("| Lon: {0:7.2f}; Lat: {1:6.2f} |".format( + sta.longitude, sta.latitude)) + print("| Start time: {0:19s} |".format( + sta.startdate.strftime("%Y-%m-%d %H:%M:%S"))) + print("| End time: {0:19s} |".format( + sta.enddate.strftime("%Y-%m-%d %H:%M:%S"))) + print("|-----------------------------------------------|") + print("| Searching Possible events: |") + print("| Start: {0:19s} |".format( + tstart.strftime("%Y-%m-%d %H:%M:%S"))) + print("| End: {0:19s} |".format( + tend.strftime("%Y-%m-%d %H:%M:%S"))) + if args.maxmag is None: + print("| Mag: >{0:3.1f}", format(args.minmag) + + " |") + else: + msg = "| Mag: {0:3.1f}".format(args.minmag) + \ + " - {0:3.1f}".format(args.maxmag) + \ + " |" + print(msg) + + print("| ... |") + + # Get catalogue using deployment start and end + cat = event_client.get_events( + starttime=tstart, endtime=tend, + minmagnitude=args.minmag, maxmagnitude=args.maxmag) + + # Total number of events in Catalogue + nevK = 0 + nevtT = len(cat) + print( + "| Found {0:5d}".format(nevtT) + + " possible events |") + ievs = range(0, nevtT) + + # Get Local Data Availabilty + if len(args.localdata) > 0: + print("|-----------------------------------------------|") + print("| Cataloging Local Data... |") + if args.useNet: + stalcllist = utils.list_local_data_stn( + lcldrs=args.localdata, sta=sta.station, + net=sta.network, dtype=args.dtype, altnet=sta.altnet) + print("| {0:>2s}.{1:5s}: {2:6d}".format( + sta.network, sta.station, len(stalcllist)) + + " files |") + #print(stalcllist[0:10]) + else: + stalcllist = utils.list_local_data_stn( + lcldrs=args.localdata, sta=sta.station, dtype=args.dtype) + print("| {0:5s}: {1:6d} files " + + " |".format( + sta.station, len(stalcllist))) + else: + stalcllist = [] + print("|===============================================|") + + # Select order of processing + if args.reverse: + ievs = range(0, nevtT) + else: + ievs = range(nevtT-1, -1, -1) + + # Read through catalogue + for iev in ievs: + + # Extract event + ev = cat[iev] + + # Initialize RF object with station info + rfdata = RFData(sta) + + # Add event to rfdata object + accept = rfdata.add_event( + ev, gacmin=args.mindist, gacmax=args.maxdist, + phase=args.phase, returned=True) + + # Define time stamp + yr = str(rfdata.meta.time.year).zfill(4) + jd = str(rfdata.meta.time.julday).zfill(3) + hr = str(rfdata.meta.time.hour).zfill(2) + + # If event is accepted (data exists) + if accept: + + # Display Event Info + nevK = nevK + 1 + if args.reverse: + inum = iev + 1 + else: + inum = nevtT - iev + 1 + print(" ") + print("*"*50) + print("* #{0:d} ({1:d}/{2:d}): {3:13s} {4}".format( + nevK, inum, nevtT, rfdata.meta.time.strftime( + "%Y%m%d_%H%M%S"), stkey)) + if args.verb: + print("* Phase: {}".format(args.phase)) + print("* Origin Time: " + + rfdata.meta.time.strftime("%Y-%m-%d %H:%M:%S")) + print( + "* Lat: {0:6.2f}; Lon: {1:7.2f}".format( + rfdata.meta.lat, rfdata.meta.lon)) + print( + "* Dep: {0:6.2f} km; Mag: {1:3.1f}".format( + rfdata.meta.dep, rfdata.meta.mag)) + print( + "* Dist: {0:7.2f} km;".format(rfdata.meta.epi_dist) + + " Epi dist: {0:6.2f} deg\n".format(rfdata.meta.gac) + + "* Baz: {0:6.2f} deg;".format(rfdata.meta.baz) + + " Az: {0:6.2f} deg".format(rfdata.meta.az)) + + # Event Folder + timekey = rfdata.meta.time.strftime("%Y%m%d_%H%M%S") + evtdir = datapath / timekey + ZNEfile = evtdir / 'ZNE_Data.pkl' + metafile = evtdir / 'Meta_Data.pkl' + stafile = evtdir / 'Station_Data.pkl' + + # Check if RF data already exist and overwrite has been set + if evtdir.exists(): + if RFfile.exists(): + if not args.ovr: + continue + + # Get data + has_data = rfdata.download_data( + client=data_client, dts=args.dts, stdata=stalcllist, + ndval=args.ndval, dtype=args.dtype, new_sr=args.new_sampling_rate, + returned=True, verbose=args.verb) + + if not has_data: + continue + + # Create Folder if it doesn't exist + if not evtdir.exists(): + evtdir.mkdir(parents=True) + + # Save ZNE Traces + pickle.dump(rfdata.data, open(ZNEfile, "wb")) + + # Save Z12 if components exist + if hasattr(rfdata, "dataZ12"): + Z12file = evtdir / 'Z12_Data.pkl' + pickle.dump(rfdata.dataZ12, open(Z12file, "wb")) + + # Save event meta data + pickle.dump(rfdata.meta, open(metafile, "wb")) + + # Save Station Data + pickle.dump(rfdata.sta, open(stafile, "wb")) + + # Update + if args.verb: + print("* Wrote Output Files to: ") + print("* "+str(evtdir)) + print("*"*50) + + +if __name__ == "__main__": + + # Run main program + main() diff --git a/setup.py b/setup.py index c9f78dc..89e94d5 100644 --- a/setup.py +++ b/setup.py @@ -36,8 +36,10 @@ def find_version(*paths): python_requires='>=3.6', packages=setuptools.find_packages(), entry_points={'console_scripts':[ - 'rfpy_calc=rfpy.scripts.rfpy_calc:main', - 'rfpy_recalc=rfpy.scripts.rfpy_recalc:main', + 'rfpy_get_data=rfpy.scripts.rfpy_get_data:main', + 'rfpy_calc_simdec=rfpy.scripts.rfpy_calc_simdec:main', + # 'rfpy_recalc=rfpy.scripts.rfpy_recalc:main', + 'rfpy_calc_sse=rfpy.scripts.rfpy_recalc_sse:main', 'rfpy_plot=rfpy.scripts.rfpy_plot:main', 'rfpy_harmonics=rfpy.scripts.rfpy_harmonics:main', 'rfpy_hk=rfpy.scripts.rfpy_hk:main', From f3542471d50988c7864ddfb5ba0a0442d571995d Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 17 Nov 2021 12:29:40 -0500 Subject: [PATCH 19/54] re-inserted skipped commits --- rfpy/rfdata.py | 38 ++++++++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 9c393ec..743ea84 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -338,6 +338,16 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, """ + def _resample(tr, new_sr): + + if tr.stats.sampling_rate > new_sr: + # only filter when downsampling + tr.filter('lowpass', freq=0.5*new_sr, corners=2, zerophase=True) + + tr.resample(new_sr, no_filter=True) + return tr + + if self.meta is None: raise(Exception("Requires event data as attribute - aborting")) @@ -368,13 +378,19 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, trZ = stream.select(component='Z')[0] self.data = Stream(traces=[trZ, trN, trE]) - # Filter Traces and resample - self.data.filter('lowpass', freq=0.5*new_sr, - corners=2, zerophase=True) + # Filter Traces when downsampling + if self.data[0].stats.sampling_rate > new_sr: + self.data.filter('lowpass', freq=0.5*new_sr, + corners=2, zerophase=True) + # Resample self.data.resample(new_sr, no_filter=True) # If there is no ZNE, perhaps there is Z12? - except: + except Exception as e: + if verbose: + print('* Met exception: ') + print('* ' + str(e)) + print('* Trying to access Z12 data.') try: tr1 = stream.select(component='1')[0] @@ -382,9 +398,11 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, trZ = stream.select(component='Z')[0] self.data = Stream(traces=[trZ, tr1, tr2]) - # Filter Traces and resample - self.data.filter('lowpass', freq=0.5*new_sr, - corners=2, zerophase=True) + # Filter Traces when downsampling + if self.data[0].stats.sampling_rate > new_sr: + self.data.filter('lowpass', freq=0.5*new_sr, + corners=2, zerophase=True) + # Resample self.data.resample(new_sr, no_filter=True) # Save Z12 components in case it's necessary for later @@ -393,7 +411,11 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, # Rotate from Z12 to ZNE using StDb azcorr attribute self.rotate(align='ZNE') - except: + except Exception as e: + if verbose: + print('* Met exception: ') + print('* ' + str(e)) + print('* Cannot process') self.meta.accept = False if returned: From 699cbbbb3d5ed7cf38faf0e6bd98c7187f5b7c6b Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Thu, 18 Nov 2021 13:59:18 -0500 Subject: [PATCH 20/54] next big commit --- rfpy/plotting.py | 72 ++++ rfpy/rfdata.py | 9 +- rfpy/scripts/rfpy_calc.py | 682 +++++++------------------------ rfpy/scripts/rfpy_calc_simdec.py | 93 ++++- rfpy/scripts/rfpy_calc_sse.py | 412 ------------------- rfpy/scripts/rfpy_plot.py | 28 +- 6 files changed, 335 insertions(+), 961 deletions(-) delete mode 100644 rfpy/scripts/rfpy_calc_sse.py diff --git a/rfpy/plotting.py b/rfpy/plotting.py index 11ac6c0..3e9dafd 100644 --- a/rfpy/plotting.py +++ b/rfpy/plotting.py @@ -447,6 +447,78 @@ def wiggle_single_event(rfdata, filt=None, pre_filt=None, trange=None): plt.show() +def panel(stream1, stream2, sort=None, tmin=0., tmax=30, normalize=True, + save=False, title=None, form='png'): + + taxis = stream1[0].stats.taxis + tselect = (taxis > tmin) & (taxis < tmax) + nt = len(taxis[tselect]) + # print(nt) + # print(taxis[tselect]) + # print(len(stream1[0].data)) + # print(len(stream1[0].data[tselect])) + + baz_array = np.array([[i, tr.stats.baz] for i, tr in enumerate(stream1)]) + slow_array = np.array([[i, tr.stats.slow] for i, tr in enumerate(stream1)]) + + # build panel + panelV = np.zeros((len(stream1), nt)) + for i, tr in enumerate(stream1): + panelV[i,:] = tr.data[tselect] + panelH = np.zeros((len(stream2), nt)) + for i, tr in enumerate(stream2): + panelH[i,:] = tr.data[tselect] + + if normalize: + vmin = np.mean(panelV) - 3.*np.std(panelV) + vmax = np.mean(panelV) + 3.*np.std(panelV) + else: + vmin = -0.5 + vmax = 0.5 + + plt.figure(figsize=(8,5)) + + plt.subplot(311) + plt.pcolor(panelV.T, cmap='bwr', vmin=vmin, vmax=vmax) + plt.xlim((0,len(stream1))) + plt.ylim((0,nt)) + plt.gca().invert_yaxis() + plt.ylabel('Time (s)') + plt.gca().set_yticks(np.arange(0., 600.05, 200.)) + plt.gca().set_yticklabels([str(0), str(10), str(20), str(30)]) + plt.gca().set_xticklabels([]) + + plt.subplot(312) + plt.pcolor(panelH.T, cmap='bwr', vmin=vmin, vmax=vmax) + plt.xlim((0,len(stream1))) + plt.ylim((0,nt)) + plt.gca().invert_yaxis() + plt.ylabel('Time (s)') + plt.gca().set_yticks(np.arange(0., 600.05, 200.)) + plt.gca().set_yticklabels([str(0), str(10), str(20), str(30)]) + plt.gca().set_xticklabels([]) + + plt.subplot(313) + plt.plot(baz_array[:,0], baz_array[:,1]) + plt.xlim((0,len(stream1))) + plt.ylim((0,360)) + plt.ylabel('Back-azimuth (deg)') + plt.xlabel('Receiver function #') + plt.twinx() + plt.plot(slow_array[:,0], slow_array[:,1], '.') + plt.ylabel('Slowness (s/km)') + plt.grid() + + if title: + plt.suptitle(title) + + if save: + plt.savefig('RF_PLOTS/' + stream1[0].stats.station + + '.' + title + '.' + form, format=form) + else: + plt.show() + + def event_dist(stream, phase='P', save=False, title=None, form='png'): import cartopy.crs as ccrs diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 743ea84..9e37d1f 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -1133,14 +1133,15 @@ def _add_stats(trace, store): trace.stats.vs = self.meta.vs trace.stats.phase = self.meta.phase trace.stats.align = self.meta.align - nn = self.data[0].stats.npts - sr = self.data[0].stats.sampling_rate - trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn) - if store == 'rf': trace.stats.is_rf = True + nn = self.rf[0].stats.npts + sr = self.rf[0].stats.sampling_rate elif store == 'specs': trace.stats.is_specs = True + nn = self.specs[0].stats.npts + sr = self.specs[0].stats.sampling_rate + trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn) return trace diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 251fde5..74cb57a 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -24,19 +24,16 @@ # Import modules and functions +from argparse import ArgumentParser +from os.path import exists as exist import numpy as np import pickle import stdb -from obspy.clients.fdsn import Client -from rfpy import utils, RFData +from rfpy import RFData from pathlib import Path -from argparse import ArgumentParser -from os.path import exists as exist -from obspy import UTCDateTime -from numpy import nan -def get_calc_arguments(argv=None): +def get_recalc_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. @@ -46,13 +43,12 @@ def get_calc_arguments(argv=None): parser = ArgumentParser( usage="%(prog)s [arguments] ", - description="Script used to download and pre-process " + - "three-component ('Z', 'N', and 'E'), seismograms for individual " + - "events and calculate teleseismic P-wave receiver functions" + - "This version requests data on the fly for a given date " + - "range. Data are requested from the internet using the " + - "client services framework. The stations are processed one " + - "by one and the data are stored to disk.") + description="Script used to re-calculate receiver functions " + + "that already exist on disk, but using different " + + "processing options. The stations are processed one " + + "by one and the data are stored to disk. " + + " \n" + + "Note: The sampling rate cannot be changed to a new rate") # General Settings parser.add_argument( @@ -77,13 +73,6 @@ def get_calc_arguments(argv=None): dest="verb", default=False, help="Specify to increase verbosity.") - parser.add_argument( - "-O", "--overwrite", - action="store_true", - dest="ovr", - default=False, - help="Force the overwriting of pre-existing data. " + - "[Default False]") parser.add_argument( "-L", "--long-name", action="store_true", @@ -94,183 +83,34 @@ def get_calc_arguments(argv=None): "names, regardless of the key type of the database." ) - # Server Settings - ServerGroup = parser.add_argument_group( - title="Server Settings", - description="Settings associated with which " - "datacenter to log into.") - ServerGroup.add_argument( - "-S", "--Server", - action="store", - type=str, - dest="Server", - default="IRIS", - help="Specify the server to connect to. Options include: " + - "BGR, ETH, GEONET, GFZ, INGV, IPGP, IRIS, KOERI, " + - "LMU, NCEDC, NEIP, NERIES, ODC, ORFEUS, RESIF, SCEDC, USGS, USP. " + - "[Default IRIS]") - ServerGroup.add_argument( - "-U", "--User-Auth", - action="store", - type=str, - dest="UserAuth", - default="", - help="Enter your IRIS Authentification Username and Password " + - "(--User-Auth='username:authpassword') to " + - "access and download restricted data. " + - "[Default no user and password]") - - # Database Settings - DataGroup = parser.add_argument_group( - title="Local Data Settings", - description="Settings associated with defining " + - "and using a local data base of pre-downloaded " + - "day-long SAC or MSEED files.") - DataGroup.add_argument( - "--local-data", - action="store", - type=str, - dest="localdata", - default=None, - help="Specify a comma separated list of paths containing " + - "day-long sac or mseed files of data already downloaded. " + - "If data exists for a seismogram is already present on disk, " + - "it is selected preferentially over downloading " + - "the data using the Client interface") - DataGroup.add_argument( - "--dtype", - action="store", - type=str, - dest="dtype", - default='SAC', - help="Specify the data archive file type, either SAC " + - " or MSEED. Note the default behaviour is to search for " + - "SAC files. Local archive files must have extensions of '.SAC' "+ - " or '.MSEED. These are case dependent, so specify the correct case"+ - "here.") - DataGroup.add_argument( - "--no-data-zero", - action="store_true", - dest="ndval", - default=False, - help="Specify to force missing data to be set as zero, rather " + - "than default behaviour which sets to nan. Note this is applied " + - "only to the SAC data") - DataGroup.add_argument( - "--no-local-net", - action="store_false", - dest="useNet", - default=True, - help="Specify to prevent using the Network code in the " + - "search for local data (sometimes for CN stations " + - "the dictionary name for a station may disagree with that " + - "in the filename. [Default Network used]") - DataGroup.add_argument( - "--save-Z12", - action="store_true", - dest="saveZ12", - default=False, - help="Specify to save Z12 (un-rotated) components. [Default " + - "False]") - - # Event Selection Criteria - EventGroup = parser.add_argument_group( - title="Event Settings", - description="Settings associated with refining " + - "the events to include in matching event-station pairs") - EventGroup.add_argument( - "--start", - action="store", - type=str, - dest="startT", - default="", - help="Specify a UTCDateTime compatible string representing " + - "the start time for the event search. This will override any " + - "station start times. [Default start date of station]") - EventGroup.add_argument( - "--end", - action="store", - type=str, - dest="endT", - default="", - help="Specify a UTCDateTime compatible string representing " + - "the end time for the event search. This will override any " + - "station end times [Default end date of station]") - EventGroup.add_argument( - "--reverse", "-R", + # Constants Settings + ConstGroup = parser.add_argument_group( + title='Parameter Settings', + description="Miscellaneous default values and settings") + ConstGroup.add_argument( + "--Z12", action="store_true", - dest="reverse", + dest="Z12", default=False, - help="Reverse order of events. Default behaviour starts at " + - "oldest event and works towards most recent. Specify reverse " + - "order and instead the program will start with the most recent " + - "events and work towards older") - EventGroup.add_argument( - "--minmag", - action="store", - type=float, - dest="minmag", - default=6.0, - help="Specify the minimum magnitude of event for which to search. " + - "[Default 6.0]") - EventGroup.add_argument( - "--maxmag", - action="store", - type=float, - dest="maxmag", - default=9.0, - help="Specify the maximum magnitude of event for which to search. " + - "[Default None, i.e. no limit]") - - # Geometry Settings - PhaseGroup = parser.add_argument_group( - title="Geometry Settings", - description="Settings associatd with the " - "event-station geometries for the specified phase") - PhaseGroup.add_argument( + help="Use Z12 data if available. [Default uses ZNE data]") + ConstGroup.add_argument( "--phase", action="store", type=str, dest="phase", - default='P', + default='allP', help="Specify the phase name to use. Be careful with the distance. " + - "setting. Options are 'P' or 'PP'. [Default 'P']") - PhaseGroup.add_argument( - "--mindist", - action="store", - type=float, - dest="mindist", - default=None, - help="Specify the minimum great circle distance (degrees) between " + - "the station and event. [Default depends on phase]") - PhaseGroup.add_argument( - "--maxdist", - action="store", - type=float, - dest="maxdist", - default=None, - help="Specify the maximum great circle distance (degrees) between " + - "the station and event. [Default depends on phase]") - - # Constants Settings - ConstGroup = parser.add_argument_group( - title='Parameter Settings', - description="Miscellaneous default values and settings") + "setting. Options are 'P', 'PP', 'allP', 'S', 'SKS' or 'allS'. " + + "[Default 'allP']") ConstGroup.add_argument( - "--sampling-rate", + "--resample", action="store", type=float, - dest="new_sampling_rate", - default=10., - help="Specify new sampling rate in Hz. [Default 10.]") - ConstGroup.add_argument( - "--dts", - action="store", - type=float, - dest="dts", - default=150., - help="Specify the window length in sec (symmetric about arrival " + - "time). [Default 150.]") + dest="resample", + default=None, + help="Specify the new sampling-rate for the receiver functions. " + + "Note the sampling rate of the original data (ZNE or Z12) stored " + + "on disk is unchanged. [Default None]") ConstGroup.add_argument( "--align", action="store", @@ -317,7 +157,7 @@ def get_calc_arguments(argv=None): type=float, dest="fmin", default=0.05, - help="Specify the minimum frequency corner for SNR and CC " + + help="Specify the minimum frequency corner for SNR " + "filter (Hz). [Default 0.05]") ConstGroup.add_argument( "--fmax", @@ -325,10 +165,10 @@ def get_calc_arguments(argv=None): type=float, dest="fmax", default=1.0, - help="Specify the maximum frequency corner for SNR and CC " + + help="Specify the maximum frequency corner for SNR " + "filter (Hz). [Default 1.0]") - # Constants Settings + # Deconvolution Settings DeconGroup = parser.add_argument_group( title='Deconvolution Settings', description="Parameters for deconvolution") @@ -367,109 +207,16 @@ def get_calc_arguments(argv=None): if len(args.stkeys) > 0: args.stkeys = args.stkeys.split(',') - # construct start time - if len(args.startT) > 0: - try: - args.startT = UTCDateTime(args.startT) - except: - parser.error( - "Cannot construct UTCDateTime from start time: " + - args.startT) - else: - args.startT = None - - # construct end time - if len(args.endT) > 0: - try: - args.endT = UTCDateTime(args.endT) - except: - parser.error( - "Cannot construct UTCDateTime from end time: " + - args.endT) - else: - args.endT = None - - # Parse User Authentification - if not len(args.UserAuth) == 0: - tt = args.UserAuth.split(':') - if not len(tt) == 2: - parser.error( - "Error: Incorrect Username and Password " + - "Strings for User Authentification") - else: - args.UserAuth = tt - else: - args.UserAuth = [] - - # Parse Local Data directories - if args.localdata is not None: - args.localdata = args.localdata.split(',') - else: - args.localdata = [] - - # Check Datatype specification - if (not args.dtype.upper() == 'MSEED') and (not args.dtype.upper() == 'SAC'): - parser.error( - "Error: Local Data Archive must be of types 'SAC'" + - "or MSEED. These must match the file extensions for " + - " the archived data.") - - # Check NoData Value - if args.ndval: - args.ndval = 0.0 - else: - args.ndval = nan - - # Check distances for selected phase - if args.phase not in ['P', 'PP', 'S', 'SKS']: + if args.phase not in ['P', 'PP', 'allP', 'S', 'SKS', 'allS']: parser.error( - "Error: choose between 'P', 'PP', 'S' and 'SKS'.") - if args.phase == 'P': - if not args.mindist: - args.mindist = 30. - if not args.maxdist: - args.maxdist = 100. - if args.mindist < 30. or args.maxdist > 100.: - parser.error( - "Distances should be between 30 and 100 deg. for " + - "teleseismic 'P' waves.") - elif args.phase == 'PP': - if not args.mindist: - args.mindist = 100. - if not args.maxdist: - args.maxdist = 180. - if args.mindist < 100. or args.maxdist > 180.: - parser.error( - "Distances should be between 100 and 180 deg. for " + - "teleseismic 'PP' waves.") - elif args.phase == 'S': - if not args.mindist: - args.mindist = 55. - if not args.maxdist: - args.maxdist = 85. - if args.mindist < 55. or args.maxdist > 85.: - parser.error( - "Distances should be between 55 and 85 deg. for " + - "teleseismic 'S' waves.") - elif args.phase == 'SKS': - if not args.mindist: - args.mindist = 85. - if not args.maxdist: - args.maxdist = 115. - if args.mindist < 85. or args.maxdist > 115.: - parser.error( - "Distances should be between 85 and 115 deg. for " + - "teleseismic 'SKS' waves.") - - if args.pre_filt is not None: - args.pre_filt = [float(val) for val in args.pre_filt.split(',')] - args.pre_filt = sorted(args.pre_filt) - if (len(args.pre_filt)) != 2: - parser.error( - "Error: --pre-filt should contain 2 " + - "comma-separated floats") + "Error: choose between 'P', 'PP', 'allP', 'S', 'SKS' and 'allS'.") + if args.phase == 'allP': + args.listphase = ['P', 'PP'] + elif args.phase == 'allS': + args.listphase = ['S', 'SKS'] + else: + args.listphase = [args.phase] - # Check alignment arguments if args.align is None: args.align = 'ZRT' elif args.align not in ['ZRT', 'LQT', 'PVH']: @@ -477,39 +224,46 @@ def get_calc_arguments(argv=None): "Error: Incorrect alignment specifier. Should be " + "either 'ZRT', 'LQT', or 'PVH'.") - if args.dt_snr > args.dts: - args.dt_snr = args.dts - 10. - print("SNR window > data window. Defaulting to data " + - "window minus 10 sec.") - if args.method not in ['wiener', 'water', 'multitaper']: parser.error( "Error: 'method' should be either 'wiener', 'water' or " + "'multitaper'") + if args.pre_filt is not None: + args.pre_filt = [float(val) for val in args.pre_filt.split(',')] + args.pre_filt = sorted(args.pre_filt) + if (len(args.pre_filt)) != 2: + parser.error( + "Error: --pre-filt should contain 2 " + + "comma-separated floats") + return args def main(): print() - print("##############################################") - print("# __ _ #") - print("# _ __ / _|_ __ _ _ ___ __ _| | ___ #") - print("# | '__| |_| '_ \| | | | / __/ _` | |/ __| #") - print("# | | | _| |_) | |_| | | (_| (_| | | (__ #") - print("# |_| |_| | .__/ \__, |___\___\__,_|_|\___| #") - print("# |_| |___/_____| #") - print("# #") - print("##############################################") + print("########################################################") + print("# #") + print("# __ _ #") + print("# _ __ / _|_ __ _ _ _ __ ___ ___ __ _| | ___ #") + print("# | '__| |_| '_ \| | | | | '__/ _ \/ __/ _` | |/ __| #") + print("# | | | _| |_) | |_| | | | | __/ (_| (_| | | (__ #") + print("# |_| |_| | .__/ \__, |___|_| \___|\___\__,_|_|\___| #") + print("# |_| |___/_____| #") + print("# #") + print("########################################################") print() # Run Input Parser - args = get_calc_arguments() + args = get_recalc_arguments() # Load Database db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Track processed folders + procfold = [] + # Loop over station keys for stkey in list(stkeys): @@ -522,36 +276,12 @@ def main(): stfld = stkey.split('.')[0]+"."+stkey.split('.')[1] # Define path to see if it exists - if args.phase in ['P', 'PP']: + if args.phase in ['P', 'PP', 'allP']: datapath = Path('P_DATA') / stfld - elif args.phase in ['S', 'SKS']: + elif args.phase in ['S', 'SKS', 'allS']: datapath = Path('S_DATA') / stfld - if not datapath.exists(): - print('Path to '+str(datapath)+' doesn`t exist - creating it') - datapath.mkdir(parents=True) - - # Establish client for data - if len(args.UserAuth) == 0: - data_client = Client(args.Server) - else: - data_client = Client(args.Server, user=args.UserAuth[0], - password=args.UserAuth[1]) - - # Establish client for events - event_client = Client() - - # Get catalogue search start time - if args.startT is None: - tstart = sta.startdate - else: - tstart = args.startT - - # Get catalogue search end time - if args.endT is None: - tend = sta.enddate - else: - tend = args.endT - if tstart > sta.enddate or tend < sta.startdate: + if not datapath.is_dir(): + print('Path to ' + str(datapath) + ' doesn`t exist - continuing') continue # Temporary print locations @@ -567,213 +297,113 @@ def main(): print(" ") print(" ") print("|===============================================|") - print("|===============================================|") print("| {0:>8s} |".format( sta.station)) print("|===============================================|") - print("|===============================================|") print("| Station: {0:>2s}.{1:5s} |".format( sta.network, sta.station)) print("| Channel: {0:2s}; Locations: {1:15s} |".format( sta.channel, ",".join(tlocs))) print("| Lon: {0:7.2f}; Lat: {1:6.2f} |".format( sta.longitude, sta.latitude)) - print("| Start time: {0:19s} |".format( - sta.startdate.strftime("%Y-%m-%d %H:%M:%S"))) - print("| End time: {0:19s} |".format( - sta.enddate.strftime("%Y-%m-%d %H:%M:%S"))) print("|-----------------------------------------------|") - print("| Searching Possible events: |") - print("| Start: {0:19s} |".format( - tstart.strftime("%Y-%m-%d %H:%M:%S"))) - print("| End: {0:19s} |".format( - tend.strftime("%Y-%m-%d %H:%M:%S"))) - if args.maxmag is None: - print("| Mag: >{0:3.1f}", format(args.minmag) + - " |") - else: - msg = "| Mag: {0:3.1f}".format(args.minmag) + \ - " - {0:3.1f}".format(args.maxmag) + \ - " |" - print(msg) - - print("| ... |") - - # Get catalogue using deployment start and end - cat = event_client.get_events( - starttime=tstart, endtime=tend, - minmagnitude=args.minmag, maxmagnitude=args.maxmag) - - # Total number of events in Catalogue - nevK = 0 - nevtT = len(cat) - print( - "| Found {0:5d}".format(nevtT) + - " possible events |") - ievs = range(0, nevtT) - - # Get Local Data Availabilty - if len(args.localdata) > 0: - print("|-----------------------------------------------|") - print("| Cataloging Local Data... |") - if args.useNet: - stalcllist = utils.list_local_data_stn( - lcldrs=args.localdata, sta=sta.station, - net=sta.network, dtype=args.dtype, altnet=sta.altnet) - print("| {0:>2s}.{1:5s}: {2:6d}".format( - sta.network, sta.station, len(stalcllist)) + - " files |") - #print(stalcllist[0:10]) - else: - stalcllist = utils.list_local_data_stn( - lcldrs=args.localdata, sta=sta.station, dtype=args.dtype) - print("| {0:5s}: {1:6d} files " + - " |".format( - sta.station, len(stalcllist))) - else: - stalcllist = [] - print("|===============================================|") - # Select order of processing - if args.reverse: - ievs = range(0, nevtT) - else: - ievs = range(nevtT-1, -1, -1) + # Check for folder already processed + if stfld in procfold: + print(' {0} already processed...skipping '.format(stfld)) + continue - # Read through catalogue - for iev in ievs: + datafiles = [x for x in datapath.iterdir() if x.is_dir()] + for folder in datafiles: - # Extract event - ev = cat[iev] + # Skip hidden folders + if folder.name.startswith('.'): + continue - # Initialize RF object with station info + # Re-initialize RFData object rfdata = RFData(sta) - # Add event to rfdata object - accept = rfdata.add_event( - ev, gacmin=args.mindist, gacmax=args.maxdist, - phase=args.phase, returned=True) - - # Define time stamp - yr = str(rfdata.meta.time.year).zfill(4) - jd = str(rfdata.meta.time.julday).zfill(3) - hr = str(rfdata.meta.time.hour).zfill(2) - - # If event is accepted (data exists) - if accept: - - # Display Event Info - nevK = nevK + 1 - if args.reverse: - inum = iev + 1 - else: - inum = nevtT - iev + 1 - print(" ") - print("*"*50) - print("* #{0:d} ({1:d}/{2:d}): {3:13s} {4}".format( - nevK, inum, nevtT, rfdata.meta.time.strftime( - "%Y%m%d_%H%M%S"), stkey)) - if args.verb: - print("* Phase: {}".format(args.phase)) - print("* Origin Time: " + - rfdata.meta.time.strftime("%Y-%m-%d %H:%M:%S")) - print( - "* Lat: {0:6.2f}; Lon: {1:7.2f}".format( - rfdata.meta.lat, rfdata.meta.lon)) - print( - "* Dep: {0:6.2f} km; Mag: {1:3.1f}".format( - rfdata.meta.dep, rfdata.meta.mag)) - print( - "* Dist: {0:7.2f} km;".format(rfdata.meta.epi_dist) + - " Epi dist: {0:6.2f} deg\n".format(rfdata.meta.gac) + - "* Baz: {0:6.2f} deg;".format(rfdata.meta.baz) + - " Az: {0:6.2f} deg".format(rfdata.meta.az)) - - # Event Folder - timekey = rfdata.meta.time.strftime("%Y%m%d_%H%M%S") - evtdir = datapath / timekey - RFfile = evtdir / 'RF_Data.pkl' - ZNEfile = evtdir / 'ZNE_Data.pkl' - metafile = evtdir / 'Meta_Data.pkl' - stafile = evtdir / 'Station_Data.pkl' - - # Check if RF data already exist and overwrite has been set - if evtdir.exists(): - if RFfile.exists(): - if not args.ovr: - continue - - # Get data - has_data = rfdata.download_data( - client=data_client, dts=args.dts, stdata=stalcllist, - ndval=args.ndval, dtype=args.dtype, new_sr=args.new_sampling_rate, - returned=True, verbose=args.verb) - - if not has_data: - continue - - # Create Folder if it doesn't exist - if not evtdir.exists(): - evtdir.mkdir(parents=True) - - # Save ZNE Traces - pickle.dump(rfdata.data, open(ZNEfile, "wb")) - - # Save Z12 if components exist - if hasattr(rfdata, "dataZ12"): - Z12file = evtdir / 'Z12_Data.pkl' - pickle.dump(rfdata.dataZ12, open(Z12file, "wb")) - - # Rotate from ZNE to 'align' ('ZRT', 'LQT', or 'PVH') - rfdata.rotate(vp=args.vp, vs=args.vs, align=args.align) - - # Calculate snr over dt_snr seconds - rfdata.calc_snr( - dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax) - if args.verb: - print("* SNR: {}".format(rfdata.meta.snr)) - - # Make sure no processing happens for NaNs - if np.isnan(rfdata.meta.snr): - if args.verb: - print("* SNR NaN...Skipping") - print("*"*50) - continue - - # Calculate spectra - rfdata.calc_spectra( - vp=args.vp, vs=args.vs, - align=args.align, method=args.method, - wavelet='complete', envelope_threshold=0.05, - time=5) - - # Deconvolve data - rfdata.deconvolve( - method=args.method, gfilt=args.gfilt, - wlevel=args.wlevel) - - # Get cross-correlation QC - rfdata.calc_cc() - if args.verb: - print("* CC: {}".format(rfdata.meta.cc)) - - # Convert to Stream - rfstream = rfdata.to_stream(store='rf') - - # Save event meta data - pickle.dump(rfdata.meta, open(metafile, "wb")) - - # Save Station Data - pickle.dump(rfdata.sta, open(stafile, "wb")) - - # Save RF Traces - pickle.dump(rfstream, open(RFfile, "wb")) - - # Update - if args.verb: - print("* Wrote Output Files to: ") - print("* "+str(evtdir)) - print("*"*50) + # Load meta data + metafile = folder / "Meta_Data.pkl" + if not metafile.is_file(): + continue + rfdata.meta = pickle.load(open(metafile, 'rb')) + + # Skip data not in list of phases + if rfdata.meta.phase not in args.listphase: + continue + + if args.verb: + print("* Station: {0}; folder: {1}".format(stkey, folder)) + + if args.Z12: + try: + # Try loading Z12_Data + Z12file = folder / "Z12_Data.pkl" + rfdata.data = pickle.load(open(Z12file, 'rb')) + # Remove rotated flag and snr flag + rfdata.meta.rotated = False + rfdata.rotate(align='ZNE') + except: + + print("Z12_Data.pkl not available - using ZNE_Data.pkl") + # Load ZNE data + ZNEfile = folder / "ZNE_Data.pkl" + rfdata.data = pickle.load(open(ZNEfile, 'rb')) + else: + # Load ZNE data + ZNEfile = folder / "ZNE_Data.pkl" + rfdata.data = pickle.load(open(ZNEfile, 'rb')) + + # Resample if requested + if args.resample: + rfdata.data.resample( + args.resample, no_filter=False) + + # Remove rotated flag and snr flag + rfdata.meta.rotated = False + rfdata.meta.snr = None + + # Rotate from ZNE to 'align' ('ZRT', 'LQT', or 'PVH') + rfdata.rotate(vp=args.vp, vs=args.vs, align=args.align) + + # Calculate SNR + rfdata.calc_snr(dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax) + + if args.verb: + print("* SNR: {}".format(rfdata.meta.snr)) + + # Deconvolve data + rfdata.deconvolve( + vp=args.vp, vs=args.vs, + align=args.align, method=args.method, + gfilt=args.gfilt, wlevel=args.wlevel, + pre_filt=args.pre_filt) + + # Get cross-correlation QC + rfdata.calc_cc() + + if args.verb: + print("* CC: {}".format(rfdata.meta.cc)) + + # Convert to Stream + rfstream = rfdata.to_stream() + + # Save RF Traces + RFfile = folder / "RF_Data.pkl" + pickle.dump(rfstream, open(RFfile, "wb")) + + # Save Meta data + metafile = folder / "Meta_Data.pkl" + pickle.dump(rfdata.meta, open(metafile, 'wb')) + + # Update + if args.verb: + print("* Output files written") + print("**************************************************") + + # Update processed folders + procfold.append(stfld) if __name__ == "__main__": diff --git a/rfpy/scripts/rfpy_calc_simdec.py b/rfpy/scripts/rfpy_calc_simdec.py index b9bf65f..eb4c69f 100644 --- a/rfpy/scripts/rfpy_calc_simdec.py +++ b/rfpy/scripts/rfpy_calc_simdec.py @@ -229,12 +229,60 @@ def get_simdec_arguments(argv=None): action="store", dest="nslow", type=int, - default=16, + default=20, help="Specify integer number of slowness bins to consider " + "(range of slowness values is 0.04 to 0.08). " + "If not None, the plot will show receiver " + - "functions sorted by slowness values. [Default 16]") - + "functions sorted by slowness values. [Default 20]") + PreGroup.add_argument( + "--bp", + action="store", + type=str, + dest="bp", + default=None, + help="Specify the corner frequencies for the bandpass filter. " + + "[Default no filtering]") + + PlotGroup = parser.add_argument_group( + title='Plot Settings', + description="Options for plot format") + PlotGroup.add_argument( + "--normalize", + action="store_true", + dest="norm", + default=False, + help="Set this option to produce receiver functions normalized " + + "by the max amplitude of stacked RFs. [Default False]") + PlotGroup.add_argument( + "--trange", + action="store", + default=None, + type=str, + dest="trange", + help="Specify the time range for the x-axis (sec). Negative times " + + "are allowed [Default 0., 30.]") + PlotGroup.add_argument( + "--save", + action="store_true", + dest="saveplot", + default=False, + help="Set this option if you wish to save the figure. [Default " + + "does not save figure]") + PlotGroup.add_argument( + "--title", + action="store", + dest="titleplot", + type=str, + default='', + help="Specify title of figure. [Default None]") + PlotGroup.add_argument( + "--format", + action="store", + type=str, + dest="form", + default="png", + help="Specify format of figure. Can be any one of the valid" + + "matplotlib formats: 'png', 'jpg', 'eps', 'pdf'. [Default 'png']") args = parser.parse_args(argv) @@ -268,6 +316,24 @@ def get_simdec_arguments(argv=None): "Error: 'method' should be either 'wiener', 'water' or " + "'multitaper'") + if args.bp is not None: + args.bp = [float(val) for val in args.bp.split(',')] + args.bp = sorted(args.bp) + if (len(args.bp)) != 2: + parser.error( + "Error: --bp should contain 2 " + + "comma-separated floats") + + if args.trange is None: + args.trange = [0., 30.] + else: + args.trange = [float(val) for val in args.trange.split(',')] + args.trange = sorted(args.trange) + if (len(args.trange)) != 2: + parser.error( + "Error: --trange should contain 2 " + + "comma-separated floats") + if args.pre_filt is not None: args.pre_filt = [float(val) for val in args.pre_filt.split(',')] args.pre_filt = sorted(args.pre_filt) @@ -437,15 +503,14 @@ def main(): NNstreams.append(rfstream[3]) # Bin into baz and slowness bins - LL_tmp, = binning.bin_baz_slow(LLstreams, nbaz=args.nbaz, nslow=args.nslow, + LL_tmp, = binning.bin_baz_slow(LLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, pws=False, include_empty=False) - QL_tmp, = binning.bin_baz_slow(QLstreams, nbaz=args.nbaz, nslow=args.nslow, + QL_tmp, = binning.bin_baz_slow(QLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, pws=False, include_empty=False) - TL_tmp, = binning.bin_baz_slow(TLstreams, nbaz=args.nbaz, nslow=args.nslow, + TL_tmp, = binning.bin_baz_slow(TLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, pws=False, include_empty=False) - NN_tmp, = binning.bin_baz_slow(NNstreams, nbaz=args.nbaz, nslow=args.nslow, + NN_tmp, = binning.bin_baz_slow(NNstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, pws=False, include_empty=False) - print(QL_tmp) # Initialize new empty streams RFL = Stream() @@ -471,6 +536,18 @@ def main(): RFQ.append(rfdata.rf[1]) RFT.append(rfdata.rf[2]) + # Filter + if args.bp: + RFQ.filter('bandpass', freqmin=args.bp[0], + freqmax=args.bp[1], corners=2, + zerophase=True) + RFT.filter('bandpass', freqmin=args.bp[0], + freqmax=args.bp[1], corners=2, + zerophase=True) + + plotting.panel(RFQ, RFT, tmin=args.trange[0], tmax=args.trange[1], + normalize=True, save=args.saveplot, + title=args.titleplot, form=args.form) # # Update processed folders # procfold.append(stfld) diff --git a/rfpy/scripts/rfpy_calc_sse.py b/rfpy/scripts/rfpy_calc_sse.py deleted file mode 100644 index 74cb57a..0000000 --- a/rfpy/scripts/rfpy_calc_sse.py +++ /dev/null @@ -1,412 +0,0 @@ -#!/usr/bin/env python - -# Copyright 2019 Pascal Audet -# -# This file is part of RfPy. -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in -# all copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. - - -# Import modules and functions -from argparse import ArgumentParser -from os.path import exists as exist -import numpy as np -import pickle -import stdb -from rfpy import RFData -from pathlib import Path - - -def get_recalc_arguments(argv=None): - """ - Get Options from :class:`~optparse.OptionParser` objects. - - This function is used for data processing on-the-fly (requires web connection) - - """ - - parser = ArgumentParser( - usage="%(prog)s [arguments] ", - description="Script used to re-calculate receiver functions " + - "that already exist on disk, but using different " + - "processing options. The stations are processed one " + - "by one and the data are stored to disk. " + - " \n" + - "Note: The sampling rate cannot be changed to a new rate") - - # General Settings - parser.add_argument( - "indb", - help="Station Database to process from.", - type=str) - parser.add_argument( - "--keys", - action="store", - type=str, - dest="stkeys", - default="", - help="Specify a comma separated list of station keys for " + - "which to perform the analysis. These must be " + - "contained within the station database. Partial keys will " + - "be used to match against those in the dictionary. For " + - "instance, providing IU will match with all stations in " + - "the IU network [Default processes all stations in the database]") - parser.add_argument( - "-v", "-V", "--verbose", - action="store_true", - dest="verb", - default=False, - help="Specify to increase verbosity.") - parser.add_argument( - "-L", "--long-name", - action="store_true", - dest="lkey", - default=False, - help="Force folder names to use long-key form (NET.STN.CHN). " + - "Default behaviour uses short key form (NET.STN) for the folder " + - "names, regardless of the key type of the database." - ) - - # Constants Settings - ConstGroup = parser.add_argument_group( - title='Parameter Settings', - description="Miscellaneous default values and settings") - ConstGroup.add_argument( - "--Z12", - action="store_true", - dest="Z12", - default=False, - help="Use Z12 data if available. [Default uses ZNE data]") - ConstGroup.add_argument( - "--phase", - action="store", - type=str, - dest="phase", - default='allP', - help="Specify the phase name to use. Be careful with the distance. " + - "setting. Options are 'P', 'PP', 'allP', 'S', 'SKS' or 'allS'. " + - "[Default 'allP']") - ConstGroup.add_argument( - "--resample", - action="store", - type=float, - dest="resample", - default=None, - help="Specify the new sampling-rate for the receiver functions. " + - "Note the sampling rate of the original data (ZNE or Z12) stored " + - "on disk is unchanged. [Default None]") - ConstGroup.add_argument( - "--align", - action="store", - type=str, - dest="align", - default=None, - help="Specify component alignment key. Can be either " + - "ZRT, LQT, or PVH. [Default ZRT]") - ConstGroup.add_argument( - "--vp", - action="store", - type=float, - dest="vp", - default=6.0, - help="Specify near-surface Vp to use with --align=PVH (km/s). " + - "[Default 6.0]") - ConstGroup.add_argument( - "--vs", - action="store", - type=float, - dest="vs", - default=3.5, - help="Specify near-surface Vs to use with --align=PVH (km/s). " + - "[Default 3.5]") - ConstGroup.add_argument( - "--dt-snr", - action="store", - type=float, - dest="dt_snr", - default=30., - help="Specify the window length over which to calculate " + - "the SNR in sec. [Default 30.]") - ConstGroup.add_argument( - "--pre-filt", - action="store", - type=str, - dest="pre_filt", - default=None, - help="Specify two floats with low and high frequency corners for " + - "pre-filter (before deconvolution). [Default None]") - ConstGroup.add_argument( - "--fmin", - action="store", - type=float, - dest="fmin", - default=0.05, - help="Specify the minimum frequency corner for SNR " + - "filter (Hz). [Default 0.05]") - ConstGroup.add_argument( - "--fmax", - action="store", - type=float, - dest="fmax", - default=1.0, - help="Specify the maximum frequency corner for SNR " + - "filter (Hz). [Default 1.0]") - - # Deconvolution Settings - DeconGroup = parser.add_argument_group( - title='Deconvolution Settings', - description="Parameters for deconvolution") - DeconGroup.add_argument( - "--method", - action="store", - dest="method", - type=str, - default="wiener", - help="Specify the deconvolution method. Available methods " + - "include 'wiener', 'water' and 'multitaper'. [Default 'wiener']") - DeconGroup.add_argument( - "--gfilt", - action="store", - dest="gfilt", - type=float, - default=None, - help="Specify the Gaussian filter width in Hz. " + - "[Default None]") - DeconGroup.add_argument( - "--wlevel", - action="store", - dest="wlevel", - type=float, - default=0.01, - help="Specify the water level, used in the 'water' method. " + - "[Default 0.01]") - - args = parser.parse_args(argv) - - # Check inputs - if not exist(args.indb): - parser.error("Input file " + args.indb + " does not exist") - - # create station key list - if len(args.stkeys) > 0: - args.stkeys = args.stkeys.split(',') - - if args.phase not in ['P', 'PP', 'allP', 'S', 'SKS', 'allS']: - parser.error( - "Error: choose between 'P', 'PP', 'allP', 'S', 'SKS' and 'allS'.") - if args.phase == 'allP': - args.listphase = ['P', 'PP'] - elif args.phase == 'allS': - args.listphase = ['S', 'SKS'] - else: - args.listphase = [args.phase] - - if args.align is None: - args.align = 'ZRT' - elif args.align not in ['ZRT', 'LQT', 'PVH']: - parser.error( - "Error: Incorrect alignment specifier. Should be " + - "either 'ZRT', 'LQT', or 'PVH'.") - - if args.method not in ['wiener', 'water', 'multitaper']: - parser.error( - "Error: 'method' should be either 'wiener', 'water' or " + - "'multitaper'") - - if args.pre_filt is not None: - args.pre_filt = [float(val) for val in args.pre_filt.split(',')] - args.pre_filt = sorted(args.pre_filt) - if (len(args.pre_filt)) != 2: - parser.error( - "Error: --pre-filt should contain 2 " + - "comma-separated floats") - - return args - - -def main(): - - print() - print("########################################################") - print("# #") - print("# __ _ #") - print("# _ __ / _|_ __ _ _ _ __ ___ ___ __ _| | ___ #") - print("# | '__| |_| '_ \| | | | | '__/ _ \/ __/ _` | |/ __| #") - print("# | | | _| |_) | |_| | | | | __/ (_| (_| | | (__ #") - print("# |_| |_| | .__/ \__, |___|_| \___|\___\__,_|_|\___| #") - print("# |_| |___/_____| #") - print("# #") - print("########################################################") - print() - - # Run Input Parser - args = get_recalc_arguments() - - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) - - # Track processed folders - procfold = [] - - # Loop over station keys - for stkey in list(stkeys): - - # Extract station information from dictionary - sta = db[stkey] - - # Construct Folder Name - stfld = stkey - if not args.lkey: - stfld = stkey.split('.')[0]+"."+stkey.split('.')[1] - - # Define path to see if it exists - if args.phase in ['P', 'PP', 'allP']: - datapath = Path('P_DATA') / stfld - elif args.phase in ['S', 'SKS', 'allS']: - datapath = Path('S_DATA') / stfld - if not datapath.is_dir(): - print('Path to ' + str(datapath) + ' doesn`t exist - continuing') - continue - - # Temporary print locations - tlocs = sta.location - if len(tlocs) == 0: - tlocs = [''] - for il in range(0, len(tlocs)): - if len(tlocs[il]) == 0: - tlocs[il] = "--" - sta.location = tlocs - - # Update Display - print(" ") - print(" ") - print("|===============================================|") - print("| {0:>8s} |".format( - sta.station)) - print("|===============================================|") - print("| Station: {0:>2s}.{1:5s} |".format( - sta.network, sta.station)) - print("| Channel: {0:2s}; Locations: {1:15s} |".format( - sta.channel, ",".join(tlocs))) - print("| Lon: {0:7.2f}; Lat: {1:6.2f} |".format( - sta.longitude, sta.latitude)) - print("|-----------------------------------------------|") - - # Check for folder already processed - if stfld in procfold: - print(' {0} already processed...skipping '.format(stfld)) - continue - - datafiles = [x for x in datapath.iterdir() if x.is_dir()] - for folder in datafiles: - - # Skip hidden folders - if folder.name.startswith('.'): - continue - - # Re-initialize RFData object - rfdata = RFData(sta) - - # Load meta data - metafile = folder / "Meta_Data.pkl" - if not metafile.is_file(): - continue - rfdata.meta = pickle.load(open(metafile, 'rb')) - - # Skip data not in list of phases - if rfdata.meta.phase not in args.listphase: - continue - - if args.verb: - print("* Station: {0}; folder: {1}".format(stkey, folder)) - - if args.Z12: - try: - # Try loading Z12_Data - Z12file = folder / "Z12_Data.pkl" - rfdata.data = pickle.load(open(Z12file, 'rb')) - # Remove rotated flag and snr flag - rfdata.meta.rotated = False - rfdata.rotate(align='ZNE') - except: - - print("Z12_Data.pkl not available - using ZNE_Data.pkl") - # Load ZNE data - ZNEfile = folder / "ZNE_Data.pkl" - rfdata.data = pickle.load(open(ZNEfile, 'rb')) - else: - # Load ZNE data - ZNEfile = folder / "ZNE_Data.pkl" - rfdata.data = pickle.load(open(ZNEfile, 'rb')) - - # Resample if requested - if args.resample: - rfdata.data.resample( - args.resample, no_filter=False) - - # Remove rotated flag and snr flag - rfdata.meta.rotated = False - rfdata.meta.snr = None - - # Rotate from ZNE to 'align' ('ZRT', 'LQT', or 'PVH') - rfdata.rotate(vp=args.vp, vs=args.vs, align=args.align) - - # Calculate SNR - rfdata.calc_snr(dt=args.dt_snr, fmin=args.fmin, fmax=args.fmax) - - if args.verb: - print("* SNR: {}".format(rfdata.meta.snr)) - - # Deconvolve data - rfdata.deconvolve( - vp=args.vp, vs=args.vs, - align=args.align, method=args.method, - gfilt=args.gfilt, wlevel=args.wlevel, - pre_filt=args.pre_filt) - - # Get cross-correlation QC - rfdata.calc_cc() - - if args.verb: - print("* CC: {}".format(rfdata.meta.cc)) - - # Convert to Stream - rfstream = rfdata.to_stream() - - # Save RF Traces - RFfile = folder / "RF_Data.pkl" - pickle.dump(rfstream, open(RFfile, "wb")) - - # Save Meta data - metafile = folder / "Meta_Data.pkl" - pickle.dump(rfdata.meta, open(metafile, 'wb')) - - # Update - if args.verb: - print("* Output files written") - print("**************************************************") - - # Update processed folders - procfold.append(stfld) - - -if __name__ == "__main__": - - # Run main program - main() diff --git a/rfpy/scripts/rfpy_plot.py b/rfpy/scripts/rfpy_plot.py index 4715c64..4dcf040 100644 --- a/rfpy/scripts/rfpy_plot.py +++ b/rfpy/scripts/rfpy_plot.py @@ -149,10 +149,10 @@ def get_plot_arguments(argv=None): action="store", dest="nbaz", type=int, - default=None, + default=36, help="Specify integer number of back-azimuth bins to consider " + "(typically 36 or 72). If not None, the plot will show receiver " + - "functions sorted by back-azimuth values. [Default None]") + "functions sorted by back-azimuth values. [Default 36]") PreGroup.add_argument( "--nslow", action="store", @@ -301,9 +301,6 @@ def get_plot_arguments(argv=None): args.nbaz = 72 print("'nbaz' or 'nslow' not specified - plotting using " + "'nbaz=72'") - elif args.nbaz is not None and args.nslow is not None: - parser.error( - "Error: Cannot specify both 'nbaz' and 'nslow'") if args.trange is None: args.trange = [0., 30.] @@ -321,8 +318,6 @@ def get_plot_arguments(argv=None): if args.nbaz is None and args.nslow is None: parser.error("Specify at least one of --nbaz or --nslow") - elif args.nbaz is not None and args.nslow is not None: - parser.error("Specify only one of --nbaz or --nslow") return args @@ -507,18 +502,24 @@ def main(): print("Number of transverse RF data: " + str(len(rfTstream))) print('') - if args.nbaz: + if args.nbaz and args.nslow is None: # Bin according to BAZ rf_tmp = binning.bin(rfRstream, rfTstream, typ='baz', nbin=args.nbaz+1, pws=args.pws) - elif args.nslow: + elif args.nslow and args.nbaz is None: # Bin according to slowness rf_tmp = binning.bin(rfRstream, rfTstream, typ='slow', nbin=args.nslow+1, pws=args.pws) + elif args.nbaz is not None and args.nslow is not None: + # Bin into baz and slowness bins + rf_tmp = binning.bin_baz_slow(rfRstream, rfTstream, + nbaz=args.nbaz+1, nslow=args.nslow+1, + pws=False) + # Check bin counts: for tr in rf_tmp[0]: if (tr.stats.nbin < args.binlim): @@ -559,19 +560,24 @@ def main(): tr2 = None # Now plot - if args.nbaz: + if args.nbaz and args.nslow is None: plotting.wiggle_bins(rf_tmp[0], rf_tmp[1], tr1=tr1, tr2=tr2, btyp='baz', scale=args.scale, tmin=args.trange[0], tmax=args.trange[1], norm=norm, save=args.saveplot, title=args.titleplot, form=args.form) - elif args.nslow: + elif args.nslow and args.nbaz is None: plotting.wiggle_bins(rf_tmp[0], rf_tmp[1], tr1=tr1, tr2=tr2, btyp='slow', scale=args.scale, tmin=args.trange[0], tmax=args.trange[1], norm=norm, save=args.saveplot, title=args.titleplot, form=args.form) + elif args.nbaz is not None and args.nslow is not None: + plotting.panel(rf_tmp[0], rf_tmp[1], tmin=args.trange[0], tmax=args.trange[1], + normalize=True, save=args.saveplot, + title=args.titleplot, form=args.form) + # Update processed folders procfold.append(stfld) From 5cbe0d8c40eb89e63872b9beab41f384a6e834b7 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Thu, 18 Nov 2021 15:04:16 -0500 Subject: [PATCH 21/54] third large commit --- rfpy/plotting.py | 26 +++++++++----------------- rfpy/rfdata.py | 2 +- rfpy/scripts/rfpy_calc.py | 14 +++++++++----- setup.py | 3 +-- 4 files changed, 20 insertions(+), 25 deletions(-) diff --git a/rfpy/plotting.py b/rfpy/plotting.py index 3e9dafd..87123b1 100644 --- a/rfpy/plotting.py +++ b/rfpy/plotting.py @@ -447,16 +447,12 @@ def wiggle_single_event(rfdata, filt=None, pre_filt=None, trange=None): plt.show() -def panel(stream1, stream2, sort=None, tmin=0., tmax=30, normalize=True, +def panel(stream1, stream2, sort=None, tmin=0., tmax=30., normalize=True, save=False, title=None, form='png'): taxis = stream1[0].stats.taxis tselect = (taxis > tmin) & (taxis < tmax) nt = len(taxis[tselect]) - # print(nt) - # print(taxis[tselect]) - # print(len(stream1[0].data)) - # print(len(stream1[0].data[tselect])) baz_array = np.array([[i, tr.stats.baz] for i, tr in enumerate(stream1)]) slow_array = np.array([[i, tr.stats.slow] for i, tr in enumerate(stream1)]) @@ -470,32 +466,28 @@ def panel(stream1, stream2, sort=None, tmin=0., tmax=30, normalize=True, panelH[i,:] = tr.data[tselect] if normalize: - vmin = np.mean(panelV) - 3.*np.std(panelV) - vmax = np.mean(panelV) + 3.*np.std(panelV) + vmin = np.mean(panelV) - 4.*np.std(panelV) + vmax = np.mean(panelV) + 4.*np.std(panelV) else: vmin = -0.5 vmax = 0.5 + x, y = np.meshgrid(np.arange(len(stream1)), taxis[tselect]) + plt.figure(figsize=(8,5)) plt.subplot(311) - plt.pcolor(panelV.T, cmap='bwr', vmin=vmin, vmax=vmax) - plt.xlim((0,len(stream1))) - plt.ylim((0,nt)) + plt.pcolor(x, y, panelV.T, shading='auto', cmap='bwr', vmin=vmin, vmax=vmax) + plt.ylim(tmin, tmax) plt.gca().invert_yaxis() plt.ylabel('Time (s)') - plt.gca().set_yticks(np.arange(0., 600.05, 200.)) - plt.gca().set_yticklabels([str(0), str(10), str(20), str(30)]) plt.gca().set_xticklabels([]) plt.subplot(312) - plt.pcolor(panelH.T, cmap='bwr', vmin=vmin, vmax=vmax) - plt.xlim((0,len(stream1))) - plt.ylim((0,nt)) + plt.pcolor(x, y, panelH.T, shading='auto', cmap='bwr', vmin=vmin, vmax=vmax) + plt.ylim(tmin, tmax) plt.gca().invert_yaxis() plt.ylabel('Time (s)') - plt.gca().set_yticks(np.arange(0., 600.05, 200.)) - plt.gca().set_yticklabels([str(0), str(10), str(20), str(30)]) plt.gca().set_xticklabels([]) plt.subplot(313) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 9e37d1f..ba2e21f 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -793,7 +793,7 @@ def _Pwavelet(parent, method='complete', overhang=5., def _calc_specs(parent, daughter1, daughter2, noise, nn, method): - # Get length, zero padding parameters and frequencies + # Sampling interval dt = parent.stats.delta # Wiener or Water level deconvolution diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 74cb57a..3ed9ad6 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -373,12 +373,16 @@ def main(): if args.verb: print("* SNR: {}".format(rfdata.meta.snr)) + rfdata.calc_spectra( + vp=args.vp, vs=args.vs, + align=args.align, method=args.method, + wavelet='complete', envelope_threshold=0.05, + time=5) + # Deconvolve data rfdata.deconvolve( - vp=args.vp, vs=args.vs, - align=args.align, method=args.method, - gfilt=args.gfilt, wlevel=args.wlevel, - pre_filt=args.pre_filt) + align=args.align, method=args.method, + gfilt=args.gfilt, wlevel=args.wlevel) # Get cross-correlation QC rfdata.calc_cc() @@ -387,7 +391,7 @@ def main(): print("* CC: {}".format(rfdata.meta.cc)) # Convert to Stream - rfstream = rfdata.to_stream() + rfstream = rfdata.to_stream(store='rf') # Save RF Traces RFfile = folder / "RF_Data.pkl" diff --git a/setup.py b/setup.py index 89e94d5..40ad9c1 100644 --- a/setup.py +++ b/setup.py @@ -38,8 +38,7 @@ def find_version(*paths): entry_points={'console_scripts':[ 'rfpy_get_data=rfpy.scripts.rfpy_get_data:main', 'rfpy_calc_simdec=rfpy.scripts.rfpy_calc_simdec:main', - # 'rfpy_recalc=rfpy.scripts.rfpy_recalc:main', - 'rfpy_calc_sse=rfpy.scripts.rfpy_recalc_sse:main', + 'rfpy_calc=rfpy.scripts.rfpy_calc:main', 'rfpy_plot=rfpy.scripts.rfpy_plot:main', 'rfpy_harmonics=rfpy.scripts.rfpy_harmonics:main', 'rfpy_hk=rfpy.scripts.rfpy_hk:main', From a6a0d87f5cbf21c9a1def8f427b88bd88e42492f Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Thu, 18 Nov 2021 15:37:21 -0500 Subject: [PATCH 22/54] a bit of cleanup --- rfpy/rfdata.py | 8 +- rfpy/rfstreams.py | 560 ------------------------------- rfpy/scripts/rfpy_calc_simdec.py | 2 +- 3 files changed, 8 insertions(+), 562 deletions(-) delete mode 100644 rfpy/rfstreams.py diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index ba2e21f..087995b 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -646,7 +646,7 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.): def calc_spectra(self, phase='P', vp=None, vs=None, align=None, method='wiener', wavelet='complete', envelope_threshold=0.05, time=5., pre_filt=None, - writeto=False): + writeto=False, norm=False): """ Calculate the numerators and denominator of the spectral division for receiver function calculation using one component as the source wavelet. @@ -954,6 +954,12 @@ def _calc_specs(parent, daughter1, daughter2, noise, nn, method): corners=2, zerophase=True) for tr in [trL, trQ, trT, trNl, trNq]] + if norm: + # Standardize all trace data - useful in simdec + meanN = np.mean(trNl.data) + stdN = np.std(trNl.data) + [(tr.data - meanN)/stdN for tr in [trL, trQ, trT, trNl, trNq]] + if writeto: with open(writeto, 'wb') as f: pickle.dump(Stream(traces=[trL, trQ, trT]), f) diff --git a/rfpy/rfstreams.py b/rfpy/rfstreams.py deleted file mode 100644 index 22249ba..0000000 --- a/rfpy/rfstreams.py +++ /dev/null @@ -1,560 +0,0 @@ -def RFStreams(object): - - def _init__(self, rfdata=None): - - self.rfdata = [] - - if isinstance(rfdata, RFData): - rfdata = [rfdata] - if rfdata: - self.rfdata.extend(rfdata) - - def __add__(self, other): - """ - Add two `:class:`~plateflex.classes.Grid` objects or a - :class:`~plateflex.classes.Project` object with a single grid. - - """ - if isinstance(other, RFData): - other = RFStreams([other]) - if not isinstance(other, RFStreams): - raise TypeError - rfdata = self.rfdata + other.rfdata - return self.__class__(rfdata=rfdata) - - def __iter__(self): - """ - Return a robust iterator for :class:`~plateflex.classes.Grid` - objects - - """ - return list(self.rfdata).__iter__() - - def append(self, rfdata): - """ - Append a single :class:`~plateflex.classes.Grid` object to the - current `:class:`~plateflex.classes.Project` object. - - :type grid: :class:`~plateflex.classes.Grid` - :param grid: object to append to project - - .. rubric:: Example - - >>> import numpy as np - >>> from plateflex import Grid, Project - >>> nn = 200; dd = 10. - >>> grid = Grid(np.random.randn(nn, nn), dd, dd) - >>> project = Project() - >>> project.append(grid) - """ - - if isinstance(rfdata, RFData): - self.rfdata.append(rfdata) - else: - msg = 'Append only supports a single RFData object as an argument.' - raise TypeError(msg) - - return self - - def extend(self, rfdata_list): - """ - Extend the current Project object with a list of Grid objects. - - :param trace_list: list of :class:`~plateflex.classes.Grid` objects or - :class:`~plateflex.classes.Project`. - - .. rubric:: Example - - >>> import numpy as np - >>> from plateflex import Grid, Project - >>> nn = 200; dd = 10. - >>> grid1 = Grid(np.random.randn(nn, nn), dd, dd) - >>> grid2 = Grid(np.random.randn(nn, nn), dd, dd) - >>> project = Project() - >>> project.extend(grids=[grid1, grid2]) - - """ - if isinstance(rfdata_list, list): - for _i in rfdata_list: - # Make sure each item in the list is a Grid object. - if not isinstance(_i, RFData): - msg = 'Extend only accepts a list of RFData objects.' - raise TypeError(msg) - self.rfdata.extend(rfdata_list) - elif isinstance(rfdata_list, RFStreams): - self.rfdata.extend(rfdata_list.rfdata) - else: - msg = 'Extend only supports a list of RFData objects as argument.' - raise TypeError(msg) - return self - - - def calc_spectra(self, wavelet='complete',): - - def _npow2(x): - return 1 if x == 0 else 2**(x-1).bit_length() - - def _pad(array, n): - tmp = np.zeros(n) - tmp[:array.shape[0]] = array - return tmp - - def _Pwavelet(parent, method='complete', overhang=5, - envelope_threshold=0.05, time=5): - - """ - Select wavelet from the parent function for deconvolution using method. - parent: obspy.Trace - wavefrom to extract the wavelet from - method: str - 'complete' use complete parent signal after P arrival (current - implementation) - 'envelope' use only the part of the parent signal after the - P arrival where - envelope > envelope_threshold*max(envelope) - fall back to 'complete' if condition not reached - 'time' use only this many seonds after P arrival` - fall back to 'complete' if longer than parent - overhang: float - seconds before start and after end of wavelet to be used for - tapering - envelope_threshold: float - fraction of the envelope that defines wavelet (for - method='envelope') - time: float - window (seconds) that defines the wavelet (for method='time') - minimum time (seconds) of the wavelet (for method='envelope') - - Return: - left, right: (obspy.UTCDateTime) start and end time of wavelet - """ - - import obspy.signal.filter as osf - - if method not in ['complete', 'envelope', 'time', 'noise']: - msg = 'Unknow method for wavelet extraction: ' + method - raise NotImplementedError(msg) - - errflag = False - - if method == 'envelope': - split = int((self.meta.time + self.meta.ttime + - time - parent.stats.starttime ) * - parent.stats.sampling_rate) - env = osf.envelope(parent.data) - env /= max(env) # normalize - env = env[split:] # only look after P + time - try: - i = np.nonzero( - np.diff( - np.array( - env > envelope_threshold, dtype=int))==-1)[0][0] - except IndexError: - i = len(parent.data)-1 - dts = i * parent.stats.delta + time - left = self.meta.time + self.meta.ttime - overhang - right = self.meta.time + self.meta.ttime + dts + overhang - if left < parent.stats.starttime or right > parent.stats.endtime: - print('Envelope wavelet longer than trace.') - print('Falling back to "complete" wavelet') - errflag = True - - if method == 'time': - left = self.meta.time + self.meta.ttime - overhang - right = self.meta.time + self.meta.ttime + time + overhang - if left < parent.stats.starttime or right > parent.stats.endtime: - print('Time wavelet longer than trace.') - print('Falling back to "complete" wavelet') - errflag = True - - if method == 'complete' or errflag: - dts = len(parent.data)*parent.stats.delta/2. - left = self.meta.time + self.meta.ttime - overhang - right = self.meta.time + self.meta.ttime + dts - 2*overhang - - if method == 'noise': - dts = len(parent.data)*parent.stats.delta/2. - left = self.meta.time + self.meta.ttime - dts - right = self.meta.time + self.meta.ttime - overhang - - return left, right - - def _specs(trL, trQ, trT, trN, nn, method): - - # Get length, zero padding parameters and frequencies - dt = parent.stats.delta - - # npad = _npow2(nn*2) - npad = nn - freqs = np.fft.fftfreq(npad, d=dt) - - # Fourier transform - Fl = np.fft.fft(trL.data, n=npad) - Fq = np.fft.fft(trQ.data, n=npad) - Ft = np.fft.fft(trT.data, n=npad) - Fn = np.fft.fft(trN.data, n=npad) - - # Copy traces - Sll = trL.copy() - Sql = trQ.copy() - Stl = trT.copy() - Snn = trN.copy() - - # Auto and cross spectra - Sll.data = np.real(Fl*np.conjugate(Fl)) - Sql.data = Fq*np.conjugate(Fl) - Stl.data = Ft*np.conjugate(Fl) - Snn.data = np.real(Fn*np.conjugate(Fn)) - - return Sll, Sql, Stl, Snn - - - for rfdata in self.rfdata: - - # Get the name of components (order is critical here) - cL = stream.meta.align[0] - cQ = stream.meta.align[1] - cT = stream.meta.align[2] - - # Define signal and noise - trL = stream.data.select(component=cL)[0].copy() - trQ = stream.data.select(component=cQ)[0].copy() - trT = stream.data.select(component=cT)[0].copy() - trN = stream.data.select(component=cL)[0].copy() - - # Standardize all trace data - meanN = np.mean(trN.data) - stdN = np.std(trN.data) - [tr.data = (tr.data - meanN)/stdN for tr in [trL, trQ, trT, trN]] - - trL.stats.channel = 'WV' + stream.meta.align[0] - trQ.stats.channel = 'WV' + stream.meta.align[1] - trT.stats.channel = 'WV' + stream.meta.align[2] - trN.stats.channel = 'WV' + stream.meta.align[0] - - if phase == 'P' or 'PP': - - # Get signal length (i.e., seismogram to deconvolve) from trace length - over = 5 - dtsqt = len(trL.data)*trL.stats.delta/2. - - # Traces will be zero-paded to this length (samples) - nn = int(round((dtsqt+over)*trL.stats.sampling_rate)) + 1 - - sig_left, sig_right = _Pwavelet(trL, method=wavelet, - envelope_threshold=envelope_threshold, time=time, overhang=over) - - # Trim wavelet - trL.trim(sig_left, sig_right, nearest_sample=False, pad=True, - fill_value=0.) - - # Signal window (-5. to dtsqt-10 sec) - sig_left, sig_right = _Pwavelet(trQ, method='complete', overhang=over) - - # Trim signal traces - [tr.trim(sig_left, sig_right, nearest_sample=False, - pad=True, fill_value=0.) for tr in [trQ, trT]] - - # Noise window (-dtsqt to -5. sec) - noise_left, noise_right = _Pwavelet(trQ, method='noise', overhang=over) - - # Trim noise trace - trN.trim(noise_left, noise_right, nearest_sample=False, - pad=True, fill_value=0.) - - elif phase == 'S' or 'SKS': - - # Get signal length (i.e., seismogram to deconvolve) from trace length - dts = len(trL.data)*trL.stats.delta/2. - - # Trim signal traces (-5. to dts-10 sec) - trL.trim(self.meta.time+self.meta.ttime+25.-dts/2., - self.meta.time+self.meta.ttime+25.) - trQ.trim(self.meta.time+self.meta.ttime+25.-dts/2., - self.meta.time+self.meta.ttime+25.) - trT.trim(self.meta.time+self.meta.ttime+25.-dts/2., - self.meta.time+self.meta.ttime+25.) - - # Trim noise traces (-dts to -5 sec) - trN.trim(self.meta.time+self.meta.ttime-dts, - self.meta.time+self.meta.ttime-dts/2.) - - # Taper traces - only necessary processing after trimming - # TODO: What does this to the multitaper method - [tr.taper(max_percentage=0.05, max_length=2.) - for tr in [trL, trQ, trT, trN]] - - # Spectra - if phase == 'P' or 'PP': - sll, sql, stl, snn = _specs(trL, trQ, trT, trN, nn) - - rfdata.specs = Stream(traces=[sll, sql, stl, snn]) - - - # elif phase == 'S' or 'SKS': - # rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method) - - - def bin_baz_slow_specs(self, nbaz=72, nslow=1): - - """ - Function to stack auto- and cross-component spectra into back-azimuth - and slowness bins. - - Parameters - ---------- - nbaz : int - Number of back-azimuth samples in bins - nslow : int - Number of slowness samples in bins - include_empty : bool - Return empty bins as null traces (default omits them) - - Returns - ------- - stack : :class:`~obspy.core.Stream` - Stream containing one or two stacked traces, - depending on the number of input streams - - """ - - # Define back-azimuth and slowness bins - baz_bins = np.linspace(0, 360, nbaz) - slow_bins = np.linspace(0.04, 0.08, nslow) - - # Extract baz and slowness - baz = [stream[i].stats.baz for i in range(len(stream))] - slow = [stream[i].stats.slow for i in range(len(stream))] - - # Digitize baz and slowness - ibaz = np.digitize(baz, baz_bins) - islow = np.digitize(slow, slow_bins) - - final_stream = [] - - for stream in self.streams: - try: - # Define empty streams - binned_stream = Stream() - - # Loop through baz_bins - for i in range(nbaz): - for j in range(nslow): - - nbin = 0 - array = np.zeros(len(stream[0].specs.data)) - - # Loop all traces - for k, tr in enumerate(stream): - - # If index of baz_bins is equal to ibaz - if i == ibaz[k] and j == islow[k]: - - nbin += 1 - array += tr.data - - continue - - if nbin > 0 or include_empty: - - # Average and update stats - array /= nbin - weight = np.real(abs(weight/nbin)) - - trace = Trace(header=stream[0].stats) - trace.stats.baz = baz_bins[i] - trace.stats.slow = slow_bins[j] - trace.stats.nbin = nbin - - if not pws: - weight = np.ones(len(stream[0].data)) - trace.data = weight*array - binned_stream.append(trace) - - final_stream.append(binned_stream) - - except: - continue - - return final_stream - - - def deconvolve(self, phase='P', vp=None, vs=None, - align=None, method='wiener', - envelope_threshold=0.05, time=5, pre_filt=None, - gfilt=None, wlevel=0.01, writeto=None): - """ - Deconvolves three-component data using one component as the source wavelet. - The source component is always taken as the dominant compressional - component, which can be either 'Z', 'L', or 'P'. - - Parameters - ---------- - vp : float - P-wave velocity at surface (km/s) - vs : float - S-wave velocity at surface (km/s) - align : str - Alignment of coordinate system for rotation - ('ZRT', 'LQT', or 'PVH') - method : str - Method for deconvolution. Options are 'wiener', 'water' or - 'multitaper' - wavelet : str - Type of wavelet for deconvolution. Options are 'complete', 'time' or - 'envelope' - envelope_threshold : float - Threshold [0-1] used in ``wavelet='envelope'``. - time : float - Window length used in ``wavelet='time'``. - Minimum window length for ``wavelet='envelope'``. - pre_filt : list of 2 floats - Low and High frequency corners of bandpass filter applied - before deconvolution - gfilt : float - Center frequency of Gaussian filter (Hz). - wlevel : float - Water level used in ``method='water'``. - writeto : str or None - Write wavelets for deconvolution to file. - - Attributes - ---------- - rf : :class:`~obspy.core.Stream` - Stream containing the receiver function traces - - """ - - if not self.meta.accept: - return - - - def _gauss_filt(dt, nft, f0): - df = 1./(nft*dt) - nft21 = int(0.5*nft + 1) - f = df*np.arange(nft21) - w = 2.*np.pi*f - gauss = np.zeros(nft) - gauss[:nft21] = np.exp(-0.25*(w/f0)**2.)/dt - gauss[nft21:] = np.flip(gauss[1:nft21-1]) - return gauss - - - def _decon(parent, daughter1, daughter2, noise, nn, method): - - - # Final processing depends on method - if method == 'wiener': - Sdenom = Spp + Snn - elif method == 'water': - phi = np.amax(Spp)*wlevel - Sdenom = Spp - Sdenom[Sdenom < phi] = phi - - # Apply Gaussian filter? - if gfilt: - gauss = _gauss_filt(dt, npad, gfilt) - gnorm = np.sum(gauss)*(freqs[1]-freqs[0])*dt - else: - gauss = np.ones(npad) - gnorm = 1. - - # Is this correct? - #parent de-noised = np.fft.ifftshift(np.real(np.fft.ifft(gauss*Sdenom))/gnorm) - #daughter1 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd1p))/gnorm) - #daughter2 de-noised = np.fft.ifftshift(np.real(np.fft.ifft( gauss*Sd2p))/gnorm) - - # Copy traces - rfp = parent.copy() - rfd1 = daughter1.copy() - rfd2 = daughter2.copy() - - # Spectral division and inverse transform - rfp.data = np.fft.ifftshift(np.real(np.fft.ifft( - gauss*Spp/Sdenom))/gnorm) - rfd1.data = np.fft.ifftshift(np.real(np.fft.ifft( - gauss*Sd1p/Sdenom))/gnorm) - rfd2.data = np.fft.ifftshift(np.real(np.fft.ifft( - gauss*Sd2p/Sdenom))/gnorm) - - return rfp, rfd1, rfd2 - - if not self.meta.rotated: - print("Warning: Data have not been rotated yet - rotating now") - self.rotate(vp=vp, vs=vs, align=align) - - # v--True if None v--True if nan, error if None - if not self.meta.snr or not np.isfinite(self.meta.snr): - print("Warning: SNR has not been calculated - " + - "calculating now using default") - self.calc_snr() - - if hasattr(self, 'rf'): - print("Warning: Data have been deconvolved already - passing") - return - - - # Pre-filter waveforms before deconvolution - if pre_filt: - [tr.filter('bandpass', freqmin=pre_filt[0], freqmax=pre_filt[1], - corners=2, zerophase=True) - for tr in [trL, trQ, trT, trNl, trNq]] - - if writeto: - with open(writeto, 'wb') as f: - pickle.dump(Stream(traces=[trL, trQ, trT]), f) - - # Deconvolve - if phase == 'P' or 'PP': - rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method) - - elif phase == 'S' or 'SKS': - rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method) - - # Update stats of streams - rfL.stats.channel = 'RF' + self.meta.align[0] - rfQ.stats.channel = 'RF' + self.meta.align[1] - rfT.stats.channel = 'RF' + self.meta.align[2] - - self.rf = Stream(traces=[rfL, rfQ, rfT]) - - - def to_stream(self): - """ - Method to switch from RFData object to Stream object. - This allows easier manipulation of the receiver functions - for post-processing. - - """ - - if not self.meta.accept: - return - - def _add_specstats(trace): - trace.stats.snr = self.meta.snr - trace.stats.snrh = self.meta.snrh - trace.stats.cc = self.meta.cc - trace.stats.slow = self.meta.slow - trace.stats.baz = self.meta.baz - trace.stats.gac = self.meta.gac - trace.stats.stlo = self.sta.longitude - trace.stats.stla = self.sta.latitude - trace.stats.evlo = self.meta.lon - trace.stats.evla = self.meta.lat - trace.stats.vp = self.meta.vp - trace.stats.vs = self.meta.vs - trace.stats.phase = self.meta.phase - trace.stats.is_rf = True - nn = self.rf[0].stats.npts - sr = self.rf[0].stats.sampling_rate - trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn) - return trace - - if not hasattr(self, 'rf'): - raise(Exception("Warning: Receiver functions are not available")) - - stream = self.rf - for tr in stream: - tr = _add_rfstats(tr) - - return stream diff --git a/rfpy/scripts/rfpy_calc_simdec.py b/rfpy/scripts/rfpy_calc_simdec.py index eb4c69f..aec351f 100644 --- a/rfpy/scripts/rfpy_calc_simdec.py +++ b/rfpy/scripts/rfpy_calc_simdec.py @@ -491,7 +491,7 @@ def main(): vp=args.vp, vs=args.vs, align=args.align, method=args.method, wavelet='complete', envelope_threshold=0.05, - time=5) + time=5, norm=True) # Convert to Stream rfstream = rfdata.to_stream(store='specs') From 015b265f8bb7d771c5371d5944dd12511cdb853d Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 22 Nov 2021 16:09:41 -0500 Subject: [PATCH 23/54] fixed type of data trace during binning --- rfpy/binning.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 0aa6530..4e8dc79 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -209,7 +209,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, for j in range(nslow): nbin = 0 - array = np.zeros(len(stream[0].data), dtype=type(stream[0].data)) + array = np.zeros(len(stream[0].data), dtype=type(stream[0].data[0])) weight = np.zeros(len(stream[0].data), dtype=complex) # Loop all traces From bd0a09db1aac663553e056042dea18d67e0af56d Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 22 Nov 2021 16:10:24 -0500 Subject: [PATCH 24/54] terminal strings --- rfpy/scripts/rfpy_calc.py | 20 ++++++++++---------- rfpy/scripts/rfpy_calc_simdec.py | 20 ++++++++++---------- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 3ed9ad6..6171f44 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -243,16 +243,16 @@ def get_recalc_arguments(argv=None): def main(): print() - print("########################################################") - print("# #") - print("# __ _ #") - print("# _ __ / _|_ __ _ _ _ __ ___ ___ __ _| | ___ #") - print("# | '__| |_| '_ \| | | | | '__/ _ \/ __/ _` | |/ __| #") - print("# | | | _| |_) | |_| | | | | __/ (_| (_| | | (__ #") - print("# |_| |_| | .__/ \__, |___|_| \___|\___\__,_|_|\___| #") - print("# |_| |___/_____| #") - print("# #") - print("########################################################") + print("##############################################") + print("# #") + print("# __ _ #") + print("# _ __ / _|_ __ _ _ ___ __ _| | ___ #") + print("# | '__| |_| '_ \| | | | / __/ _` | |/ __| #") + print("# | | | _| |_) | |_| | | (_| (_| | | (__ #") + print("# |_| |_| | .__/ \__, |___\___\__,_|_|\___| #") + print("# |_| |___/_____| #") + print("# #") + print("##############################################") print() # Run Input Parser diff --git a/rfpy/scripts/rfpy_calc_simdec.py b/rfpy/scripts/rfpy_calc_simdec.py index aec351f..6667454 100644 --- a/rfpy/scripts/rfpy_calc_simdec.py +++ b/rfpy/scripts/rfpy_calc_simdec.py @@ -348,16 +348,16 @@ def get_simdec_arguments(argv=None): def main(): print() - print("########################################################") - print("# #") - print("# __ _ #") - print("# _ __ / _|_ __ _ _ _ __ ___ ___ __ _| | ___ #") - print("# | '__| |_| '_ \| | | | | '__/ _ \/ __/ _` | |/ __| #") - print("# | | | _| |_) | |_| | | | | __/ (_| (_| | | (__ #") - print("# |_| |_| | .__/ \__, |___|_| \___|\___\__,_|_|\___| #") - print("# |_| |___/_____| #") - print("# #") - print("########################################################") + print("##################################################################################") + print("# #") + print("# __ _ _ _ #") + print("# _ __ / _|_ __ _ _ ___ __ _| | ___ ___(_)_ __ ___ __| | ___ ___ #") + print("# | '__| |_| '_ \| | | | / __/ _` | |/ __| / __| | '_ ` _ \ / _` |/ _ \/ __| #") + print("# | | | _| |_) | |_| | | (_| (_| | | (__ \__ \ | | | | | | (_| | __/ (__ #") + print("# |_| |_| | .__/ \__, |___\___\__,_|_|\___|___|___/_|_| |_| |_|\__,_|\___|\___| #") + print("# |_| |___/_____| |_____| #") + print("# #") + print("##################################################################################") print() # Run Input Parser From 85e497ad9ca23a010480dec9b37a22dcc5161bd2 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Fri, 10 Dec 2021 15:03:11 -0800 Subject: [PATCH 25/54] RFdata: handle missing align argument and to_stream store='data' --- rfpy/rfdata.py | 33 +++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 087995b..e4884b4 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -693,6 +693,14 @@ def calc_spectra(self, phase='P', vp=None, vs=None, if not self.meta.accept: return + if not align: + try: + align = self.meta.align + except AttributeError: + msg = "No alignment found. " + msg += "Please supply via align keyword." + raise ValueError(msg) + def _npow2(x): return 1 if x == 0 else 2**(x-1).bit_length() @@ -992,7 +1000,7 @@ def _calc_specs(parent, daughter1, daughter2, noise, nn, method): def deconvolve(self, align=None, method='wiener', - gfilt=None, wlevel=0.01, writeto=None): + gfilt=None, wlevel=0.01): """ Deconvolves three-component data using one component as the source wavelet. The source component is always taken as the dominant compressional @@ -1010,8 +1018,6 @@ def deconvolve(self, align=None, method='wiener', Center frequency of Gaussian filter (Hz). wlevel : float Water level used in ``method='water'``. - writeto : str or None - Write wavelets for deconvolution to file. Attributes ---------- @@ -1020,10 +1026,18 @@ def deconvolve(self, align=None, method='wiener', """ - try: - len(self.specs) > 0 - except: - Exception("spectra have not been calculated") + if not hasattr(self, 'specs'): + msg = "Spectra have not been calculated." + msg += "Call RFData.calc_spectra() first." + raise Exception(msg) + + if not align: + try: + align = self.meta.align + except AttributeError: + msg = "No alignment found. " + msg += "Please supply via align keyword." + raise ValueError(msg) # Make everything explicit @@ -1147,11 +1161,14 @@ def _add_stats(trace, store): trace.stats.is_specs = True nn = self.specs[0].stats.npts sr = self.specs[0].stats.sampling_rate + else: + nn = self.data[0].stats.npts + sr = self.data[0].stats.sampling_rate trace.stats.taxis = np.fft.fftshift(np.fft.fftfreq(nn, sr)*nn) return trace - if store is None: + if store is None or store is 'data': stream = self.data for tr in stream: tr = _add_stats(tr, store) From 3e3ccbeaf7e8cbaa2bbcccac8a81de12557142fd Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 6 Jan 2022 11:33:28 -0800 Subject: [PATCH 26/54] rfdata: Allow RFData to be initialized without station --- rfpy/rfdata.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index e4884b4..35c3bf4 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -164,6 +164,8 @@ class RFData(object): ---------- sta : object Object containing station information - from :mod:`~stdb` database. + When RFData is initialized without sta, functions requiring station info + coordinate will not be working. meta : :class:`~rfpy.rfdata.Meta` Object of metadata information for single event (initially set to None) data : :class:`~obspy.core.Stream` @@ -172,7 +174,7 @@ class RFData(object): """ - def __init__(self, sta): + def __init__(self, sta=None): # Load example data if initializing empty object if sta == 'demo' or sta == 'Demo': @@ -184,7 +186,10 @@ def __init__(self, sta): "examples/data", "MMPY.pkl"), 'rb'))['NY.MMPY'] # Attributes from parameters - self.sta = sta + if sta: + self.sta = sta + else: + self.sta = None # Initialize meta and data objects as None self.meta = None From 65b175e047237430c367b77803e3768ef675d4b5 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 24 Feb 2022 15:18:47 -0800 Subject: [PATCH 27/54] binning: only compute hilbert transform when pws=True --- rfpy/binning.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 4e8dc79..3700232 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -110,9 +110,10 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, nb += 1 array += tr.data - hilb = hilbert(tr.data) - phase = np.arctan2(hilb.imag, hilb.real) - weight += np.exp(1j*phase) + if pwd: + hilb = hilbert(tr.data) + phase = np.arctan2(hilb.imag, hilb.real) + weight += np.exp(1j*phase) continue @@ -220,9 +221,10 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, nbin += 1 array += tr.data - hilb = hilbert(np.real(tr.data)) - phase = np.arctan2(hilb.imag, hilb.real) - weight += np.exp(1j*phase) + if pws: + hilb = hilbert(np.real(tr.data)) + phase = np.arctan2(hilb.imag, hilb.real) + weight += np.exp(1j*phase) continue @@ -289,9 +291,10 @@ def bin_all(stream1, stream2=None, pws=False): # Get phase weights for tr in stream: array += tr.data - hilb = hilbert(tr.data) - phase = np.arctan2(hilb.imag, hilb.real) - pweight += np.exp(1j*phase) + if pws: + hilb = hilbert(tr.data) + phase = np.arctan2(hilb.imag, hilb.real) + pweight += np.exp(1j*phase) # Normalize array = array/len(stream) From 248d569051a513091b8eaf81b3ae4578c865124d Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 24 Feb 2022 17:15:11 -0800 Subject: [PATCH 28/54] rfdata: Fixed broken computation of waterlevel --- rfpy/rfdata.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 35c3bf4..4c43b66 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -1062,7 +1062,7 @@ def deconvolve(self, align=None, method='wiener', Sdenom = SLL.data + SNN.data elif method == 'water': - phi = np.amax(SLL)*wlevel + phi = np.amax(SLL.data)*wlevel Sdenom = SLL.data Sdenom[Sdenom < phi] = phi From 6124fd69caf81ffdba63749b63d2602814be4e63 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 7 Mar 2022 13:58:54 -0800 Subject: [PATCH 29/54] binning: new stack attribute: index_list --- rfpy/binning.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 3700232..8da810d 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -101,6 +101,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, nb = 0 array = np.zeros(len(stream[0].data)) weight = np.zeros(len(stream[0].data), dtype=complex) + index_list = [] # Loop through stat for j, tr in enumerate(stream): @@ -110,6 +111,8 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, nb += 1 array += tr.data + index_list.append(j) + if pwd: hilb = hilbert(tr.data) phase = np.arctan2(hilb.imag, hilb.real) @@ -125,19 +128,18 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, trace = Trace(header=stream[0].stats) trace.stats.nbin = nb + trace.stats.index_list = index_list + if typ == 'baz': trace.stats.baz = bins[i] trace.stats.slow = None - trace.stats.nbin = nb elif typ == 'slow': trace.stats.slow = bins[i] trace.stats.baz = None - trace.stats.nbin = nb elif typ == 'dist': trace.stats.dist = bins[i] trace.stats.slow = None trace.stats.baz = None - trace.stats.nbin = nb if not pws: weight = np.ones(len(stream[0].data)) trace.data = weight*array @@ -180,6 +182,12 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, Stream containing one or two stacked traces, depending on the number of input streams + Note + ---- + + Sets the following attributes of the stack: + nbin: Number of bins + index_list: Indices constituent traces in source stream """ # Define back-azimuth and slowness bins @@ -210,6 +218,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, for j in range(nslow): nbin = 0 + index_list = [] array = np.zeros(len(stream[0].data), dtype=type(stream[0].data[0])) weight = np.zeros(len(stream[0].data), dtype=complex) @@ -221,6 +230,8 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, nbin += 1 array += tr.data + index_list.append(k) + if pws: hilb = hilbert(np.real(tr.data)) phase = np.arctan2(hilb.imag, hilb.real) @@ -238,6 +249,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, trace.stats.baz = baz_bins[i] trace.stats.slow = slow_bins[j] trace.stats.nbin = nbin + trace.stats.index_list = index_list if not pws: weight = np.ones(len(stream[0].data)) From c1e06377a0a30d50b7903b8e01e228093e59451e Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Mon, 7 Mar 2022 14:06:01 -0800 Subject: [PATCH 30/54] rfdata: deconvolution method 'water2' --- rfpy/rfdata.py | 106 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 101 insertions(+), 5 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 4c43b66..a34bf7f 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -668,7 +668,7 @@ def calc_spectra(self, phase='P', vp=None, vs=None, Alignment of coordinate system for rotation ('ZRT', 'LQT', or 'PVH') method : str - Method for deconvolution. Options are 'wiener', 'water' or + Method for deconvolution. Options are 'wiener', 'water', 'water2' or 'multitaper' wavelet : str Type of wavelet for deconvolution. Options are 'complete', 'time' or @@ -684,7 +684,7 @@ def calc_spectra(self, phase='P', vp=None, vs=None, gfilt : float Center frequency of Gaussian filter (Hz). wlevel : float - Water level used in ``method='water'``. + Water level used in ``method='water'`` and ``method='water2'`` writeto : str or None Write wavelets for deconvolution to file. @@ -810,7 +810,7 @@ def _calc_specs(parent, daughter1, daughter2, noise, nn, method): dt = parent.stats.delta # Wiener or Water level deconvolution - if method == 'wiener' or method == 'water': + if method == 'wiener' or method == 'water' or method == 'water2': # npad = _npow2(nn*2) npad = nn @@ -1031,6 +1031,92 @@ def deconvolve(self, align=None, method='wiener', """ + def _gcv_beta(trL, trQ, betas, npad): + try: + # TODO: Conserve individual spectra in the binning process + nt = self.specs[0].stats.nbin + except AttributeError: + # This condition is currently always met + nt = 1 + + # Fourier transform + wft = np.fft.fft(trL, n=npad) + vft = np.fft.fft(trQ, n=npad) + + # Auto and cross spectra + wwft = wft*np.conjugate(wft) + vwft = vft*np.conjugate(wft) + + # Number of frequencies + nf = len(wft) + + misfits = np.zeros_like(betas) + modnorm = np.zeros_like(betas) + gcvf = np.zeros_like(betas) + for ib, beta in enumerate(betas): + # Define operator W W* / (W W* + B) and deconvolve to get + # impulse response in frequency domain. + wwft2 = wwft + beta + rft = vwft / wwft2 + xft = wwft / wwft2 + + # Compute model norm. + modnorm[ib] = np.linalg.norm(rft)**2 + + # Compute data misfit. Note misfit is numerator of GCV function + # Note also earlier mistake where norm(nft)^2 was norm(nft). + if nt == 1: + nft = vft - wft*rft # this should be SP + misfits[ib] = np.linalg.norm(nft)**2 + else: + raise NotImplementedError + # TODO: criterion never met + # TODO: iterate over individual spectra + for it in range(nt): + nft = Sdp[it] - Spp[it]*rft + misfits[ib] = misfits[ib] + np.linalg.norm(nft)**2 + + # Compute denominator and GCV function. + den = (nf*nt - np.real(np.sum(xft)))**2 + gcvf[ib] = misfits[ib]/den + + # Compute best beta. + ibest = np.argmin(gcvf) + + # If minimum not found inform user. + if ibest == 0 or ibest == len(betas): + print('WARNING: No minimum found for GCV') + print('change search limits') + print('index at minimum and no of spectra') + print(ibest, nt) + print('Using median value in range:') + ibest = len(betas)//2 + print(betas[ibest]) + + + if False: + import matplotlib.pyplot as mp + + fig, axs = mp.subplots(ncols=2) + ax = axs[0] + ax.plot(modnorm, misfits, color='black') + ax.plot(modnorm, misfits, marker='+', color='red') + ax.plot(modnorm[ibest], misfits[ibest], marker='o', color='green') + ax.set_xlabel('Model Norm') + ax.set_ylabel('Data Misfit') + ax = axs[1] + ax.plot(betas, gcvf) + ax.plot(betas, gcvf, marker='+', color='red') + ax.plot(betas[ibest], gcvf[ibest], color='green', marker='o') + ax.set_xscale('log') + ax.set_xlabel('Regularization Parameter') + ax.set_ylabel('GCV Function') + fig.show() + input('Press key to continue') + mp.close(fig) + + return betas[ibest] + if not hasattr(self, 'specs'): msg = "Spectra have not been calculated." msg += "Call RFData.calc_spectra() first." @@ -1044,7 +1130,6 @@ def deconvolve(self, align=None, method='wiener', msg += "Please supply via align keyword." raise ValueError(msg) - # Make everything explicit SLL = self.specs[0].copy() SQL = self.specs[1].copy() @@ -1057,7 +1142,6 @@ def deconvolve(self, align=None, method='wiener', # Wiener or Water level deconvolution if method in ['wiener', 'multitaper']: - # Denominator (Spp + Snn) Sdenom = SLL.data + SNN.data @@ -1066,6 +1150,18 @@ def deconvolve(self, align=None, method='wiener', Sdenom = SLL.data Sdenom[Sdenom < phi] = phi + elif method == 'water2': + # TODO: try to determine waterlevel automatically with: + # TODO: May only make sense in overdetermined case + # (i.e. simultanous deconvolution) + # beta = _gcv_beta(self.data[0].data, self.data[1].data, wlevel, npad) + + # David Gubbins + # Time Series Analysis and Inverse Therory for Geophysicists + # "Wiener Filter", Eq. 10.21 + beta = np.amax(SLL.data)*wlevel + Sdenom = SLL.data + beta + # Apply Gaussian filter? if gfilt: gauss = _gauss_filt(dt, npad, gfilt) From 7858460026569dfcb9463504dce42b9557267395 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Tue, 8 Mar 2022 10:56:07 -0800 Subject: [PATCH 31/54] utils: Catch None in station list --- rfpy/utils.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rfpy/utils.py b/rfpy/utils.py index b77f709..8a18677 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -80,6 +80,8 @@ def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='SAC', altnet=[]) fpathmatch = [] # Loop over all local data directories for lcldr in lcldrs: + if not lcldr: + continue # Recursiely walk through directory for root, dirnames, filenames in walk(lcldr): # Keep paths only for those matching the station From d951ce787ead4481ecf70c4b772d4b80bddcf89f Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Tue, 8 Mar 2022 10:58:34 -0800 Subject: [PATCH 32/54] binning: doc cleanup, added function baz_slow_bin_indices() --- rfpy/binning.py | 63 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index 8da810d..c04c27b 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -64,6 +64,11 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, Stream containing one or two stacked traces, depending on the number of input streams + Note + ---- + Sets the following attributes of the stack: + nbin: Number of bins + index_list: Indices of constituent traces in source stream """ if not typ in ['baz', 'slow', 'dist']: @@ -184,10 +189,9 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, Note ---- - Sets the following attributes of the stack: nbin: Number of bins - index_list: Indices constituent traces in source stream + index_list: Indices of constituent traces in source stream """ # Define back-azimuth and slowness bins @@ -265,8 +269,61 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, return final_stream -def bin_all(stream1, stream2=None, pws=False): +def baz_slow_bin_indices(stream, nbaz=36+1, nslow=20+1, include_empty=False): """ + Sort traces of streams into backazimuth (nbaz) and slowness (nslow) bins. + + Parameters + ---------- + stream : :class:`~obspy.core.Stream` + Stream of seismograms to be sorted in bins + bazs : int + Number of back-azimuth samples in bins + slows : int + Number of slowness samples in bins + include_empty : bool + Include empty bins (default omits them) + + Returns + ------- + index_lists : list of lists + Indices that sorts traces of stream into bins defined by nbaz and nslow + bazslow_list : list of 2*tuple + Backazimuth and slowness values of each element of index list + """ + + # Define back-azimuth and slowness bins + baz_bins = np.linspace(0, 360, nbaz) + slow_bins = np.linspace(0.04, 0.08, nslow) + + # Extract baz and slowness + baz = [s.stats.baz for s in stream] + slow = [s.stats.slow for s in stream] + + # Digitize baz and slowness + ibaz = np.digitize(baz, baz_bins) + islow = np.digitize(slow, slow_bins) + + index_lists = [] + bazslow_list = [] + + for i in range(nbaz): + for j in range(nslow): + index_list = [] + + for k, tr in enumerate(stream): + if i == ibaz[k] and j == islow[k]: + index_list.append(k) + + if len(index_list) > 0 or include_empty: + index_lists.append(index_list) + bazslow_list.append((baz_bins[i], slow_bins[j])) + + return index_lists, bazslow_list + + +def bin_all(stream1, stream2=None, pws=False): + """ Function to bin all streams into a single trace. This can be done using a linear stack (i.e., simple mean), or using phase-weighted stacking. From e3322e99010753817995806b567e51906728fb11 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 10 Mar 2022 13:35:12 -0800 Subject: [PATCH 33/54] binning: * New method stack() * baz_slow_bin_indices() now roks with list of RFData --- rfpy/binning.py | 76 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 70 insertions(+), 6 deletions(-) diff --git a/rfpy/binning.py b/rfpy/binning.py index c04c27b..1f38554 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -269,14 +269,78 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, return final_stream -def baz_slow_bin_indices(stream, nbaz=36+1, nslow=20+1, include_empty=False): +def stack(rfdatas, which='rf', pws=False, + attributes={'bin': None, 'slow': None, 'index_list': None}): + """ + Function to stack receiver functions stored in a RFData list. This can + be done using a linear stack (i.e., simple mean), or using phase-weighted + stacking. + + Parameters + ---------- + rfdatas : list of :class:`rfpy.rfdata.RFdata` + List of RFData objects to be stacked + which : str + 'rf' stack receiver functions + 'specs' stack spectra + pws : bool + Whether or not to perform phase-weighted stacking (which=rf only) + attributes : dict + Attributes passed to the traces of stack + + Returns + ------- + stack : :class:`~obspy.core.Stream` + Stream containing stacked traces + """ + + nbin = len(rfdatas) + + template = rfdatas[0].__dict__[which] + typ = float + if which == 'specs': + typ = complex + + array = np.zeros((len(template), len(template[0].data)), dtype=typ) + weight = np.zeros((len(template), len(template[0].data)), dtype=complex) + + for rfdata in rfdatas: + stream = rfdata.__dict__[which] + + for n, tr in enumerate(stream): + array[n, :] += tr.data + + if pws: + hilb = hilbert(np.real(tr.data)) + phase = np.arctan2(hilb.imag, hilb.real) + weight[n, :] += np.exp(1j*phase) + + if pws: + weight = np.real(abs(weight/nbin)) + array *= weight + array /= nbin + + stacks = Stream() + for n in range(len(template)): + stack = Trace(header=template[n].stats) + stack.data = array[n, :] + stack.stats.nbin = nbin + for att in attributes: + val = attributes[att] + stack.stats.__dict__[att] = val + stacks.append(stack) + + return stacks + + +def baz_slow_bin_indices(rfs, nbaz=36+1, nslow=20+1, include_empty=False): """ Sort traces of streams into backazimuth (nbaz) and slowness (nslow) bins. Parameters ---------- - stream : :class:`~obspy.core.Stream` - Stream of seismograms to be sorted in bins + rfs : list of :class:`~rfpy.RFData` + List of receiver functions to be sorted in bins bazs : int Number of back-azimuth samples in bins slows : int @@ -297,8 +361,8 @@ def baz_slow_bin_indices(stream, nbaz=36+1, nslow=20+1, include_empty=False): slow_bins = np.linspace(0.04, 0.08, nslow) # Extract baz and slowness - baz = [s.stats.baz for s in stream] - slow = [s.stats.slow for s in stream] + baz = [rf.meta.baz for rf in rfs] + slow = [rf.meta.slow for rf in rfs] # Digitize baz and slowness ibaz = np.digitize(baz, baz_bins) @@ -311,7 +375,7 @@ def baz_slow_bin_indices(stream, nbaz=36+1, nslow=20+1, include_empty=False): for j in range(nslow): index_list = [] - for k, tr in enumerate(stream): + for k in range(len(rfs)): if i == ibaz[k] and j == islow[k]: index_list.append(k) From ebaac7e04260151b36f222dbc7dcdfced22854c4 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 10 Mar 2022 13:37:11 -0800 Subject: [PATCH 34/54] Automatic optimal_wlevel() now part of utils --- rfpy/utils.py | 105 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) diff --git a/rfpy/utils.py b/rfpy/utils.py index 8a18677..4e83ef4 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -611,3 +611,108 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, else: print("* Waveforms Retrieved...") return False, st + + +def optimal_wlevel(rfdatas, wlevels=10**(np.arange(-5, 0.25, 0.25)), iplot=False): + """ + Optimal water level for 'water2' deconvolution using Generalized Cross + Validation + + Accepts + ------- + + rfdatas : list of rfpy.RFData + rfdata elements hold rotated seismograms in rfdata.data attribute + wlevels : iterable + Range of water levels (fraction of maximum stacked amplitude) to try. + Default: 1e-5 to 10 + + Returns: + -------- + wlevel : float + Optimal water level (fraction of maximum stacked amplitude) + """ + + nt = len(rfdatas) + nf = len(rfdatas[0].specs[0].data) + + wft = np.zeros((nt, nf), dtype=complex) + vft = np.zeros((nt, nf), dtype=complex) + + for it, rfdata in enumerate(rfdatas): + + if not rfdata.meta.rotated: + msg = 'Element {:} of rfdatas must be rotated.'.format(it) + raise ValueError(msg) + + wft[it, :] = np.fft.fft(rfdata.data[0].data, n=nf) + vft[it, :] = np.fft.fft(rfdata.data[1].data, n=nf) + + # Auto and cross spectra + wwft = wft*np.conjugate(wft) + vwft = vft*np.conjugate(wft) + + if nt > 0: + wwft = np.sum(wwft) + vwft = np.sum(vwft) + + betas = wlevels*np.amax(wwft) + misfits = np.zeros_like(betas) + modnorm = np.zeros_like(betas) + gcvf = np.zeros_like(betas) + + for ib, beta in enumerate(betas): + # Define operator W W* / (W W* + B) and deconvolve to get + # impulse response in frequency domain. + wwft2 = wwft + beta + rft = vwft / wwft2 + xft = wwft / wwft2 + + # Compute model norm. + modnorm[ib] = np.linalg.norm(rft)**2 + + # Compute data misfit. Note misfit is numerator of GCV function + # Note also earlier mistake where norm(nft)^2 was norm(nft). + if nt == 1: + nft = vft - wft*rft # this should be SP + misfits[ib] = np.linalg.norm(nft)**2 + else: + for it in range(nt): + nft = vft[it, :] - wft[it, :]*rft + misfits[ib] = misfits[ib] + np.linalg.norm(nft)**2 + + # Compute denominator and GCV function. + den = (nf*nt - np.real(np.sum(xft)))**2 + gcvf[ib] = misfits[ib]/den + + # Compute best beta. + ibest = np.argmin(gcvf) + + if iplot: + import matplotlib.pyplot as mp + + fig, axs = mp.subplots(ncols=2) + ax = axs[0] + ax.plot(modnorm, misfits, color='black') + ax.plot(modnorm, misfits, marker='+', color='red') + ax.plot(modnorm[ibest], misfits[ibest], marker='o', color='green') + ax.set_xlabel('Model Norm') + ax.set_ylabel('Data Misfit') + ax = axs[1] + ax.plot(betas, gcvf) + ax.plot(betas, gcvf, marker='+', color='red') + ax.plot(betas[ibest], gcvf[ibest], color='green', marker='o') + ax.set_xscale('log') + ax.set_xlabel('Regularization Parameter') + ax.set_ylabel('GCV Function') + fig.show() + input('Press key to continue') + mp.close(fig) + + if ibest == 0 or ibest == len(betas): + # Minimum not found + msg = 'No minimum found. Try extending wlevels.' + raise ValueError(msg) + + wlevel = float(betas[ibest]/np.amax(wwft)) + return wlevel From 887481a1f62f272c46508f54bd269e0720550c89 Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 10 Mar 2022 13:43:59 -0800 Subject: [PATCH 35/54] Removed optimal_wlevel() method from RFData. Now part of utils --- rfpy/rfdata.py | 97 ++++---------------------------------------------- 1 file changed, 6 insertions(+), 91 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index a34bf7f..e85abcc 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -171,7 +171,12 @@ class RFData(object): data : :class:`~obspy.core.Stream` Stream object containing the three-component seismograms (either un-rotated or rotated by the method :func:`~rfpy.rfdata.rotate`) - + rf : :class:`~obspy.core.Stream` + Stream object containing the receiver functions + (created when executing :func:`~rfpy.rfdata.deconvlove`) + specs : :class:`~obspy.core.Stream` + Stream object containing the spectra of the data + (created when executing :func:`~rfpy.rfdata.calc_spectra`) """ def __init__(self, sta=None): @@ -1031,91 +1036,6 @@ def deconvolve(self, align=None, method='wiener', """ - def _gcv_beta(trL, trQ, betas, npad): - try: - # TODO: Conserve individual spectra in the binning process - nt = self.specs[0].stats.nbin - except AttributeError: - # This condition is currently always met - nt = 1 - - # Fourier transform - wft = np.fft.fft(trL, n=npad) - vft = np.fft.fft(trQ, n=npad) - - # Auto and cross spectra - wwft = wft*np.conjugate(wft) - vwft = vft*np.conjugate(wft) - - # Number of frequencies - nf = len(wft) - - misfits = np.zeros_like(betas) - modnorm = np.zeros_like(betas) - gcvf = np.zeros_like(betas) - for ib, beta in enumerate(betas): - # Define operator W W* / (W W* + B) and deconvolve to get - # impulse response in frequency domain. - wwft2 = wwft + beta - rft = vwft / wwft2 - xft = wwft / wwft2 - - # Compute model norm. - modnorm[ib] = np.linalg.norm(rft)**2 - - # Compute data misfit. Note misfit is numerator of GCV function - # Note also earlier mistake where norm(nft)^2 was norm(nft). - if nt == 1: - nft = vft - wft*rft # this should be SP - misfits[ib] = np.linalg.norm(nft)**2 - else: - raise NotImplementedError - # TODO: criterion never met - # TODO: iterate over individual spectra - for it in range(nt): - nft = Sdp[it] - Spp[it]*rft - misfits[ib] = misfits[ib] + np.linalg.norm(nft)**2 - - # Compute denominator and GCV function. - den = (nf*nt - np.real(np.sum(xft)))**2 - gcvf[ib] = misfits[ib]/den - - # Compute best beta. - ibest = np.argmin(gcvf) - - # If minimum not found inform user. - if ibest == 0 or ibest == len(betas): - print('WARNING: No minimum found for GCV') - print('change search limits') - print('index at minimum and no of spectra') - print(ibest, nt) - print('Using median value in range:') - ibest = len(betas)//2 - print(betas[ibest]) - - - if False: - import matplotlib.pyplot as mp - - fig, axs = mp.subplots(ncols=2) - ax = axs[0] - ax.plot(modnorm, misfits, color='black') - ax.plot(modnorm, misfits, marker='+', color='red') - ax.plot(modnorm[ibest], misfits[ibest], marker='o', color='green') - ax.set_xlabel('Model Norm') - ax.set_ylabel('Data Misfit') - ax = axs[1] - ax.plot(betas, gcvf) - ax.plot(betas, gcvf, marker='+', color='red') - ax.plot(betas[ibest], gcvf[ibest], color='green', marker='o') - ax.set_xscale('log') - ax.set_xlabel('Regularization Parameter') - ax.set_ylabel('GCV Function') - fig.show() - input('Press key to continue') - mp.close(fig) - - return betas[ibest] if not hasattr(self, 'specs'): msg = "Spectra have not been calculated." @@ -1151,11 +1071,6 @@ def _gcv_beta(trL, trQ, betas, npad): Sdenom[Sdenom < phi] = phi elif method == 'water2': - # TODO: try to determine waterlevel automatically with: - # TODO: May only make sense in overdetermined case - # (i.e. simultanous deconvolution) - # beta = _gcv_beta(self.data[0].data, self.data[1].data, wlevel, npad) - # David Gubbins # Time Series Analysis and Inverse Therory for Geophysicists # "Wiener Filter", Eq. 10.21 From c6d56c4497b21df9e6aed8cb116ae421759f7a1b Mon Sep 17 00:00:00 2001 From: Wasja Bloch Date: Thu, 10 Mar 2022 13:44:55 -0800 Subject: [PATCH 36/54] rfpy_calc_simdec: * Support for optimal waterlevel * Compacter and easier binning and deconvolution routines --- rfpy/scripts/rfpy_calc_simdec.py | 69 +++++++++++++++++++++----------- 1 file changed, 46 insertions(+), 23 deletions(-) diff --git a/rfpy/scripts/rfpy_calc_simdec.py b/rfpy/scripts/rfpy_calc_simdec.py index 6667454..0d7a494 100644 --- a/rfpy/scripts/rfpy_calc_simdec.py +++ b/rfpy/scripts/rfpy_calc_simdec.py @@ -30,7 +30,7 @@ import pickle import stdb from obspy import Stream -from rfpy import RFData, binning, plotting +from rfpy import RFData, binning, plotting, utils from pathlib import Path @@ -194,6 +194,11 @@ def get_simdec_arguments(argv=None): default=0.01, help="Specify the water level, used in the 'water' method. " + "[Default 0.01]") + DeconGroup.add_argument( + "--try-auto-wlevel", + action="store_true", + dest="auto_wlevel", + help="Try to determine waterlevel using Generalized Cross Validation.") PreGroup = parser.add_argument_group( title='Pre-processing Settings', @@ -418,10 +423,11 @@ def main(): # print(' {0} already processed...skipping '.format(stfld)) # continue - LLstreams = Stream() - QLstreams = Stream() - TLstreams = Stream() - NNstreams = Stream() + #LLstreams = Stream() + #QLstreams = Stream() + #TLstreams = Stream() + #NNstreams = Stream() + RFs = np.array([]) datafiles = [x for x in datapath.iterdir() if x.is_dir()] for folder in datafiles: @@ -493,24 +499,26 @@ def main(): wavelet='complete', envelope_threshold=0.05, time=5, norm=True) - # Convert to Stream - rfstream = rfdata.to_stream(store='specs') + ## Convert to Stream + #rfstream = rfdata.to_stream(store='specs') - # Add to Stream objects - LLstreams.append(rfstream[0]) - QLstreams.append(rfstream[1]) - TLstreams.append(rfstream[2]) - NNstreams.append(rfstream[3]) + ## Add to Stream objects + #LLstreams.append(rfstream[0]) + #QLstreams.append(rfstream[1]) + #TLstreams.append(rfstream[2]) + #NNstreams.append(rfstream[3]) + RFs = np.append(RFs, RF) # Bin into baz and slowness bins - LL_tmp, = binning.bin_baz_slow(LLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, - pws=False, include_empty=False) - QL_tmp, = binning.bin_baz_slow(QLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, - pws=False, include_empty=False) - TL_tmp, = binning.bin_baz_slow(TLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, - pws=False, include_empty=False) - NN_tmp, = binning.bin_baz_slow(NNstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, - pws=False, include_empty=False) + #LL_tmp, = binning.bin_baz_slow(LLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, + #pws=False, include_empty=False) + #QL_tmp, = binning.bin_baz_slow(QLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, + #pws=False, include_empty=False) + #TL_tmp, = binning.bin_baz_slow(TLstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, + #pws=False, include_empty=False) + #NN_tmp, = binning.bin_baz_slow(NNstreams, nbaz=args.nbaz+1, nslow=args.nslow+1, + #pws=False, include_empty=False) + iibin, baz_slows = binning.baz_slow_bin_indices(RFs, nbaz, nslow) # Initialize new empty streams RFL = Stream() @@ -518,18 +526,33 @@ def main(): RFT = Stream() # Cycle through binned spectra - for i in range(len(LL_tmp)): + #for i in range(len(LL_tmp)): + for n, ibin in enumerate(iibin): + + # Receiver functions in bin + RFset = RFs[ibin] # Initialize new RFData object with empty Meta header rfdata = RFData(sta) # Extract the specs from the binned streams - rfdata.specs = Stream(traces=[LL_tmp[i], QL_tmp[i], TL_tmp[i], NN_tmp[i]]) + #rfdata.specs = Stream(traces=[LL_tmp[i], QL_tmp[i], TL_tmp[i], NN_tmp[i]]) + rfdata.specs = binning.stack(RFset, which='specs') + + # Set waterlevel + wlevel = args.wlevel + if args.auto_wlevel: + try: + wlevel = utils.optimal_wlevel(RFset) + except ValueError: + print('Water level GCV function had no minimum.') + print('Falling back to user supplied value.') + wlevel = args.wlevel # Deconvolve rfdata.deconvolve( align=args.align, method=args.method, - gfilt=args.gfilt, wlevel=args.wlevel) + gfilt=args.gfilt, wlevel=wlevel) # Append RFs to streams RFL.append(rfdata.rf[0]) From da8239b43e777137d7b20e44e794fb8b9eb8243e Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 14 Mar 2022 15:18:33 -0400 Subject: [PATCH 37/54] add example data in package --- setup.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/setup.py b/setup.py index 40ad9c1..25d6287 100644 --- a/setup.py +++ b/setup.py @@ -35,11 +35,14 @@ def find_version(*paths): install_requires=['numpy', 'obspy', 'stdb>=0.2.0'], python_requires='>=3.6', packages=setuptools.find_packages(), - entry_points={'console_scripts':[ - 'rfpy_get_data=rfpy.scripts.rfpy_get_data:main', - 'rfpy_calc_simdec=rfpy.scripts.rfpy_calc_simdec:main', - 'rfpy_calc=rfpy.scripts.rfpy_calc:main', - 'rfpy_plot=rfpy.scripts.rfpy_plot:main', - 'rfpy_harmonics=rfpy.scripts.rfpy_harmonics:main', - 'rfpy_hk=rfpy.scripts.rfpy_hk:main', - 'rfpy_ccp=rfpy.scripts.rfpy_ccp:main']}) + include_package_data=True, + package_data={'': ['examples/data/*.pkl']}, + entry_points={ + 'console_scripts': + ['rfpy_get_data=rfpy.scripts.rfpy_get_data:main', + 'rfpy_calc_simdec=rfpy.scripts.rfpy_calc_simdec:main', + 'rfpy_calc=rfpy.scripts.rfpy_calc:main', + 'rfpy_plot=rfpy.scripts.rfpy_plot:main', + 'rfpy_harmonics=rfpy.scripts.rfpy_harmonics:main', + 'rfpy_hk=rfpy.scripts.rfpy_hk:main', + 'rfpy_ccp=rfpy.scripts.rfpy_ccp:main']}) From 61a3621e8f18bde2a900e2cc377d4b489c00f2c6 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 16 Mar 2022 14:42:07 -0400 Subject: [PATCH 38/54] added Wasja --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3182aaf..74cb084 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ RfPy is a software to calculate single event-station receiver functions from the Installation, Usage, API documentation and scripts are described at https://paudetseis.github.io/RfPy/. -Authors: [`Pascal Audet`](https://www.uogeophysics.com/authors/admin/) (Developer and Maintainer) and [`Jeremy Gosselin`](https://www.uogeophysics.com/authors/gosselin/) (Contributor) +Authors: [Pascal Audet](https://www.uogeophysics.com/authors/admin/) and [Wasja Bloch](https://www.eoas.ubc.ca/people/wasjabloch)