Plot seismogram (SAC file), events, stations in Python (Part 1)

Here is an example of plotting SAC files in Python. The sample SAC files can be downloaded here and the Jupyter notebook can be downloaded here.

First, import some useful packages, including obspy, pandas, numpy and Basemap. By the way, they are all great packages (obspy is amazing for anyone who uses seismic data)

from obspy import read
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

#Ignore warnings due to python 2 and 3 conflict
import warnings

Let’s read the sample Z component using read from obspy

stream = read("2015.0529.0700/*Z.sac")

In here, we use the header from SAC file using tr.stats.sac.(SAC_header)

# setup mercator map projection.
m = Basemap(lon_0=180,projection='hammer')
evlat = stream[0].stats.sac.evla; evlon = stream[0].stats.sac.evlo

#Plot the event
xx,yy = m(evlon,evlat)
m.scatter(xx, yy, marker = "*" ,s=150, c="r" , edgecolors = "k", alpha = 1)

for tr in stream:
    stlat = tr.stats.sac.stla; stlon = tr.stats.sac.stlo 
    xx,yy = m(stlon,stlat)
    m.scatter(xx, yy, marker = "^" ,s=150, c="g" , edgecolors = "k", alpha = 1)
plt.title("Event-station map")

I used a simple trick to plot the seismogram with distance by make the y:
y = data + dist*weight_factor
with data is the amplitude of seismic trace, dist: distance in km (SAC header) and weight_factor = 0.01

The red line indicate the predicted P arrival time that I have calculated and store in SAC header t3

for tr in stream:
    dist = tr.stats.sac.dist
plt.ylabel("x100 km")    
Plot the same seismogram but using filled colors (which is more suitable to plot other kind of seismic traces)
for tr in stream:
    dist = tr.stats.sac.dist*0.01
    x = tr.times()
    y =
    plt.fill_between(x,y, dist, y > dist, color='r', alpha = 0.8)
    plt.fill_between(x,y, dist, y < dist, color='b', alpha = 0.8)
plt.ylabel("x100 km")    
Nguyen Cong Nghia

Time-series Analysis using Python I


Time-series analysis is essential in most fields of science including geophysics, economics, etc. Most of the geophysical data comes in a time-series format including the seismic recordings. In this part of the series of tutorial, we will see how we can quickly load the data, and visualize it.


This tutorial does not require the reader to have any basic understanding of Python or any programming language. But we expect the reader to have installed the jupyter notebook on their system. If the reader has not installed it yet, then they can follow the previous post where we went through the steps involved in getting started with Python.

What is Time-series?

Time-series is a collection of data at fixed time intervals. This can be analyzed to obtain long-term trends, statistics, and many other sorts of inferences depending on the subject.


We also need some data to undergo the analysis. We demonstrate the analysis using our GPS data. It can be downloaded from here.

Let’s get started

The first step is always to start the Python interpreter. In our case, we will use the jupyter notebook.

Jupyter notebook can be started using the terminal. Firstly, navigate to your directory containing the data and the type “jupyter notebook” on your terminal.

jupyter notebook

Next, we create a new Python 3 notebook, rename it as pythontut1. Then, we need to import some of the libraries:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib.pyplot import rcParams
rcParams['figure.figsize'] = 15, 6

Loading the data

Now, we load the data using the pandas library functions. Here, we use the function read_csv. But, before that let’s observe the format of the data:

!head 0498.COR

The prefix “!” can be used to execute any Linux command in the notebook.

We can see that the data has no header information, and 8 columns. The columns are namely, “year”, “latitude”, “longitude”, “Height”, “dN”, “dE”, “dU”.

So, now we read the data and set the above names to the different columns.

df=pd.read_csv("0498.COR", header=None, delimiter='\s+', names=['Year',"Lat","Long","Hgt","dN","dE","dU","nav"])

It is essential to understand the above command. We gave the argument of the filename, header (default is the first line), delimiter (default is a comma) and the names of each column, respectively. Then we output the first 5 lines of the data using the df.head() command.

Our data is now loaded, and if we want to extract any section of the data, we can easily do that.


Here, we have used the column names to extract the two columns only. We can also use the index values.


When we use “.loc” method to extract the section of the data, then we need to use the column name whereas when we use the “.iloc” method then we use the index values. Here, df.iloc[:,3] extracts all the rows of the 3rd column (“Hgt”).


