.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/DataCollectionSAC.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_DataCollectionSAC.py: SAC Database ============ In this Tutorial we are going to get all good receiver functions for the year 2018 for station IU-HRV ([Adam Dziewonski Observatory](http://www.seismology.harvard.edu/hrv.html)). To Compute the receiver function we need to download and organize the observed data. PyGLImER will automatically take care of these things for you in the background, but do check out the 'Database Docs' if you aren't familiar. Downloading Event & Station Metadata ------------------------------------ To download the data we use the :class:`pyglimer.waveform.request.Request` class. The first method from this class that we are going to use is the download event catalog public method :func:`pyglimer.waveform.request.Request.download_evtcat`, to get a set of events that contains all wanted earthquakes. This method is launched automatically upon initialization. To initialize said `class` we setup a parameter dictionary, with all the needed information. Let's look at the expected information: .. GENERATED FROM PYTHON SOURCE LINES 28-31 .. code-block:: default # sphinx_gallery_thumbnail_number = 1 # sphinx_gallery_dummy_images = 1 .. GENERATED FROM PYTHON SOURCE LINES 32-33 First let's get a path where to create the data. .. GENERATED FROM PYTHON SOURCE LINES 33-80 .. code-block:: default # Some needed Imports import os from obspy import UTCDateTime from pyglimer.waveform.request import Request # Get notebook path for future reference of the database: try: db_base_path = ipynb_path except NameError: try: db_base_path = os.path.dirname(os.path.realpath(__file__)) except NameError: db_base_path = os.getcwd() # Define file locations proj_dir = os.path.join(db_base_path, 'tmp', 'database_sac') request_dict = { # Necessary arguments 'proj_dir': proj_dir, 'raw_subdir': os.path.join('waveforms', 'raw'),# Directory of the waveforms 'prepro_subdir': os.path.join('waveforms', 'preprocessed'), # Directory of the preprocessed waveforms 'rf_subdir': os.path.join('waveforms', 'RF'), # Directory of the receiver functions 'statloc_subdir': 'stations', # Directory stations 'evt_subdir': 'events', # Directory of the events 'log_subdir': 'log', # Directory for the logs 'loglvl': 'WARNING', # logging level 'format': 'sac', # Format to save database in "phase": "P", # 'P' or 'S' receiver functions "rot": "RTZ", # Coordinate system to rotate to "deconmeth": "waterlevel", # Deconvolution method "starttime": UTCDateTime(2021, 1, 1, 0, 0, 0), # Starttime of database. # Here, starttime of HRV "endtime": UTCDateTime(2021, 7, 1, 0, 0, 0), # Endtimetime of database # kwargs below "pol": 'v', # Source wavelet polaristion. Def. "v" --> SV "minmag": 5.5, # Earthquake minimum magnitude. Def. 5.5 "event_coords": None, # Specific event?. Def. None "network": "IU", # Restricts networks. Def. None "station": "HRV", # Restricts stations. Def. None "waveform_client": ["IRIS"], # FDSN server client (s. obspy). Def. None "evtcat": None, # If you have already downloaded a set of # events previously, you can use them here } .. GENERATED FROM PYTHON SOURCE LINES 81-83 Now that all parameters are in place, let's initialize the :class:`pyglimer.waveform.request.Request` .. GENERATED FROM PYTHON SOURCE LINES 83-87 .. code-block:: default # Initializing the Request class and downloading the data R = Request(**request_dict) .. GENERATED FROM PYTHON SOURCE LINES 88-90 The initialization will look for all events for which data is available. To see whether the events make sense we plot a map of the events: .. GENERATED FROM PYTHON SOURCE LINES 90-101 .. code-block:: default import matplotlib.pyplot as plt from pyglimer.plot.plot_utils import plot_catalog from pyglimer.plot.plot_utils import set_mpl_params # Setting plotting parameters set_mpl_params() # Plotting the catalog plot_catalog(R.evtcat) .. image-sg:: /tutorials/images/sphx_glr_DataCollectionSAC_001.png :alt: DataCollectionSAC :srcset: /tutorials/images/sphx_glr_DataCollectionSAC_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 102-103 We can also quickly check how many events we gathered. .. GENERATED FROM PYTHON SOURCE LINES 103-106 .. code-block:: default print(f"There are {len(R.evtcat)} available events") .. rst-class:: sphx-glr-script-out .. code-block:: none There are 289 available events .. GENERATED FROM PYTHON SOURCE LINES 107-128 Downloading waveform data and station information ------------------------------------------------- The next step on our journey to receiver functions is retrieving the corresponding waveform data. To retrieve the waveform data, we can two different public methods of `Request`: `download_waveforms()` or, as in this case `download_waveforms_small_db()`. Both methods use the station and event locations to get viable records of receiver function depending on epicentral distance and traveltimes. `download_waveforms()` relies on obspy's massdownloader and is best suited for extremely large databases (i.e., from many different networks and stations). `download_waveforms_small_db()` is faster for downloads from few stations and, additionally, has the advantage that the desired networks and stations can be defined as lists and not only as strings. Note that this method requires you to define the channels, you'd like to download (as a string, wildcards allowed). ***NOTE:*** This might take a while. .. GENERATED FROM PYTHON SOURCE LINES 128-131 .. code-block:: default R.download_waveforms_small_db(channel='BH?') .. rst-class:: sphx-glr-script-out .. code-block:: none --> Enter TOA event loop for station IU.HRV .. GENERATED FROM PYTHON SOURCE LINES 132-133 Let's have a quick look at how many miniseeds we actually downloaded. .. GENERATED FROM PYTHON SOURCE LINES 133-143 .. code-block:: default from glob import glob # Path to the where the miniseeds are stored data_storage = os.path.join( proj_dir, 'waveforms','raw','P','**', '*.mseed') # Print output print(f"Number of downloaded waveforms: {len(glob(data_storage))}") .. rst-class:: sphx-glr-script-out .. code-block:: none Number of downloaded waveforms: 69 .. GENERATED FROM PYTHON SOURCE LINES 144-171 The final step to get you receiver function data is the preprocessing. Although it is hidden in a single function, which is :func:`pyglimer.waveform.request.Request.preprocess` A lot of decisions are being made: Processing steps: 1. Clips waveform to the right length (tz before and ta after theorethical arrival.) 2. Demean & Detrend 3. Tapering 4. Remove Instrument response, convert to velocity & simulate havard station 5. Rotation to NEZ and, subsequently, to RTZ. 6. Compute SNR for highpass filtered waveforms (highpass f defined in qc.lowco) If SNR lower than in qc.SNR_criteria for all filters, rejects waveform. 7. Write finished and filtered waveforms to folder specified in qc.outputloc. 8. Write info file with shelf containing station, event and waveform information. 9. (Optional) If we had chosen a different coordinate system in ``rot`` than RTZ, it would now cast the preprocessed waveforms information that very coordinate system. 10. Deconvolution with method ``deconmeth`` from our dict is perfomed. It again uses the request class to perform this. The ``if __name__ ...`` expression is needed for running this examples .. GENERATED FROM PYTHON SOURCE LINES 171-174 .. code-block:: default R.preprocess(hc_filt=1.5, client='single') .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/3 [00:00 .. GENERATED FROM PYTHON SOURCE LINES 233-234 Let's zoom into the first 20 seconds (~200km) .. GENERATED FROM PYTHON SOURCE LINES 234-237 .. code-block:: default rftrace.plot(lim=[0,20]) .. image-sg:: /tutorials/images/sphx_glr_DataCollectionSAC_003.png :alt: DataCollectionSAC :srcset: /tutorials/images/sphx_glr_DataCollectionSAC_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 238-243 Plot RF section +++++++++++++++ Since we have an entire stream of receiver functions at hand, we can plot a section .. GENERATED FROM PYTHON SOURCE LINES 243-246 .. code-block:: default rfstream.plot(scalingfactor=1) .. image-sg:: /tutorials/images/sphx_glr_DataCollectionSAC_004.png :alt: PRF component :srcset: /tutorials/images/sphx_glr_DataCollectionSAC_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 247-249 Similar to the single RF plot we can provide time and epicentral distance limits: .. GENERATED FROM PYTHON SOURCE LINES 249-256 .. code-block:: default timelimits = (0, 20) # seconds epilimits = (32, 36) # epicentral distance rfstream.plot( scalingfactor=0.25, lim=timelimits, epilimits=epilimits, linewidth=0.75) .. image-sg:: /tutorials/images/sphx_glr_DataCollectionSAC_005.png :alt: PRF component :srcset: /tutorials/images/sphx_glr_DataCollectionSAC_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 257-259 By increasing the scaling factor and removing the plotted lines, we can already see trends: .. GENERATED FROM PYTHON SOURCE LINES 259-265 .. code-block:: default rfstream.plot( scalingfactor=0.5, lim=timelimits, epilimits=epilimits, line=False) .. image-sg:: /tutorials/images/sphx_glr_DataCollectionSAC_006.png :alt: PRF component :srcset: /tutorials/images/sphx_glr_DataCollectionSAC_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 266-267 As simple as that you can create your own receiver functions with just a single smalle script. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 1 minutes 12.934 seconds) .. _sphx_glr_download_tutorials_DataCollectionSAC.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: DataCollectionSAC.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: DataCollectionSAC.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_