Monthly Archives: February 2017

Searching the USGS Earthquake Catalog with Mathematica

Posted on by

This note documents a simple in-class activity that I use in my GEOSC 109H earthquakes and society class. In class I describe the Gutenberg-Richter Relation that represents the distribution of earthquakes as a function of magnitude (a nice paper-based activity is available from SCEC). I show the students a plot from the GCMT catalog in class, in this activity, they use Mathematica to download data for a few randomly chosen years and see if the pattern I showed in class holds up.

The Mathematica notebook begins with a few utility functions that allow students to use web services to acquire seismicity data and a few functions to list some events (which is partly a check to see that things are working). Below is a listing of the search script (an example usage is provided below). 

comcatRectangularSearchPath[start_, end_, minlat_, maxlat_, minlon_, maxlon_] := Module[{},
  "http://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=" <> 
   DateString[start, "ISODateTime"] <> "&endtime=" <> 
   DateString[end, "ISODateTime"] <> 
   "&minmagnitude=4.0&minlatitude=" <> 
   ToString[Round[minlat, 0.001]] <>
   "&maxlatitude=" <> 
   ToString[Round[maxlat, 0.001]] <>
   "&minlongitude=" <> 
   ToString[Round[minlon, 0.001]] <> "&maxlongitude=" <> 
   ToString[Round[maxlon, 0.001]] <> "&orderby=time-asc"
  ]

comcatRectangularSearch[start_, end_, minlat_: - 90, maxlat_: 90, 
  minlon_: - 180, maxlon_: 180] := Module[{path},
   path = comcatRectangularSearchPath[start, end, minlat, maxlat, minlon, maxlon];
  {path, Import[path, "CSV"]}
  ]

The scripts are simplistic, and to avoid errors for some of the more active years, I limit the magnitude range of the analysis to events greater than or equal to magnitude 4.0. The code below is a function to list the largest events in the search results.

Continue reading →

A Day-of-Year Calculator with Python

Posted on by

Many seismic data are stored using a numerical day of the year in place of a month and day. This is efficient for the data storage, but it leads to the need to convert back and forth from command month-day to day-of-year values. Over the years, I have seen printed tables of values that some seismologists kept nearby to do the conversion back and forth. Now, of course, we work at computers much of the day and we generally have an interactive terminal available, so I created a simple python utility to perform the conversion.

Modern python has extensive packages to deal with time quantities, but I use a simple, lightweight package that I wrote to handle the conversions (back in the mid 90’s when python had many fewer packages and most everyone was using PERL).

The date to day-of-year conversion is relatively straight forward, so there’s not much code to write – here is what I use.

doy.py

!/usr/bin/python
import sys,string
import pytutil
#
if len(sys.argv) == 4:
	
	month = string.atoi(sys.argv[1])
	day   = string.atoi(sys.argv[2])
	year  = string.atoi(sys.argv[3])
	
	print "   %02d/%02d/%04d is the %03d of the year." \
		  % (month,day,year,cja_pytutil.doy(year,month,day))
elif len(sys.argv) == 3:
	doy   = string.atoi(sys.argv[1])
	year  = string.atoi(sys.argv[2])
	output = cja_pytutil.monthdate(year,doy)
	print "   Day %03d of %4d is %02d/%02d/%04d " \
		  % (doy,year,output[0],output[1],year,)
else:
	print " "
	print "  Usage: month day year"
	print "      or"
	print "  Usage: doy year"
	print " "

Continue reading →

Exploring Seismic Motions with an Arduino

Posted on by

In my Earthquakes and Society course I use an arduino and MEMs sensor to introduce students to seismograms and to help them develop an intuitive understanding of seismogram analysis. The idea is to have them build their own sensors and explore the motions they can produce with the simple system (while it is connected to the computer in the classroom). The picture is pretty much all the information they need to connect the hardware, but I provide a little background on a breadboard, etc. Within a few minutes, non-science students are quite comfortable with the system.

Arduino accelerometer system. The arduino is the larger board, the small breakout contains an accelerometer.

Continue reading →

Seismometer Responses with Mathematica

Posted on by

Seismometers, Poles, & Zeros

