From ba5282339374dc840944a4287ec3b1b3c34c9b6e Mon Sep 17 00:00:00 2001 From: Andrew Schaeffer Date: Thu, 27 Mar 2025 17:21:20 -0700 Subject: [PATCH] A Schaeffer. Adapt OrientPy to be compatible with NFSI data. (1) - allow selectable vertical component identifier. NFSI Data is 123, not 12Z (2) - split the waveform and catalogue clients to separate services, and allow for base_urls to be specified Signed-off-by: Andrew Schaeffer 2025-03-27 17:20 PDT --- orientpy/classes.py | 73 ++++++++++++----------- orientpy/io.py | 96 +++++++++++++++++-------------- orientpy/plotting.py | 6 +- orientpy/scripts/bng_calc_auto.py | 64 ++++++++++++++++++--- orientpy/scripts/dl_average.py | 37 +++++++++--- orientpy/scripts/dl_calc.py | 69 +++++++++++++++++----- orientpy/utils.py | 4 +- 7 files changed, 239 insertions(+), 110 deletions(-) diff --git a/orientpy/classes.py b/orientpy/classes.py index e6de2f2..c9cfb4d 100644 --- a/orientpy/classes.py +++ b/orientpy/classes.py @@ -119,7 +119,7 @@ class Orient(object): """ - def __init__(self, sta): + def __init__(self, sta, zcomp='Z'): # Attributes from parameters self.sta = sta @@ -127,6 +127,7 @@ def __init__(self, sta): # Initialize meta and data objects as None self.meta = None self.data = None + self.zcomp = zcomp def add_event(self, event, gacmin=None, gacmax=None, depmax=1000., returned=False): @@ -193,6 +194,9 @@ def download_data(self, client, stdata=[], ndval=np.nan, new_sr=None, :class:`~obspy.core.UTCDateTime` object. returned : bool Whether or not to return the ``accept`` attribute + zcomp : str + Character representing the vertical component. Should be + a single character :str:. Default Z Returns ------- @@ -228,14 +232,19 @@ def download_data(self, client, stdata=[], ndval=np.nan, new_sr=None, # Download data err, stream = io.download_data( client=client, sta=self.sta, start=tstart, end=tend, - stdata=stdata, ndval=ndval, new_sr=new_sr, + stdata=stdata, ndval=ndval, new_sr=new_sr, zcomp=self.zcomp, verbose=verbose) + if err or stream is None: + print("* Waveform Request Failed...Skipping") + self.meta.accept = False + return self.meta.accept + # Store as attributes with traces in dictionary try: trE = stream.select(component='E')[0] trN = stream.select(component='N')[0] - trZ = stream.select(component='Z')[0] + trZ = stream.select(component=self.zcomp)[0] self.data = Stream(traces=[trZ, trN, trE]) # Filter Traces and resample @@ -246,26 +255,24 @@ def download_data(self, client, stdata=[], ndval=np.nan, new_sr=None, # If there is no ZNE, perhaps there is Z12? except: + tr1 = stream.select(component='1')[0] + tr2 = stream.select(component='2')[0] + trZ = stream.select(component=self.zcomp)[0] + self.data = Stream(traces=[trZ, tr1, tr2]) - try: - tr1 = stream.select(component='1')[0] - tr2 = stream.select(component='2')[0] - trZ = stream.select(component='Z')[0] - self.data = Stream(traces=[trZ, tr1, tr2]) - - # Filter Traces and resample - if new_sr: - self.data.filter('lowpass', freq=0.5*new_sr, - corners=2, zerophase=True) - self.data.resample(new_sr) - - if not np.sum([np.std(tr.data) for tr in self.data]) > 0.: - print('Error: Data are all zeros') - self.meta.accept = False + # Filter Traces and resample + if new_sr: + self.data.filter('lowpass', freq=0.5*new_sr, + corners=2, zerophase=True) + self.data.resample(new_sr) - except: + if not np.sum([np.std(tr.data) for tr in self.data]) > 0.: + print('Error: Data are all zeros') self.meta.accept = False + # except: + # self.meta.accept = False + if returned: return self.meta.accept @@ -311,9 +318,9 @@ class BNG(Orient): """ - def __init__(self, sta): + def __init__(self, sta, zcomp='Z'): - Orient.__init__(self, sta) + Orient.__init__(self, sta, zcomp=zcomp) def calc(self, dphi, dts, tt, bp=None, showplot=False): @@ -366,7 +373,8 @@ def calc(self, dphi, dts, tt, bp=None, showplot=False): test_sets = { 'ZNE':{'Z', 'N', 'E'}, 'Z12':{'Z', '1', '2'}, - 'Z23':{'Z', '2', '3'}, + 'Z23':{'Z', '2', '3'}, + '312':{'3', '1', '2'}, '123':{'1', '2', '3'} # probably should raise an exception if this is the case, } # as no correction is estimated for the vertical component for test_key in test_sets: @@ -374,13 +382,14 @@ def calc(self, dphi, dts, tt, bp=None, showplot=False): if test_set.issubset(set(comps_id)): # use sets to avoid sorting complications comps_codes = list(test_key) break - + #-- temporarily modify channel codes, assuming that N/E are not oriented properly channel_code_prefix = stream[0].stats.channel[:2] # prefix should be the same for all # 3 components by now stream.select(component=comps_codes[1])[0].stats.channel = channel_code_prefix + '1' stream.select(component=comps_codes[2])[0].stats.channel = channel_code_prefix + '2' - stream.select(component=comps_codes[0])[0].stats.channel = channel_code_prefix + 'Z' + stream.select(component=comps_codes[0])[0].stats.channel = channel_code_prefix + self.zcomp.upper() + # Filter if specified if bp: @@ -396,8 +405,8 @@ def calc(self, dphi, dts, tt, bp=None, showplot=False): # Define signal and noise tr1 = stdata.select(component='1')[0].copy() tr2 = stdata.select(component='2')[0].copy() - trZ = stdata.select(component='Z')[0].copy() - ntrZ = stnoise.select(component='Z')[0].copy() + trZ = stdata.select(component=self.zcomp.upper())[0].copy() + ntrZ = stnoise.select(component=self.zcomp.upper())[0].copy() # Calculate and store SNR as attribute self.meta.snr = 10.*np.log10( @@ -517,9 +526,9 @@ class DL(Orient): """ - def __init__(self, sta): + def __init__(self, sta, zcomp='Z'): - Orient.__init__(self, sta) + Orient.__init__(self, sta, zcomp=zcomp) def calc(self, showplot=False): @@ -609,15 +618,15 @@ def calc(self, showplot=False): # R1 path R1phi[k], R1cc[k] = utils.DLcalc( - stream.copy(), item[0], item[1], + stream, item[0], item[1], item[2], self.meta.epi_dist, self.meta.baz, Ray1, - winlen=item[3], ptype=0) + winlen=item[3], ptype=0, zcomp=self.zcomp) # R2 path R2phi[k], R2cc[k] = utils.DLcalc( - stream.copy(), item[0], item[1], + stream, item[0], item[1], item[2], dist2, baz2, Ray2, - winlen=item[4], ptype=0) + winlen=item[4], ptype=0, zcomp=self.zcomp) # Store azimuths and CC values as attributes self.meta.R1phi = R1phi diff --git a/orientpy/io.py b/orientpy/io.py index 6b10028..be31968 100644 --- a/orientpy/io.py +++ b/orientpy/io.py @@ -400,7 +400,8 @@ def parse_localdata_for_comp(comp='Z', stdata=[], sta=None, def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, - stdata=[], ndval=nan, new_sr=0., verbose=False): + stdata=[], ndval=nan, new_sr=0., verbose=False, zcomp='Z', + coordsys=2): """ 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 @@ -445,8 +446,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, from math import floor # Output - print(("* {0:s}.{1:2s} - ZNE:".format(sta.station, - sta.channel.upper()))) + print("* {0:s}.{1:2s}:".format(sta.station, sta.channel.upper())) # Set Error Default to True erd = True @@ -456,7 +456,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Only a single day: Search for local data # Get Z localdata errZ, stZ = parse_localdata_for_comp( - comp='Z', stdata=stdata, sta=sta, start=start, end=end, + comp=zcomp, stdata=stdata, sta=sta, start=start, end=end, ndval=ndval) # Get N localdata errN, stN = parse_localdata_for_comp( @@ -481,53 +481,65 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # Construct location name if len(tloc) == 0: tloc = "--" - # Construct Channel List - channelsZNE = sta.channel.upper() + 'Z,' + sta.channel.upper() + \ - 'N,' + sta.channel.upper() + 'E' - print(("* {1:2s}[ZNE].{2:2s} - Checking Network".format( - sta.station, sta.channel.upper(), tloc))) - + # Get waveforms, with extra 1 second to avoid # traces cropped too short - traces are trimmed later + st = None + try: - st = client.get_waveforms( - network=sta.network, - station=sta.station, location=loc, - channel=channelsZNE, starttime=start, - endtime=end+1., attach_response=False) - if len(st) == 3: - print("* - ZNE Data Downloaded") - - # It's possible if len(st)==1 that data is Z12 - else: - # Construct Channel List - channelsZ12 = sta.channel.upper() + 'Z,' + \ - sta.channel.upper() + '1,' + \ - sta.channel.upper() + '2' - msg = "* {1:2s}[Z12].{2:2s} - Checking Network".format( - sta.station, sta.channel.upper(), tloc) - print(msg) - try: - st = client.get_waveforms( - network=sta.network, - station=sta.station, location=loc, - channel=channelsZ12, starttime=start, - endtime=end, attach_response=False) - if len(st) == 3: - print("* - Z12 Data Downloaded") - else: - print("* Stream has less than 3 components") - st = None - except: - st = None + # Construct Channel List + channelsZNE = sta.channel.upper() + zcomp.upper() + ',' + sta.channel.upper() + 'N,' + sta.channel.upper() + 'E' + print("* {0:2s}[{1:1s}NE].{2:2s} - Checking FDSNWS".format( + sta.channel.upper(), zcomp.upper(), tloc)) + st = client.get_waveforms(network=sta.network,station=sta.station, location=loc,channel=channelsZNE, starttime=start,endtime=end+1., attach_response=False) except: + print("* - No Data Available") + + #-- check if download got all needed data + if st is not None and len(st.select(component=zcomp)) >= 1 and len(st.select(component="N")) >= 1 and len(st.select(component='E')) >= 1: + print("* - "+zcomp.upper() + "NE Data Downloaded") + break + + + else: + # There was no data for above, so try other channels. + print("* - "+zcomp.upper() + "NE missing data") + + # Construct Channel List + channelsZ12 = sta.channel.upper() + zcomp.upper() +',' + sta.channel.upper() + '1,' + sta.channel.upper() + '2' + print("* {0:2s}[{1:1s}12].{2:2s} - Checking Network".format( + sta.channel.upper(), zcomp.upper(), tloc)) + try: + st = client.get_waveforms( + network=sta.network, + station=sta.station, location=loc, + channel=channelsZ12, starttime=start, + endtime=end, attach_response=False) + except: + print("* - No Data Available") + st = None + erd = True + + #-- check if download got all needed data + if st is not None and len(st.select(component=zcomp)) >= 1 and len(st.select(component="1")) >= 1 and len(st.select(component='2')) >= 1: + print("* - "+zcomp.upper() + "12 Data Downloaded") + break + else: + print("* - "+zcomp.upper() + "12 missing data") + + + if st is None: + print("* - Stream is missing components") st = None + erd = True + # Break if we successfully obtained 3 components in st if not erd: - break + #print (st) + # Check the correct 3 components exist if st is None: print("* Error retrieving waveforms") @@ -558,7 +570,7 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, str(tr.stats.starttime)+" " + str(tr.stats.endtime)) for tr in st] print("* True start: "+str(start)) - print("* -> Shifting traces to true 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)]) diff --git a/orientpy/plotting.py b/orientpy/plotting.py index be1af9f..2f3346a 100644 --- a/orientpy/plotting.py +++ b/orientpy/plotting.py @@ -93,7 +93,7 @@ def plot_bng_waveforms(bng, stream, dts, tt): """ fig, ax = plt.subplots(5, 1, sharey=True) - cmpts = ['Z', 'R', 'T', '1', '2'] + cmpts = [bng.zcomp.upper(), 'R', 'T', '1', '2'] taxis = np.linspace(-dts, dts, stream[0].stats.npts) for item in list(zip(ax, cmpts)): @@ -405,9 +405,9 @@ def plot_dl_results(stkey, R1phi, R1cc, R2phi, R2cc, ind, val, # Add text as title text = "Station "+stkey + \ - ": $\phi$ = {0:.1f} $\pm$ {1:.1f}".format(val, err) + ": $\phi$ = {0:.1f} $\pm$ {1:.1f} (n={2:.0f}, cc={3:.2f})".format(val, err, sum(ind),cc0) plt.suptitle(text, fontsize=12) gs.tight_layout(f, rect=[0, 0, 1, 0.95]) return plt - \ No newline at end of file + diff --git a/orientpy/scripts/bng_calc_auto.py b/orientpy/scripts/bng_calc_auto.py index 4950df0..a6c0bd8 100755 --- a/orientpy/scripts/bng_calc_auto.py +++ b/orientpy/scripts/bng_calc_auto.py @@ -6,7 +6,7 @@ import numpy as np from obspy.clients.fdsn import Client from obspy.taup import TauPyModel -from orientpy import BNG +from orientpy import BNG, utils from pathlib import Path from argparse import ArgumentParser @@ -105,6 +105,15 @@ def get_bng_calc_arguments(argv=None): "Options include: BGR, ETH, GEONET, GFZ, INGV, IPGP, IRIS, KOERI, " + "LMU, NCEDC, NEIP, NERIES, ODC, ORFEUS, RESIF, SCEDC, USGS, USP. " + "[Default IRIS]") + Svparm.add_argument( + "--cat_url", + action="store", + type=str, + dest="cat_url", + default=None, + help="Specify the obspy base_url server address (and port if needed) " + + "to open for the catalogue client. Overrides any settings to cat_client. " + + "[Default None]") Svparm.add_argument( "--waveform-source", action="store", @@ -115,6 +124,15 @@ def get_bng_calc_arguments(argv=None): "Options include: BGR, ETH, GEONET, GFZ, INGV, IPGP, IRIS, KOERI, " + "LMU, NCEDC, NEIP, NERIES, ODC, ORFEUS, RESIF, SCEDC, USGS, USP. " + "[Default IRIS]") + Svparm.add_argument( + "--wf_url", + action="store", + type=str, + dest="wf_url", + default=None, + help="Specify the obspy base_url server address (and port if needed) " + + "to open for the waveform client. Overrides any settings to wf_client. " + + "[Default None]") Svparm.add_argument( "-U", "--User-Auth", @@ -136,15 +154,25 @@ def get_bng_calc_arguments(argv=None): type=str, default="", help="Specify list of Station Keys in the database to process.") + stparm.add_argument( + "--Zcomp", + dest="zcomp", + type=str, + default='Z', + help="Specify the Vertical Component Channel Identifier. "+ + "[Default Z].") stparm.add_argument( "-c", "--coord-system", dest="nameconv", type=int, default=2, help="Coordinate system specification of instrument. " + - "(0) Attempt Autodetect between 1 and 2; (1) HZ, HN, HE; " + - "(2) Left Handed: HZ, H2 90 CW H1; (3) Right Handed: HZ, H2 90 CCW " + - "H1. [Default 2]") + "(0) Attempt Autodetect between 1 and 2; " + + "(1) HZ, HN, HE; " + + "(2) Left Handed: HZ, H2 90 CW H1; " + + "(3) Right Handed: HZ, H2 90 CCW H1 " + + "(4) Left Handed Numeric: H3, H2 90 CW H1 " + + "[Default 2]") #-- Timing Tmparm = parser.add_argument_group( @@ -358,7 +386,10 @@ def main(args=None): # Establish client for catalogue if args.verb > 1: print(" Establishing Catalogue Client...") - cat_client = Client(args.cat_client) + if args.cat_url is not None: + cat_client = Client(base_url=args.cat_url) + else: + cat_client = Client(args.cat_client) if args.verb > 1: print(" Done") @@ -366,9 +397,17 @@ def main(args=None): if args.verb > 1: print(" Establishing Waveform Client...") if len(args.UserAuth) == 0: - wf_client = Client(args.wf_client) + if args.wf_url is not None: + wf_client = Client(base_url=args.wf_url) + else: + wf_client = Client(args.wf_client) else: - wf_client = Client(args.wf_client, + if args.wf_url is not None: + wf_client = Client(base_url=args.wf_url, + user=args.UserAuth[0], + password=args.UserAuth[1]) + else: + wf_client = Client(args.wf_client, user=args.UserAuth[0], password=args.UserAuth[1]) if args.verb > 1: @@ -437,17 +476,24 @@ def main(args=None): cat = cat_client.get_events( starttime=tstart, endtime=tend, minmagnitude=args.minmag, maxmagnitude=args.maxmag) + + # get index of repeat events, save for later + reps = np.unique(utils.catclean(cat)) + except: raise(Exception(" Fatal Error: Cannot download Catalogue")) if args.verb > 1: - print("| Retrieved {0} events ".format(len(cat.events))) + print("| Retrieved {0} unique events of {1}".format(len(cat.events)-len(reps),len(cat.events))) print() for i, ev in enumerate(cat): + if i in reps: + continue + # Initialize BNGData object with station info - bngdata = BNG(sta) + bngdata = BNG(sta,zcomp=args.zcomp) # Add event to object accept = bngdata.add_event( diff --git a/orientpy/scripts/dl_average.py b/orientpy/scripts/dl_average.py index f1dd426..70b9235 100644 --- a/orientpy/scripts/dl_average.py +++ b/orientpy/scripts/dl_average.py @@ -72,6 +72,12 @@ def get_dl_average_arguments(argv=None): type=float, dest="cc", help="Cross-correlation threshold for final estimate. [Default 0.8]") + parser.add_argument( + "--min-mag", + default=5.5, + type=float, + dest="minmag", + help="Specify default minimum magnitude to include in average. [Default 5.5]") # Station Selection Parameters stparm = parser.add_argument_group( @@ -155,10 +161,11 @@ def main(args=None): continue meta = pickle.load(open(filename, 'rb')) - R1phi.append(meta.R1phi) - R2phi.append(meta.R2phi) - R1cc.append(meta.R1cc) - R2cc.append(meta.R2cc) + if meta.mag > args.minmag: + R1phi.append(meta.R1phi) + R2phi.append(meta.R2phi) + R1cc.append(meta.R1cc) + R2cc.append(meta.R2cc) R1phi = np.array(R1phi).flatten() R1cc = np.array(R1cc).flatten() @@ -168,18 +175,32 @@ def main(args=None): phi = np.concatenate((R1phi, R2phi), axis=None) cc = np.concatenate((R1cc, R2cc), axis=None) ind = cc > args.cc - val, err = utils.estimate(phi, ind) # output results to termianl - print("| D-L mean, error, data included: " + + print("| D-L mean, error, # robust: " + "{0:.2f}, {1:.2f}, {2}".format(val, err, np.sum(ind))) - print("| D-L CC level: {0:.1f}".format(args.cc)) + print("| D-L CC level: {0:.2f}".format(args.cc)) print() + if np.sum(np.isnan(np.array([val, err])))>0: continue + + #-- Save Text results + outfilename = indir / "{0:2s}.{1:s}.dat".format(sta.network,sta.station) + if not outfilename.is_file(): + fid=open(outfilename,'w') + fid.writelines(",