Methods for using Python and ObsPy: The first examples

Methods for using Python and ObsPy: The first examples

You will need to open a terminal window.  I personally use the XQuartz App in MacOS, which allows me X-window graphics for my own applications, but a standard Terminal window should be OK.  Python handles its own graphics.

start with switching to the BASH shell, then invoking Python

bash
python

you get a new prompt with Python, meaning that you have entered the world of Python syntax.  Fear not.  To conclude your Python session, you can CTRL-D or enter the command "exit()" to leave.

It is also possible to invoke the fancier IPython package in a similar way:

bash
ipython

try some simple commands that are easy to verify

q=2
q==2
q==3
q
q+1
b=2*q
print(b)

Then some syntax tricks: multiple commands on a line separated by semicolons

a=1;b=2;c=3
a
b
c==3

command split into two lines with a backslash, but only at an operation or equal sign

xval=3.14
xval
yval=\
3.14
yval
zval= yval +\
xval
zval

You must import modules in order to use them

import numpy

You can import modules and abbreviate their names

import numpy as np

Python syntax requires that you attach any function to the library it comes from.  A simple command like

d=sqrt(2)

will fail.  The square root function sqrt is in numpy, and to use it you need to specify where the function comes from

d=numpy.sqrt(2)
d

This is why abbreviations of Python modules are useful

ee=np.sqrt(3)
ee

dont want 16-digit precision?  Use the round command to limit the answer to 4 or 6 digits:

round(ee,4)

round(numpy.sqrt(2),6)

You can import all the functions in numpy if you like, and never need to use the above syntax. This could cause confusion if two functions in two separate modules have the same name.

import single functions:

from numpy import sqrt
sqrt(6)

to import all numpy functions:

from numpy import *

In IPython, it is possible to clear the memory of the python package, which erases the values of all the variables, with the command

%reset

In regular Python the %reset command wont work.

Unlike high-level languages like Fortran and C, the variables in Python are assigned their variable-type by context. The commands

aval = 1.0
bval = "BigNose"

will create aval as a floating-point variable, and bval as a string variable

L= [1,2,3]

creates an integer list of length 3. Try

L.reverse()
L
L.pop(2)
L

A list can have dissimilar objects, e.g.,

LL = [1,"D",3]

A tuple is a list that behaves more like an array in Fortran or C, with objects that are all the same.

M=(1,2,3)

Integer tuples are used to define arrays of numbers

amat = numpy.zeros((3,5))

sets amat to be a 3 by 5 matrix of zeros, with 3 rows and 5 columns

Python does loops differently from what you are accustomed to.  The loops are declared with the word "for" and a colon. The loops are not closed by a statement, but rather by a change in indentation.  Try

b,c,d=3,4,5
b
c
d

for a in np.arange(c,d,0.1):
   x = (-b + np.sqrt(b**2 + 4*c*a))/(2.0*a)
   print(a,b,c,x)

note that the function arange(x1,x2,x3) steps in equal increments of x3 from x1 to x2

The loop executes differently if the last line is not indented.  The indentation defines the loop.  Nested loops are indicated by deeper indentations. To print the results at the end of the loop rather than within it, add a blank line to close the loop.

for a in np.arange(c,d,0.1):
   x = (-b + np.sqrt(b**2 + 4*c*a))/(2.0*a)

print(a,b,c,x)

Loops and branching commands "while" and "if" and "else" follow similar syntax.  Python requires a colon to tell it where the loop-definition ends and the loop-operations begin

Vector operations can be created in place of loops.  The code

b,c,d=3,4,5
a = np.arange(c,d,0.1)
x = (-b + np.sqrt(b**2 + 4*c*a))/(2.0*a)
print(a)
print(x)
print(a,x)

generates a vector a to hold the output of the arange function.  Placing the vector a into the equation line that follows vectorizes the equation, so that the variable x is also a vector.  The print commands will print vectors as well as scalar variables.

you will want to write scriptfiles with common command sequences, e.g., importing the usual libraries.  A scriptfile in Python conventionally uses the suffix ".py"

