Physical Acoustics Lab

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

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

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)

In GEOPHYS330, we use python to download, plot and do a regression on a hundred years of sea level measurements in the Auckland harbour:

# these are the modules that contain functions we will use in the script:
import csv # this is to read a "comma separate values" file 
from urllib.request import urlopen # to read csv data from a website
import matplotlib.pyplot as plt # plotting functions
import numpy as np # numerical tools
from scipy.stats import linregress # linear regression function
import pandas as pd

# you can find and load the sea level data for Auckland here:
url="http://www.psmsl.org/data/obtaining/rlr.monthly.data/150.rlrdata"
df = pd.read_csv(url,delimiter=';')

# "bad" values are defined by a large negative number, so we only
# accept positive values:
df = df[df.iloc[:,1] >= 0]

# convert the pandas dataframe to numpy matrix:
npdat =df.as_matrix()
time= npdat[:,0]
height = npdat[:,1]

slope, intercept, r_value, p_value, std_err = linregress(time,height)
print("r-squared: %f", r_value**2)

# plot the sea height data:
plt.plot(time,height,'ro')
plt.plot(time,intercept+slope*time,'k')
plt.grid()
plt.xlabel('Date (year)')
plt.ylabel('Sea height (mm)')
plt.title('Sea height at Princes Wharf, Auckland (New Zealand)')
plt.legend(['measurements','best linear fit'],loc=4)
plt.axis('tight')
#plt.savefig('sealevel_auckland.png')
plt.show()