From 519e50b3df3b8c574eca02a4889a4aba284b3815 Mon Sep 17 00:00:00 2001 From: pasansherath Date: Wed, 13 Sep 2023 13:18:13 +1200 Subject: [PATCH 01/22] Added deconv method in Audet 2010 BSSA paper --- rfpy/rfdata.py | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 4caf39a..56b27b5 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -89,6 +89,9 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'): if self.dep is not None: if self.dep > 1000.: self.dep = self.dep/1000. + if self.dep < 0.: + # Some events have negative depths in km + self.dep = np.abs(self.dep) else: self.dep = 10. @@ -647,7 +650,6 @@ def deconvolve(self, phase='P', vp=None, vs=None, Stream containing the receiver function traces """ - if not self.meta.accept: return @@ -669,13 +671,12 @@ def _gauss_filt(dt, nft, f0): gauss[nft21:] = np.flip(gauss[1:nft21-1]) return gauss - def _decon(parent, daughter1, daughter2, noise, nn, method): - + def _decon(parent, daughter1, daughter2, noise_parent, noise_daughter1, nn, method): # Get length, zero padding parameters and frequencies dt = parent.stats.delta # Wiener or Water level deconvolution - if method == 'wiener' or method == 'water': + if method == 'wiener' or method == 'water' or method == "wiener_audet_bssa2010": # npad = _npow2(nn*2) npad = nn @@ -685,17 +686,23 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): Fp = np.fft.fft(parent.data, n=npad) Fd1 = np.fft.fft(daughter1.data, n=npad) Fd2 = np.fft.fft(daughter2.data, n=npad) - Fn = np.fft.fft(noise.data, n=npad) + Fpn = np.fft.fft(noise_parent.data, n=npad) + Fd1n = np.fft.fft(noise_daughter1.data, n=npad) # Auto and cross spectra Spp = np.real(Fp*np.conjugate(Fp)) Sd1p = Fd1*np.conjugate(Fp) Sd2p = Fd2*np.conjugate(Fp) - Snn = np.real(Fn*np.conjugate(Fn)) + Snpp = np.real(Fpn*np.conjugate(Fpn)) + Snd1d1 = np.real(Fd1n*np.conjugate(Fd1n)) + Snpd1 = np.abs(Fd1n*np.conjugate(Fpn)) # Final processing depends on method if method == 'wiener': - Sdenom = Spp + Snn + Sdenom = Spp + Snpp + # Wiener ad-hoc deconvolution in Audet 2010 - BSSA + elif method == "wiener_audet_bssa2010": + Sdenom = Spp + 0.25*Snpp + 0.25*Snd1d1 + 0.5*Snpd1 elif method == 'water': phi = np.amax(Spp)*wlevel Sdenom = Spp @@ -712,11 +719,11 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): [tr.stats.npts for tr in [parent, daughter1, daughter2, - noise]], npad): + noise_parent]], npad): parent.data = _pad(parent.data, npad) daughter1.data = _pad(daughter1.data, npad) daughter2.data = _pad(daughter2.data, npad) - noise.data = _pad(noise.data, npad) + noise_parent.data = _pad(noise_parent.data, npad) freqs = np.fft.fftfreq(npad, d=dt) @@ -732,7 +739,7 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): Fd2 = np.fft.fft(np.multiply(tapers.transpose(), daughter2.data)) Fn = np.fft.fft(np.multiply(tapers.transpose(), - noise.data)) + noise_parent.data)) # Auto and cross spectra Spp = np.sum(np.real(Fp*np.conjugate(Fp)), axis=0) @@ -848,10 +855,10 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): # Deconvolve if phase == 'P' or 'PP': - rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method) + rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, trNq, nn, method) elif phase == 'S' or 'SKS': - rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, nn, method) + rfQ, rfL, rfT = _decon(trQ, trL, trT, trNq, trNl, nn, method) # Update stats of streams rfL.stats.channel = 'RF' + self.meta.align[0] From 3b52c630c8fd7ff7c30b9b66a03464ac025ffefa Mon Sep 17 00:00:00 2001 From: pasansherath Date: Wed, 13 Sep 2023 13:20:37 +1200 Subject: [PATCH 02/22] Added Audet 2010 BSSA deconv method option to rfpy_calc.py and rfpy_recalc.py scripts --- rfpy/scripts/rfpy_calc.py | 7 ++++--- rfpy/scripts/rfpy_recalc.py | 7 ++++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 362588f..ab2ab3c 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -484,10 +484,11 @@ def get_calc_arguments(argv=None): print("SNR window > data window. Defaulting to data " + "window minus 10 sec.") - if args.method not in ['wiener', 'water', 'multitaper']: + if args.method not in ['wiener', 'water', 'multitaper', \ + 'wiener_audet_bssa2010']: parser.error( "Error: 'method' should be either 'wiener', 'water' or " + - "'multitaper'") + "'multitaper', or 'wiener_audet_bssa2010'") return args @@ -655,7 +656,7 @@ def main(): print("| {0:>2s}.{1:5s}: {2:6d}".format( sta.network, sta.station, len(stalcllist)) + " files |") - #print(stalcllist[0:10]) + print(stalcllist[0:10]) else: stalcllist = utils.list_local_data_stn( lcldrs=args.localdata, sta=sta.station, dtype=args.dtype) diff --git a/rfpy/scripts/rfpy_recalc.py b/rfpy/scripts/rfpy_recalc.py index 74cb57a..83e3d83 100644 --- a/rfpy/scripts/rfpy_recalc.py +++ b/rfpy/scripts/rfpy_recalc.py @@ -224,10 +224,11 @@ def get_recalc_arguments(argv=None): "Error: Incorrect alignment specifier. Should be " + "either 'ZRT', 'LQT', or 'PVH'.") - if args.method not in ['wiener', 'water', 'multitaper']: + if args.method not in ['wiener', 'water', 'multitaper',\ + 'wiener_audet_bssa2010']: parser.error( - "Error: 'method' should be either 'wiener', 'water' or " + - "'multitaper'") + "Error: 'method' should be either 'wiener', 'water', " + + "'multitaper', or 'wiener_audet_bssa2010'") if args.pre_filt is not None: args.pre_filt = [float(val) for val in args.pre_filt.split(',')] From 5016c220de687282f008587cb3b1025bc9e57073 Mon Sep 17 00:00:00 2001 From: pasansherath Date: Wed, 13 Sep 2023 13:23:05 +1200 Subject: [PATCH 03/22] Included Audet 2010 BSSA paper deconv method in the scripts.rst doc --- docs/scripts.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/scripts.rst b/docs/scripts.rst index b353306..7a3d227 100644 --- a/docs/scripts.rst +++ b/docs/scripts.rst @@ -153,8 +153,8 @@ Usage Parameters for deconvolution --method METHOD Specify the deconvolution method. Available methods - include 'wiener', 'water' and 'multitaper'. [Default - 'wiener'] + include 'wiener', 'water', 'multitaper', and + 'wiener_audet_bssa2010'. [Default 'wiener'] --gfilt GFILT Specify the Gaussian filter width in Hz. [Default None] --wlevel WLEVEL Specify the water level, used in the 'water' method. @@ -240,8 +240,8 @@ Usage Parameters for deconvolution --method METHOD Specify the deconvolution method. Available methods - include 'wiener', 'water' and 'multitaper'. [Default - 'wiener'] + include 'wiener', 'water', 'multitaper' and + 'wiener_audet_bssa2010'. [Default 'wiener'] --gfilt GFILT Specify the Gaussian filter width in Hz. [Default None] --wlevel WLEVEL Specify the water level, used in the 'water' method. [Default 0.01] @@ -262,7 +262,7 @@ Usage .. code-block:: - $ rfpy_recalc -h + $ rfpy_plot -h ################################################# # __ _ _ # From a86ef54c74e3cd5b67179031bfd28842fc79de3c Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 12 May 2025 20:25:07 +1200 Subject: [PATCH 04/22] initializing v0.2.1 --- meson.build | 2 +- pyproject.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/meson.build b/meson.build index 00981a9..18b2c98 100644 --- a/meson.build +++ b/meson.build @@ -1,5 +1,5 @@ project('rfpy', 'c', - version : '0.2.0', + version : '0.2.1', license: 'MIT', meson_version: '>=0.64.0', ) diff --git a/pyproject.toml b/pyproject.toml index dec229e..083738b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ requires = ["meson-python>0.15.0", "numpy >= 1.25.0"] [project] name = "rfpy" -version = "0.2.0" +version = "0.2.1" description = "Python package for calculating and post-processing receiver functions" authors = [ { name = "Pascal Audet", email = "pascal.audet@uottawa.ca" } From 45f382ab5d85dfce8a8d0fb569c999f5ebf4f83e Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 12 May 2025 20:25:32 +1200 Subject: [PATCH 05/22] __init__ to 0.2.1 --- rfpy/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rfpy/__init__.py b/rfpy/__init__.py index 93019d7..eaf7b06 100644 --- a/rfpy/__init__.py +++ b/rfpy/__init__.py @@ -20,7 +20,7 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. -__version__ = '0.2.0' +__version__ = '0.2.1' __author__ = 'Pascal Audet' From 9025dfc0c477db93959dba6fedc4f36fa2da6910 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 12 May 2025 20:26:28 +1200 Subject: [PATCH 06/22] Adding SDS client for local data --- rfpy/rfdata.py | 64 ++- rfpy/scripts/rfpy_calc.py | 135 +++--- rfpy/utils.py | 902 ++++++++++++++++++-------------------- 3 files changed, 523 insertions(+), 578 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index c3d1ed7..5fcbd3a 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -54,13 +54,13 @@ class Meta(object): az : float Azimuth - pointing to station from earthquake (degrees) ttime : float - Predicted arrival time (sec) + Predicted arrival time (sec) ph : str - Phase name + Phase name slow : float - Horizontal slowness of phase + Horizontal slowness of phase inc : float - Incidence angle of phase at surface + Incidence angle of phase at surface vp : float P-wave velocity at surface (km/s) vs : float @@ -72,8 +72,8 @@ class Meta(object): Whether or not data have been rotated to ``align`` coordinate system zcomp : str - Vertical Component Identifier. Should be a single character. - This is different then 'Z' only for fully unknown component + Vertical Component Identifier. Should be a single character. + This is different then 'Z' only for fully unknown component orientation (i.e., components are 1, 2, 3) """ @@ -155,8 +155,8 @@ def __init__(self, sta, event, gacmin=30., gacmax=90., phase='P'): 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 + station information with a single event (i.e., earthquake) + metadata, corresponding raw and rotated seismograms and receiver functions. Note @@ -198,8 +198,8 @@ def __init__(self, sta, zcomp='Z'): def add_event(self, event, gacmin=30., gacmax=90., phase='P', returned=False): """ - Adds event metadata to RFData object, including travel time info - of P wave. + Adds event metadata to RFData object, including travel time info + of P wave. Parameters ---------- @@ -281,7 +281,7 @@ def add_data(self, stream, returned=False, new_sr=5.): print(stream) if not isinstance(stream, Stream): - raise(Exception("Event has incorrect type")) + raise Exception("Event has incorrect type") try: self.data = stream @@ -296,17 +296,17 @@ def add_data(self, stream, returned=False, new_sr=5.): corners=2, zerophase=True) self.data.resample(new_sr, no_filter=True) - except: + except Exception as e: print("Error: Not all channels are available") self.meta.accept = False if returned: return self.meta.accept - 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): + def download_data(self, client, new_sr=5., dts=120., + remove_response=False, returned=False, + verbose=False): + """ Downloads seismograms based on event origin time and P phase arrival. @@ -315,14 +315,10 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, ---------- client : :class:`~obspy.client.fdsn.Client` Client object - ndval : float - Fill in value for missing data new_sr : float New sampling rate (Hz) dts : float 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() @@ -344,7 +340,7 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, """ if self.meta is None: - raise(Exception("Requires event data as attribute - aborting")) + raise Exception("Requires event data as attribute - aborting") if not self.meta.accept: return @@ -353,18 +349,16 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, tstart = self.meta.time + self.meta.ttime - dts 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")) + if verbose: + 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( client=client, sta=self.sta, start=tstart, end=tend, - stdata=stdata, dtype=dtype, ndval=ndval, new_sr=new_sr, - remove_response=remove_response, - local_response_dir=local_response_dir, - verbose=verbose, zcomp=self.zcomp) + new_sr=new_sr, verbose=verbose, remove_response=remove_response, + zcomp=self.zcomp) # Store as attributes with traces in dictionary try: @@ -380,7 +374,7 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, self.data.resample(new_sr, no_filter=True) # If there is no ZNE, perhaps there is Z12 (or zcomp12)? - except: + except Exception as e: try: tr1 = stream.select(component='1')[0] @@ -404,7 +398,7 @@ 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: self.meta.accept = False if returned: @@ -413,11 +407,11 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, def rotate(self, vp=None, vs=None, align=None): """ Rotates 3-component seismograms from vertical (Z), - east (E) and north (N) to longitudinal (L), + east (E) and north (N) to longitudinal (L), radial (Q) and tangential (T) components of motion. Note that the method 'rotate' from ``obspy.core.stream.Stream`` is used for the rotation ``'ZNE->ZRT'`` and ``'ZNE->LQT'``. - Rotation ``'ZNE->PVH'`` is implemented separately here + Rotation ``'ZNE->PVH'`` is implemented separately here due to different conventions. Can also rotate Z12 to ZNE. @@ -545,8 +539,7 @@ def rotate(self, vp=None, vs=None, align=None): self.meta.rotated = True else: - raise(Exception("incorrect 'align' argument")) - + raise Exception("incorrect 'align' argument") def calc_snr(self, dt=30., fmin=0.05, fmax=1.): """ @@ -634,7 +627,6 @@ 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, align=None, method='wiener', wavelet='complete', envelope_threshold=0.05, time=5, pre_filt=None, diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 09b3394..8e206f4 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -27,7 +27,8 @@ import numpy as np import pickle import stdb -from obspy.clients.fdsn import Client +from obspy.clients.fdsn import Client as FDSN_Client +from obspy.clients.filesystem.sds import Client as SDS_Client from obspy import Catalog, UTCDateTime from http.client import IncompleteRead from rfpy import utils, RFData @@ -155,22 +156,22 @@ def get_calc_arguments(argv=None): type=str, dest="localdata", default=None, - help="Specify a comma separated list of paths containing " + + help="Specify path 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") + "the data using the FDSN 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.") + "--dtype", + action="store", + type=str, + dest="dtype", + default='MSEED', + 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", @@ -419,26 +420,27 @@ def get_calc_arguments(argv=None): if args.userauth is not None: tt = args.userauth.split(':') if not len(tt) == 2: - msg = ("Error: Incorrect Username and Password Strings for User " - "Authentification") + msg = ( + "Error: Incorrect Username and Password Strings " + + "for User Authentification") parser.error(msg) else: args.userauth = tt else: args.userauth = [None, None] - # Parse Local Data directories - if args.localdata is not None: - args.localdata = args.localdata.split(',') - else: - args.localdata = [] + # # 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.") + if args.dtype.upper() not in ['MSEED', '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: @@ -557,14 +559,19 @@ def main(): datapath.mkdir(parents=True) # Establish client - data_client = Client( - base_url=args.server, - user=args.userauth[0], - password=args.userauth[1], - eida_token=args.tokenfile) + if args.localdata is None: + data_client = FDSN_Client( + base_url=args.server, + user=args.userauth[0], + password=args.userauth[1], + eida_token=args.tokenfile) + else: + data_client = SDS_Client( + args.localdata, + format=args.dtype) # Establish client for events - event_client = Client() + event_client = FDSN_Client() # Get catalogue search start time if args.startT is None: @@ -581,13 +588,10 @@ def main(): 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 + tlocs = [] + for il in range(0, len(sta.location)): + if len(sta.location[il]) == 0: + tlocs.append("--") # Update Display print(" ") @@ -630,15 +634,16 @@ def main(): cat = event_client.get_events( starttime=tstart, endtime=tend, minmagnitude=args.minmag, maxmagnitude=args.maxmag) + except IncompleteRead: # See http.client.IncompleteRead # https://github.com/obspy/obspy/issues/3028#issue-1179808237 - + # Read yearly chunks print('| Warning: Unable to get entire catalog at once |') print('| Trying to get one year at a time |') print('| |') - + chunk = 365 * 86400 # Query length in seconds cat = Catalog() tend = min(tend, UTCDateTime.now()) @@ -668,27 +673,27 @@ def main(): " 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("|===============================================|") + # # 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: @@ -761,9 +766,11 @@ def main(): # 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) + client=data_client, + dts=args.dts, + new_sr=args.new_sampling_rate, + returned=True, + verbose=args.verb) if not has_data: continue diff --git a/rfpy/utils.py b/rfpy/utils.py index 6a7d718..c500145 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -39,378 +39,374 @@ def traceshift(trace, tt): return rtrace -def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='SAC', altnet=[]): - """ - Function to take the list of local directories and recursively - find all data that matches the station name - - Parameters - ---------- - lcldrs : List - List of local directories - sta : Dict - Station metadata from :mod:`~StDb` - net : str - Network name - altnet : List - List of alternative networks - - Returns - ------- - fpathmatch : List - Sorted list of matched directories - - """ - from fnmatch import filter - from os import walk - from os.path import join - - if sta is None: - return [] - else: - if net is None: - sstrings = ['*.{0:s}.*.{1:s}'.format(sta, dtype)] - else: - sstrings = ['*.{0:s}.{1:s}.*.{2:s}'.format(net, sta, dtype)] - if len(altnet) > 0: - for anet in altnet: - sstrings.append( - '*.{0:s}.{1:s}.*.{2:s}'.format(anet, sta, dtype)) - - fpathmatch = [] - # Loop over all local data directories - for lcldr in lcldrs: - # Recursiely walk through directory - for root, dirnames, filenames in walk(lcldr): - # Keep paths only for those matching the station - for sstring in sstrings: - for filename in filter(filenames, sstring): - # No hidden directories - if not '/.' in root: - fpathmatch.append(join(root, filename)) - - fpathmatch.sort() - - return fpathmatch - - -def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, - start=UTCDateTime, end=UTCDateTime, ndval=nan): - """ - Function to determine the path to data for a given component and alternate network - - Parameters - ---------- - comp : str - Channel for seismogram (one letter only) - stdata : List - Station list - sta : Dict - Station metadata from :mod:`~StDb` data base - start : :class:`~obspy.core.utcdatetime.UTCDateTime` - Start time for request - end : :class:`~obspy.core.utcdatetime.UTCDateTime` - End time for request - ndval : float or nan - Default value for missing data - - Returns - ------- - err : bool - Boolean for error handling (`False` is associated with success) - st : :class:`~obspy.core.Stream` - Stream containing North, East and Vertical components of motion - - """ - - from fnmatch import filter - - # Get start and end parameters - styr = start.strftime("%Y") - stjd = start.strftime("%j") - edyr = end.strftime("%Y") - edjd = end.strftime("%j") - - # Intialize to default positive error - erd = True - - print( - ("* {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: - - 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: - print("* - Data Unavailable") - return erd, None - - # Process the local Files - for sacfile in lclfiles: - # Read File - st = read(sacfile) - # st = read(sacfile, format="SAC") - - if dtype.upper() == 'MSEED': - if len(st) > 1: - st.merge(method=1, interpolation_samples=- - 1, fill_value=-123456789) - - # Should only be one component, otherwise keep reading If more - # than 1 component, error - if len(st) != 1: - pass - - else: - # Check start/end times in range - if (st[0].stats.starttime <= start and - st[0].stats.endtime >= end): - st.trim(starttime=start, endtime=end) - - eddt = False - # Check for NoData and convert to NaN if a SAC file - if dtype.upper() == 'SAC': - 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 - - # Check for Nan in stream for SAC - if True in isnan(st[0].data): - print( - "* !!! Missing Data Present !!! " + - "Skipping (NaNs)") - # Check for ND Val in stream for MSEED - elif -123456789 in st[0].data: - print( - "* !!! Missing Data Present !!! " + - "Skipping (MSEED fill)") - else: - if eddt and (ndval == 0.0): - if any(st[0].data == 0.0): - print( - "* !!! Missing Data Present " + - "!!! (Set to Zero)") - - st[0].stats.update() - tloc = st[0].stats.location - if len(tloc) == 0: - tloc = "--" - - # Processed succesfully...Finish - print(("* {1:3s}.{2:2s} - From Disk".format( - st[0].stats.station, st[0].stats.channel.upper(), - tloc))) - return False, st - - # Time Window spans Multiple days - else: - 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: - print("* - Data Unavailable") - return erd, None - - # Now try to merge the two separate day files - if len(lclfiles1) > 0 and len(lclfiles2) > 0: - # Loop over first day file options - for sacf1 in lclfiles1: - st1 = read(sacf1) - if dtype.upper() == 'MSEED': - if len(st1) > 1: - st1.merge(method=1, interpolation_samples=-1, - fill_value=-123456789) - - # Loop over second day file options - for sacf2 in lclfiles2: - st2 = read(sacf2) - if dtype.upper() == 'MSEED': - if len(st2) > 1: - st2.merge( - method=1, interpolation_samples=-1, - fill_value=-123456789) - - # Check time overlap of the two files. - if st1[0].stats.endtime >= \ - st2[0].stats.starttime-st2[0].stats.delta: - # eddt1 = False - # eddt2 = False - # if dtype.upper() == 'SAC': - # # Check for NoData and convert to NaN - # st1nd = st1[0].stats.sac['user9'] - # st2nd = st2[0].stats.sac['user9'] - # if (not st1nd == 0.0) and (not st1nd == -12345.0): - # st1[0].data[st1[0].data == st1nd] = ndval - # eddt1 = True - # if (not st2nd == 0.0) and (not st2nd == -12345.0): - # st2[0].data[st2[0].data == st2nd] = ndval - # eddt2 = True - - st = st1 + st2 - # Need to work on this HERE (AJS OCT 2015). - # If Calibration factors are different, - # then the traces cannot be merged. - try: - st.merge(method=1, interpolation_samples=- - 1, fill_value=-123456789) - - # Should only be one component, otherwise keep - # reading If more than 1 component, error - if len(st) != 1: - print(st) - print("merge failed?") - - else: - if (st[0].stats.starttime <= start and - st[0].stats.endtime >= end): - st.trim(starttime=start, endtime=end) - - eddt = False - # Check for NoData and convert to NaN if a SAC file - if dtype.upper() == 'SAC': - 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 - - # Check for Nan in stream for SAC - if True in isnan(st[0].data): - print( - "* !!! Missing Data " + - "Present !!! Skipping (NaNs)") - # Check for ND Val in stream for MSEED - elif -123456789 in st[0].data: - print( - "* !!! Missing Data Present !!! " + - "Skipping (MSEED fill)") - else: - if (eddt1 or eddt2) and (ndval == 0.0): - if any(st[0].data == 0.0): - print( - "* !!! Missing " + - "Data Present !!! (Set " + - "to Zero)") - - st[0].stats.update() - tloc = st[0].stats.location - if len(tloc) == 0: - tloc = "--" - - # Processed succesfully...Finish - print(("* {1:3s}.{2:2s} - " + - "From Disk".format( - st[0].stats.station, - st[0].stats.channel.upper(), - tloc))) - return False, st - - except: - pass - else: - st2ot = st2[0].stats.endtime-st2[0].stats.delta - 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. - print("* - Data Unavailable") - 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, local_response_dir='', zcomp='Z'): +# def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='MSEED', altnet=[]): +# """ +# Function to take the list of local directories and recursively +# find all data that matches the station name + +# Parameters +# ---------- +# lcldrs : List +# List of local directories +# sta : Dict +# Station metadata from :mod:`~StDb` +# net : str +# Network name +# altnet : List +# List of alternative networks + +# Returns +# ------- +# fpathmatch : List +# Sorted list of matched directories + +# """ +# from fnmatch import filter +# from os import walk +# from os.path import join + +# if sta is None: +# return [] +# else: +# if net is None: +# sstrings = ['*.{0:s}.*.{1:s}'.format(sta, dtype)] +# else: +# sstrings = ['*.{0:s}.{1:s}.*.{2:s}'.format(net, sta, dtype)] +# if len(altnet) > 0: +# for anet in altnet: +# sstrings.append( +# '*.{0:s}.{1:s}.*.{2:s}'.format(anet, sta, dtype)) + +# fpathmatch = [] +# # Loop over all local data directories +# for lcldr in lcldrs: +# # Recursiely walk through directory +# for root, dirnames, filenames in walk(lcldr): +# # Keep paths only for those matching the station +# for sstring in sstrings: +# for filename in filter(filenames, sstring): +# # No hidden directories +# if not '/.' in root: +# fpathmatch.append(join(root, filename)) + +# fpathmatch.sort() + +# return fpathmatch + + +# def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, +# start=UTCDateTime, end=UTCDateTime, ndval=nan): +# """ +# Function to determine the path to data for a given component and alternate network + +# Parameters +# ---------- +# comp : str +# Channel for seismogram (one letter only) +# stdata : List +# Station list +# sta : Dict +# Station metadata from :mod:`~StDb` data base +# start : :class:`~obspy.core.utcdatetime.UTCDateTime` +# Start time for request +# end : :class:`~obspy.core.utcdatetime.UTCDateTime` +# End time for request +# ndval : float or nan +# Default value for missing data + +# Returns +# ------- +# err : bool +# Boolean for error handling (`False` is associated with success) +# st : :class:`~obspy.core.Stream` +# Stream containing North, East and Vertical components of motion + +# """ + +# from fnmatch import filter + +# # Get start and end parameters +# styr = start.strftime("%Y") +# stjd = start.strftime("%j") +# edyr = end.strftime("%Y") +# edjd = end.strftime("%j") + +# # Intialize to default positive error +# erd = True + +# print( +# ("* {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: + +# 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: +# print("* - Data Unavailable") +# return erd, None + +# # Process the local Files +# for sacfile in lclfiles: +# # Read File +# st = read(sacfile) +# # st = read(sacfile, format="SAC") + +# if dtype.upper() == 'MSEED': +# if len(st) > 1: +# st.merge(method=1, interpolation_samples=- +# 1, fill_value=-123456789) + +# # Should only be one component, otherwise keep reading If more +# # than 1 component, error +# if len(st) != 1: +# pass + +# else: +# # Check start/end times in range +# if (st[0].stats.starttime <= start and +# st[0].stats.endtime >= end): +# st.trim(starttime=start, endtime=end) + +# eddt = False +# # Check for NoData and convert to NaN if a SAC file +# if dtype.upper() == 'SAC': +# 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 + +# # Check for Nan in stream for SAC +# if True in isnan(st[0].data): +# print( +# "* !!! Missing Data Present !!! " + +# "Skipping (NaNs)") +# # Check for ND Val in stream for MSEED +# elif -123456789 in st[0].data: +# print( +# "* !!! Missing Data Present !!! " + +# "Skipping (MSEED fill)") +# else: +# if eddt and (ndval == 0.0): +# if any(st[0].data == 0.0): +# print( +# "* !!! Missing Data Present " + +# "!!! (Set to Zero)") + +# st[0].stats.update() +# tloc = st[0].stats.location +# if len(tloc) == 0: +# tloc = "--" + +# # Processed succesfully...Finish +# print(("* {1:3s}.{2:2s} - From Disk".format( +# st[0].stats.station, st[0].stats.channel.upper(), +# tloc))) +# return False, st + +# # Time Window spans Multiple days +# else: +# 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: +# print("* - Data Unavailable") +# return erd, None + +# # Now try to merge the two separate day files +# if len(lclfiles1) > 0 and len(lclfiles2) > 0: +# # Loop over first day file options +# for sacf1 in lclfiles1: +# st1 = read(sacf1) +# if dtype.upper() == 'MSEED': +# if len(st1) > 1: +# st1.merge(method=1, interpolation_samples=-1, +# fill_value=-123456789) + +# # Loop over second day file options +# for sacf2 in lclfiles2: +# st2 = read(sacf2) +# if dtype.upper() == 'MSEED': +# if len(st2) > 1: +# st2.merge( +# method=1, interpolation_samples=-1, +# fill_value=-123456789) + +# # Check time overlap of the two files. +# if st1[0].stats.endtime >= \ +# st2[0].stats.starttime-st2[0].stats.delta: +# # eddt1 = False +# # eddt2 = False +# # if dtype.upper() == 'SAC': +# # # Check for NoData and convert to NaN +# # st1nd = st1[0].stats.sac['user9'] +# # st2nd = st2[0].stats.sac['user9'] +# # if (not st1nd == 0.0) and (not st1nd == -12345.0): +# # st1[0].data[st1[0].data == st1nd] = ndval +# # eddt1 = True +# # if (not st2nd == 0.0) and (not st2nd == -12345.0): +# # st2[0].data[st2[0].data == st2nd] = ndval +# # eddt2 = True + +# st = st1 + st2 +# # Need to work on this HERE (AJS OCT 2015). +# # If Calibration factors are different, +# # then the traces cannot be merged. +# try: +# st.merge(method=1, interpolation_samples=- +# 1, fill_value=-123456789) + +# # Should only be one component, otherwise keep +# # reading If more than 1 component, error +# if len(st) != 1: +# print(st) +# print("merge failed?") + +# else: +# if (st[0].stats.starttime <= start and +# st[0].stats.endtime >= end): +# st.trim(starttime=start, endtime=end) + +# eddt = False +# # Check for NoData and convert to NaN if a SAC file +# if dtype.upper() == 'SAC': +# 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 + +# # Check for Nan in stream for SAC +# if True in isnan(st[0].data): +# print( +# "* !!! Missing Data " + +# "Present !!! Skipping (NaNs)") +# # Check for ND Val in stream for MSEED +# elif -123456789 in st[0].data: +# print( +# "* !!! Missing Data Present !!! " + +# "Skipping (MSEED fill)") +# else: +# if (eddt1 or eddt2) and (ndval == 0.0): +# if any(st[0].data == 0.0): +# print( +# "* !!! Missing " + +# "Data Present !!! (Set " + +# "to Zero)") + +# st[0].stats.update() +# tloc = st[0].stats.location +# if len(tloc) == 0: +# tloc = "--" + +# # Processed succesfully...Finish +# print(("* {1:3s}.{2:2s} - " + +# "From Disk".format( +# st[0].stats.station, +# st[0].stats.channel.upper(), +# tloc))) +# return False, st + +# except: +# pass +# else: +# st2ot = st2[0].stats.endtime-st2[0].stats.delta +# 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. +# print("* - Data Unavailable") +# 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(), new_sr=0., verbose=False, + remove_response=False, zcomp='Z'): """ 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 given data is already available locally. - Note - ---- - Currently only supports NEZ Components! - Parameters ---------- client : :class:`~obspy.client.fdsn.Client` @@ -421,20 +417,12 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, Start time for request end : :class:`~obspy.core.utcdatetime.UTCDateTime` End time for request - stdata : 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``. zcomp: str - Vertical Component Identifier. Should be a single character. - This is different then 'Z' only for fully unknown component + Vertical Component Identifier. Should be a single character. + This is different then 'Z' only for fully unknown component orientation (i.e., components are 1, 2, 3) Returns @@ -456,106 +444,61 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, from numpy import any from math import floor - # Output - print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, - sta.channel.upper()))) - - # Set Error Default to True - erd = True - - # Check if there is local data - if len(stdata) > 0: - # Only a single day: Search for local data - # Get Z localdata - errZ, stZ = parse_localdata_for_comp( - comp='Z', stdata=stdata, dtype=dtype, sta=sta, start=start, end=end, - ndval=ndval) - # Get N localdata - errN, stN = parse_localdata_for_comp( - comp='N', stdata=stdata, dtype=dtype, sta=sta, start=start, end=end, - ndval=ndval) - # Get E localdata - errE, stE = parse_localdata_for_comp( - comp='E', stdata=stdata, dtype=dtype, sta=sta, start=start, end=end, - ndval=ndval) - # Retreived Succesfully? - erd = errZ or errN or errE - 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: - erd = False - - for loc in sta.location: - tloc = loc - # Construct location name - if len(tloc) == 0: - tloc = "--" - - # Construct Channel List - 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) - - # Use 'zcomp' only for 1, 2, 3 components, where 3 is likely Z - chaZ12 = (sta.channel.upper() + zcomp + ',' + - 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: - if len(st) == 3: - # It's possible if len(st)==1 that data is Z12 - print("* - Data Downloaded") - break - - # Break if we successfully obtained 3 components in st - if not erd: - break + # # Output + # print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, + # sta.channel.upper()))) + + for loc in sta.location: + + # Construct location name + if len(loc) == 0: + tloc = "--" + + # Construct Channel List + cha = sta.channel.upper() + '?' + msg = "* {0:s}.{1:2s}?.{2:2s} - Checking Network".format( + sta.station, sta.channel.upper(), tloc) + 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=cha, + starttime=start, + endtime=end+1.) + except Exception as e: + if verbose: + print("* Met exception:") + print("* " + e.__repr__()) + st = None + else: + if len(st) == 3: + # It's possible if len(st)==1 that data is Z12 + print("* - Data Downloaded") + # break # Check the correct 3 components exist - if st is None: + if st is None or len(st) < 3: print("* Error retrieving waveforms") print("**************************************************") return True, None # Three components successfully retrieved if remove_response: - st.remove_response() - print("*") - print("* Restituted stream to true ground velocity.") - print("*") + try: + st.remove_response() + if verbose: + print("*") + print("* Restituted stream to true ground velocity.") + print("*") + except Exception as e: + print("*") + print('* Cannot remove response, moving on.') + print("*") # Detrend and apply taper st.detrend('demean').detrend('linear').taper( @@ -563,12 +506,13 @@ 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+" " + - str(tr.stats.starttime)+" " + - str(tr.stats.endtime)) for tr in st] - print("* True start: "+str(start)) - print("* -> Shifting traces to true start") + if verbose: + 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)]) @@ -578,17 +522,19 @@ 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") + if verbose: + 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: + except Exception as e: print("* Unable to trim") print("* -> Skipping") print("**************************************************") + return True, None # Check final lengths - they should all be equal if start times From c95519492bdc7a4830d74744e4366d6a468dd415 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Mon, 12 May 2025 20:45:05 +1200 Subject: [PATCH 07/22] some clean up --- rfpy/rfdata.py | 2 +- rfpy/utils.py | 367 +------------------------------------------------ 2 files changed, 5 insertions(+), 364 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 5fcbd3a..96888e5 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -309,7 +309,7 @@ def download_data(self, client, new_sr=5., dts=120., """ Downloads seismograms based on event origin time and - P phase arrival. + P phase arrival and resamples the waveforms. Parameters ---------- diff --git a/rfpy/utils.py b/rfpy/utils.py index c500145..b0a2231 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -39,373 +39,14 @@ def traceshift(trace, tt): return rtrace -# def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='MSEED', altnet=[]): -# """ -# Function to take the list of local directories and recursively -# find all data that matches the station name - -# Parameters -# ---------- -# lcldrs : List -# List of local directories -# sta : Dict -# Station metadata from :mod:`~StDb` -# net : str -# Network name -# altnet : List -# List of alternative networks - -# Returns -# ------- -# fpathmatch : List -# Sorted list of matched directories - -# """ -# from fnmatch import filter -# from os import walk -# from os.path import join - -# if sta is None: -# return [] -# else: -# if net is None: -# sstrings = ['*.{0:s}.*.{1:s}'.format(sta, dtype)] -# else: -# sstrings = ['*.{0:s}.{1:s}.*.{2:s}'.format(net, sta, dtype)] -# if len(altnet) > 0: -# for anet in altnet: -# sstrings.append( -# '*.{0:s}.{1:s}.*.{2:s}'.format(anet, sta, dtype)) - -# fpathmatch = [] -# # Loop over all local data directories -# for lcldr in lcldrs: -# # Recursiely walk through directory -# for root, dirnames, filenames in walk(lcldr): -# # Keep paths only for those matching the station -# for sstring in sstrings: -# for filename in filter(filenames, sstring): -# # No hidden directories -# if not '/.' in root: -# fpathmatch.append(join(root, filename)) - -# fpathmatch.sort() - -# return fpathmatch - - -# def parse_localdata_for_comp(comp='Z', stdata=[], dtype='SAC', sta=None, -# start=UTCDateTime, end=UTCDateTime, ndval=nan): -# """ -# Function to determine the path to data for a given component and alternate network - -# Parameters -# ---------- -# comp : str -# Channel for seismogram (one letter only) -# stdata : List -# Station list -# sta : Dict -# Station metadata from :mod:`~StDb` data base -# start : :class:`~obspy.core.utcdatetime.UTCDateTime` -# Start time for request -# end : :class:`~obspy.core.utcdatetime.UTCDateTime` -# End time for request -# ndval : float or nan -# Default value for missing data - -# Returns -# ------- -# err : bool -# Boolean for error handling (`False` is associated with success) -# st : :class:`~obspy.core.Stream` -# Stream containing North, East and Vertical components of motion - -# """ - -# from fnmatch import filter - -# # Get start and end parameters -# styr = start.strftime("%Y") -# stjd = start.strftime("%j") -# edyr = end.strftime("%Y") -# edjd = end.strftime("%j") - -# # Intialize to default positive error -# erd = True - -# print( -# ("* {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: - -# 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: -# print("* - Data Unavailable") -# return erd, None - -# # Process the local Files -# for sacfile in lclfiles: -# # Read File -# st = read(sacfile) -# # st = read(sacfile, format="SAC") - -# if dtype.upper() == 'MSEED': -# if len(st) > 1: -# st.merge(method=1, interpolation_samples=- -# 1, fill_value=-123456789) - -# # Should only be one component, otherwise keep reading If more -# # than 1 component, error -# if len(st) != 1: -# pass - -# else: -# # Check start/end times in range -# if (st[0].stats.starttime <= start and -# st[0].stats.endtime >= end): -# st.trim(starttime=start, endtime=end) - -# eddt = False -# # Check for NoData and convert to NaN if a SAC file -# if dtype.upper() == 'SAC': -# 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 - -# # Check for Nan in stream for SAC -# if True in isnan(st[0].data): -# print( -# "* !!! Missing Data Present !!! " + -# "Skipping (NaNs)") -# # Check for ND Val in stream for MSEED -# elif -123456789 in st[0].data: -# print( -# "* !!! Missing Data Present !!! " + -# "Skipping (MSEED fill)") -# else: -# if eddt and (ndval == 0.0): -# if any(st[0].data == 0.0): -# print( -# "* !!! Missing Data Present " + -# "!!! (Set to Zero)") - -# st[0].stats.update() -# tloc = st[0].stats.location -# if len(tloc) == 0: -# tloc = "--" - -# # Processed succesfully...Finish -# print(("* {1:3s}.{2:2s} - From Disk".format( -# st[0].stats.station, st[0].stats.channel.upper(), -# tloc))) -# return False, st - -# # Time Window spans Multiple days -# else: -# 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: -# print("* - Data Unavailable") -# return erd, None - -# # Now try to merge the two separate day files -# if len(lclfiles1) > 0 and len(lclfiles2) > 0: -# # Loop over first day file options -# for sacf1 in lclfiles1: -# st1 = read(sacf1) -# if dtype.upper() == 'MSEED': -# if len(st1) > 1: -# st1.merge(method=1, interpolation_samples=-1, -# fill_value=-123456789) - -# # Loop over second day file options -# for sacf2 in lclfiles2: -# st2 = read(sacf2) -# if dtype.upper() == 'MSEED': -# if len(st2) > 1: -# st2.merge( -# method=1, interpolation_samples=-1, -# fill_value=-123456789) - -# # Check time overlap of the two files. -# if st1[0].stats.endtime >= \ -# st2[0].stats.starttime-st2[0].stats.delta: -# # eddt1 = False -# # eddt2 = False -# # if dtype.upper() == 'SAC': -# # # Check for NoData and convert to NaN -# # st1nd = st1[0].stats.sac['user9'] -# # st2nd = st2[0].stats.sac['user9'] -# # if (not st1nd == 0.0) and (not st1nd == -12345.0): -# # st1[0].data[st1[0].data == st1nd] = ndval -# # eddt1 = True -# # if (not st2nd == 0.0) and (not st2nd == -12345.0): -# # st2[0].data[st2[0].data == st2nd] = ndval -# # eddt2 = True - -# st = st1 + st2 -# # Need to work on this HERE (AJS OCT 2015). -# # If Calibration factors are different, -# # then the traces cannot be merged. -# try: -# st.merge(method=1, interpolation_samples=- -# 1, fill_value=-123456789) - -# # Should only be one component, otherwise keep -# # reading If more than 1 component, error -# if len(st) != 1: -# print(st) -# print("merge failed?") - -# else: -# if (st[0].stats.starttime <= start and -# st[0].stats.endtime >= end): -# st.trim(starttime=start, endtime=end) - -# eddt = False -# # Check for NoData and convert to NaN if a SAC file -# if dtype.upper() == 'SAC': -# 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 - -# # Check for Nan in stream for SAC -# if True in isnan(st[0].data): -# print( -# "* !!! Missing Data " + -# "Present !!! Skipping (NaNs)") -# # Check for ND Val in stream for MSEED -# elif -123456789 in st[0].data: -# print( -# "* !!! Missing Data Present !!! " + -# "Skipping (MSEED fill)") -# else: -# if (eddt1 or eddt2) and (ndval == 0.0): -# if any(st[0].data == 0.0): -# print( -# "* !!! Missing " + -# "Data Present !!! (Set " + -# "to Zero)") - -# st[0].stats.update() -# tloc = st[0].stats.location -# if len(tloc) == 0: -# tloc = "--" - -# # Processed succesfully...Finish -# print(("* {1:3s}.{2:2s} - " + -# "From Disk".format( -# st[0].stats.station, -# st[0].stats.channel.upper(), -# tloc))) -# return False, st - -# except: -# pass -# else: -# st2ot = st2[0].stats.endtime-st2[0].stats.delta -# 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. -# print("* - Data Unavailable") -# 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(), new_sr=0., verbose=False, remove_response=False, zcomp='Z'): """ - 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 - given data is already available locally. + Function to build a stream object for a seismogram in a given time window + by getting data from a client object, either from a local SDS archive or + from an FDSN web-service. The function performs sanity checks for + the start times, sampling rates and window lengths. Parameters ---------- From d5aabdf6838ecc41b204bc926707053af35362bc Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 14 May 2025 11:26:23 +1200 Subject: [PATCH 08/22] non-bare except statements --- rfpy/scripts/rfpy_calc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 8e206f4..e9b7e25 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -395,7 +395,7 @@ def get_calc_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Cannot construct UTCDateTime from start time: " + args.startT) @@ -406,7 +406,7 @@ def get_calc_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Cannot construct UTCDateTime from end time: " + args.endT) From 7896b8835356010da633f270f12410fedfd85e53 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 14 May 2025 11:26:51 +1200 Subject: [PATCH 09/22] deep copy for tlocs --- rfpy/utils.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/rfpy/utils.py b/rfpy/utils.py index b0a2231..6307dac 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -2,6 +2,7 @@ from obspy import UTCDateTime from numpy import nan, isnan, abs import numpy as np +import copy from obspy import Stream, Inventory from obspy import read, read_inventory @@ -58,6 +59,10 @@ def download_data(client=None, sta=None, start=UTCDateTime(), Start time for request end : :class:`~obspy.core.utcdatetime.UTCDateTime` End time for request + new_st : float + New sampling rate (Hz) + verbose : bool + Whether or not to print messages to screen during run-time remove_response : bool Remove instrument response from seismogram and resitute to true ground velocity (m/s) using obspy.core.trace.Trace.remove_response() @@ -94,6 +99,8 @@ def download_data(client=None, sta=None, start=UTCDateTime(), # Construct location name if len(loc) == 0: tloc = "--" + else: + tloc = copy.copy(loc) # Construct Channel List cha = sta.channel.upper() + '?' From c151cbed7f4770a3101eef11c6c1230fd5cd1ff0 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 14 May 2025 13:44:48 +1200 Subject: [PATCH 10/22] better managing of locs --- rfpy/scripts/rfpy_calc.py | 11 ++++++----- rfpy/scripts/rfpy_harmonics.py | 6 +++--- rfpy/scripts/rfpy_hk.py | 6 +++--- rfpy/scripts/rfpy_plot.py | 6 +++--- rfpy/scripts/rfpy_recalc.py | 8 +++----- rfpy/utils.py | 8 ++++---- 6 files changed, 22 insertions(+), 23 deletions(-) diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index e9b7e25..4c4494a 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -27,6 +27,7 @@ import numpy as np import pickle import stdb +import copy from obspy.clients.fdsn import Client as FDSN_Client from obspy.clients.filesystem.sds import Client as SDS_Client from obspy import Catalog, UTCDateTime @@ -43,8 +44,6 @@ def get_calc_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( @@ -588,9 +587,11 @@ def main(): continue # Temporary print locations - tlocs = [] - for il in range(0, len(sta.location)): - if len(sta.location[il]) == 0: + tlocs = copy.copy(sta.location) + if len(tlocs) == 0: + tlocs = [''] + for il in range(0, len(tlocs)): + if len(tlocs[il]) == 0: tlocs.append("--") # Update Display diff --git a/rfpy/scripts/rfpy_harmonics.py b/rfpy/scripts/rfpy_harmonics.py index 273cb4d..98bbafe 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -27,6 +27,7 @@ import numpy as np import pickle import stdb +import copy from obspy.clients.fdsn import Client from obspy.core import Stream, UTCDateTime from rfpy import binning, plotting, Harmonics @@ -391,13 +392,12 @@ def main(): continue # Temporary print locations - tlocs = sta.location + tlocs = copy.copy(sta.location) if len(tlocs) == 0: tlocs = [''] for il in range(0, len(tlocs)): if len(tlocs[il]) == 0: - tlocs[il] = "--" - sta.location = tlocs + tlocs.append("--") # Update Display print(" ") diff --git a/rfpy/scripts/rfpy_hk.py b/rfpy/scripts/rfpy_hk.py index 5741818..66a2adc 100644 --- a/rfpy/scripts/rfpy_hk.py +++ b/rfpy/scripts/rfpy_hk.py @@ -27,6 +27,7 @@ import numpy as np import pickle import stdb +import copy from obspy.clients.fdsn import Client from obspy.core import Stream, UTCDateTime from rfpy import binning, plotting, HkStack @@ -542,13 +543,12 @@ def main(): continue # Temporary print locations - tlocs = sta.location + tlocs = copy.copy(sta.location) if len(tlocs) == 0: tlocs = [''] for il in range(0, len(tlocs)): if len(tlocs[il]) == 0: - tlocs[il] = "--" - sta.location = tlocs + tlocs.append("--") # Update Display print(" ") diff --git a/rfpy/scripts/rfpy_plot.py b/rfpy/scripts/rfpy_plot.py index 310d782..6c69754 100644 --- a/rfpy/scripts/rfpy_plot.py +++ b/rfpy/scripts/rfpy_plot.py @@ -27,6 +27,7 @@ import numpy as np import pickle import stdb +import copy from obspy import Stream, UTCDateTime from rfpy import binning, plotting from pathlib import Path @@ -368,13 +369,12 @@ def main(): continue # Temporary print locations - tlocs = sta.location + tlocs = copy.copy(sta.location) if len(tlocs) == 0: tlocs = [''] for il in range(0, len(tlocs)): if len(tlocs[il]) == 0: - tlocs[il] = "--" - sta.location = tlocs + tlocs.append("--") # Update Display print(" ") diff --git a/rfpy/scripts/rfpy_recalc.py b/rfpy/scripts/rfpy_recalc.py index 74cb57a..73b58bd 100644 --- a/rfpy/scripts/rfpy_recalc.py +++ b/rfpy/scripts/rfpy_recalc.py @@ -29,6 +29,7 @@ import numpy as np import pickle import stdb +import copy from rfpy import RFData from pathlib import Path @@ -37,8 +38,6 @@ 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( @@ -285,13 +284,12 @@ def main(): continue # Temporary print locations - tlocs = sta.location + tlocs = copy.copy(sta.location) if len(tlocs) == 0: tlocs = [''] for il in range(0, len(tlocs)): if len(tlocs[il]) == 0: - tlocs[il] = "--" - sta.location = tlocs + tlocs.append("--") # Update Display print(" ") diff --git a/rfpy/utils.py b/rfpy/utils.py index 6307dac..1dcba51 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -97,15 +97,15 @@ def download_data(client=None, sta=None, start=UTCDateTime(), for loc in sta.location: # Construct location name - if len(loc) == 0: - tloc = "--" + if loc == "--": + tloc = "" else: tloc = copy.copy(loc) # Construct Channel List cha = sta.channel.upper() + '?' msg = "* {0:s}.{1:2s}?.{2:2s} - Checking Network".format( - sta.station, sta.channel.upper(), tloc) + sta.station, sta.channel.upper(), loc) print(msg) # Get waveforms, with extra 1 second to avoid @@ -114,7 +114,7 @@ def download_data(client=None, sta=None, start=UTCDateTime(), st = client.get_waveforms( network=sta.network, station=sta.station, - location=loc, + location=tloc, channel=cha, starttime=start, endtime=end+1.) From 5b30da84659c83336271377a11accc9e5e69026b Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 14 May 2025 14:27:40 +1200 Subject: [PATCH 11/22] add SDS in documentation --- docs/rfpy.rst | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/docs/rfpy.rst b/docs/rfpy.rst index da3095e..5011b6f 100644 --- a/docs/rfpy.rst +++ b/docs/rfpy.rst @@ -83,6 +83,48 @@ Installing from source pip install . +Using local data +================ + +The main script packaged with ``RfPy`` uses FDSN web services through and ``ObsPy`` `Client` to load waveform data. For waveform data locally stored on your hard drive, the scripts can use a `Client` that reads a `SeisComP Data Structure `_ archive containing SAC or miniSEED waveform data. Check out the scripts ``rfpy_calc`` below and the argument ``--local-data`` and ``--dtype`` for more details. + +Station Metadata +---------------- + +If you have data stored locally on your drive, it is likely you also have a station `XML `_ file containing the metadata. The corresponding ObsPy documentation is `here `_. + +To convert the station `XML` file to an input that can be read by ``OrientPy``, you run the command ``gen_stdb station.xml`` (only available on StDb version 0.2.7), which will create the file ``station.pkl``. If you don't have a station `XML` file but you have a dataless SEED file, you can convert it first to `XML` using `this tools `_. + +Waveform Data +------------- + +The SDS folder containing the waveform data has the structure: + +.. code-block:: python + + archive + + year + + network code + + station code + + channel code + type + + one file per day and location, e.g. NET.STA.LOC.CHAN.TYPE.YEAR.DOY + + +For example: + +.. code-block:: python + + SDS/ + 2014/ + YH/ + LOBS3/ + HH1.D/ + YH.LOBS3..CH1.D.2014.332 + ... + + +Note, the filename does not include the extension (`.MSEED` or `.SAC`), and the characters `.D` (for type Data) that appear in both the channel code and the filename. Note also the two dots (`..`). If there is a location code, it should appear between those dots (e.g., for a location code `10`, the corresponding filename should be `YH.LOBS3.10.HH1.D.2014.332`). There is no location code for the YH.LOBS3 data, and this field is simply absent from the filenames. Finally, the day-of-year (DOY) field must be zero-padded to be exactly 3 characters. + Basic Usage =========== From d756197da562c55c8e00ec860592f778c82f89b2 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 14 May 2025 14:29:37 +1200 Subject: [PATCH 12/22] remove docs from deploy actions --- .github/workflows/deploy.yml | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/.github/workflows/deploy.yml b/.github/workflows/deploy.yml index 3c5096c..4653f12 100644 --- a/.github/workflows/deploy.yml +++ b/.github/workflows/deploy.yml @@ -52,22 +52,22 @@ jobs: # pytest -v --cov=orientpy ../orientpy/tests/ # bash <(curl -s https://codecov.io/bash) - - name: Make docs - if: matrix.python-version == '3.12' - shell: bash -l {0} - run: | - cd docs - conda install sphinx - pip install sphinx_rtd_theme - make html - touch _build/html/.nojekyll - cd .. + # - name: Make docs + # if: matrix.python-version == '3.12' + # shell: bash -l {0} + # run: | + # cd docs + # conda install sphinx + # pip install sphinx_rtd_theme + # make html + # touch _build/html/.nojekyll + # cd .. - - name: Deploy 🚀 - if: matrix.python-version == '3.12' - uses: JamesIves/github-pages-deploy-action@3.7.1 - with: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - BRANCH: gh-pages # The branch the action should deploy to. - FOLDER: docs/_build/html # The folder the action should deploy. - CLEAN: true # Automatically remove deleted files from the deploy branch + # - name: Deploy 🚀 + # if: matrix.python-version == '3.12' + # uses: JamesIves/github-pages-deploy-action@3.7.1 + # with: + # GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + # BRANCH: gh-pages # The branch the action should deploy to. + # FOLDER: docs/_build/html # The folder the action should deploy. + # CLEAN: true # Automatically remove deleted files from the deploy branch From 728d8e25b4f2a6ac1b0e8e8f47daf0855aa7a0de Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 16 May 2025 09:11:40 +1200 Subject: [PATCH 13/22] some updates --- docs/index.rst | 5 ----- rfpy/utils.py | 11 ++--------- 2 files changed, 2 insertions(+), 14 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index c67bfae..e6fc605 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -17,8 +17,6 @@ used in command-line scripts. :target: https://doi.org/10.5281/zenodo.3905414 .. image:: https://travis-ci.com/paudetseis/RfPy.svg?branch=master :target: https://travis-ci.com/paudetseis/RfPy -.. image:: https://codecov.io/gh/paudetseis/RfPy/branch/master/graph/badge.svg - :target: https://codecov.io/gh/paudetseis/RfPy .. note:: @@ -33,9 +31,6 @@ used in command-line scripts. * ``RfData`` objects are used to calculate single-station and single-event receiver functions, whereas ``rf`` can handle multiple stations at once. -.. warning:: - ``RfPy`` was recently updated to fix a problem when running the scripts under Windows OS. The consequence is that version ``0.1.1`` will throw an error if the extension .py is specified when calling the scripts. The accompanying documentation also uses version ``0.2.1`` of ``StDb`` in the Tutorials section. - .. toctree:: :maxdepth: 1 :caption: Quick Links diff --git a/rfpy/utils.py b/rfpy/utils.py index 1dcba51..38ee25c 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -3,8 +3,7 @@ from numpy import nan, isnan, abs import numpy as np import copy -from obspy import Stream, Inventory -from obspy import read, read_inventory +from obspy import Stream, read def floor_decimal(n, decimals=0): @@ -59,7 +58,7 @@ def download_data(client=None, sta=None, start=UTCDateTime(), Start time for request end : :class:`~obspy.core.utcdatetime.UTCDateTime` End time for request - new_st : float + new_sr : float New sampling rate (Hz) verbose : bool Whether or not to print messages to screen during run-time @@ -84,12 +83,6 @@ def download_data(client=None, sta=None, start=UTCDateTime(), """ - from fnmatch import filter - from obspy import read, Stream - from os.path import dirname, join, exists - from numpy import any - from math import floor - # # Output # print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, # sta.channel.upper()))) From 4332aae96b1cbdf730b7eb855986bd5e52fd5346 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 16 May 2025 11:33:22 +1200 Subject: [PATCH 14/22] updated --local-data arg description --- rfpy/rfdata.py | 2 +- rfpy/scripts/rfpy_calc.py | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 96888e5..448ef30 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -255,7 +255,7 @@ def add_data(self, stream, returned=False, new_sr=5.): Attributes ---------- - zne_data : :class:`~obspy.core.Stream` + data : :class:`~obspy.core.Stream` Stream container for NEZ seismograms Returns diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 4c4494a..cd8685c 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -155,11 +155,12 @@ def get_calc_arguments(argv=None): type=str, dest="localdata", default=None, - help="Specify path 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 FDSN Client interface") + help="Specify absolute path to a SeisComP Data Structure (SDS) " + + "archive containing day-long SAC or MSEED files" + + "(e.g., --local-data=/Home/username/Data/SDS). " + + "See https://www.seiscomp.de/seiscomp3/doc/applications/slarchive/SDS.html " + + "for details on the SDS format. If this option is used, it takes " + + "precedence over the --server settings.") DataGroup.add_argument( "--dtype", action="store", From a219eaf813df43b0d185c3233ed80091ec9f68ed Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 16 May 2025 13:03:28 +1200 Subject: [PATCH 15/22] esthetic change --- rfpy/scripts/rfpy_calc.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index cd8685c..0d73f97 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -634,8 +634,10 @@ def main(): # Get catalogue using deployment start and end try: cat = event_client.get_events( - starttime=tstart, endtime=tend, - minmagnitude=args.minmag, maxmagnitude=args.maxmag) + starttime=tstart, + endtime=tend, + minmagnitude=args.minmag, + maxmagnitude=args.maxmag) except IncompleteRead: # See http.client.IncompleteRead @@ -653,8 +655,10 @@ def main(): print("| Start: {0:19s} |".format( tstart.strftime("%Y-%m-%d %H:%M:%S"))) cat += event_client.get_events( - starttime=tstart, endtime=tstart + chunk, - minmagnitude=args.minmag, maxmagnitude=args.maxmag) + starttime=tstart, + endtime=tstart + chunk, + minmagnitude=args.minmag, + maxmagnitude=args.maxmag) # Make sure that we go all the way to tend if tstart + chunk > tend: From baca92af01937206a2f3a041fa742cd0cce986e1 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Tue, 20 May 2025 20:30:31 +1200 Subject: [PATCH 16/22] view source code in docs --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index a72cb0e..58bb415 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -29,7 +29,7 @@ # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon'] +extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon', 'sphinx.ext.viewcode'] autodoc_member_order = 'bysource' From 69ea18af62cb4a065a9c3fc0510f108ca04108a2 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Tue, 20 May 2025 21:16:28 +1200 Subject: [PATCH 17/22] cosmetic changes --- docs/scripts.rst | 12 +- docs/tutorials.rst | 10 ++ rfpy/binning.py | 23 ++- rfpy/ccp.py | 258 +++++++++++++++++++-------------- rfpy/harmonics.py | 19 ++- rfpy/hk.py | 86 +++++------ rfpy/plotting.py | 65 +-------- rfpy/rfdata.py | 170 ++++++++++++++-------- rfpy/scripts/rfpy_calc.py | 42 +----- rfpy/scripts/rfpy_ccp.py | 45 +++--- rfpy/scripts/rfpy_harmonics.py | 20 +-- rfpy/scripts/rfpy_hk.py | 16 +- rfpy/scripts/rfpy_recalc.py | 9 +- rfpy/utils.py | 31 +++- 14 files changed, 409 insertions(+), 397 deletions(-) diff --git a/docs/scripts.rst b/docs/scripts.rst index d955b94..afc060d 100644 --- a/docs/scripts.rst +++ b/docs/scripts.rst @@ -53,7 +53,7 @@ Usage instance, providing IU will match with all stations in the IU network [Default processes all stations in the database] - -v, -V, --verbose Specify to increase verbosity. + -V, --verbose Specify to increase verbosity. -O, --overwrite Force the overwriting of pre-existing data. [Default False] --zcomp ZCOMP Specify the Vertical Component Channel Identifier. @@ -228,7 +228,7 @@ Usage providing IU will match with all stations in the IU network [Default processes all stations in the database] - -v, -V, --verbose Specify to increase verbosity. + -V, --verbose Specify to increase verbosity. -L, --long-name 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 @@ -315,7 +315,7 @@ Usage instance, providing IU will match with all stations in the IU network [Default processes all stations in the database] - -v, -V, --verbose Specify to increase verbosity. + -V, --verbose Specify to increase verbosity. -O, --overwrite Force the overwriting of pre-existing figures. [Default False] -L, --long-name Force folder names to use long-key form (NET.STN.CHN). @@ -426,7 +426,7 @@ Usage instance, providing IU will match with all stations in the IU network [Default processes all stations in the database] - -v, -V, --verbose Specify to increase verbosity. + -V, --verbose Specify to increase verbosity. -O, --overwrite Force the overwriting of pre-existing data. [Default False] -L, --long-name Force folder names to use long-key form (NET.STN.CHN). @@ -565,7 +565,7 @@ Usage 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] - -v, -V, --verbose Specify to increase verbosity. + -V, --verbose Specify to increase verbosity. -O, --overwrite Force the overwriting of pre-existing data. [Default False] -L, --long-name Force folder names to use long-key form (NET.STN.CHN). @@ -689,7 +689,7 @@ Usage providing IU will match with all stations in the IU network [Default processes all stations in the database] - -v, -V, --verbose Specify to increase verbosity. + -V, --verbose Specify to increase verbosity. -O, --overwrite Force the overwriting of pre-existing data. [Default False] -L, --long-name Force folder names to use long-key form (NET.STN.CHN). diff --git a/docs/tutorials.rst b/docs/tutorials.rst index bbbc06c..6d3854e 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -68,6 +68,16 @@ An example log printed on the terminal will look like: .. code-block:: + ############################################## + # __ _ # + # _ __ / _|_ __ _ _ ___ __ _| | ___ # + # | '__| |_| '_ \| | | | / __/ _` | |/ __| # + # | | | _| |_) | |_| | | (_| (_| | | (__ # + # |_| |_| | .__/ \__, |___\___\__,_|_|\___| # + # |_| |___/_____| # + # # + ############################################## + |===============================================| |===============================================| | MMPY | diff --git a/rfpy/binning.py b/rfpy/binning.py index ea87ee0..f783f25 100644 --- a/rfpy/binning.py +++ b/rfpy/binning.py @@ -34,7 +34,7 @@ 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 mean), or using phase-weighted stacking. @@ -66,8 +66,8 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, """ - if not typ in ['baz', 'slow', 'dist']: - raise(Exception("Type has to be 'baz' or 'slow' or 'dist'")) + if typ not in ['baz', 'slow', 'dist']: + raise Exception("Type has to be 'baz' or 'slow' or 'dist'") if typ == 'baz': stat = [stream1[i].stats.baz for i in range(len(stream1))] @@ -116,7 +116,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, hilb = hilbert(tr.data) phase = np.arctan2(hilb.imag, hilb.real) pweight += np.exp(1j*phase) - + continue if nb > 0 or include_empty: @@ -131,7 +131,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, trace.stats.baz = bins[i] trace.stats.slow = alt trace.stats.nbin = nb - elif typ == 'slow': + elif typ == 'slow': trace.stats.slow = bins[i] trace.stats.baz = alt trace.stats.nbin = nb @@ -148,7 +148,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+1, pws=False, final_stream.append(binned_stream) - except: + except Exception: continue return final_stream @@ -156,7 +156,7 @@ def bin(stream1, stream2=None, typ='baz', nbin=36+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 mean), or using phase-weighted stacking. @@ -223,7 +223,7 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, hilb = hilbert(tr.data) phase = np.arctan2(hilb.imag, hilb.real) pweight += np.exp(1j*phase) - + continue if nbin > 0 or include_empty: @@ -244,14 +244,14 @@ def bin_baz_slow(stream1, stream2=None, nbaz=36+1, nslow=20+1, pws=False, final_stream.append(binned_stream) - except: + except Exception: continue return final_stream 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. @@ -302,8 +302,7 @@ def bin_all(stream1, stream2=None, pws=False): # Put back into traces stack.append(Trace(data=weight*array, header=stats)) - except: + except Exception: continue return stack - diff --git a/rfpy/ccp.py b/rfpy/ccp.py index bc49894..c50ad50 100644 --- a/rfpy/ccp.py +++ b/rfpy/ccp.py @@ -43,41 +43,41 @@ class CCPimage(object): Common Conversion Point (CCP) stacks for each of the main three main phases (Ps, Pps and Pss) using radial-component receiver functions. The object is used to project the stacks along - a linear profile, specified by start and end geographical coordinate - locations, which are subsequently averaged to produce a final - CCP image. The averaging can be done using a linear weighted sum, + a linear profile, specified by start and end geographical coordinate + locations, which are subsequently averaged to produce a final + CCP image. The averaging can be done using a linear weighted sum, or a phase-weighted sum. Methods should be used in the appropriate sequence (see ``rfpy_ccp.py`` for details). Note ---- - By default, the object initializes with un-defined coordinate locations for + By default, the object initializes with un-defined coordinate locations for the profile. If not specified during initialization, make sure they are specified later, before the other methods are used, e.g. - ``ccpimage.xs_lat1 = 10.; ccpimage.xs_lon1 = 110.``, etc. Note also that the - default 1D velocity model may not be applicable to your region of + ``ccpimage.xs_lat1 = 10.; ccpimage.xs_lon1 = 110.``, etc. Note also that + the default 1D velocity model may not be applicable to your region of interest and a different model can be implemented during initialization - or later during processing. + or later during processing. Parameters ---------- coord_start : list List of two floats corresponding to the (latitude, longitude) - pair for the start point of the profile + pair for the start point of the profile coord_end : list List of two floats corresponding to the (latitude, longitude) - pair for the end point of the profile + pair for the end point of the profile weights : list List of three floats with corresponding weights for the Ps, Pps and Pss phases used during linear, weighted averaging dep : :class:`~numpy.ndarray` - Array of depth values defining the 1D background seismic velocity model. - Note that the maximum depth defined here sets the maximum depth + Array of depth values defining the 1D background seismic velocity + model. Note that the maximum depth defined here sets the maximum depth in each of the CCP stacks and the final CCP image. vp : :class:`~numpy.ndarray` Array of Vp values defining the 1D background seismic velocity model vpvs : float - Constant Vp/Vs ratio for the 1D model. + Constant Vp/Vs ratio for the 1D model. Other Parameters ---------------- @@ -129,8 +129,11 @@ def __init__(self, coord_start=[None, None], coord_end=[None, None], self.dx = dx # Get total length of grid from end points - xlength = haversine(self.xs_lat1, self.xs_lon1, - self.xs_lat2, self.xs_lon2) + xlength = haversine( + self.xs_lat1, + self.xs_lon1, + self.xs_lat2, + self.xs_lon2) # number of cells laterally and vertically self.nx = int(np.rint(xlength/self.dx)) @@ -144,7 +147,7 @@ def __init__(self, coord_start=[None, None], coord_end=[None, None], vs = vp/vpvs else: if not vp.shape == vs.shape: - raise(Exception("vp and vs arrays have a different shape")) + raise Exception("vp and vs arrays have a different shape") self.vp = sp.interpolate.interp1d(dep, vp, kind='linear')(self.zarray) self.vs = sp.interpolate.interp1d(dep, vs, kind='linear')(self.zarray) @@ -154,7 +157,7 @@ def add_rfstream(self, rfstream): Method to add a :class:`~obspy.core.Stream` object to the list ``radialRF``. With at least one stream in ``radialRF``, the object is ready for the ``prep_data`` method, and the corresponding flag is - updated. + updated. Parameters ---------- @@ -163,8 +166,8 @@ def add_rfstream(self, rfstream): """ if len(rfstream) > 0: - # fftshift if the time axis starts at negative lags - if rfstream[0].stats.taxis[0]<0.: + # fftshift if the time axis starts at negative lags + if rfstream[0].stats.taxis[0] < 0.: for tr in rfstream: tr.data = np.fft.fftshift(tr.data) @@ -174,13 +177,13 @@ def add_rfstream(self, rfstream): def prep_data(self, f1=0.05, f2ps=0.5, f2pps=0.25, f2pss=0.2, nbaz=36+1, nslow=40+1): """ - Method to pre-process the data and calculate the CCP points for each + Method to pre-process the data and calculate the CCP points for each of the receiver functions. Pre-processing includes the binning to back-azimuth and slowness bins to reduce processing time and generate cleaner stacks, as well as filtering to emphasize the energy of the - various phases as a function of depth. As a general rule, the high - frequency corner of the 2 reverberated phases (Pps and Pss) should be - approximately half of the high frequency corner of the direct + various phases as a function of depth. As a general rule, the high + frequency corner of the 2 reverberated phases (Pps and Pss) should be + approximately half of the high frequency corner of the direct (Ps) phase. At the end of this step, the object is updated with the amplitude at the lat, lon and depth values corresponding to the raypath for each receiver function. The object is now ready for the @@ -189,13 +192,17 @@ def prep_data(self, f1=0.05, f2ps=0.5, f2pps=0.25, f2pss=0.2, Parameters ---------- f1 : float - Low-frequency corner of the bandpass filter used for all phases (Hz) + Low-frequency corner of the bandpass filter used for all phases + (Hz) f2ps : float - High-frequency corner of the bandpass filter used for the Ps phase (Hz) + High-frequency corner of the bandpass filter used for the Ps phase + (Hz) f2pps : float - High-frequency corner of the bandpass filter used for the Pps phase (Hz) + High-frequency corner of the bandpass filter used for the Pps phase + (Hz) f2pss : float - High-frequency corner of the bandpass filter used for the Pss phase (Hz) + High-frequency corner of the bandpass filter used for the Pss phase + (Hz) nbaz : int Number of increments in the back-azimuth bins nslow : int @@ -212,7 +219,8 @@ def prep_data(self, f1=0.05, f2ps=0.5, f2pps=0.25, f2pss=0.2, amp_pss_depth : :class:`numpy.ndarray` 2D array of amplitudes as a function of depth for the Pss phase lon_depth : :class:`numpy.ndarray` - 2D array of longitude as a function of depth (i.e., piercing points) + 2D array of longitude as a function of depth (i.e., piercing + points) lat_depth : :class:`numpy.ndarray` 2D array of latitude as a function of depth (i.e., piercing points) is_ready_for_presstack : boolean @@ -223,7 +231,7 @@ def prep_data(self, f1=0.05, f2ps=0.5, f2pps=0.25, f2pss=0.2, """ if not self.is_ready_for_prep: - raise(Exception("CCPimage not ready for pre-prep")) + raise Exception("CCPimage not ready for pre-prep") ikey = 0 total_traces = 0 @@ -250,22 +258,34 @@ def prep_data(self, f1=0.05, f2ps=0.5, f2pps=0.25, f2pss=0.2, # Filter Ps, Pps and Pss st_ps.filter( - 'bandpass', freqmin=f1, freqmax=f2ps, - corners=4, zerophase=True) + 'bandpass', + freqmin=f1, + freqmax=f2ps, + corners=4, + zerophase=True) st_pps.filter( - 'bandpass', freqmin=f1, freqmax=f2pps, - corners=4, zerophase=True) + 'bandpass', + freqmin=f1, + freqmax=f2pps, + corners=4, + zerophase=True) st_pss.filter( - 'bandpass', freqmin=f1, freqmax=f2pss, - corners=4, zerophase=True) + 'bandpass', + freqmin=f1, + freqmax=f2pss, + corners=4, + zerophase=True) del RFbin print("Station: "+st_ps[0].stats.station) for itr in _progressbar(range(len(st_ps)), '', 25): # Get raypath and travel time for all phases - tt_ps, tt_pps, tt_pss, plon, plat = \ - raypath(st_ps[itr], dep=self.zarray, vp=self.vp, vs=self.vs) + tt_ps, tt_pps, tt_pss, plon, plat = raypath( + st_ps[itr], + dep=self.zarray, + vp=self.vp, + vs=self.vs) # Now get amplitude of RF at corresponding travel # time along the raypath @@ -327,13 +347,13 @@ def prep_data(self, f1=0.05, f2ps=0.5, f2pps=0.25, f2pss=0.2, def prestack(self): """ - Method to project the raypaths onto the 2D profile for each of the three - phases. The final grid is defined here, using the parameter ``dx`` in km. - The horizontal extent is pre-determined from the start and end points of - the profile. At the end of this step, the object contains the set of - amplitudes at each of the 2D grid points, for each of the three phases. - The object is now ready for the methods ``ccp`` and/or ``gccp``, with - the corresponding flag updated. + Method to project the raypaths onto the 2D profile for each of the + three phases. The final grid is defined here, using the parameter + ``dx`` in km. The horizontal extent is pre-determined from the start + and end points of the profile. At the end of this step, the object + contains the set of amplitudes at each of the 2D grid points, for + each of the three phases. The object is now ready for the methods + ``ccp`` and/or ``gccp``, with the corresponding flag updated. The following attributes are added to the object: @@ -356,7 +376,7 @@ def prestack(self): """ if not self.is_ready_for_prestack: - raise(Exception("CCPimage not ready for prestack")) + raise Exception("CCPimage not ready for prestack") xs_latitudes = np.asarray( np.linspace(self.xs_lat1, self.xs_lat2, self.nx)) @@ -383,7 +403,7 @@ def prestack(self): minimum_distance = np.amin(distance_tests) ix = np.where(distance_tests == - np.amin(distance_tests))[0][0] + np.amin(distance_tests))[0][0] nonzero_count = np.count_nonzero( xs_amps_ps[iz, ix, :]) @@ -423,10 +443,10 @@ def prestack(self): def ccp(self): """ - Method to average the amplitudes at each grid point to produce 2D images - for each of the three phases. At the end of this step, the object - contains the three 2D arrays that can be further averaged into a single - final image. + Method to average the amplitudes at each grid point to produce 2D + images for each of the three phases. At the end of this step, the + object contains the three 2D arrays that can be further averaged + into a single final image. The following attributes are added to the object: @@ -442,7 +462,7 @@ def ccp(self): """ if not self.is_ready_for_ccp: - raise(Exception("CCPimage not ready for ccp")) + raise Exception("CCPimage not ready for ccp") xs_ps_avg = np.zeros((self.nz, self.nx)) xs_pps_avg = np.zeros((self.nz, self.nx)) @@ -477,12 +497,12 @@ def ccp(self): def gccp(self, wlen=15.): """ - Method to average the amplitudes at each grid point to produce 2D images - for each of the three phases. In this method, the grid points are further - smoothed in the horizontal direction using a Gaussian function to simulate - P-wave sensitivity kernels. At the end of this step, the object - contains the three 2D smoothed arrays that can be further averaged into a - single final image. + Method to average the amplitudes at each grid point to produce 2D + images for each of the three phases. In this method, the grid points + are further smoothed in the horizontal direction using a Gaussian + function to simulate P-wave sensitivity kernels. At the end of this + step, the object contains the three 2D smoothed arrays that can be + further averaged into a single final image. Parameters ---------- @@ -494,16 +514,19 @@ def gccp(self, wlen=15.): Other Parameters ---------------- xs_gauss_ps : :class:`numpy.ndarray` - 2D array of stacked and Gaussian-filtered amplitudes for the Ps phase + 2D array of stacked and Gaussian-filtered amplitudes for the + Ps phase xs_gauss_pps : :class:`numpy.ndarray` - 2D array of stacked and Gaussian-filtered amplitudes for the Pps phase + 2D array of stacked and Gaussian-filtered amplitudes for the + Pps phase xs_gauss_pss : :class:`numpy.ndarray` - 2D array of stacked and Gaussian-filtered amplitudes for the Pss phase + 2D array of stacked and Gaussian-filtered amplitudes for the + Pss phase """ if not self.is_ready_for_gccp: - raise(Exception("CCPimage not ready for gccp")) + raise Exception("CCPimage not ready for gccp") if not hasattr(self, 'xs_ps_avg'): self.ccp() @@ -533,19 +556,20 @@ def linear_stack(self, typ='ccp'): Other Parameters ---------------- tot_trace : :class:`numpy.ndarray` - 2D array of amplitudes for the linearly combined Ps, Pps and Pss phases + 2D array of amplitudes for the linearly combined + Ps, Pps and Pss phases """ tot_trace = np.zeros((self.nz, self.nx)) - if typ=='ccp': + if typ == 'ccp': if not hasattr(self, "xs_ps_avg"): self.ccp() xs_ps = self.xs_ps_avg*self.weights[0] xs_pps = self.xs_pps_avg*self.weights[1] xs_pss = self.xs_pss_avg*self.weights[2] - elif typ=='gccp': + elif typ == 'gccp': if not hasattr(self, 'xs_gauss_ps'): self.gccp() xs_ps = self.xs_gauss_ps*self.weights[0] @@ -562,7 +586,7 @@ def linear_stack(self, typ='ccp'): def phase_weighted_stack(self, typ='gccp'): """ - Method to average the three 2D smoothed images into a final, + Method to average the three 2D smoothed images into a final, phase-weighted CCP image. Parameters @@ -575,20 +599,20 @@ def phase_weighted_stack(self, typ='gccp'): Other Parameters ---------------- tot_trace : :class:`numpy.ndarray` - 2D array of amplitudes for the phase-weighted, combined + 2D array of amplitudes for the phase-weighted, combined Ps, Pps and Pss phases """ tot_trace = np.zeros((self.nz, self.nx)) - if typ=='ccp': + if typ == 'ccp': if not hasattr(self, "xs_ps_avg"): self.ccp() xs_ps = self.xs_ps_avg*self.weights[0] xs_pps = self.xs_pps_avg*self.weights[1] xs_pss = self.xs_pss_avg*self.weights[2] - elif typ=='gccp': + elif typ == 'gccp': if not hasattr(self, 'xs_gauss_ps'): self.gccp() xs_ps = self.xs_gauss_ps*self.weights[0] @@ -619,7 +643,6 @@ def phase_weighted_stack(self, typ='gccp'): self.tot_trace = tot_trace - def save(self, title): """ Method to save the `ccpimage` object to file. @@ -635,16 +658,15 @@ def save(self, title): if title is None: title = "CCP_image.pkl" - if not ".pkl" in title: + if ".pkl" not in title: file = open(title+".pkl", "wb") else: file = open(title, "wb") pickle.dump(self, file) file.close() - - def plot_ccp(self, vmin=-0.05, vmax=0.05, - save=False, title='', fmt='png'): + def plot_ccp(self, vmin=-0.05, vmax=0.05, + save=False, title='', fmt='png'): """ Method to plot the final CCP stacks along the line. @@ -669,10 +691,13 @@ def plot_ccp(self, vmin=-0.05, vmax=0.05, 4, 1, figsize=(xm, 8)) # plt.pcolormesh(xarray,zarray,xs_ps_avg,cmap=cm.coolwarm,vmin=vmin,vmax=vmax) - im1 = ax1.pcolormesh(self.xarray, self.zarray, - self.xs_ps_avg*self.weights[0], - cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im1 = ax1.pcolormesh( + self.xarray, + self.zarray, + self.xs_ps_avg*self.weights[0], + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im1, ax=ax1) ax1.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -684,10 +709,13 @@ def plot_ccp(self, vmin=-0.05, vmax=0.05, ax1.invert_yaxis() # plt.pcolormesh(xarray,zarray,xs_pps_avg,cmap=cm.coolwarm,vmin=vmin,vmax=vmax) - im2 = ax2.pcolormesh(self.xarray, self.zarray, - self.xs_pps_avg*self.weights[1], - cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im2 = ax2.pcolormesh( + self.xarray, + self.zarray, + self.xs_pps_avg*self.weights[1], + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im2, ax=ax2) ax2.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -698,10 +726,13 @@ def plot_ccp(self, vmin=-0.05, vmax=0.05, ax2.set_title('Pps CCP image', size=10) ax2.invert_yaxis() - im3 = ax3.pcolormesh(self.xarray, self.zarray, - self.xs_pss_avg*self.weights[2], - cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im3 = ax3.pcolormesh( + self.xarray, + self.zarray, + self.xs_pss_avg*self.weights[2], + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im3, ax=ax3) ax3.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -712,9 +743,13 @@ def plot_ccp(self, vmin=-0.05, vmax=0.05, ax3.set_title('Pss CCP image', size=10) ax3.invert_yaxis() - im4 = ax4.pcolormesh(self.xarray, self.zarray, - self.tot_trace, cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im4 = ax4.pcolormesh( + self.xarray, + self.zarray, + self.tot_trace, + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im4, ax=ax4) ax4.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -734,7 +769,7 @@ def plot_ccp(self, vmin=-0.05, vmax=0.05, plt.show() def plot_gccp(self, vmin=-0.015, vmax=0.015, - save=False, title='', fmt='png'): + save=False, title='', fmt='png'): """ Method to plot the final GCCP stacks along the line. @@ -759,10 +794,13 @@ def plot_gccp(self, vmin=-0.015, vmax=0.015, 4, 1, figsize=(xm, 8)) # plt.pcolormesh(xarray,zarray,xs_ps_avg,cmap=cm.coolwarm,vmin=vmin,vmax=vmax) - im1 = ax1.pcolormesh(self.xarray, self.zarray, - self.xs_gauss_ps*self.weights[0], - cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im1 = ax1.pcolormesh( + self.xarray, + self.zarray, + self.xs_gauss_ps*self.weights[0], + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im1, ax=ax1) ax1.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -774,10 +812,13 @@ def plot_gccp(self, vmin=-0.015, vmax=0.015, ax1.invert_yaxis() # plt.pcolormesh(xarray,zarray,xs_pps_avg,cmap=cm.coolwarm,vmin=vmin,vmax=vmax) - im2 = ax2.pcolormesh(self.xarray, self.zarray, - self.xs_gauss_pps*self.weights[1], - cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im2 = ax2.pcolormesh( + self.xarray, + self.zarray, + self.xs_gauss_pps*self.weights[1], + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im2, ax=ax2) ax2.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -788,10 +829,13 @@ def plot_gccp(self, vmin=-0.015, vmax=0.015, ax2.set_title('Pps GCCP image', size=10) ax2.invert_yaxis() - im3 = ax3.pcolormesh(self.xarray, self.zarray, - self.xs_gauss_pss*self.weights[2], - cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im3 = ax3.pcolormesh( + self.xarray, + self.zarray, + self.xs_gauss_pss*self.weights[2], + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im3, ax=ax3) ax3.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -806,9 +850,13 @@ def plot_gccp(self, vmin=-0.015, vmax=0.015, self.tot_trace = ndimage.filters.gaussian_filter( self.tot_trace, sigma=(3, 0), order=0) - im4 = ax4.pcolormesh(self.xarray, self.zarray, - self.tot_trace, cmap=cm.RdBu_r, - vmin=vmin, vmax=vmax) + im4 = ax4.pcolormesh( + self.xarray, + self.zarray, + self.tot_trace, + cmap=cm.RdBu_r, + vmin=vmin, + vmax=vmax) bar = plt.colorbar(im4, ax=ax4) ax4.set_xlim((min(self.xarray)), (max(self.xarray))) @@ -904,10 +952,10 @@ def ttime(tr, delta_z, vp, vs, phase=None): # Calculate travel time for phase if phase == 'Ps': tt = delta_z*(np.sqrt((1./vs)**2 - slow**2) - - np.sqrt((1./vp)**2 - slow**2)) + np.sqrt((1./vp)**2 - slow**2)) elif phase == 'Pps': tt = delta_z*(np.sqrt((1./vs)**2 - slow**2) + - np.sqrt((1./vp)**2 - slow**2)) + np.sqrt((1./vp)**2 - slow**2)) elif phase == 'Pss': tt = 2.*delta_z*(np.sqrt((1./vs)**2 - slow**2)) else: diff --git a/rfpy/harmonics.py b/rfpy/harmonics.py index 4c7e3bc..866cc1f 100644 --- a/rfpy/harmonics.py +++ b/rfpy/harmonics.py @@ -91,10 +91,11 @@ def __init__(self, radialRF, transvRF=None, azim=0, xmin=0., xmax=10.): file.close() if not transvRF: - raise TypeError("__init__() missing 1 required positional argument: 'transvRF'") + raise TypeError("__init__() missing 1 required positional " + + "argument: 'transvRF'") # fftshift if the time axis starts at negative lags - if radialRF[0].stats.taxis[0]<0.: + if radialRF[0].stats.taxis[0] < 0.: for tr in radialRF: tr.data = np.fft.fftshift(tr.data) for tr in transvRF: @@ -257,8 +258,8 @@ def dcomp_fix_azim(self, azim=None): else: self.azim = azim - print('Decomposing receiver functions into baz harmonics for azimuth = ', - azim) + print("Decomposing receiver functions into baz harmonics for " + + "azimuth = ", azim) # Some integers nbin = len(self.radialRF) @@ -355,7 +356,7 @@ def forward(self, baz_list=None): """ if not hasattr(self, 'hstream'): - raise(Exception("Decomposition has not been performed yet")) + raise Exception("Decomposition has not been performed yet") if not baz_list: print("Warning: no BAZ specified - using all baz from " + @@ -415,7 +416,6 @@ def forward(self, baz_list=None): self.radial_forward.append(trR) self.transv_forward.append(trT) - def plot(self, ymax=30., scale=10., save=False, title=None, form='png'): """ Method to plot the 5 harmonic components. @@ -475,8 +475,11 @@ def plot(self, ymax=30., scale=10., save=False, title=None, form='png'): ax1.grid() if save: - plt.savefig('FIGURES/'+sta+'.'+title+'.'+form, dpi=300, - bbox_inches='tight', format=form) + plt.savefig( + 'FIGURES/'+sta+'.'+title+'.'+form, + dpi=300, + bbox_inches='tight', + format=form) plt.show() def save(self, file): diff --git a/rfpy/hk.py b/rfpy/hk.py index 5cbf6a5..c9f1910 100644 --- a/rfpy/hk.py +++ b/rfpy/hk.py @@ -24,7 +24,8 @@ Functions to calculate thickness (H) and Vp/Vs ratio (R) of the crust based on moveout times of direct Ps and reverberated Pps and Pss phases. -The stacks are obtained from the median weighted by the phase of individual signals. +The stacks are obtained from the median weighted by the phase of individual +signals. """ @@ -39,13 +40,13 @@ class HkStack(object): """ A HkStack object contains attributes and methods to stack radial - receiver functions along moveout curves for point measurements + receiver functions along moveout curves for point measurements of the depth to Moho (H) and P-to-S velocity ratio (k) beneath a seismic stations. The object is initialized with at least one :class:`~obspy.core.Stream` containing observed (or synthetised) radial receiver functions. The methods available can produce linear weighted stacks or product stacks, and can be used in the presence - of flat or dipping Moho (with known strike and dip). + of flat or dipping Moho (with known strike and dip). Note ---- @@ -59,7 +60,7 @@ class HkStack(object): ---------- rfV1 : :class:`~obspy.core.Stream` Stream object containing the radial-component receiver function - seismograms + seismograms rfV2 : :class:`~obspy.core.Stream` Stream object containing the radial-component receiver function seismograms (typically filtered at lower frequencies) @@ -77,13 +78,15 @@ class HkStack(object): dk : float Spacing between adjacent Vp/Vs search values hbound : list - List of 2 floats that determine the range of Moho depth values to search + List of 2 floats that determine the range of Moho depth values + to search dh : float Spacing between adjacent Moho depth search values weights : list List of 3 floats that determine the relative weights of the individual - phase stacks to be used in the final stack. The third weight is negative - since Pss amplitudes are opposite to those of the Ps and Pps phases. + phase stacks to be used in the final stack. The third weight is + negative since Pss amplitudes are opposite to those of the Ps and + Pps phases. phases : list List of 3 strings ('ps', 'pps', 'pss') corresponding to the thre phases of interest (`do not modify this attribute`) @@ -113,7 +116,7 @@ def __init__(self, rfV1, rfV2=None, strike=None, dip=None, vp=6.0): for tr in rfV2: tr.data = np.fft.fftshift(tr.data)[0:int(nn/2)] tr.stats.taxis = np.fft.fftshift(tr.data)[0:int(nn/2)] - except: + except Exception: pass self.rfV1 = rfV1 @@ -132,7 +135,7 @@ def stack(self, vp=None): """ Method to calculate Hk stacks from radial receiver functions. The stacks are calculated using phase-weighted stacking for - individual phases and take the median of the weighted stack + individual phases and take the median of the weighted stack to avoid bias by outliers. Note @@ -140,14 +143,14 @@ def stack(self, vp=None): If two streams are available as attributes, the method will assume that the second stream should be used for stacking along the Pps and Pss move out curves (e.g., if the second stream contains - lower frequency signals). Furthermore, If the ``vp`` argument is - not specified, the method will use the - value set during initialization (``vp=6.0`` km/s) + lower frequency signals). Furthermore, If the ``vp`` argument is + not specified, the method will use the value set during + initialization (``vp=6.0`` km/s) Parameters ---------- vp : float - Mean crust P-wave velocity (km/s). + Mean crust P-wave velocity (km/s). Attributes ---------- @@ -164,7 +167,7 @@ def stack(self, vp=None): if not vp: try: vp = self.rfV1[0].stats.vp - except: + except Exception: vp = self.vp # Station name @@ -199,14 +202,6 @@ def stack(self, vp=None): weight += np.exp(1j*tphase[0]) amp[i] = trace[0] - # ### Attempt at speeding things up - # ind = (np.abs(rfV.stats.taxis - tt)).argmin() - # trace = rfV.copy() - # thilb = hilbert(trace.data) - # tphase = np.arctan2(thilb.imag, thilb.real) - # weight += np.exp(1j*tphase[ind]) - # amp[i] = trace.data[ind] - weight = abs(weight/len(self.rfV1))**4 sig[ih, ik, ip] = np.var(amp)*np.real(weight) pws[ih, ik, ip] = np.median(amp)*np.real(weight) @@ -219,7 +214,7 @@ def stack_dip(self, vp=None, strike=None, dip=None): Method to calculate Hk stacks from radial receiver functions using known stike and dip angles of the Moho. The stacks are calculated using phase-weighted stacking for - individual phases and take the median of the weighted stack + individual phases and take the median of the weighted stack to avoid bias by outliers. Note @@ -227,17 +222,18 @@ def stack_dip(self, vp=None, strike=None, dip=None): If two streams are available as attributes, the method will assume that the second stream should be used for stacking along the Pps and Pss move out curves (e.g., if the second stream contains - lower frequency signals). Furthermore, - If the arguments are not specified, the method will use the - values set during initialization (``vp=6.0`` km/s, + lower frequency signals). Furthermore, + If the arguments are not specified, the method will use the + values set during initialization (``vp=6.0`` km/s, ``strike=0.``, ``dip=0.``) Parameters ---------- vp : float - Mean crust P-wave velocity (km/s). + Mean crust P-wave velocity (km/s). strike : float - Strike angle of dipping Moho (has to be known or estimated a priori) + Strike angle of dipping Moho (has to be known or estimated + a priori) dip : float Dip angle of Moho (has to be known or estimated a priori) @@ -265,7 +261,7 @@ def stack_dip(self, vp=None, strike=None, dip=None): if not vp: try: vp = self.rfV1[0].stats.vp - except: + except Exception: vp = 6.0 sta = self.rfV1[0].stats.station @@ -310,20 +306,20 @@ def stack_dip(self, vp=None, strike=None, dip=None): def average(self, typ='sum', q=0.05, err_method='amp'): """ Method to combine the phase-weighted stacks to produce a final - stack, from which to estimate the H and k parameters and their + stack, from which to estimate the H and k parameters and their associated errors. Parameters ---------- typ : str How the phase-weigthed stacks should be combined to produce - a final stack. Available options are: weighted sum (``typ=sum``) + a final stack. Available options are: weighted sum (``typ=sum``) or product (``typ=product``). q : float Confidence level for the error estimate err_method : str How errors should be estimated. Options are ``err_method='amp'`` - to estimate errors from amplitude, or ``err_method='stats'`` to + to estimate errors from amplitude, or ``err_method='stats'`` to use a statistical F test from the residuals. """ @@ -336,11 +332,11 @@ def average(self, typ='sum', q=0.05, err_method='amp'): ps = self.pws[:, :, 0]*self.weights[0] try: pps = self.pws[:, :, 1]*self.weights[1] - except: + except Exception: pps = None try: pss = self.pws[:, :, 2]*self.weights[2] - except: + except Exception: pss = None # Get stacks @@ -359,7 +355,7 @@ def average(self, typ='sum', q=0.05, err_method='amp'): pss = 1. stack = ps*pps*pss else: - raise(Exception("'typ' must be either 'sum' or 'product'")) + raise Exception("'typ' must be either 'sum' or 'product'") self.typ = typ @@ -372,7 +368,7 @@ def average(self, typ='sum', q=0.05, err_method='amp'): try: self.error() - except: + except Exception: self.err_k0 = 0. self.err_h0 = 0. @@ -388,7 +384,7 @@ def error(self, q=0.05, err_method='amp'): Confidence level for the error estimate err_method : str How errors should be estimated. Options are ``err_method='amp'`` - to estimate errors from amplitude, or ``err_method='stats'`` to + to estimate errors from amplitude, or ``err_method='stats'`` to use a statistical F test from the residuals. """ @@ -431,7 +427,7 @@ def error(self, q=0.05, err_method='amp'): err = np.where(msf > self.err_contour) else: - raise(Exception("'err_method' must be either 'stats' or 'amp'")) + raise Exception("'err_method' must be either 'stats' or 'amp'") self.err_method = err_method # Estimate uncertainty (q confidence interval) @@ -470,11 +466,11 @@ def plot(self, save=False, title=None, form='png'): ps[ps < 0] = 0. try: pps[pps < 0] = 0. - except: + except Exception: pass try: pss[pss < 0] = 0. - except: + except Exception: pass # Set up figure @@ -514,8 +510,8 @@ def plot(self, save=False, title=None, form='png'): ax4.set_title('Stack') ax4.set_xlabel('Thickness (km)') - #cbar = fig.colorbar(im, ticks=[-vmax, 0, vmax]) - #cbar.ax.set_yticklabels(['min', '0', 'max']) + # cbar = fig.colorbar(im, ticks=[-vmax, 0, vmax]) + # cbar.ax.set_yticklabels(['min', '0', 'max']) # Get confidence intervals if hasattr(self, 'err_contour'): @@ -536,7 +532,7 @@ def plot(self, save=False, title=None, form='png'): # Add star showing best fit try: ax4.scatter(self.h0, self.k0, 60, marker='*', color='white') - except: + except Exception: print("'h0' and 'k0' are not available") if title: @@ -552,9 +548,7 @@ def plot(self, save=False, title=None, form='png'): plt.close() -## JMG ## def save(self, file): - ## JMG ## """ Saves HkStack object to file @@ -571,7 +565,7 @@ def save(self, file): output.close() def _residuals(self): - """ + """ Internal method to obtain residuals between observed and predicted receiver functions given the Moho depth and Vp/Vs obtained from the Hk stack. diff --git a/rfpy/plotting.py b/rfpy/plotting.py index 4311152..716f770 100644 --- a/rfpy/plotting.py +++ b/rfpy/plotting.py @@ -71,7 +71,7 @@ def wiggle(stream1, stream2=None, sort=None, tmin=0., tmax=30, normalize=True, if sort is not None: try: stream1.traces.sort(key=lambda x: x.stats[sort], reverse=False) - except: + except Exception: print("Warning: stats attribute " + sort + " is not available in stats") pass @@ -188,11 +188,11 @@ def wiggle_bins(stream1, stream2=None, tr1=None, tr2=None, st2 = stream2.copy() if not (btyp == 'baz' or btyp == 'slow' or btyp == 'dist'): - raise(Exception("Type has to be 'baz' or 'slow' or 'dist'")) + raise Exception("Type has to be 'baz' or 'slow' or 'dist'") if not (xtyp == 'time' or xtyp == 'depth'): - raise(Exception("Type has to be 'time' or 'depth'")) + raise Exception("Type has to be 'time' or 'depth'") if btyp == 'slow' and xtyp == 'depth': - raise(Exception("Cannot plot by slowness if data is migrated")) + raise Exception("Cannot plot by slowness if data is migrated") # Figure out scaling here if norm is not None: @@ -456,60 +456,3 @@ def wiggle_single_event(rfdata, filt=None, pre_filt=None, trange=None): rfdata.meta.snr, rfdata.meta.cc)) plt.show() - - -def event_dist(stream, phase='P', save=False, title=None, form='png'): - - import cartopy.crs as ccrs - - # Specify azimuthal equidistant projection - aeqd = ccrs.AzimuthalEquidistant(central_longitude=stream[0].stats.stlo, - central_latitude=stream[0].stats.stla, - globe=ccrs.Globe()) - - # Extract latitude and longitude - evlat = [] - evlon = [] - phlist = [] - for tr in stream: - evlat.append(tr.stats.evla) - evlon.append(tr.stats.evlo) - phlist.append(tr.stats.phase) - - # Turn into arrays - evlat = np.array(evlat) - evlon = np.array(evlon) - phlist = np.array(phlist) - - # Select phase to plot - phP = phlist == 'P' - phPP = phlist == 'PP' - - # Now plot - fig = plt.figure(figsize=(3, 3)) - ax = fig.add_subplot(1, 1, 1, projection=aeqd) - if np.sum(phP) > 0: - # ax.scatter(evlon[phP], evlat[phP], c='royalblue', label='P') - ax.scatter(evlon[phP], evlat[phP], c='royalblue', label='P', - transform=ccrs.Geodetic()) - if np.sum(phPP) > 0: - ind = phase == 'P' - ax.scatter(evlon[phPP], evlat[phPP], c='coral', label='PP', - transform=ccrs.Geodetic()) - - ax.scatter(stream[0].stats.stlo, stream[0].stats.stla, c='grey', - marker='v', transform=ccrs.Geodetic()) - ax.coastlines() - - if title: - plt.suptitle(title) - - if save: - plt.savefig('RF_PLOTS/' + stream[0].stats.station + - '.' + title + '.event_dist.' + form, format=form) - else: - plt.show() - - plt.close() - - # ax.gridlines() diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 448ef30..bc6b71d 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -232,7 +232,7 @@ def add_event(self, event, gacmin=30., gacmax=90., phase='P', print(event.short_str()) if not isinstance(event, Event): - raise(Exception("Event has incorrect type")) + raise Exception("Event has incorrect type") # Store as object attributes self.meta = Meta(sta=self.sta, event=event, @@ -266,7 +266,7 @@ def add_data(self, stream, returned=False, new_sr=5.): """ if not self.meta: - raise(Exception("No meta data available - aborting")) + raise Exception("No meta data available - aborting") if not self.meta.accept: return @@ -287,7 +287,8 @@ def add_data(self, stream, returned=False, new_sr=5.): self.data = stream if not np.allclose( - [tr.stats.npts for tr in stream[1:]], stream[0].stats.npts): + [tr.stats.npts for tr in stream[1:]], + stream[0].stats.npts): self.meta.accept = False # Filter Traces @@ -320,8 +321,8 @@ def download_data(self, client, new_sr=5., dts=120., dts : float Time duration (sec) remove_response : bool - Remove instrument response from seismogram and resitute to true ground - velocity (m/s) using obspy.core.trace.Trace.remove_response() + Remove instrument response from seismogram and resitute to true + ground velocity (m/s) returned : bool Whether or not to return the ``accept`` attribute verbose : bool @@ -356,8 +357,13 @@ def download_data(self, client, new_sr=5., dts=120., # Download data err, stream = utils.download_data( - client=client, sta=self.sta, start=tstart, end=tend, - new_sr=new_sr, verbose=verbose, remove_response=remove_response, + client=client, + sta=self.sta, + start=tstart, + end=tend, + new_sr=new_sr, + verbose=verbose, + remove_response=remove_response, zcomp=self.zcomp) # Store as attributes with traces in dictionary @@ -458,8 +464,6 @@ def rotate(self, vp=None, vs=None, align=None): Z, N, E = rotate2zne(trZ.data, 0., -90., trN.data, azim, 0., trE.data, azim+90., 0.) - # Z, N, E = rotate2zne(trZ.data, 0., -90., trN.data, - # azim, 0., trE.data, azim+90., 0.) trN.data = N trE.data = E @@ -632,9 +636,9 @@ def deconvolve(self, phase='P', vp=None, vs=None, 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'. + 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 ---------- @@ -646,11 +650,11 @@ def deconvolve(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' or 'multitaper' wavelet : str - Type of wavelet for deconvolution. Options are 'complete', 'time' or - 'envelope' + Type of wavelet for deconvolution. Options are 'complete', 'time' + or 'envelope' envelope_threshold : float Threshold [0-1] used in ``wavelet='envelope'``. time : float @@ -660,7 +664,7 @@ def deconvolve(self, phase='P', vp=None, vs=None, Low and High frequency corners of bandpass filter applied before deconvolution gfilt : float - Center frequency of Gaussian filter (Hz). + Center frequency of Gaussian filter (Hz). wlevel : float Water level used in ``method='water'``. writeto : str or None @@ -695,33 +699,39 @@ def _gauss_filt(dt, nft, f0): return gauss def _Pwavelet(parent, method='complete', overhang=5, - envelope_threshold=0.05, time=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 + Select wavelet from the parent function for deconvolution using + method. + + Parameters + ---------- + parent : obspy.Trace + Wavefrom to extract the wavelet from + method : str + 'complete' uses complete parent signal after P arrival + (current implementation) + 'envelope' uses 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 + 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 + envelope_threshold : float + Fraction of the envelope that defines wavelet (for method='envelope') - time: float - window (seconds) that defines the wavelet (for method='time') + 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 + ------- + left, right : (obspy.UTCDateTime) + Start and end time of wavelet """ import obspy.signal.filter as osf @@ -734,8 +744,8 @@ def _Pwavelet(parent, method='complete', overhang=5, if method == 'envelope': split = int((self.meta.time + self.meta.ttime + - time - parent.stats.starttime ) * - parent.stats.sampling_rate) + 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 @@ -743,7 +753,8 @@ def _Pwavelet(parent, method='complete', overhang=5, i = np.nonzero( np.diff( np.array( - env > envelope_threshold, dtype=int))==-1)[0][0] + env > envelope_threshold, + dtype=int)) == -1)[0][0] except IndexError: i = len(parent.data)-1 dts = i * parent.stats.delta + time @@ -766,14 +777,14 @@ def _Pwavelet(parent, method='complete', overhang=5, 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 @@ -861,9 +872,9 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): 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) + # 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() @@ -885,7 +896,7 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): 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): + if self.meta.snr is None or not np.isfinite(self.meta.snr): print("Warning: SNR has not been calculated - " + "calculating now using default") self.calc_snr() @@ -914,37 +925,53 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): if phase == 'P' or 'PP': - # Get signal length (i.e., seismogram to deconvolve) from trace length + # 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) + 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, + 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) + 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]] + [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) + noise_left, noise_right = _Pwavelet( + trQ, + method='noise', + overhang=over) # Trim noise traces - [tr.trim(noise_left, noise_right, nearest_sample=False, - pad=True, fill_value=0.) for tr in [trNl, trNq]] + [tr.trim(noise_left, noise_right, nearest_sample=False, + pad=True, fill_value=0.) for tr in [trNl, trNq]] elif phase == 'S' or 'SKS': - # Get signal length (i.e., seismogram to deconvolve) from trace length + # 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) @@ -954,7 +981,7 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): 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.) @@ -963,13 +990,13 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): # Taper traces - only necessary processing after trimming # TODO: What does this to the multitaper method - [tr.taper(max_percentage=0.05, max_length=2.) + [tr.taper(max_percentage=0.05, max_length=2.) for tr in [trL, trQ, trT, trNl, trNq]] # Pre-filter waveforms before deconvolution if pre_filt: [tr.filter('bandpass', freqmin=pre_filt[0], freqmax=pre_filt[1], - corners=2, zerophase=True) + corners=2, zerophase=True) for tr in [trL, trQ, trT, trNl, trNq]] if writeto: @@ -990,14 +1017,13 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): self.rf = Stream(traces=[rfL, rfQ, rfT]) - def calc_cc(self): if not self.meta.accept: return if not hasattr(self, 'rf'): - raise(Exception("Warning: Receiver functions are not available")) + raise Exception("Warning: Receiver functions are not available") obs_L = self.data[0].copy() obs_Q = self.data[1].copy() @@ -1005,11 +1031,30 @@ def calc_cc(self): sr = obs_L.stats.sampling_rate # Filter using SNR bandpass - obs_L.detrend().taper(max_percentage=0.05, max_length=2.) - obs_Q.detrend().taper(max_percentage=0.05, max_length=2.) - obs_L.filter('bandpass', freqmin=0.05, freqmax=1., corners=2, zerophase=True) - obs_Q.filter('bandpass', freqmin=0.05, freqmax=1., corners=2, zerophase=True) - obs_rfQ.filter('bandpass', freqmin=0.05, freqmax=1., corners=2, zerophase=True) + obs_L.detrend().taper( + max_percentage=0.05, + max_length=2.) + obs_Q.detrend().taper( + max_percentage=0.05, + max_length=2.) + obs_L.filter( + 'bandpass', + freqmin=0.05, + freqmax=1., + corners=2, + zerophase=True) + obs_Q.filter( + 'bandpass', + freqmin=0.05, + freqmax=1., + corners=2, + zerophase=True) + obs_rfQ.filter( + 'bandpass', + freqmin=0.05, + freqmax=1., + corners=2, + zerophase=True) # Convolve L with rfQ to obtain predicted Q pred_Q = obs_Q.copy() @@ -1027,7 +1072,6 @@ def calc_cc(self): # Get cross correlation coefficient between observed and predicted Q self.meta.cc = np.corrcoef(obs_Q.data, pred_Q.data)[0][1] - def to_stream(self): """ Method to switch from RFData object to Stream object. @@ -1060,7 +1104,7 @@ def _add_rfstats(trace): return trace if not hasattr(self, 'rf'): - raise(Exception("Warning: Receiver functions are not available")) + raise Exception("Warning: Receiver functions are not available") stream = self.rf for tr in stream: diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 0d73f97..1aa0219 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -41,10 +41,6 @@ def get_calc_arguments(argv=None): - """ - Get Options from :class:`~optparse.OptionParser` objects. - - """ parser = ArgumentParser( usage="%(prog)s [arguments] ", @@ -74,7 +70,7 @@ def get_calc_arguments(argv=None): "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", + "-V", "--verbose", action="store_true", dest="verb", default=False, @@ -131,11 +127,11 @@ def get_calc_arguments(argv=None): "waveform server (--user-auth='username:authpassword') to access " + "and download restricted data. [Default no user and password]") ServerGroup.add_argument( - "--eida-token", - action="store", + "--eida-token", + action="store", type=str, - dest="tokenfile", - default=None, + dest="tokenfile", + default=None, help="Token for EIDA authentication mechanism, see " + "http://geofon.gfz-potsdam.de/waveform/archive/auth/index.php. " "If a token is provided, argument --user-auth will be ignored. " @@ -429,12 +425,6 @@ def get_calc_arguments(argv=None): else: args.userauth = [None, None] - # # Parse Local Data directories - # if args.localdata is not None: - # args.localdata = args.localdata.split(',') - # else: - # args.localdata = [] - # Check Datatype specification if args.dtype.upper() not in ['MSEED', 'SAC']: parser.error( @@ -679,28 +669,6 @@ def main(): " 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) diff --git a/rfpy/scripts/rfpy_ccp.py b/rfpy/scripts/rfpy_ccp.py index d367684..ff28722 100644 --- a/rfpy/scripts/rfpy_ccp.py +++ b/rfpy/scripts/rfpy_ccp.py @@ -37,12 +37,6 @@ def get_ccp_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] ", @@ -86,8 +80,7 @@ def get_ccp_arguments(argv=None): 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." - ) + "names, regardless of the key type of the database.") LineGroup = parser.add_argument_group( title='Line Geometry Settings', @@ -469,9 +462,10 @@ def main(): print() # Initialize CCPimage object - ccpimage = CCPimage(coord_start=args.coord_start, - coord_end=args.coord_end, - dz=args.dz, dx=args.dx) + ccpimage = CCPimage( + coord_start=args.coord_start, + coord_end=args.coord_end, + dz=args.dz, dx=args.dx) # Loop over station keys for stkey in list(stkeys): @@ -616,9 +610,12 @@ def main(): ccpfile = open(load_file, "rb") ccpimage = pickle.load(ccpfile) ccpfile.close() - ccpimage.prep_data(f1=args.f1, f2ps=args.f2ps, - f2pps=args.f2pps, f2pss=args.f2pss, - nbaz=args.nbaz, nslow=args.nslow) + ccpimage.prep_data( + f1=args.f1, + f2ps=args.f2ps, + f2pps=args.f2pps, + f2pss=args.f2pss, + nbaz=args.nbaz, nslow=args.nslow) ccpimage.is_ready_for_prestack = True ccpimage.save(prep_file) print() @@ -666,7 +663,7 @@ def main(): else: prestack_file = Path("CCP_prestack.pkl") if not prestack_file.is_file(): - raise(Exception("No CCP_prestack.pkl file available - aborting")) + raise Exception("No CCP_prestack.pkl file available - aborting") else: if args.linear: print() @@ -698,8 +695,12 @@ def main(): print("CCPimage saved to {0}".format(str(ccp_file))) if args.ccp_figure: - ccpimage.plot_ccp(save=args.save_figure, fmt=args.fmt, - vmin=-1.*args.cbound, vmax=args.cbound, title=args.title) + ccpimage.plot_ccp( + save=args.save_figure, + fmt=args.fmt, + vmin=-1.*args.cbound, + vmax=args.cbound, + title=args.title) else: pass @@ -714,7 +715,7 @@ def main(): else: prestack_file = Path("CCP_prestack.pkl") if not prestack_file.is_file(): - raise(Exception("No CCP_prestack.pkl file available - aborting")) + raise Exception("No CCP_prestack.pkl file available - aborting") else: if args.linear: print() @@ -746,8 +747,12 @@ def main(): print("CCPimage saved to {0}".format(str(gccp_file))) if args.ccp_figure: - ccpimage.plot_gccp(save=args.save_figure, fmt=args.fmt, - vmin=-1.*args.cbound, vmax=args.cbound, title=args.title) + ccpimage.plot_gccp( + save=args.save_figure, + fmt=args.fmt, + vmin=-1.*args.cbound, + vmax=args.cbound, + title=args.title) else: pass diff --git a/rfpy/scripts/rfpy_harmonics.py b/rfpy/scripts/rfpy_harmonics.py index 98bbafe..ae07f8f 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -38,12 +38,6 @@ def get_harmonics_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] ", @@ -68,7 +62,7 @@ def get_harmonics_arguments(argv=None): "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", + "-V", "--verbose", action="store_true", dest="verb", default=False, @@ -87,8 +81,7 @@ def get_harmonics_arguments(argv=None): 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." - ) + "names, regardless of the key type of the database.") # Event Selection Criteria TimeGroup = parser.add_argument_group( @@ -273,7 +266,7 @@ def get_harmonics_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Cannot construct UTCDateTime from start time: " + args.startT) @@ -284,7 +277,7 @@ def get_harmonics_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Cannot construct UTCDateTime from end time: " + args.endT) @@ -338,8 +331,7 @@ def main(): print("# __ _ _ #") print("# _ __ / _|_ __ _ _ | |__ __ _ _ __ _ __ ___ ___ _ __ (_) ___ ___ #") print("# | '__| |_| '_ \| | | | | '_ \ / _` | '__| '_ ` _ \ / _ \| '_ \| |/ __/ __| #") - print( - "# | | | _| |_) | |_| | | | | | (_| | | | | | | | | (_) | | | | | (__\__ \ #") + print("# | | | _| |_) | |_| | | | | | (_| | | | | | | | | (_) | | | | | (__\__ \ #") print("# |_| |_| | .__/ \__, |___|_| |_|\__,_|_| |_| |_| |_|\___/|_| |_|_|\___|___/ #") print("# |_| |___/_____| #") print("# #") @@ -510,7 +502,7 @@ def main(): harmonics.dcomp_fix_azim(azim=args.azim) if args.save_plot and not Path('FIGURES').is_dir(): - Path('FIGURES').mkdir(parents=True) + Path('FIGURES').mkdir(parents=True) if args.plot: harmonics.plot(args.ymax, args.scale, diff --git a/rfpy/scripts/rfpy_hk.py b/rfpy/scripts/rfpy_hk.py index 66a2adc..e975836 100644 --- a/rfpy/scripts/rfpy_hk.py +++ b/rfpy/scripts/rfpy_hk.py @@ -38,12 +38,6 @@ def get_hk_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] ", @@ -360,7 +354,7 @@ def get_hk_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Cannot construct UTCDateTime from start time: " + args.startT) @@ -371,7 +365,7 @@ def get_hk_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Cannot construct UTCDateTime from end time: " + args.endT) @@ -397,7 +391,7 @@ def get_hk_arguments(argv=None): "Error: --bp should contain 2 " + "comma-separated floats") -## JMG ## + # JMG # if args.slowbound is None: args.slowbound = [0.04, 0.08] else: @@ -417,7 +411,7 @@ def get_hk_arguments(argv=None): parser.error( "Error: --bazbound should contain 2 " + "comma-separated floats") -## JMG ## + # JMG # if args.phase not in ['P', 'PP', 'allP', 'S', 'SKS', 'allS']: parser.error( @@ -700,7 +694,7 @@ def main(): try: hkstack = HkStack(rfRstream, rfV2=rfRstream_copy, strike=args.strike, dip=args.dip, vp=args.vp) - except: + except Exception: hkstack = HkStack(rfRstream, strike=args.strike, dip=args.dip, vp=args.vp) diff --git a/rfpy/scripts/rfpy_recalc.py b/rfpy/scripts/rfpy_recalc.py index 73b58bd..4aff0ac 100644 --- a/rfpy/scripts/rfpy_recalc.py +++ b/rfpy/scripts/rfpy_recalc.py @@ -35,10 +35,6 @@ def get_recalc_arguments(argv=None): - """ - Get Options from :class:`~optparse.OptionParser` objects. - - """ parser = ArgumentParser( usage="%(prog)s [arguments] ", @@ -67,7 +63,7 @@ def get_recalc_arguments(argv=None): "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", + "-V", "--verbose", action="store_true", dest="verb", default=False, @@ -342,8 +338,7 @@ def main(): # Remove rotated flag and snr flag rfdata.meta.rotated = False rfdata.rotate(align='ZNE') - except: - + except Exception: print("Z12_Data.pkl not available - using ZNE_Data.pkl") # Load ZNE data ZNEfile = folder / "ZNE_Data.pkl" diff --git a/rfpy/utils.py b/rfpy/utils.py index 38ee25c..d7b3731 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -1,9 +1,30 @@ +# 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. + +# -*- coding: utf-8 -*- import math -from obspy import UTCDateTime -from numpy import nan, isnan, abs +from obspy import UTCDateTime, Stream, read import numpy as np import copy -from obspy import Stream, read def floor_decimal(n, decimals=0): @@ -83,10 +104,6 @@ def download_data(client=None, sta=None, start=UTCDateTime(), """ - # # Output - # print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, - # sta.channel.upper()))) - for loc in sta.location: # Construct location name From e16ab5fe1884dd68b76c3df805c9f978482f9964 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 13 Jun 2025 09:15:13 +1200 Subject: [PATCH 18/22] added stationXML support --- docs/rfpy.rst | 4 +-- docs/scripts.rst | 5 ++- rfpy/scripts/rfpy_calc.py | 28 +++++++++++----- rfpy/scripts/rfpy_ccp.py | 22 ++++++++++--- rfpy/scripts/rfpy_harmonics.py | 20 +++++++++--- rfpy/scripts/rfpy_hk.py | 22 ++++++++++--- rfpy/scripts/rfpy_plot.py | 22 ++++++++++--- rfpy/scripts/rfpy_recalc.py | 26 +++++++++++---- rfpy/utils.py | 60 ++++++++++++++++++++++++++++++++++ 9 files changed, 170 insertions(+), 39 deletions(-) diff --git a/docs/rfpy.rst b/docs/rfpy.rst index 5011b6f..c881c37 100644 --- a/docs/rfpy.rst +++ b/docs/rfpy.rst @@ -58,7 +58,7 @@ Install the `StDb` dependency using ``pip`` inside the ``rfpy`` environment: .. sourcecode:: bash - pip install stdb + pip install stdb@git+https://github.com/schaefferaj/stdb Installing from GitHub development branch ----------------------------------------- @@ -86,7 +86,7 @@ Installing from source Using local data ================ -The main script packaged with ``RfPy`` uses FDSN web services through and ``ObsPy`` `Client` to load waveform data. For waveform data locally stored on your hard drive, the scripts can use a `Client` that reads a `SeisComP Data Structure `_ archive containing SAC or miniSEED waveform data. Check out the scripts ``rfpy_calc`` below and the argument ``--local-data`` and ``--dtype`` for more details. +The main script packaged with ``RfPy`` uses FDSN web services through and ``ObsPy`` `Client` to load waveform data. For waveform data locally stored on your hard drive, the scripts can use a `Client` that reads a `SeisComP Data Structure `_ archive containing SAC or miniSEED waveform data. Check out the scripts ``rfpy_calc`` below and the argument ``--SDS-path`` and ``--dtype`` for more details. Station Metadata ---------------- diff --git a/docs/scripts.rst b/docs/scripts.rst index 6d741fa..ce29df7 100644 --- a/docs/scripts.rst +++ b/docs/scripts.rst @@ -85,10 +85,9 @@ Usage Settings associated with defining and using a local data base of pre-downloaded day-long SAC or MSEED files. - --local-data LOCALDATA + --SDS-path LOCALDATA Specify absolute path to a SeisComP Data Structure (SDS) archive - containing day-long SAC or MSEED files(e.g., --local- - data=/Home/username/Data/SDS). See + containing day-long SAC or MSEED files(e.g., --SDS-path=/Home/username/Data/SDS). See https://www.seiscomp.de/seiscomp3/doc/applications/slarchive/SDS.html for details on the SDS format. If this option is used, it takes precedence over the --server settings. diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 0b2a9a8..827db0f 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -30,7 +30,7 @@ import copy from obspy.clients.fdsn import Client as FDSN_Client from obspy.clients.filesystem.sds import Client as SDS_Client -from obspy import Catalog, UTCDateTime +from obspy import Catalog, UTCDateTime, read_inventory from http.client import IncompleteRead from rfpy import utils, RFData from pathlib import Path @@ -146,14 +146,14 @@ def get_calc_arguments(argv=None): "and using a local data base of pre-downloaded " + "day-long SAC or MSEED files.") DataGroup.add_argument( - "--local-data", + "--SDS-path", action="store", type=str, - dest="localdata", + dest="sds_path", default=None, help="Specify absolute path to a SeisComP Data Structure (SDS) " + "archive containing day-long SAC or MSEED files" + - "(e.g., --local-data=/Home/username/Data/SDS). " + + "(e.g., --SDS-path=/Home/username/Data/SDS). " + "See https://www.seiscomp.de/seiscomp3/doc/applications/slarchive/SDS.html " + "for details on the SDS format. If this option is used, it takes " + "precedence over the --server settings.") @@ -503,8 +503,20 @@ def main(): # Run Input Parser args = get_calc_arguments() - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml']: + print( + "Error: Must supply a station list in .pkl or .xml format ") + exit() + + if ext == 'pkl': + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + elif ext == 'xml': + inv = read_inventory(args.indb) + db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) # Loop over station keys for stkey in list(stkeys): @@ -527,7 +539,7 @@ def main(): datapath.mkdir(parents=True) # Establish client - if args.localdata is None: + if args.sds_path is None: data_client = FDSN_Client( base_url=args.server, user=args.userauth[0], @@ -535,7 +547,7 @@ def main(): eida_token=args.tokenfile) else: data_client = SDS_Client( - args.localdata, + args.sds_path, format=args.dtype) # Establish client for events diff --git a/rfpy/scripts/rfpy_ccp.py b/rfpy/scripts/rfpy_ccp.py index ff28722..23ede63 100644 --- a/rfpy/scripts/rfpy_ccp.py +++ b/rfpy/scripts/rfpy_ccp.py @@ -28,8 +28,8 @@ import pickle import stdb from obspy.clients.fdsn import Client -from obspy.core import Stream, UTCDateTime -from rfpy import binning, plotting, CCPimage +from obspy import Stream, UTCDateTime, read_inventory +from rfpy import binning, plotting, CCPimage, utils from pathlib import Path from argparse import ArgumentParser from os.path import exists as exist @@ -61,7 +61,7 @@ def get_ccp_arguments(argv=None): "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", + "-V", "--verbose", action="store_true", dest="verb", default=False, @@ -431,8 +431,20 @@ def main(): # Run Input Parser args = get_ccp_arguments() - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml']: + print( + "Error: Must supply a station list in .pkl or .xml format ") + exit() + + if ext == 'pkl': + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + elif ext == 'xml': + inv = read_inventory(args.indb) + db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) # Track processed folders procfold = [] diff --git a/rfpy/scripts/rfpy_harmonics.py b/rfpy/scripts/rfpy_harmonics.py index ae07f8f..c4ef5c7 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -29,8 +29,8 @@ import stdb import copy from obspy.clients.fdsn import Client -from obspy.core import Stream, UTCDateTime -from rfpy import binning, plotting, Harmonics +from obspy import Stream, UTCDateTime, read_inventory +from rfpy import binning, plotting, Harmonics, utils from pathlib import Path from argparse import ArgumentParser from os.path import exists as exist @@ -341,8 +341,20 @@ def main(): # Run Input Parser args = get_harmonics_arguments() - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml']: + print( + "Error: Must supply a station list in .pkl or .xml format ") + exit() + + if ext == 'pkl': + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + elif ext == 'xml': + inv = read_inventory(args.indb) + db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) # Track processed folders procfold = [] diff --git a/rfpy/scripts/rfpy_hk.py b/rfpy/scripts/rfpy_hk.py index e975836..05e2c82 100644 --- a/rfpy/scripts/rfpy_hk.py +++ b/rfpy/scripts/rfpy_hk.py @@ -29,8 +29,8 @@ import stdb import copy from obspy.clients.fdsn import Client -from obspy.core import Stream, UTCDateTime -from rfpy import binning, plotting, HkStack +from obspy import Stream, UTCDateTime, read_inventory +from rfpy import binning, plotting, HkStack, utils from pathlib import Path from argparse import ArgumentParser from os.path import exists as exist @@ -62,7 +62,7 @@ def get_hk_arguments(argv=None): "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", + "-V", "--verbose", action="store_true", dest="verb", default=False, @@ -488,8 +488,20 @@ def main(): # Run Input Parser args = get_hk_arguments() - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml']: + print( + "Error: Must supply a station list in .pkl or .xml format ") + exit() + + if ext == 'pkl': + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + elif ext == 'xml': + inv = read_inventory(args.indb) + db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) # Track processed folders procfold = [] diff --git a/rfpy/scripts/rfpy_plot.py b/rfpy/scripts/rfpy_plot.py index 6c69754..aa0503b 100644 --- a/rfpy/scripts/rfpy_plot.py +++ b/rfpy/scripts/rfpy_plot.py @@ -28,8 +28,8 @@ import pickle import stdb import copy -from obspy import Stream, UTCDateTime -from rfpy import binning, plotting +from obspy import Stream, UTCDateTime, read_inventory +from rfpy import binning, plotting, utils from pathlib import Path from argparse import ArgumentParser from os.path import exists as exist @@ -64,7 +64,7 @@ def get_plot_arguments(argv=None): "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", + "-V", "--verbose", action="store_true", dest="verb", default=False, @@ -342,8 +342,20 @@ def main(): # Run Input Parser args = get_plot_arguments() - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml']: + print( + "Error: Must supply a station list in .pkl or .xml format ") + exit() + + if ext == 'pkl': + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + elif ext == 'xml': + inv = read_inventory(args.indb) + db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) # Track processed folders procfold = [] diff --git a/rfpy/scripts/rfpy_recalc.py b/rfpy/scripts/rfpy_recalc.py index 23cb033..5de76bd 100644 --- a/rfpy/scripts/rfpy_recalc.py +++ b/rfpy/scripts/rfpy_recalc.py @@ -30,8 +30,9 @@ import pickle import stdb import copy -from rfpy import RFData +from rfpy import RFData, utils from pathlib import Path +from obspy import read_inventory def get_recalc_arguments(argv=None): @@ -220,11 +221,10 @@ def get_recalc_arguments(argv=None): "Error: Incorrect alignment specifier. Should be " + "either 'ZRT', 'LQT', or 'PVH'.") - if args.method not in ['wiener', 'water', 'multitaper',\ - 'wiener_audet_bssa2010']: + if args.method not in ['wiener', 'wiener-mod', 'water', 'multitaper']: parser.error( - "Error: 'method' should be either 'wiener', 'water', " + - "'multitaper', or 'wiener_audet_bssa2010'") + "Error: 'method' should be either 'wiener', 'wiener-mod', " + + "'water', or 'multitaper'") if args.pre_filt is not None: args.pre_filt = [float(val) for val in args.pre_filt.split(',')] @@ -255,8 +255,20 @@ def main(): # Run Input Parser args = get_recalc_arguments() - # Load Database - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml']: + print( + "Error: Must supply a station list in .pkl or .xml format ") + exit() + + if ext == 'pkl': + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) + + elif ext == 'xml': + inv = read_inventory(args.indb) + db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) # Track processed folders procfold = [] diff --git a/rfpy/utils.py b/rfpy/utils.py index d7b3731..671ea8b 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -25,6 +25,66 @@ from obspy import UTCDateTime, Stream, read import numpy as np import copy +import re +from stdb import StDbElement + + +def get_stkeys(inventory, keys=[]): + + allkeys = [] + station_list = inventory.get_contents()['stations'] + allkeys = [s.split(' ')[0] for s in station_list] + + if len(keys) > 0: + # Extract key subset + + stkeys = [] + for key in keys: + # Convert the pattern to a regex pattern + # Replace '.' with '\.' to match literal dots + # Replace '*' with '.*' to match any sequence of characters + # Replace '?' with '.' to match any single character + pattern = key.replace('.', r'\.').replace('*', '.*').replace('?', '.') + + # Compile the regex pattern + regex = re.compile(f'^.*{pattern}.*$') + + # Filter allkeys based on the compiled regex + stkeys.extend([key for key in allkeys if regex.match(key)]) + + else: + stkeys = allkeys + + return stkeys + + +def inv2stdb(inventory, keys=[]): + + stkeys = get_stkeys(inventory, keys) + + stations = {} + for key in stkeys: + net = key.split('.')[0] + sta = key.split('.')[1] + cha = '?H?' + inv = inventory.select(network=net, station=sta, channel=cha) + seed_id = inv.get_contents()['channels'][0] + coords = inv.get_coordinates(seed_id) + + stdb_element = StDbElement( + station=sta, + network=net, + channel=seed_id.split('.')[3][0:2], + location=seed_id.split('.')[2], + latitude=coords['latitude'], + longitude=coords['longitude'], + elevation=coords['elevation'], + startdate=inv[0].stations[0].start_date, + enddate=inv[0].stations[0].end_date + ) + stations[key] = stdb_element + + return stations, stkeys def floor_decimal(n, decimals=0): From 6b36315c3a290871e962d04457280ef32e0ed7fc Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 13 Jun 2025 09:47:05 +1200 Subject: [PATCH 19/22] fixed locs in inv2stdb --- rfpy/utils.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/rfpy/utils.py b/rfpy/utils.py index 671ea8b..2f64bd3 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -70,12 +70,16 @@ def inv2stdb(inventory, keys=[]): inv = inventory.select(network=net, station=sta, channel=cha) seed_id = inv.get_contents()['channels'][0] coords = inv.get_coordinates(seed_id) + channel = seed_id.split('.')[3][0:2] + location = list(seed_id.split('.')[2]) + if len(location) == 0: + location.append('--') stdb_element = StDbElement( station=sta, network=net, - channel=seed_id.split('.')[3][0:2], - location=seed_id.split('.')[2], + channel=channel, + location=list(set(location)), latitude=coords['latitude'], longitude=coords['longitude'], elevation=coords['elevation'], @@ -164,6 +168,7 @@ def download_data(client=None, sta=None, start=UTCDateTime(), """ + print(sta) for loc in sta.location: # Construct location name @@ -197,10 +202,13 @@ def download_data(client=None, sta=None, start=UTCDateTime(), if len(st) == 3: # It's possible if len(st)==1 that data is Z12 print("* - Data Downloaded") - # break + elif len(st) < 3: + print("* Error retrieving waveforms") + print("**************************************************") + return True, None # Check the correct 3 components exist - if st is None or len(st) < 3: + if st is None: print("* Error retrieving waveforms") print("**************************************************") return True, None From 5cd61975d7c1e2821dabf93a9482d31683cf4c3b Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 13 Jun 2025 09:53:25 +1200 Subject: [PATCH 20/22] cleanup of utils --- rfpy/utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/rfpy/utils.py b/rfpy/utils.py index 2f64bd3..0d8375c 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -168,7 +168,6 @@ def download_data(client=None, sta=None, start=UTCDateTime(), """ - print(sta) for loc in sta.location: # Construct location name @@ -201,7 +200,7 @@ def download_data(client=None, sta=None, start=UTCDateTime(), else: if len(st) == 3: # It's possible if len(st)==1 that data is Z12 - print("* - Data Downloaded") + print("* - Data Downloaded") elif len(st) < 3: print("* Error retrieving waveforms") print("**************************************************") From 992f11949c0235dd453ba805572f709d4f6ccea3 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Fri, 13 Jun 2025 10:04:23 +1200 Subject: [PATCH 21/22] fix docstring to add 'wiener-mod' --- rfpy/rfdata.py | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index 70e7adb..e69e470 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -256,11 +256,6 @@ def add_data(self, stream, returned=False, new_sr=5.): returned : bool Whether or not to return the ``accept`` attribute - Attributes - ---------- - data : :class:`~obspy.core.Stream` - Stream container for NEZ seismograms - Returns ------- accept : bool @@ -336,11 +331,6 @@ def download_data(self, client, new_sr=5., dts=120., accept : bool Whether or not the object is accepted for further analysis - Attributes - ---------- - data : :class:`~obspy.core.Stream` - Stream containing :class:`~obspy.core.Trace` objects - """ if self.meta is None: @@ -653,8 +643,8 @@ def deconvolve(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 - 'multitaper' + Method for deconvolution. Options are 'wiener', 'wiener-mod', + 'water', or 'multitaper' wavelet : str Type of wavelet for deconvolution. Options are 'complete', 'time' or 'envelope' From 8d74fdb96e52e50257aeeb4f9355508b0c2642a5 Mon Sep 17 00:00:00 2001 From: Pascal Audet Date: Wed, 9 Jul 2025 12:43:50 -0400 Subject: [PATCH 22/22] clean up docs and stationXML import from StDb --- docs/rfpy.rst | 15 ++++++- docs/scripts.rst | 12 ++--- rfpy/rfdata.py | 93 +++++++++++++++++++++++---------------- rfpy/scripts/rfpy_calc.py | 33 ++++++-------- rfpy/utils.py | 63 -------------------------- 5 files changed, 88 insertions(+), 128 deletions(-) diff --git a/docs/rfpy.rst b/docs/rfpy.rst index c881c37..73aa40a 100644 --- a/docs/rfpy.rst +++ b/docs/rfpy.rst @@ -91,9 +91,20 @@ The main script packaged with ``RfPy`` uses FDSN web services through and ``ObsP Station Metadata ---------------- -If you have data stored locally on your drive, it is likely you also have a station `XML `_ file containing the metadata. The corresponding ObsPy documentation is `here `_. +If you have data stored locally on your drive, it is likely you also +have a station `XML `_ file +containing the metadata. The corresponding ObsPy documentation is +`here `_. -To convert the station `XML` file to an input that can be read by ``OrientPy``, you run the command ``gen_stdb station.xml`` (only available on StDb version 0.2.7), which will create the file ``station.pkl``. If you don't have a station `XML` file but you have a dataless SEED file, you can convert it first to `XML` using `this tools `_. +You can now use a stationXML (`.xml`) file instead of the StDb `.pkl` format. +Alternatively, you can convert the stationXML file to an StDb `.pkl` file +by running the command ``gen_stdb station.xml`` (these options are only +available on StDb version 0.2.7. If you don't have a station `XML` file but you have +a dataless SEED file, you can convert it first to XML using `this tools `_. + +.. note:: + Please note that using the stationXML directly as input means you cannot + correct the orientation of H1 and H2 components using the azimuth correction term stored as ``azcorr`` in the StDb file, as this information is not stored in the stationXML file. Waveform Data ------------- diff --git a/docs/scripts.rst b/docs/scripts.rst index ce29df7..dcef120 100644 --- a/docs/scripts.rst +++ b/docs/scripts.rst @@ -44,7 +44,9 @@ Usage services framework. The stations are processed one by one and the data are stored to disk. positional arguments: - indb Station Database to process from. + indb Station Database to process from. Available + formats are: StDb (.pkl or .csv) or stationXML + (.xml) options: -h, --help show this help message and exit @@ -61,7 +63,7 @@ Usage regardless of the key type of the database. Server Settings: - Settings associated with which datacenter to log into. + Settings associated with FDSN datacenters for archived data. --server SERVER Base URL of FDSN web service compatible server (e.g. “http://service.iris.edu”) or key string for recognized server (one @@ -82,10 +84,10 @@ Usage PGP message in it. [Default None] Local Data Settings: - Settings associated with defining and using a local data base of pre-downloaded day-long - SAC or MSEED files. + Settings associated with a SeisComP database for locally archived + data. - --SDS-path LOCALDATA + --SDS-path SDS_PATH Specify absolute path to a SeisComP Data Structure (SDS) archive containing day-long SAC or MSEED files(e.g., --SDS-path=/Home/username/Data/SDS). See https://www.seiscomp.de/seiscomp3/doc/applications/slarchive/SDS.html diff --git a/rfpy/rfdata.py b/rfpy/rfdata.py index e69e470..931b0ed 100644 --- a/rfpy/rfdata.py +++ b/rfpy/rfdata.py @@ -625,8 +625,8 @@ def calc_snr(self, dt=30., fmin=0.05, fmax=1.): self.meta.snrh = 10*np.log10(srms*srms/nrms/nrms) def deconvolve(self, phase='P', vp=None, vs=None, - align=None, method='wiener', wavelet='complete', - envelope_threshold=0.05, time=5, pre_filt=None, + align=None, method='wiener', wavelet='envelope', + envelope_threshold=0.05, time=30, pre_filt=None, gfilt=None, wlevel=0.01, writeto=None): """ Deconvolves three-component data using one component as the source @@ -925,46 +925,63 @@ def _decon(parent, daughter1, daughter2, noise_parent, noise_daughter1, nn, meth # Get signal length (i.e., seismogram to deconvolve) from # trace length - over = 5 - dtsqt = len(trL.data)*trL.stats.delta/2. + # over = 5 + dts = 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) + nn = int(round((dts - 5.)*trL.stats.sampling_rate)) + 1 + + sig_left = self.meta.time + self.meta.ttime - 5. + sig_right = self.meta.time + self.meta.ttime + dts - 10. # Trim signal traces - [tr.trim(sig_left, sig_right, nearest_sample=False, - pad=True, fill_value=0.) for tr in [trQ, trT]] + [tr.trim(sig_left, sig_right, nearest_sample=False, + pad=nn, fill_value=0.) for tr in [trL, trQ, trT]] - # Noise window (-dtsqt to -5. sec) - noise_left, noise_right = _Pwavelet( - trQ, - method='noise', - overhang=over) + # Noise window (-dts to -5. sec) + noise_left = self.meta.time + self.meta.ttime - dts + noise_right = self.meta.time + self.meta.ttime - 5. # Trim noise traces - [tr.trim(noise_left, noise_right, nearest_sample=False, - pad=True, fill_value=0.) for tr in [trNl, trNq]] + [tr.trim(noise_left, noise_right, nearest_sample=False, + pad=nn, fill_value=0.) for tr in [trNl, trNq]] + + # dtsqt = len(trL.data)*trL.stats.delta/2. + # nn = int(round((dtsqt-5.)*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 traces + # [tr.trim(noise_left, noise_right, nearest_sample=False, + # pad=True, fill_value=0.) for tr in [trNl, trNq]] elif phase == 'S' or 'SKS': @@ -997,9 +1014,9 @@ def _decon(parent, daughter1, daughter2, noise_parent, noise_daughter1, nn, meth 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) + # if writeto: + # with open(writeto, 'wb') as f: + # pickle.dump(Stream(traces=[trL, trQ, trT]), f) # Deconvolve if phase == 'P' or 'PP': diff --git a/rfpy/scripts/rfpy_calc.py b/rfpy/scripts/rfpy_calc.py index 827db0f..6c73ff7 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -55,7 +55,8 @@ def get_calc_arguments(argv=None): # General Settings parser.add_argument( "indb", - help="Station Database to process from.", + help="Station Database to process from. Available formats are: " + + "StDb (.pkl or .csv) or stationXML (.xml)", type=str) parser.add_argument( "--keys", @@ -102,8 +103,7 @@ def get_calc_arguments(argv=None): # Server Settings ServerGroup = parser.add_argument_group( title="Server Settings", - description="Settings associated with which " - "datacenter to log into.") + description="Settings associated with FDSN datacenters for archived data.") ServerGroup.add_argument( "--server", action="store", @@ -142,9 +142,8 @@ def get_calc_arguments(argv=None): # 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.") + description="Settings associated with a SeisComP database " + + "for locally archived data.") DataGroup.add_argument( "--SDS-path", action="store", @@ -367,6 +366,12 @@ def get_calc_arguments(argv=None): if not exist(args.indb): parser.error("Input file " + args.indb + " does not exist") + # Check Extension + ext = args.indb.split('.')[-1] + + if ext not in ['pkl', 'xml', 'csv']: + parser.error("Must supply a station list in .pkl, .csv or .xml format ") + # create station key list if len(args.stkeys) > 0: args.stkeys = args.stkeys.split(',') @@ -503,20 +508,8 @@ def main(): # Run Input Parser args = get_calc_arguments() - # Check Extension - ext = args.indb.split('.')[-1] - - if ext not in ['pkl', 'xml']: - print( - "Error: Must supply a station list in .pkl or .xml format ") - exit() - - if ext == 'pkl': - db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) - - elif ext == 'xml': - inv = read_inventory(args.indb) - db, stkeys = utils.inv2stdb(inv, keys=args.stkeys) + # Get station StDb and keys + db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # Loop over station keys for stkey in list(stkeys): diff --git a/rfpy/utils.py b/rfpy/utils.py index 0d8375c..ae09dd3 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -25,72 +25,9 @@ from obspy import UTCDateTime, Stream, read import numpy as np import copy -import re from stdb import StDbElement -def get_stkeys(inventory, keys=[]): - - allkeys = [] - station_list = inventory.get_contents()['stations'] - allkeys = [s.split(' ')[0] for s in station_list] - - if len(keys) > 0: - # Extract key subset - - stkeys = [] - for key in keys: - # Convert the pattern to a regex pattern - # Replace '.' with '\.' to match literal dots - # Replace '*' with '.*' to match any sequence of characters - # Replace '?' with '.' to match any single character - pattern = key.replace('.', r'\.').replace('*', '.*').replace('?', '.') - - # Compile the regex pattern - regex = re.compile(f'^.*{pattern}.*$') - - # Filter allkeys based on the compiled regex - stkeys.extend([key for key in allkeys if regex.match(key)]) - - else: - stkeys = allkeys - - return stkeys - - -def inv2stdb(inventory, keys=[]): - - stkeys = get_stkeys(inventory, keys) - - stations = {} - for key in stkeys: - net = key.split('.')[0] - sta = key.split('.')[1] - cha = '?H?' - inv = inventory.select(network=net, station=sta, channel=cha) - seed_id = inv.get_contents()['channels'][0] - coords = inv.get_coordinates(seed_id) - channel = seed_id.split('.')[3][0:2] - location = list(seed_id.split('.')[2]) - if len(location) == 0: - location.append('--') - - stdb_element = StDbElement( - station=sta, - network=net, - channel=channel, - location=list(set(location)), - latitude=coords['latitude'], - longitude=coords['longitude'], - elevation=coords['elevation'], - startdate=inv[0].stations[0].start_date, - enddate=inv[0].stations[0].end_date - ) - stations[key] = stdb_element - - return stations, stkeys - - def floor_decimal(n, decimals=0): multiplier = 10 ** decimals return math.floor(n * multiplier) / multiplier