Now, we have the data loaded. Let’s plot the “dN”, “dE”, and “dU” values versus the year. Before doing that, let’s set the “Year” column as the index column.

df.set_index("Year", inplace=True)

We can see the “Year” column as the index of the data frame now. Plotting using Pandas is extremely easy.


We can also customize the plot easily.

plt.suptitle("GPS Data Visualization")

If we want to save the figure, then we can use the command:


This saves the figure as the pdf file named “0498Data.pdf”. The format can be set to any type “.png”, “.jpg”, ‘.eps”, etc. We set the resolution to be 300 dpi. This can be varied depending on our need. Lastly, “bbox_inches =‘tight’” crops our figure to remove all the unnecessary space.

Next Tutorial

We have loaded the data and visualized it. But we can see that our data has some trend and seasonality. In the future tutorial, we will learn how to remove that.

Python: a language for science

Python has become the most popular language for programming and its community is huge, active and ever-increasing. Before Python, MATLAB was the preferred language for the scientific community but now most of the researchers are migrating to Python. Like MATLAB, Python is quite easy to use but over that Python can be used for the system control, software developing, easy integration with other languages like FORTRAN, C , etc. This makes it extremely popular among all sorts of communities. Some P.Is want their students to work in Python because according to them, even when the students want to switch the field or simply leave the field, it could still be useful to them.

If you are familiar with other programming languages including MATLAB, then its very easy to switch to Python. Though, there are some issues the beginners and the MATLAB users face while making a move to Python. Here, let us deal with that issue and get started with Python.

Getting started with Python

Currently, there are two versions of Python – 2.7, and 3.6. But we will stick with Python 3.6. The steps for the version 2.7 is very similar.

Python can be downloaded directly from its official website

But, here, we will download Python from the Anaconda website. Steps are listed below:

  1. Download the installer for the Windows, Mac, or the Linux system. For Windows, you will get a “.exe” extension file, double-click the file and follow the installation procedure.For Mac and Linux system, the file will be with the extension “.sh” (a shell- script). You can open your favorite terminal, navigate to the directory containing the downloaded Miniconda installation file. Then, you can simply type “sh Miniconda*.sh” and follow the steps prompted. You will be asked to accept the license and enter your installation location. If you don’t have any preferred location then the installer will install Miniconda in your home directory.

  2. Now, Miniconda, as well as Python 3.6, has been installed on your machine. We need some more Python libraries to get going. We can easily install these now using Miniconda. Open your terminal and type “conda install jupyter notebook pandas numpy matplotlib”.

We now have all the necessary libraries to get started with Python. Let us open a terminal again and simply type “jupyter notebook” to launch jupyter notebook. Jupyter notebook is a very popular interface which makes using Python a fun experience. It has all the features to make things simpler and useful. You can even run the shell command in the notebook. Over that, you can write your notes using the Markdown language, which is a very straight-forward language for writing. I strongly recommend you to learn Markdown. It won’t take more than 10 minutes to learn and will be useful in so many other places.

Now, that we have all the necessary tools to use Python then let’s use python to make a very simple plot.

  1. First, start jupyter notebook by typing “jupyter notebook” on your terminal.It will open “Jupyter” in your favorite browser.
  2. Now, go to the “New” tab at the top right corner and then select “Python 3”. It will open a new page for you.
  3. Let’s import some libraries which we have installed using the conda by typing the command as directed in the figure below. To execute the command, use the “Command/Ctrl + Enter” key.
  4. We also need to execute the command “%matplotlib inline” to keep our plots inline in the notebook. This is one of the many features of jupyter notebook.
  5. Let’s define our x and y values using the numpy library which we have imported as “np” for short-form.
  6. Let’s plot our data.

Writing NetCDF4 Data using Python

For how to read a netCDF data, please refer to the previous post. Also, check the package and tools required for writing the netCDF data, check the page for reading the netCDF data.

Importing relevant libraries

import netCDF4 
import numpy as np

Screen Shot 2017-10-03 at 2.20.50 PM.png

Let us create a new empty netCDF file named “” in the “../../data” directory and open it for writing.

ncfile = netCDF4.Dataset('../../data/',mode='w',format='NETCDF4_CLASSIC') 

Screen Shot 2017-10-03 at 2.30.59 PM.png

