First Basic Exercise¶

Task: Calculate three component synthetic seismograms for the stations and events in the data/events and data/stations subdirectories and save them on disc.

Notes¶

Receiver objects can also be created from StationXML, SEED, or STATIONS files as well as obpy inventories using instaseis.Receiver.parse(); see the documentation for details.
Source objects can also be created from QuakeML, CMTSOLUTIONS, and in other ways using instaseis.Source.parse(); see the documentation for details.
The get_seismograms() method has a couple of extra arguments:
- kind: displacement, velocity, acceleration
- remove_source_shift, reconvolve_stf, dt,
... see the documentation for details.
You can use the properties of the Receiver and Source objects to create usefull filenames.

Basic lines to set up the notebook and some paths.

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

Import Instaseis and open the database:

import instaseis
db = instaseis.open_db("data/database")

1. Load Receivers¶

reminder: you can use ObsPy to load stations and plot a map:

from obspy import read_inventory

inventory = read_inventory('data/stations/all_stations.xml')
inventory.plot(projection="local", resolution="i");

This inventory can directly be used as input to generate a list of instaseis.Receiver objects:

receivers = instaseis.Receiver.parse(inventory)
for rec in receivers[:2]:
    print(rec)

Alternatively, instaseis can directly open the station xml or STATIONS file (but then you don't have the nice plot):

receivers = instaseis.Receiver.parse('data/stations/all_stations.xml')
print(receivers[0])
receivers = instaseis.Receiver.parse('data/stations/STATIONS')
print(receivers[0])

2. Load Events¶

reminder: use ObsPy to load events from a QuakeML file containing all events and plot a map:

import glob # provides iterator to loop over files

cat = obspy.core.event.Catalog()

for filename in glob.iglob('data/events/quakeml/*.xml'):
     cat += obspy.read_events(filename)
        
print(cat)
print(instaseis.Source.parse(cat.events[0]))
cat.plot();

Alternatively load QuakeML or CMTSOLUTION files directly using instaseis.Source.parse() and store the sources in a list:

sources = []

for filename in glob.iglob('data/events/quakeml/*.xml'):
    sources.append(instaseis.Source.parse(filename))
    
print(sources[0])

for filename in glob.iglob('data/events/cmtsolutions/*'):
    sources.append(instaseis.Source.parse(filename))

print(sources[0])

3. Extract Seismograms and Save to File¶

For the first solution using a ObsPy event catalog:

dt = 1.0

for event in cat:
    src = instaseis.Source.parse(event)
    srcname = '%s_Mw_%3.1f' % (src.origin_time.date, src.moment_magnitude)
    for rec in receivers:
        # create a usefull filename
        recname = '%s_%s' % (rec.network, rec.station)
        filename = '%s_%s' % (recname, srcname)
        filename = filename.replace('.', '_')
        
        # extract seismograms using instaseis
        
        
        # write to miniseed files in the data_out folder. Write as MiniSEED due to multi
        # component support.

For the second solution use a list of sources:

dt = 1.0

for src in sources:
    srcname = '%s_Mw_%3.1f' % (src.origin_time.date, src.moment_magnitude)
    for rec in receivers:
        # create a usefull filename
        recname = '%s_%s' % (rec.network, rec.station)
        filename = '%s_%s' % (recname, srcname)
        filename = filename.replace('.', '_')
        
        # extract seismograms using instaseis
        
        
        # write to miniseed files in the data_out folder. Write as MiniSEED due to multi
        # component support.

Check the data:

ls data_out/