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 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' 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/docs/rfpy.rst b/docs/rfpy.rst index da3095e..73aa40a 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 ----------------------------------------- @@ -83,6 +83,59 @@ 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 ``--SDS-path`` 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 `_. + +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 +------------- + +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 =========== diff --git a/docs/scripts.rst b/docs/scripts.rst index d955b94..dcef120 100644 --- a/docs/scripts.rst +++ b/docs/scripts.rst @@ -38,148 +38,121 @@ Usage usage: rfpy_calc [arguments] - Script used to download and pre-process three-component ('Z', 'N', and 'E'), - seismograms for individual events and calculate teleseismic P-wave receiver - functionsThis version requests data on the fly for a given date range. Data - are requested from the internet using the client services framework. The - stations are processed one by one and the data are stored to disk. + Script used to download and pre-process three-component ('Z', 'N', and 'E'), seismograms for + individual events and calculate teleseismic P-wave receiver functionsThis version requests + data on the fly for a given date range. Data are requested from the internet using the client + services framework. The stations are processed one by one and the data are stored to disk. + + positional arguments: + 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 - --keys STKEYS Specify a comma separated list of station keys for - which to perform the analysis. These must be contained - within the station database. Partial keys will be used - to match against those in the dictionary. For - instance, providing IU will match with all stations in - the IU network [Default processes all stations in the - database] - -v, -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. - [Default Z]. - -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 - database. + --keys STKEYS Specify a comma separated list of station keys for which to perform + the analysis. These must be contained within the station database. + Partial keys will be used to match against those in the dictionary. + For instance, providing IU will match with all stations in the IU + network [Default processes all stations in the database] + -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. [Default Z]. + -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 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 of 'AUSPASS', 'BGR', - 'EARTHSCOPE', 'EIDA', 'EMSC', 'ETH', 'GEOFON', - 'GEONET', 'GFZ', 'ICGC', 'IESDMC', 'INGV', 'IPGP', - 'IRIS', 'IRISPH5', 'ISC', 'KNMI', 'KOERI', 'LMU', - 'NCEDC', 'NIEP', 'NOA', 'NRCAN', 'ODC', 'ORFEUS', - 'RASPISHAKE', 'RESIF', 'RESIFPH5', 'SCEDC', 'TEXNET', - 'UIB-NORSAR', 'USGS', 'USP'). [Default 'IRIS'] - --user-auth USERAUTH Authentification Username and Password for the - waveform server (--user-auth='username:authpassword') - to access and download restricted data. [Default no - user and password] + “http://service.iris.edu”) or key string for recognized server (one + of 'AUSPASS', 'BGR', 'EARTHSCOPE', 'EIDA', 'EMSC', 'ETH', 'GEOFON', + 'GEONET', 'GFZ', 'ICGC', 'IESDMC', 'INGV', 'IPGP', 'IRIS', 'IRISPH5', + 'ISC', 'KNMI', 'KOERI', 'LMU', 'NCEDC', 'NIEP', 'NOA', 'NRCAN', + 'ODC', 'ORFEUS', 'RASPISHAKE', 'RESIF', 'RESIFPH5', 'SCEDC', + 'TEXNET', 'UIB-NORSAR', 'USGS', 'USP'). [Default 'IRIS'] + --user-auth USERAUTH Authentification Username and Password for the waveform server + (--user-auth='username:authpassword') to access and download + restricted data. [Default no user and password] --eida-token TOKENFILE - 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. - This mechanism is only available on select EIDA nodes. - The token can be provided in form of the PGP message - as a string, or the filename of a local file with the + 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. This mechanism is only + available on select EIDA nodes. The token can be provided in form of + the PGP message as a string, or the filename of a local file with the 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 files. - - --local-data LOCALDATA - Specify a comma separated list of paths containing - day-long sac 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 - --dtype DTYPE 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 casehere. - --no-data-zero Specify to force missing data to be set as zero, - rather than default behaviour which sets to nan. - --no-local-net Specify to prevent using the Network code in the - search for local data (sometimes for CN stations the - dictionary name for a station may disagree with that - in the filename. [Default Network used] - --save-Z12 Specify to save Z12 (un-rotated) components. [Default - False] + Settings associated with a SeisComP database for locally archived + data. + + --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 + for details on the SDS format. If this option is used, it takes + precedence over the --server settings. + --dtype DTYPE 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. + --save-Z12 Specify to save Z12 (un-rotated) components. [Default False] Event Settings: - Settings associated with refining the events to include in matching event- - station pairs - - --start STARTT Specify a UTCDateTime compatible string representing - the start time for the event search. This will - override any station start times. [Default start date - of station] - --end ENDT Specify a UTCDateTime compatible string representing - the end time for the event search. This will override - any station end times [Default end date of station] - --reverse, -R Reverse order of events. Default behaviour starts at - oldest event and works towards most recent. Specify - reverse order and instead the program will start with - the most recent events and work towards older - --minmag MINMAG Specify the minimum magnitude of event for which to - search. [Default 6.0] - --maxmag MAXMAG Specify the maximum magnitude of event for which to - search. [Default None, i.e. no limit] + Settings associated with refining the events to include in matching event-station pairs + + --start STARTT Specify a UTCDateTime compatible string representing the start time + for the event search. This will override any station start times. + [Default start date of station] + --end ENDT Specify a UTCDateTime compatible string representing the end time for + the event search. This will override any station end times [Default + end date of station] + --reverse Reverse order of events. Default behaviour starts at oldest event and + works towards most recent. Specify reverse order and instead the + program will start with the most recent events and work towards older + --minmag MINMAG Specify the minimum magnitude of event for which to search. [Default + 6.0] + --maxmag MAXMAG Specify the maximum magnitude of event for which to search. [Default + None, i.e. no limit] Geometry Settings: - Settings associatd with the event-station geometries for the specified - phase - - --phase PHASE Specify the phase name to use. Be careful with the - distance. setting. Options are 'P' or 'PP'. [Default - 'P'] - --mindist MINDIST Specify the minimum great circle distance (degrees) - between the station and event. [Default depends on - phase] - --maxdist MAXDIST Specify the maximum great circle distance (degrees) - between the station and event. [Default depends on - phase] + Settings associatd with the event-station geometries for the specified phase + + --phase PHASE Specify the phase name to use. Be careful with the distance. setting. + Options are 'P' or 'PP'. [Default 'P'] + --mindist MINDIST Specify the minimum great circle distance (degrees) between the + station and event. [Default depends on phase] + --maxdist MAXDIST Specify the maximum great circle distance (degrees) between the + station and event. [Default depends on phase] Parameter Settings: Miscellaneous default values and settings --sampling-rate NEW_SAMPLING_RATE Specify new sampling rate in Hz. [Default 10.] - --dts DTS Specify the window length in sec (symmetric about - arrival time). [Default 150.] - --align ALIGN Specify component alignment key. Can be either ZRT, - LQT, or PVH. [Default ZRT] - --vp VP Specify near-surface Vp to use with --align=PVH - (km/s). [Default 6.0] - --vs VS Specify near-surface Vs to use with --align=PVH - (km/s). [Default 3.5] - --dt-snr DT_SNR Specify the window length over which to calculate the - SNR in sec. [Default 30.] - --pre-filt PRE_FILT Specify two floats with low and high frequency corners - for pre-filter (before deconvolution). [Default None] - --fmin FMIN Specify the minimum frequency corner for SNR and CC - filter (Hz). [Default 0.05] - --fmax FMAX Specify the maximum frequency corner for SNR and CC - filter (Hz). [Default 1.0] + --dts DTS Specify the window length in sec (symmetric about arrival time). + [Default 150.] + --align ALIGN Specify component alignment key. Can be either ZRT, LQT, or PVH. + [Default ZRT] + --vp VP Specify near-surface Vp to use with --align=PVH (km/s). [Default 6.0] + --vs VS Specify near-surface Vs to use with --align=PVH (km/s). [Default 3.5] + --dt-snr DT_SNR Specify the window length over which to calculate the SNR in sec. + [Default 30.] + --pre-filt PRE_FILT Specify two floats with low and high frequency corners for pre-filter + (before deconvolution). [Default None] + --fmin FMIN Specify the minimum frequency corner for SNR and CC filter (Hz). + [Default 0.05] + --fmax FMAX Specify the maximum frequency corner for SNR and CC filter (Hz). + [Default 1.0] Deconvolution Settings: Parameters for deconvolution - --method METHOD Specify the deconvolution method. Available methods - include 'wiener', 'water' and 'multitaper'. [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] - + --method METHOD Specify the deconvolution method. Available methods include 'wiener', + 'wiener-mod', 'water' and 'multitaper'. [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] ``rfpy_recalc`` ++++++++++++++++ @@ -212,65 +185,54 @@ Usage usage: rfpy_recalc [arguments] - Script used to re-calculate receiver functions that already exist on disk, but - using different processing options. The stations are processed one by one and - the data are over-written to disk. + Script used to re-calculate receiver functions that already exist on disk, but using + different processing options. The stations are processed one by one and the data are stored + to disk. Note: The sampling rate cannot be changed to a new rate positional arguments: - indb Station Database to process from. - - optional arguments: - -h, --help show this help message and exit - --keys STKEYS Specify a comma separated list of station keys for - which to perform the analysis. These must be contained - within the station database. Partial keys will be used - to match against those in the dictionary. For instance, - providing IU will match with all stations in the IU - network [Default processes all stations in the - database] - -v, -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 - database. + indb Station Database to process from. + options: + -h, --help show this help message and exit + --keys STKEYS Specify a comma separated list of station keys for which to perform + the analysis. These must be contained within the station database. + Partial keys will be used to match against those in the dictionary. + For instance, providing IU will match with all stations in the IU + network [Default processes all stations in the database] + -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 database. Parameter Settings: - Miscellaneous default values and settings - - --Z12 Use Z12 data if available. [Default uses ZNE data] - --phase PHASE Specify the phase name to use. Be careful with the - distance. setting. Options are 'P', 'PP', 'allP', 'S', - 'SKS' or 'allS'. [Default 'allP'] - --resample RESAMPLE Specify the new sampling-rate for the receiver - functions. Note the sampling rate of the original data - (ZNE or Z12) stored on disk is unchanged. [Default - None] - --align ALIGN Specify component alignment key. Can be either ZRT, - LQT, or PVH. [Default ZRT] - --vp VP Specify near-surface Vp to use with --align=PVH (km/s). - [Default 6.0] - --vs VS Specify near-surface Vs to use with --align=PVH (km/s). - [Default 3.5] - --dt-snr DT_SNR Specify the window length over which to calculate the - SNR in sec. [Default 30.] - --pre-filt PRE_FILT Specify two floats with low and high frequency corners - for pre-filter (before deconvolution). [Default None] - --fmin FMIN Specify the minimum frequency corner for SNR filter - (Hz). [Default 0.05] - --fmax FMAX Specify the maximum frequency corner for SNR filter - (Hz). [Default 1.0] + Miscellaneous default values and settings + + --Z12 Use Z12 data if available. [Default uses ZNE data] + --phase PHASE Specify the phase name to use. Be careful with the distance. setting. + Options are 'P', 'PP', 'allP', 'S', 'SKS' or 'allS'. [Default 'allP'] + --resample RESAMPLE Specify the new sampling-rate for the receiver functions. Note the + sampling rate of the original data (ZNE or Z12) stored on disk is + unchanged. [Default None] + --align ALIGN Specify component alignment key. Can be either ZRT, LQT, or PVH. + [Default ZRT] + --vp VP Specify near-surface Vp to use with --align=PVH (km/s). [Default 6.0] + --vs VS Specify near-surface Vs to use with --align=PVH (km/s). [Default 3.5] + --dt-snr DT_SNR Specify the window length over which to calculate the SNR in sec. + [Default 30.] + --pre-filt PRE_FILT Specify two floats with low and high frequency corners for pre-filter + (before deconvolution). [Default None] + --fmin FMIN Specify the minimum frequency corner for SNR filter (Hz). [Default + 0.05] + --fmax FMAX Specify the maximum frequency corner for SNR filter (Hz). [Default + 1.0] Deconvolution Settings: - Parameters for deconvolution - - --method METHOD Specify the deconvolution method. Available methods - include 'wiener', 'water' and 'multitaper'. [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] + Parameters for deconvolution + --method METHOD Specify the deconvolution method. Available methods include 'wiener', + 'wiener-mod', 'water' and 'multitaper'. [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] ``rfpy_plot`` ++++++++++++++++ @@ -315,7 +277,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 +388,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 +527,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 +651,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/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" } 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' 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 c3d1ed7..931b0ed 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) """ @@ -94,6 +94,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. @@ -155,8 +158,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 +201,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 ---------- @@ -232,7 +235,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, @@ -253,11 +256,6 @@ def add_data(self, stream, returned=False, new_sr=5.): returned : bool Whether or not to return the ``accept`` attribute - Attributes - ---------- - zne_data : :class:`~obspy.core.Stream` - Stream container for NEZ seismograms - Returns ------- accept : bool @@ -266,7 +264,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 @@ -281,13 +279,14 @@ 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 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 @@ -296,36 +295,32 @@ 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. + P phase arrival and resamples the waveforms. Parameters ---------- 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() + 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 @@ -336,15 +331,10 @@ def download_data(self, client, stdata=[], dtype='SAC', ndval=np.nan, 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: - raise(Exception("Requires event data as attribute - aborting")) + raise Exception("Requires event data as attribute - aborting") if not self.meta.accept: return @@ -353,18 +343,21 @@ 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, + client=client, + sta=self.sta, + start=tstart, + end=tend, + new_sr=new_sr, + verbose=verbose, remove_response=remove_response, - local_response_dir=local_response_dir, - verbose=verbose, zcomp=self.zcomp) + zcomp=self.zcomp) # Store as attributes with traces in dictionary try: @@ -380,7 +373,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 +397,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 +406,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. @@ -464,8 +457,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 @@ -545,8 +536,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,15 +624,14 @@ 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, + 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 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 ---------- @@ -654,11 +643,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 - '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' + Type of wavelet for deconvolution. Options are 'complete', 'time' + or 'envelope' envelope_threshold : float Threshold [0-1] used in ``wavelet='envelope'``. time : float @@ -668,7 +657,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 @@ -680,7 +669,6 @@ def deconvolve(self, phase='P', vp=None, vs=None, Stream containing the receiver function traces """ - if not self.meta.accept: return @@ -703,33 +691,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 @@ -742,8 +736,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 @@ -751,7 +745,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 @@ -774,21 +769,21 @@ 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): + + 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-mod": # npad = _npow2(nn*2) npad = nn @@ -798,17 +793,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-mod": + Sdenom = Spp + 0.25*Snpp + 0.25*Snd1d1 + 0.5*Snpd1 elif method == 'water': phi = np.amax(Spp)*wlevel Sdenom = Spp @@ -825,11 +826,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) @@ -845,7 +846,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) @@ -869,9 +870,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() @@ -893,7 +894,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() @@ -922,37 +923,70 @@ def _decon(parent, daughter1, daughter2, noise, nn, method): if phase == 'P' or 'PP': - # Get signal length (i.e., seismogram to deconvolve) from trace length - over = 5 - dtsqt = len(trL.data)*trL.stats.delta/2. + # Get signal length (i.e., seismogram to deconvolve) from + # trace length + # 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) + nn = int(round((dts - 5.)*trL.stats.sampling_rate)) + 1 - # 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) + 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]] + 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]] + 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': - # 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) @@ -962,7 +996,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.) @@ -971,25 +1005,25 @@ 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: - 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': - 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] @@ -998,14 +1032,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() @@ -1013,11 +1046,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() @@ -1035,7 +1087,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. @@ -1068,7 +1119,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 09b3394..6c73ff7 100644 --- a/rfpy/scripts/rfpy_calc.py +++ b/rfpy/scripts/rfpy_calc.py @@ -27,8 +27,10 @@ import numpy as np import pickle import stdb -from obspy.clients.fdsn import Client -from obspy import Catalog, UTCDateTime +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, read_inventory from http.client import IncompleteRead from rfpy import utils, RFData from pathlib import Path @@ -39,12 +41,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( usage="%(prog)s [arguments] ", @@ -59,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", @@ -74,7 +71,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, @@ -106,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", @@ -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. " @@ -146,48 +142,31 @@ 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( - "--local-data", + "--SDS-path", action="store", type=str, - dest="localdata", + dest="sds_path", default=None, - help="Specify a comma separated list of paths containing " + - "day-long sac or mseed files of data already downloaded. " + - "If data exists for a seismogram is already present on disk, " + - "it is selected preferentially over downloading " + - "the data using the Client interface") - DataGroup.add_argument( - "--dtype", - action="store", - type=str, - dest="dtype", - default='SAC', - help="Specify the data archive file type, either SAC " + - " or MSEED. Note the default behaviour is to search for " + - "SAC files. Local archive files must have extensions of '.SAC' "+ - " or '.MSEED. These are case dependent, so specify the correct case"+ - "here.") - DataGroup.add_argument( - "--no-data-zero", - action="store_true", - dest="ndval", - default=False, - help="Specify to force missing data to be set as zero, rather " + - "than default behaviour which sets to nan. Note this is applied " + - "only to the SAC data") + help="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 " + + "for details on the SDS format. If this option is used, it takes " + + "precedence over the --server settings.") DataGroup.add_argument( - "--no-local-net", - action="store_false", - dest="useNet", - default=True, - help="Specify to prevent using the Network code in the " + - "search for local data (sometimes for CN stations " + - "the dictionary name for a station may disagree with that " + - "in the filename. [Default Network used]") + "--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( "--save-Z12", action="store_true", @@ -220,7 +199,7 @@ def get_calc_arguments(argv=None): "the end time for the event search. This will override any " + "station end times [Default end date of station]") EventGroup.add_argument( - "--reverse", "-R", + "--reverse", action="store_true", dest="reverse", default=False, @@ -362,7 +341,8 @@ def get_calc_arguments(argv=None): type=str, default="wiener", help="Specify the deconvolution method. Available methods " + - "include 'wiener', 'water' and 'multitaper'. [Default 'wiener']") + "include 'wiener', 'wiener-mod', 'water' and 'multitaper'. " + + "[Default 'wiener']") DeconGroup.add_argument( "--gfilt", action="store", @@ -386,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(',') @@ -394,7 +380,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) @@ -405,7 +391,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) @@ -419,32 +405,21 @@ 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 = [] - # Check Datatype specification - if (not args.dtype.upper() == 'MSEED') and (not args.dtype.upper() == 'SAC'): - parser.error( - "Error: Local Data Archive must be of types 'SAC'" + - "or MSEED. These must match the file extensions for " + - " the archived data.") - - # Check NoData Value - if args.ndval: - args.ndval = 0.0 - else: - args.ndval = nan + 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 distances for selected phase if args.phase not in ['P', 'PP', 'S', 'SKS']: @@ -508,10 +483,10 @@ 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', 'wiener-mod', 'water', 'multitaper']: parser.error( - "Error: 'method' should be either 'wiener', 'water' or " + - "'multitaper'") + "Error: 'method' should be either 'wiener', 'wiener-mod', " + + "'water' or 'multitaper'") return args @@ -533,7 +508,7 @@ def main(): # Run Input Parser args = get_calc_arguments() - # Load Database + # Get station StDb and keys db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # Loop over station keys @@ -557,14 +532,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.sds_path 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.sds_path, + 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 +561,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(" ") @@ -628,17 +607,20 @@ 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 # 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()) @@ -646,8 +628,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: @@ -668,28 +652,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) @@ -761,9 +723,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/scripts/rfpy_ccp.py b/rfpy/scripts/rfpy_ccp.py index d367684..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 @@ -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] ", @@ -67,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, @@ -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', @@ -438,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 = [] @@ -469,9 +474,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 +622,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 +675,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 +707,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 +727,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 +759,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 273cb4d..c4ef5c7 100644 --- a/rfpy/scripts/rfpy_harmonics.py +++ b/rfpy/scripts/rfpy_harmonics.py @@ -27,9 +27,10 @@ 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 +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 @@ -37,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] ", @@ -67,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, @@ -86,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( @@ -272,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) @@ -283,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) @@ -337,8 +331,7 @@ def main(): print("# __ _ _ #") print("# _ __ / _|_ __ _ _ | |__ __ _ _ __ _ __ ___ ___ _ __ (_) ___ ___ #") print("# | '__| |_| '_ \| | | | | '_ \ / _` | '__| '_ ` _ \ / _ \| '_ \| |/ __/ __| #") - print( - "# | | | _| |_) | |_| | | | | | (_| | | | | | | | | (_) | | | | | (__\__ \ #") + print("# | | | _| |_) | |_| | | | | | (_| | | | | | | | | (_) | | | | | (__\__ \ #") print("# |_| |_| | .__/ \__, |___|_| |_|\__,_|_| |_| |_| |_|\___/|_| |_|_|\___|___/ #") print("# |_| |___/_____| #") print("# #") @@ -348,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 = [] @@ -391,13 +396,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(" ") @@ -510,7 +514,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 5741818..05e2c82 100644 --- a/rfpy/scripts/rfpy_hk.py +++ b/rfpy/scripts/rfpy_hk.py @@ -27,9 +27,10 @@ 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 +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 @@ -37,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] ", @@ -67,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, @@ -359,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) @@ -370,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) @@ -396,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: @@ -416,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( @@ -493,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 = [] @@ -542,13 +549,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(" ") @@ -700,7 +706,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_plot.py b/rfpy/scripts/rfpy_plot.py index 310d782..aa0503b 100644 --- a/rfpy/scripts/rfpy_plot.py +++ b/rfpy/scripts/rfpy_plot.py @@ -27,8 +27,9 @@ import numpy as np import pickle import stdb -from obspy import Stream, UTCDateTime -from rfpy import binning, plotting +import copy +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 @@ -63,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, @@ -341,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 = [] @@ -368,13 +381,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..5de76bd 100644 --- a/rfpy/scripts/rfpy_recalc.py +++ b/rfpy/scripts/rfpy_recalc.py @@ -29,17 +29,13 @@ import numpy as np import pickle import stdb -from rfpy import RFData +import copy +from rfpy import RFData, utils from pathlib import Path +from obspy import read_inventory def get_recalc_arguments(argv=None): - """ - Get Options from :class:`~optparse.OptionParser` objects. - - This function is used for data processing on-the-fly (requires web connection) - - """ parser = ArgumentParser( usage="%(prog)s [arguments] ", @@ -68,7 +64,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, @@ -179,7 +175,8 @@ def get_recalc_arguments(argv=None): type=str, default="wiener", help="Specify the deconvolution method. Available methods " + - "include 'wiener', 'water' and 'multitaper'. [Default 'wiener']") + "include 'wiener', 'wiener-mod', 'water' and 'multitaper'. " + + "[Default 'wiener']") DeconGroup.add_argument( "--gfilt", action="store", @@ -224,10 +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']: + if args.method not in ['wiener', 'wiener-mod', 'water', 'multitaper']: parser.error( - "Error: 'method' should be either 'wiener', 'water' or " + - "'multitaper'") + "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(',')] @@ -258,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 = [] @@ -285,13 +294,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(" ") @@ -344,8 +352,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 6a7d718..ae09dd3 100644 --- a/rfpy/utils.py +++ b/rfpy/utils.py @@ -1,9 +1,31 @@ +# 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 -from obspy import Stream, Inventory -from obspy import read, read_inventory +import copy +from stdb import StDbElement def floor_decimal(n, decimals=0): @@ -39,377 +61,14 @@ def traceshift(trace, tt): return rtrace -def list_local_data_stn(lcldrs=list, sta=None, net=None, dtype='SAC', altnet=[]): +def download_data(client=None, sta=None, start=UTCDateTime(), + end=UTCDateTime(), new_sr=0., verbose=False, + remove_response=False, zcomp='Z'): """ - 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'): - """ - 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! + 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 ---------- @@ -421,20 +80,16 @@ 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 + new_sr : 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() - 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 @@ -450,99 +105,43 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=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()))) - - # 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 + for loc in sta.location: + + # Construct location name + 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(), loc) + 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=tloc, + 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") + elif len(st) < 3: + print("* Error retrieving waveforms") + print("**************************************************") + return True, None # Check the correct 3 components exist if st is None: @@ -552,10 +151,16 @@ def download_data(client=None, sta=None, start=UTCDateTime, end=UTCDateTime, # 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 +168,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 +184,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