Notice here that we have set the mode to be “w”, which means write mode. We can also open the data in append mode (“a”). It is safe to check whether the netCDF file has closed, using the try and except statement.

Creating Dimensions

We can now fill the netCDF files opened with the dimensions, variables, and attributes. First of all, let’s create the dimension.

lat_dim = ncfile.createDimension('lat', 73) # latitude axis
lon_dim = ncfile.createDimension('lon', 144) # longitude axis
time_dim = ncfile.createDimension('time', None) # unlimited axis (can be appended to).
for dim in ncfile.dimensions.items():

Screen Shot 2017-10-03 at 2.35.59 PM.png

Every dimension has a name and length. If we set the dimension length to be 0 or None, then it takes it as of unlimited size and can grow. Since we are following the netCDF classic format, only one dimension can be unlimited. To make more than one dimension to be unlimited follow the other format. Here, we will constrain to the classic format only as it is the simplest one.

Creating attributes

One of the nice features of netCDF data format is that we can also store the meta-data information along with the data. This information can be stored as attributes.

ncfile.title='My model data'

Screen Shot 2017-10-03 at 2.43.38 PM.png

ncfile.subtitle="My model data subtitle"
ncfile.anything="write anything"

Screen Shot 2017-10-03 at 2.45.55 PM.png

We can add as many attributes as we like.

Creating Variables

Now, let us add some variables to store some data in them. A variable has a name, a type, a shape and some data values. The shape of the variable can be stated using the tuple of the dimension names. The variable should also contain some attributes such as units to describe the data.

lat = ncfile.createVariable('lat', np.float32, ('lat',))
lat.units = 'degrees_north'
lat.long_name = 'latitude'
lon = ncfile.createVariable('lon', np.float32, ('lon',))
lon.units = 'degrees_east'
lon.long_name = 'longitude'
time = ncfile.createVariable('time', np.float64, ('time',))
time.units = 'hours since 1800-01-01'
time.long_name = 'time'
temp = ncfile.createVariable('temp',np.float64,('time','lat','lon')) # note: unlimited dimension is leftmost
temp.units = 'K' # degrees Kelvin
temp.standard_name = 'air_temperature' # this is a CF standard name

Screen Shot 2017-10-03 at 2.51.29 PM.png

Here, we create the variable using the createVariable method. This method takes 3 arguments: a variable name (string type), data types, a tuple containing the dimension. We have also added some attributes such as for the variable lat, we added the attribute of units and long_name. Also, notice the units of the time variable.

We also have defined the 3-dimensional variable “temp” which is dependent on the other variables time, lat and lon.

In addition to the custom attributes, the netCDF provides some pre-defined attributes as well.

print("-- Some pre-defined attributes for variable temp:")
print("temp.dimensions:", temp.dimensions)
print("temp.shape:", temp.shape)
print("temp.dtype:", temp.dtype)
print("temp.ndim:", temp.ndim) 

Screen Shot 2017-10-03 at 2.57.36 PM

Since no data has been added, the length of the time dimension is 0.

Writing Data

nlats = len(lat_dim); nlons = len(lon_dim); ntimes = 3
lat[:] = -90. + (180./nlats)*np.arange(nlats) # south pole to north pole
lon[:] = (180./nlats)*np.arange(nlons) # Greenwich meridian eastward
data_arr = np.random.uniform(low=280,high=330,size=(ntimes,nlats,nlons))
temp[:,:,:] = data_arr # Appends data along unlimited dimension
print("-- Wrote data, temp.shape is now ", temp.shape)
print("-- Min/Max values:", temp[:,:,:].min(), temp[:,:,:].max())

Screen Shot 2017-10-03 at 3.02.52 PM.png

The length of the lat and lon variable will be equal to its dimension. Since the length of the time variable is unlimited and is subject to grow, we can give it any size. We can treat netCDF array as a numpy array and add data to it. The above statement writes all the data at once, but we can do it iteratively as well.

Now, let’s add another time slice.

data_slice = np.random.uniform(low=280,high=330,size=(nlats,nlons))
temp[3,:,:] = data_slice 
print("-- Wrote more data, temp.shape is now ", temp.shape) 

Screen Shot 2017-10-03 at 3.10.20 PM.png

Note, that we haven’t added any data to the time variable yet.

times_arr = time[:]

