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.
Month: August 2014
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)