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__

f

__or 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

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)

## No comments:

Post a Comment