Screen Shot 2017-10-03 at 3.12.50 PM.png

The dashes indicate that there is no data available. Also, notice the 4 dashes corresponding to the four levels in of the time stacks.

Now, let us write some data to the time variable using the datetime module of Python and the date2num function of netCDF4.

import datetime as dt
from netCDF4 import date2num,num2date
dates = [dt.datetime(2014,10,1,0),dt.datetime(2014,10,2,0),dt.datetime(2014,10,3,0),dt.datetime(2014,10,4,0)]

Screen Shot 2017-10-03 at 3.17.16 PM.png

times = date2num(dates, time.units)
print(times, time.units) # numeric values

Screen Shot 2017-10-03 at 3.18.53 PM.png

Now, it’s important to close the netCDF file which has been opened previously. This flushes buffers to make sure all the data gets written. It also releases the memory resources used by the netCDF file.

# first print the Dataset object to see what we've got
# close the Dataset.
ncfile.close(); print('Dataset is closed!')

Screen Shot 2017-10-03 at 3.23.38 PM.png


Earthquake Finder

This program can be used to obtain the earthquake information from various sources. The user can search for any event for given time range, magnitude range, depth range, geographical area. It can also obtain the focal mechanism solutions for the given parameters. It gives the output file (catalog.txt by default) where all the information is stored. It also plots the output (both with the focal mechanism and without it.)
This program makes use of the obspy module of Python and retrieve data from the International Seismological Center (ISC) server


to run with default parameters


to obtain the earthquake info between today, March 2016 to the current time

python3 st=2016/3

to obtain the earthquake info between 2016, march, 29 to 2016, september, 22

python3 st=2016/3/29,et=2016/9/22

to obtain the catalog for magnitudes between 4-9 and with focal mechanism

python3 mxmag=9,mnmag=4,fm=yes

to obtain catalog for magnitude 4-7, and within radius of 10 (default) and central coordinates 22(lat),121(lon) with focal mechanism

python3 mxmag=7,mnmag=4,clat=22,clon=121,mxrad=80,fm=yes

Parameters to change (default values in the braces):

mnla(-90),mxla(90),mnlo(-180),mxlo(180),mndep(0),mxdep(700),mnmag(4),mxmag(10),mnrad(0),mxrad(10),clat(None),clon(None),st(-1 month),et (current time),fm(no)

Output format of the catalog


With Focal Mechanism -> “YEAR MONTH DAY HOUR MIN SEC LONGITUDE LATITUDE DEPTH EXP(Nm) M0 MAG Mrr Mtt Mpp Mrt Mtp Mpr Str1 Dip1 Rake1 Str2 Dip2 Rake2.”


  1. Python 3: Can be obtained from here.
  2. Extra Modules
    (a) requests: pip install requests
    (b) pandas: pip install pandas
    (c) obspy: pip install obspy
    (d) basemap: It is a great tool for creating maps using python in a simple way. The best way to obtain basemap is via Anaconda/Miniconda

    For Linux/Mac users just download the bash installer and inside the download directory, execute the following command: bash Miniconda*.sh

    For Windows user, download the exe file and install it. After the installation of Anaconda/Miniconda, basemap can be obtained by simply executing the command conda install basemap

Click here
to download the programs.


Program for extracting the earthquake informations


Often in seismology, we require the information about the earthquake such as location, origin time, depth, magnitude, etc. To obtain this information, users generally, browse the web or if he/she is experienced then just go to the predetermined sites such as global CMT catalog. Even, after reaching the correct website, we need to do some dig up to get the desired information. This program is meant to scrape the global CMT website to obtain the earthquake information.

How to obtain the program

You can download the package from here.

How to use the program

The package you download from the given site consists of two python programs. These programs require pre-installed python-3.

You can obtain python from the following websites:
1. Python software
2. The user can also download the python package anaconda

After installing the python 3, to run the program, type:
Then enter the event recognition parameter such as
And enter. It will prompt the event information.

Screen Shot 2017-08-27 at 10.07.05 PM

Using the program as a module

You can also use the package as a python module. For using it as a module, save the program in the current directory or in the computer-discoverable path.
import EQinfo_extractor as eq_ext

Download Earthquake Catalogs from Global CMT website

In seismology, we always need to download and check the event information about the events. This python script can download the event catalog from the website to the local computer for the given range of time.