A scriptfile program.py can be executed by

python program.py

but the Python package exits after executing the script.  The interactive Python version IPython can execute scriptfiles within the program using the %run command and exit within the Python environment.  You can also generate scriptfiles within IPython with the %edit command. Try

echo import numpy as np >! startup.py
echo import matplotlib.pyplot as plt >> startup.py
ipython
%run startup.py

These lines of code generate a 2-line startup script, start IPython, and execute the script.  At conclusion you will still be operating in the IPython environment. If you then execute

%edit startup.py

the IPython terminal becomes the window for the Unix vi editor.  Honestly, you might prefer to create your scriptfiles in a fancier editor, because vi is old and klunky.

In scriptfiles you can add comments, which is VERY useful.  Lines that start with a hashtag # are comments.  If a line contains """ - three doublequotes, then it declares the next group of lines to be comments, until another line with """ is found. NOTE: the MacOS TextEdit app does not create double-quotes that Python can read properly.  I use an X-window editor named "nedit" that can be downloaded from the fink archive - but you may wish to develop your own kluge for this problem.  The nedit editor is almost as old as the vi editor.

Beyreuther et al (2010) has an ObsPy coding example that reads seismic data within one time interval from one component at one station, accessed via the webdc data-aggregator.  The code in the paper must be modified to work with the 2017 implementation of ObsPy, because the syntax and function names have been updated over time.  Here is the example that generates Figure 1, and that works for my MacBook Pro.

# Figure 1 from Beyreuther etal (2010)
# commands are updated for 2017 versions of Python and ObsPy
from obspy.core import UTCDateTime
from obspy.clients.arclink import Client
# old command from Beyreuther et al (2010), function libraries and names have changed
# from obspy.signal import cornFreq2Paz, seisSim, lowpass
from obspy.signal.filter import lowpass
from obspy.signal.invsim import corn_freq_2_paz, simulate_seismometer
import numpy as np, matplotlib.pyplot as plt
#1 Retrieve Data via Arclink -- note that the "host" variable is suppressed
client = Client("webdc.eu", port=18001)
# note that ObsPy will search data archives for time intervals based on UTC time
t = UTCDateTime("2009-08-24 00:20:03")
# The data will be plotted raw and bandpassed to simulate a 1-Hz seismic sensor
one_hertz = corn_freq_2_paz(1.0) # 1Hz instrument
st = client.get_waveforms("BW", "RJOB", "", "EHZ", t, t+30)
#  PAZ means "pole and zeroes" which is one way to define a digital filter
#  zeroes in filter response arise from moving-average filters
#  poles in the filter response arise from auto-regressive filters
paz = client.get_paz("BW", "RJOB", "", "EHZ", t, t+30)
#2 Correct for frequency response of the instrument
# note that st is the stream object that holds the data, but its information is
# partitioned into data and metadata
# the data are expressed in instrument counts
# the metadata are the station parameters, such as location, instrument type and
# instrument response
res = simulate_seismometer(st[0].data.astype("float32"), st[0].stats.sampling_rate, paz, one_hertz)
# Correct for overall sensitivity, nm/s
# this syntax *= means multiply res by the number that follows
# the operator += would add a number to res
res *= 1e9 / paz["sensitivity"]
#3 Apply lowpass at 10Hz
res = lowpass(res, 10, df=st[0].stats.sampling_rate, corners=4)
#4 Plot the seismograms
# sec is the time ordinate, a vector with length equal to the data
sec = np.arange(len(res))/st[0].stats.sampling_rate
plt.subplot(211)
plt.plot(sec,st[0].data, "k")
plt.title("%s %s" % ("RJOB",t))
plt.ylabel("STS-2")
plt.subplot(212)
plt.plot(sec,res, "k")
plt.xlabel("Time [s]")
plt.ylabel("1Hz CornerFrequency")
plt.show()

Give this script a try, maybe by feeding it line by line into python or ipython on your computer.  See if you get the comparison plot that matches Figure 1 of Beyreuther et al (2010)