Merge pull request #25 from nfsi-canada/schaefferaj-master

Added ability to specify ZComp name (instead of just "Z"), to handle CH[124] data from NFSI.

Added ability to specify a local base_url for the fdsnws data_select, as well as a separate base_url for the fdsnws event selection.

Author: paudetseis

.github/workflows/deploy.yml

@@ -17,11 +17,12 @@ on:
 jobs:
   # This workflow contains a single job called "build"
   build:
+    name: ${{ matrix.os }} Python ${{ matrix.python-version }}
     # The type of runner that the job will run on
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: ["ubuntu-latest", "macos-latest"]
+        os: ["ubuntu-latest"]
         python-version: ["3.10"]
 
     steps:
@@ -38,19 +39,19 @@ jobs:
         run: |
           conda info
           conda config --add channels conda-forge
-          conda install obspy pandas pytest
+          conda install obspy pandas pytest libstdcxx-ng
           pip install stdb
           pip install geographiclib
           pip install -e .
           conda list
 
-      - name: Tests
-        shell: bash -l {0}
-        run: |
-          mkdir empty
-          cd empty
-          conda install pytest
-          pytest ../orientpy/tests/
+      # - name: Tests
+      #   shell: bash -l {0}
+      #   run: |
+      #     mkdir empty
+      #     cd empty
+      #     conda install pytest
+      #     pytest ../orientpy/tests/
       #    conda install pytest-cov
       #    pytest -v --cov=orientpy ../orientpy/tests/
       #    bash <(curl -s https://codecov.io/bash)

docs/tutorials.rst

@@ -10,13 +10,13 @@ LOBS3 and send the prompt to a logfile
 
 .. code-block::
 
-   $ query_fdsn_stdb.py -N YH -C ?H? -S LOBS3 YH_list > logfile
+   $ query_fdsn_stdb -N YH -C ?H? -S LOBS3 YH_list > logfile
 
-To check the station info for M08A, use the program ``ls_stdb.py``:
+To check the station info for M08A, use the program ``ls_stdb``:
 
 .. code-block::
 
-    $ ls_stdb.py YH_list.pkl
+    $ ls_stdb YH_list.pkl
     Listing Station Pickle: YH_list.pkl
     YH.LOBS3
     --------------------------------------------------------------------------

orientpy/classes.py

@@ -119,14 +119,15 @@ class Orient(object):
 
     """
 
-    def __init__(self, sta):
+    def __init__(self, sta, zcomp='Z'):
 
         # Attributes from parameters
         self.sta = 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,21 +373,23 @@ 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:
             test_set = test_sets[test_key]
             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

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)])

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
-    
+