Running this program is very simple:
Screen Shot 2017-08-25 at 11.47.09 PM.png

The user just needs to input the time range for the earthquakes e.g., 2000/05-2009/08.


Requirements: Python 3

To download the program, please click here.

Ray tracing through a 1-D velocity model

Refer to Chapter 4 of Shearer, Introduction to Seismology.


For a ray piercing through Earth, the ray parameter (or horizontal slowness) p is defined by several expressions:

Screenshot from 2016-12-21 16-38-07.png

where u = 1/v is the slowness, θ is the ray incidence angle, T is the travel time, X is the horizontal range and utp is the slowness at the ray turning point.

Screenshot from 2016-12-21 16-37-05.png

The vertical slowness is defined as:


and integral expressions for the surface-to-surface travel time are:

Screenshot from 2016-12-21 16-41-26.png

Screenshot from 2016-12-21 16-41-32.png

With these equations we can calculate the travel time (T) and horizontal distance (X) for a given ray with ray parameter p and velocity model v(z).

Apply these to a problem (Exercise 4.8 Shearer):

(COMPUTER) Consider MARMOD, a velocity-versus-depth model, which is typical of much of the oceanic crust (Table 4.1). Linear velocity gradients are assumed to exist at intermediate depths in the model; for example, the P velocity at 3.75 km is 6.9 km/s. Write a computer program to trace rays through this model and produce a P-wave T(X) curve, using 100 values of the ray parameter p equally spaced between 0.1236 and 0.2217 s/km. You will find it helpful to use subroutine LAYERXT (provided in Fortran in Appendix D and in the supplemental web material as a Matlab script), which gives dx and dt as a function of p for layers with linear velocity gradients. Your program will involve an outer loop over ray parameter and an inner loop over depth in the model. For each ray, set x and t to zero and then, starting with the surface layer and proceeding downward, sum the contributions, dx and dt, from LAYERXT for each layer until the ray turns. This will give x and t for the ray from the surface to the turning point. Multiply by two to obtain the total surface-to-surface values of X(p) and T(p). Now produce plots of: (a) T(X) plotted with a reduction velocity of 8 km/s, (b) X(p), and (c) τ(p). On each plot, label the prograde and retrograde branches. Where might one anticipate that the largest amplitudes will occur?

Screenshot from 2016-12-21 21-32-25.png

Using the LAYERXT subroutine (Appendix D of Shearer), the FORTRAN and MATLAB can be downloaded from here. In this example, we will use python to calculate and plot the result.


The codes can be downloaded here: main program, python LAYERXT subroutine, velocity model.

You need to install numpy, pandas and matplotlib package to run this. Normally it can be installed in the terminal by (presume python has been installed):

pip install numpy pandas matplotlib

Nguyen Cong Nghia – IESAS




Introduction to Python Part II

I. Type of objects in Python:

In Python, every object has its own class – or type of data. The in-depth tutorial can be found on the web, for example, In this tutorial, I will introduce some basic type in Python.

To check type of a variable, data, you can use function type(variable)

+ Numbers: most frequently use is float and int type. The float type uses decimal while the int rounds number. Below is the example of using these type of number. To convert to int and float type, we use int(variable) and float(variable). A number can take numeric operations like +(add), – (subtract), *(multiply), /(divide), % (modulo), ** (exponential)


+ Strings: any character information, that in between ‘ ‘ or ” “.  Strings can be joined together by using + operation.


+ List: Python often uses compound data types, used to group together with other values. The most versatile is the list, which can be written as a list of comma-separated values (items) between square brackets. Lists might contain items of different types: numbers, string or list itself (a list of lists). This type is a basic type to analyze data in Python because it makes you able to access data with ordering.

II. Function and method

In Python,methods are associated with object instances or classes; functions aren’t. When Python dispatches (calls) a method, then it binds the first parameter of that call to the appropriate object reference. A function often has the form of function(argument) while an argument can be any kind of data type (number, string, list). A method must be associated with a type of object and often have the form of object.method(argument) with an object is the suitable type to do a method.

Let’s do some practice and take the list as an example. Screen Shot 2016-12-05 at 1.29.13 PM.png

The functions used in here are len (return the length of an object – how many objects in a list) and print (display an object on screen). The methods used in here are append(add an object to a list) and reverse (reverse the order of objects in a list).