Instaseis Tutorial
Part 1: Introduction

Seismo-Live: http://seismo-live.org

Authors:

Goal: calculate seismograms on a 1D background model with Instaseis, assuming the Green's function database computed with AxiSEM is provided.

Overview:

Introduction:

  • Part 1: Very basic introduction to the API, calculate the first seismogram.

Basic Tasks:

  • Part 2: Some examples of interaction with obspy, calculate synthetics for a set of events and stations

Advanced Tasks:

  • Part 3: Plot record section
  • Part 4: Finite Source, compare to point source solution

Basic lines to set up the notebook and some paths.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os
import obspy
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (10, 8)

Basic Instaseis API Introduction

Please also have a look at our webpage - http://www.instaseis.net/ - where everything is documented extensively.

Opening a Database

To get going you have to import the package instaseis.

In [2]:
import instaseis

An Instaseis database must be opened before it can be used. These can be either a local files computed with AxiSEM, or even more easily a remote database hosted by IRIS (http://ds.iris.edu/ds/products/syngine/).

In [3]:
db = instaseis.open_db("syngine://prem_a_10s")

Some basic information about the loaded database can be reviewed by just printing it.

In [4]:
print(db)
SyngineInstaseisDB reciprocal Green's function Database (v8) generated with these parameters:
Syngine model name:      'prem_a_10s'
Syngine service version:  1.0.4
	components           : vertical and horizontal
	velocity model       : prem_ani
	attenuation          : True
	dominant period      : 10.000 s
	dump type            : displ_only
	excitation type      : monopole
	time step            : 2.441 s
	sampling rate        : 0.410 Hz
	number of samples    : 7376
	seismogram length    : 17999.6 s
	source time function : gauss_0
	source shift         : 17.084 s
	spatial order        : 4
	min/max radius       : 5671.0 - 6371.0 km
	Planet radius        : 6371.0 km
	min/max distance     : 0.0 - 180.0 deg
	time stepping scheme : symplec4
	compiler/user        : gfortran 5.2.0 by productext on dpstage
	directory/url        : http://service.iris.edu/irisws/syngine/1
	size of netCDF files : 86.0 GB
	generated by AxiSEM version v1.2-31-g69d4-dirty at 2015-11-12T09:36:01.000000Z

From this you can already glance a couple of aspects of the database used for this tutorial:

  • uses anisotropic prem as its 1D model
  • is accurate for periods down to 10 seconds
  • includes vertical and horizontal components
  • sources can have depths ranging from 0 to 700 km
  • five hour long seismograms

To avoid delays that become relevant when requesting very many seismograms, IRIS also offers the databases for download.

Receivers and Sources

Instaseis calculates seismograms for any source and receiver pair. A receiver has coordinates and optionally network and station codes. Using a reciprocal database, all receivers are assumed to be at the same depth, i.e. usually at the Earth surface.

In [5]:
rec = instaseis.Receiver(latitude=44.06238, longitude=10.59698,
                         network="IV", station="BDI")
print(rec)
Instaseis Receiver:
	Longitude :   10.6 deg
	Latitude  :   44.1 deg
	Network   : IV
	Station   : BDI
	Location  : 

Sources are naturally a bit more complex and Instaseis offers a variety of ways to define them. A straightforward way for earthquakes is to pass coordinates, moment as well as strike, dip and rake.

In [6]:
src = instaseis.Source.from_strike_dip_rake(
    latitude=27.77, longitude=85.37, depth_in_m=12000.0,
    M0=1e+21, strike=32., dip=62., rake=90.)
print(src)
Instaseis Source:
	Origin Time      : 1970-01-01T00:00:00.000000Z
	Longitude        :   85.4 deg
	Latitude         :   27.8 deg
	Depth            : 1.2e+01 km km
	Moment Magnitude :   7.93
	Scalar Moment    :   1.00e+21 Nm
	Mrr              :   8.29e+20 Nm
	Mtt              :  -2.33e+20 Nm
	Mpp              :  -5.96e+20 Nm
	Mrt              :   2.96e+20 Nm
	Mrp              :   4.74e+20 Nm
	Mtp              :  -3.73e+20 Nm

Note that origin time was not provided and hence defaults to 1970. A non double-couple source can directly be specified in terms of its moment tensor (note that instaseis uses SI units, i.e. NM, while GCMT uses dyn cm).

In [7]:
src = instaseis.Source(
    latitude=27.77, longitude=85.37, depth_in_m=12000.0,
    m_rr=8.29e+20, m_tt=-2.33e+20, m_pp=-5.96e+20,
    m_rt=2.96e+20, m_rp=4.74e+20, m_tp=-3.73e+20)
print(src)
Instaseis Source:
	Origin Time      : 1970-01-01T00:00:00.000000Z
	Longitude        :   85.4 deg
	Latitude         :   27.8 deg
	Depth            : 1.2e+01 km km
	Moment Magnitude :   7.93
	Scalar Moment    :   1.00e+21 Nm
	Mrr              :   8.29e+20 Nm
	Mtt              :  -2.33e+20 Nm
	Mpp              :  -5.96e+20 Nm
	Mrt              :   2.96e+20 Nm
	Mrp              :   4.74e+20 Nm
	Mtp              :  -3.73e+20 Nm

Sidenote: The moment tensor can be visualized using the Beachball function from obspy.imaging:

In [8]:
from obspy.imaging.beachball import beachball

mt = src.tensor / src.M0 # normalize the tensor to avoid problems in the plotting
beachball(mt, size=200, linewidth=2, facecolor='b');

Now we are ready to extract synthetic seismograms from the database. The return type is an obspy stream object, which can directly be plotted:

In [9]:
st = db.get_seismograms(source=src, receiver=rec, components='ZNERT')
print(st)
st.plot();
5 Trace(s) in Stream:
IV.BDI..LXZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T04:59:42.561567Z | 0.4 Hz, 7369 samples
IV.BDI..LXN | 1970-01-01T00:00:00.000000Z - 1970-01-01T04:59:42.561567Z | 0.4 Hz, 7369 samples
IV.BDI..LXE | 1970-01-01T00:00:00.000000Z - 1970-01-01T04:59:42.561567Z | 0.4 Hz, 7369 samples
IV.BDI..LXR | 1970-01-01T00:00:00.000000Z - 1970-01-01T04:59:42.561567Z | 0.4 Hz, 7369 samples
IV.BDI..LXT | 1970-01-01T00:00:00.000000Z - 1970-01-01T04:59:42.561567Z | 0.4 Hz, 7369 samples

Done This is all you need for a basic usage of Instaseis!