{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

\n", "\n", "\n", "

\n", "\n", "

IRIS-EarthScope Short Course

\n", "
Bloomington/IN, August 2015
\n", "
\n", "

Python/ObsPy Introduction

\n", "
\n", "\n", "
\n", "
image by Matthias Meschede
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise -- The 2008 Mt. Carmel, Illinois, Earthquake and Aftershock Series\n", "\n", "Introduction from \n", "[\"The April 18, 2008 Illinois Earthquake: An ANSS Monitoring Success\" by Robert B. Herrmann, Mitch Withers, and Harley Benz, SRL 2008](http://srl.geoscienceworld.org/content/79/6/830.extract):\n", "\n", "*\"The largest-magnitude earthquake in the past 20 years struck near **Mt. Carmel in southeastern Illinois** on Friday morning, 18 April 2008 at 09:36:59 UTC (04:37 CDT). **The Mw 5.2 earthquake was felt over an area that spanned Chicago and Atlanta**, with about 40,000 reports submitted to the U.S. Geological Survey (USGS) “Did You Feel It?” system. There were at least six felt aftershocks greater than magnitude 3 and 20 aftershocks with magnitudes greater than 2 located by regional and national seismic networks. **Portable instrumentation was deployed** by researchers of the University of Memphis and Indiana University (the first portable station was installed at about 23:00 UTC on 18 April). The portable seismographs were deployed both to capture near-source, high-frequency ground motions for significant aftershocks and to better understand structure along the active fault. [...]\"*\n", "\n", "\"World\n", "\"Felt\n", "\n", "Web page hits at USGS/NEIC during 24 hours after the earthquake:\n", "\n", " - peak rate: **3,892 hits/second**\n", " - **68 million hits** in the 24 hours after the earthquake\n", "\n", "\"USGS/NEIC\n", "\n", "Some links:\n", "\n", " - http://earthquake.usgs.gov/earthquakes/dyfi/events/us/2008qza6/us/index.html\n", " - http://earthquake.usgs.gov/earthquakes/eqinthenews/2008/us2008qza6/#summary\n", " - http://srl.geoscienceworld.org/content/79/6/830.extract\n", " - http://srl.geoscienceworld.org/content/82/5/735.short\n", " \n", "Again, please execute the following cell, it contains a few adjustments to the notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "from __future__ import print_function\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "plt.style.use('ggplot')\n", "plt.rcParams['figure.figsize'] = 12, 8" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Use the search on the [ObsPy docs](https://docs.obspy.org/) for any functionality that you do not remember/know yet..!\n", "### E.g. [searching for \"filter\"](https://docs.obspy.org/search.html?q=filter)..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Download and visualize main shock\n", "\n", "Request information on stations recording close to the event from IRIS using the [`obspy.fdsn Client`](https://docs.obspy.org/packages/obspy.fdsn.html), print the requested station information." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Download waveform data for the mainshock for one of the stations using the FDSN client (if you get an error, maybe try a different station and/or ask for help). Make the preview plot using ObsPy." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize a Spectrogram (if you have time, you can play around with the different parameters for the spectrogram). Working on a copy of the downloaded data, apply a filter, then trim the requested data to some interesting parts of the earthquake and plot the data again." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a function `plot_data(t)` that fetches waveform data for this station and shows a preview plot of 20 seconds of data starting at a given time. It should take a `UTCDateTime` object as the single argument." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test your function by calling it for the time of the main shock" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Visualize aftershock and estimate magnitude\n", "\n", "Read file \"`./data/mtcarmel.mseed`\". It contains data of stations from an aftershock network that was set up shortly after the main shock. Print the stream information and have a look at the network/station information, channel names time span of the data etc.. Make a preview plot. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The strongest aftershock you see in the given recordings is at `2008-04-21T05:38:30`. Trim the data to this aftershock and make a preview plot and a spectrogram plot of it." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a very simple approximation of the magnitude. Use the function provided below (after you execute the code box with the function you can call it anywhere in your code boxes).\n", "\n", "Demean the data and then determine the raw data's peak value of the event at one of the stations (e.g. using Python's `max` function or a NumPy method on the data array) and call the provided function for that value. (Note that this is usually done on the horizontal components.. we do it on vertical for simplicity here)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from obspy.signal import estimateMagnitude\n", "\n", "def mag_approx(peak_value, frequency, hypo_dist=20):\n", " \"\"\"\n", " Give peak value of raw data for a very crude and simple magnitude estimation.\n", "\n", " For simplicity, this is done assuming hypocentral location, peak frequency, etc.\n", " To keep it simple for now the response information is entered manually here\n", " (it is the same for all instruments used here).\n", " \"\"\"\n", " poles = [-1.48600E-01 + 1.48600E-01j,\n", " -1.48600E-01 - 1.48600E-01j,\n", " -4.14690E+02 + 0.00000E+00j,\n", " -9.99027E+02 + 9.99027E+02j,\n", " -9.99027E+02 - 9.99027E+02j]\n", " zeros = [0.0 + 0.0j,\n", " 0.0 + 0.0j,\n", " 1.1875E+03 + 0.0j]\n", " norm_factor = 7.49898E+08\n", " sensitivity = 6.97095E+05\n", " paz = {'poles': poles, 'zeros': zeros, 'gain': norm_factor,\n", " 'sensitivity': sensitivity}\n", " ml = estimateMagnitude(paz, peak_value, 0.5 / frequency, hypo_dist)\n", " return ml" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do the magnitude approximation in a for-loop for all stations in the Stream. Calculate a network magnitude as the average of all three stations." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a function `netmag(st)` that returns a network magnitude approximation. It should take a Stream object (which is assumed to be trimmed to an event) as only argument. Use the provided `mag_approx` function and calculate the mean of all traces in the stream internally." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test your function on the cut out Stream object of the large aftershock from before." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Advanced\n", "You can also download the station metadata using the FDSN client and extract poles and zeros information and directly use the `estimateMagnitude` function without using the hard-coded response information. \n", "\n", "###3. Detect aftershocks using triggering routines\n", "\n", "Read the 3-station data from file `\"./data/mtcarmel.mseed\"` again. Apply a bandpass filter, adjust it to the dominant event frequency range you have seen in the aftershock spectrogram before. Run a recursive STA/LTA trigger on the filtered data (see [ObsPy docs](https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.trigger.html)). The length of the `sta` window should be a bit longer than an arriving seismic phase, the `lta` window can be around ten to twenty times as long.\n", "\n", "Make a preview plot of the Stream object, now showing the characteristic function of the triggering. High trigger values indicate transient signals (of the frequency range of interest) that might be an event (or just a local noise burst on that station..).\n", "\n", "(play around with different values and check out the resulting characteristic function)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could now manually compare trigger values on the different stations to find small aftershocks, termed a network coincidence trigger. However, there is a convenience function in ObsPy's signal toolbox to do just that in only a few lines of code.\n", "\n", "Read the data again and apply a bandpass to the dominant frequencies of the events. Use the `coincidenceTrigger` function that returns a list of possible events (see the [ObsPy Tutorial](https://docs.obspy.org/tutorial/code_snippets/trigger_tutorial.html#network-coincidence-trigger-example) for an example of a recursive STA/LTA network coincidence trigger). Print the length of the list and adjust the trigger-on/off thresholds so that you get around 5 suspected events.\n", "\n", "Print the first trigger in the list to show information on the suspected event." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Go over the list of triggers in a for-loop. For each trigger/suspected event:\n", "\n", " - print the time of the trigger\n", " - read the waveform data, use [`starttime` and `endtime` arguments for `read()`](https://docs.obspy.org/packages/autogen/obspy.core.stream.read.html) to trim the data to the suspected event right during reading (avoiding to read the whole file again and again)\n", " - calculate and print the network magnitude using the `netmag(st)` function from earlier\n", " - make a preview plot\n", "\n", "If you're curious you can compare the crude magnitude estimates with the [table of aftershocks](http://www.seismosoc.org/publications/srl/SRL_82/srl_82-5_hamburger_et_al-esupp/Table_S2.txt) provided by the scientists that analyzed the aftershock sequence. The paper with details can be found here: [\"Aftershocks of the 2008 Mt. Carmel, Illinois, Earthquake: Evidence for Conjugate Faulting near the Termination of the Wabash Valley Fault System\" by M. W. Hamburger, K. Shoemaker, S. Horton, H. DeShon, M. Withers, G. L. Pavlis and E. Sherrill, SRL 2011](http://srl.geoscienceworld.org/content/82/5/735.short)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Advanced:\n", "You can also use event templates of some good signal-to-noise ratio aftershocks to detect more weaker aftershocks and select from weak triggers based on waveform similarities like shown in the ObsPy tutorial. \n", "\n", "## 4. Cross correlation pick alignment\n", "\n", "An unfortunate undergrad student picked his/her way through the aftershock sequence. For a high-precision relative location run (e.g. double difference, master-event relocation) you want to have more precise, cross correlation aligned phase picks.\n", "\n", "The approach by [Deichmann and Garcia-Fernandez (1992, GJI)](https://onlinelibrary.wiley.com/doi/10.1111/j.1365-246X.1992.tb02088.x/abstract;jsessionid=F73CF9485C0432A3E14EDDDD0C73E058.d03t02) is implemented in the function [**`xcorrPickCorrection`**](https://docs.obspy.org/packages/autogen/obspy.signal.cross_correlation.xcorrPickCorrection.html) in `obspy.signal.cross_correlation`. Follow the [example in the ObsPy Tutorial](https://docs.obspy.org/tutorial/code_snippets/xcorr_pick_correction.html). Get time corrections for the event picks given in `pick_times` relative to the event pick `reference_pick`.\n", "\n", "Use the data in file `\"./data/mtcarmel_100hz.mseed\"`. Before applying the correction resample to 200 Hz (the parabola fitting works better with the finer time resolution) and optionally bandpass to a relatively narrow frequency range." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reference_pick = \"2008-04-19T12:46:45.96\"\n", "\n", "pick_times = [\"2008-04-19T13:08:59.19\",\n", " \"2008-04-19T16:55:19.65\",\n", " \"2008-04-19T19:03:38.72\",\n", " \"2008-04-19T19:28:53.54\"]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }