Month: August 2014

Welcome to our sister station in Tohoku, Japan!

As of this week, we have a sister station in Tohoku Japan: THJP. The first images coming out of this station show how seismically active that region is. Of course everyone remembers the Big Tohoku Earthquake. Our colleagues at Tohoku University have been playing an important role in understanding the underground structure, and tectonic processes that shape — and shake — Japan.

firstdatafromTHJP

Click on this image  to enlarge it. You will see the signals of several local events annotated by Dr. Kentaro Emoto using Japan’s seismic network Hi-net.

Python script to plot an epicentral distance on a map

In GEOPHYS310, we learn python. One of the exercises is to plot a station (circle) on a map, with a circle representing some epicentral distance to an earthquake:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

proj = ccrs.PlateCarree(central_<wbr>longitude=180)

ax1 = plt.subplot(111, projection=proj)
ax1.set_extent([-50, 50, -40, 40], crs=proj)
ax1.stock_img()

ax1.tissot(rad_km=3200, lons=180, lats=10, alpha=0.3)
ax1.plot(0,10,'rv') # note the longitude is now with respect to the central_longitude, but NOT for the tissot circle! 

plt.show()

A python script to estimate the epicentral distance and origin time of an earthquake

In GEOPHYS330, we do a ray tracing exercise. Here is some code to perform the task of estimating the origin time and the epicentral distance of an earthquake:

from obspy import UTCDateTime
import numpy as np
import matplotlib.pyplot as plt
from obspy import read
from obspy.taup import TauPyModel

def main():
    # estimate the epicentral distance. This is one of your free parameters:
    Delta = 35 # in degrees
    # estimate the origin time of the earthquake; your other free parameter:
    t0=UTCDateTime("2014-04-13T12:25:19")

    maxamp = readandplotseismogram("https://pal.auckland.ac.nz/files/2014/08/quake.sac_.gz")
    computeandplottts(Delta,t0,maxamp)

    # tighten up the axes, and show:
    plt.axis('tight')
    plt.show()
    #return

def readandplotseismogram(seismogram):
    '''read and plot the seismogram'''
    # read the stream
    st = read(seismogram)
    starttime = UTCDateTime("2014-04-13T12:38:00")
    st.trim(starttime,starttime+1000)
    # the trace in the stream is
    tr = st[0]
    # Make a time array and an amps array:
    t1= tr.stats.starttime.timestamp
    t2= tr.stats.endtime.timestamp
    npts = tr.stats.npts
    times = np.linspace(t1,t2,npts)
    amps = tr.data
    maxamp = max(amps)
    #Plot the seismogram against a grid
    plt.figure(figsize=(12,3))
    plt.plot(times,amps,color='blue')
    #Converting tick marks to actual times:
    locs, labels = plt.xticks()
    new_label=[]
    for loc in locs:
        new_label.append(UTCDateTime(loc).strftime("%H-%M-%S"))
        plt.xticks(locs,new_label,rotation=10)
        plt.xlabel('Time (HH-MM-SS UTC)')
        plt.ylabel('Displacement')
        plt.grid()
    return maxamp

def computeandplottts(Delta,t0,maxamp):
    # compute arrival times based on a epicentral distance (and depth, but
    # I fixed depth to a constant value):
    model = TauPyModel(model="iasp91")      
    arrivals = model.get_travel_times(distance_in_degree=Delta,source_depth_in_km=10)
    #Construct vertical lines to show arrival times of predicted seismic waves
    for arrival in arrivals:
        dummy = t0+ arrival.time
        if arrival.name == "P":
            plt.vlines(dummy,-maxamp/2,maxamp/2)
            plt.text(dummy,maxamp/2+0.05*maxamp,arrival.name)
        if arrival.name == "PP":
            plt.vlines(dummy,-maxamp/2,maxamp/2)
            plt.text(dummy,maxamp/2+0.05*maxamp,arrival.name)
        if arrival.name == "S":
            plt.vlines(dummy,-maxamp/2,maxamp/2)
            plt.text(dummy,maxamp/2+0.05*maxamp,arrival.name)
        if arrival.name == "SS":
            plt.vlines(dummy,-maxamp/2,maxamp/2)
            plt.text(dummy,maxamp/2+0.05*maxamp,arrival.name)
    return()

# this will actually run the code if called stand-alone:
if __name__ == '__main__':
    main()

An implementation of raytracing software Taup in obspy

In one of our labs for GEOPHYS330, we use the ray tracing software TauP. Below is an example of how to use TauP from within python (using obspy), our computer language of choice for the lab, course, and department:

from obspy.taup.tau import plot_ray_paths
plot_ray_paths(source_depth=100,phase_list=['P'],npoints=25)