Although not an everyday calculation, seismologists often benefit from a reminder of the instrument response curves (amplitude responses). Seismometers are generally described as the ratio of two polynomials in the Laplace Transform variable, s. The polynomials are described by their roots, which we call zeros (for the numerator polynomial) and poles (for the denominator polynomial). Generally when and instrument is calibrated, the information on the calibration is communicated by the poles, zeros, and an overall gain. To compare instruments with different gains, we also often choose a normalization frequency for display.

In the following I present some Mathematica scripts that I use to examine the poles, zeros, and instrument response amplitude spectrum. Slight modifications can introduce the phase information as well, but to keep things a little short, I leave that out of the description. We’ll look at four sets of instruments, Trillium, Steckheisen, and Guralp broadband seismometers, and then some short-period instruments that are part of the IRIS Consortium instrument pool.

Nanometrics Trillium Sensor Responses

I obtained the poles and zeros and zeros for the instrument responses from the PASSCAL web site and normalized the responses to a frequency of 1 Hz.  Click here for the link.

Here is the code to set up three sets of poles and zeros as Mathematica lists – corresponding to the Trillium 120 PA, 240 (2nd Gen), and the Compact Trillium sensors.

(* 
  Poles and zeros for the Trillium 120PA from 
  http://www.passcal.nmt.edu/webfm_send/1970
*)
Zn120 = {0, 0, -90.0, -160.7, -3108};
Pn120 = {-0.03852 + 0.03658 I, -0.03852 - 0.03658 I, -178.,
         -135 + 160 I, -135 - 160 I, -671 + 1154 I, -671 - 1154 I};

(* 
  Poles and zeros for the Trillium 240 (2nd Gen) from 
  http://www.passcal.nmt.edu/webfm_send/1970 
*)
Zn240 = {0, 0, -90., -164.2, -3203.};
Pn240 = {-0.01813 + 0.01803 I, -0.01813 - 0.01803 I, -124.9 + 0 I,
         -197.5 + 256.1 I, -197.5 - 256.1 I, -569 + 1150 I, -569 + -1150 I};
(* 
  Poles and zeros for the Compact Trillium from 
  http://www.passcal.nmt.edu/webfm_send/1970 
*)
Znct = {0, 0, -434.1};
Pnct = {-0.036910 + 0.037120 I, -0.036910 + -0.037120 I, -371.2, -373.9 + 475.5 I,
   -373.9 - 475.5 I, -588. + 1508. I, -588.4 - 1508. I};

I can compute the three instrument responses using the following Mathematica functions. I could write one function and call it with the poles and zeros, but this is not that much extra work for an example. The second line in each function normalizes the response to unit amplitude at a frequency of one Hertz.

(* Instrument responses normalized at one Hz *)
Element[f, Reals];
s = I 2 \[Pi] f;
tr120resp[f_] = (Times @@ (s - Zn120)/Times @@ (s - Pn120)) / 
                (Abs[Times @@ (s - Zn120)/ Times @@ (s - Pn120)] /. f -> 1.0);

tr240resp[f_] = (Times @@ (s - Zn240)/Times @@ (s - Pn240)) / 
                (Abs[Times @@ (s - Zn240)/Times @@ (s - Pn240)] /. f -> 1.0);

ctresp[f_] = (Times @@ (s - Znct)/Times @@ (s - Pnct))/
             (Abs[Times @@ (s - Znct)/ Times @@ (s - Pnct)] /. f -> 1.0);

The variable \( s = i \omega \) is the Laplace Transform variable and has units of radians. I can display the instrument response amplitude spectrum using

p1 = LogLogPlot[{Abs[tr240resp[f]], Abs[tr120resp[f]], Abs[ctresp[f]]}, {f, 0.001, 1000 },
   Filling -> Axis, ImageSize -> 512,
   Axes -> False, Frame -> True,
   BaseStyle -> {14, FontFamily -> "Helvetica"},
   FrameLabel -> {"Frequency (Hz)", "Normalized Amplitude"},
   PlotLabel -> "Trillium Sensors" ,
   LabelStyle -> {14, FontFamily -> "Helvetica"},
   PlotLegends -> LineLegend[{"TR240", "TR120", "Compact"}]];
p1 // Print

Here is the plot.

Nanometrics Trillium seismometer response curves (amplitude).

Continue reading →