Fast and efficient computing in Python using generators

Generators don’t hold the entire result in memory. It yields one result at a time.

Ways of creating generators:

Using a function

def squares_gen(num):
        for i in num:
                yield i**2

def squares(num):
                results=[]
                for i in num:
                        results.append(i**2)
                return results 

• Elapsed time for list: 7.360722 Seconds
• Elapsed time for generators: 5.999999999950489e-06 Seconds
• Difference in time taken for the list and generators: 7.360716 Seconds for num = np.arange(1,10000000)

Using comprehension

resl = [i**2 for i in num]

resg = (i**2 for i in num)

• Elapsed time for list: 7.663468000000001 Seconds
• Elapsed time for generators: 9.999999999621423e-06 Seconds
• Difference in time taken: 7.663458000000001 Seconds for num = np.arange(1,10000000)
  

Obtaining results from the generator object:

1. Using next
   

   resg = squares_gen(num)
   		print('res of generators: ',next(resg))
   		print('res of generators: ',next(resg))
   		print('res of generators: ',next(resg))
   		

   2.Using loop:
   

   for n in resg:
   		    print(n)
   		

Advantages of using generators:

1. The generator codes are more readable.
2. Generators are much faster and uses little memory.
   

Results:

1. Using function is a faster way of creating values in Python than using loop or list comprehension for both lists and generators.
2. The difference between using list or generators is more pronounced when using a comprehension (though generators are still much faster.)
3. When we need the result of whole array at a time then the amount of time (or memory) taken to create a list or list(generators) are almost same.

Overall, generators gives a performance boost not only in execution time but with the memory as well.

Appendix

How I calculated the time taken by the process

• Calculate sum of the system and user CPU time of the current process.
  
  • time.process_time provides the system and user CPU time of the current process in seconds.
  • Use time.process_time_ns to get the result in nanoseconds
    

NOTE: The “time taken” shown in this study is subjective to different computers and varies each time depending on the state of the CPU. But each and everytime, the using generators are much faster.

Tutorials on making web app using python (Introduction)

Recently, we have joined hackathons and tried making web applications using Python. We intend to make tutorials on making these web apps. In this introduction, we present what we have done.

I. Ferryman

The website: https://ferry-man.herokuapp.com/

Github source code: https://github.com/IEScoders/ObservingEarthWithData

More information: https://www.ait.org.tw/results-of-taipei-2018-nasa-international-space-apps-challenge-hackathon/

Skills/Packages have been used:

• Handling with NASA satellite data using read, write HDF and HDF5 file (pyhdf).

• Calculations with spatial data using pandas, numpy and multiprocessing.

• Plot data with matplotlib.

• Create the website using dash.

II. Biozone

The website: https://biozone.us-east.mybluemix.net/

Github source code: https://github.com/IEScoders/BioZone

Skills/Packages have been used:

• Crawl weather data from Central Weather Bureau (Taiwan)

• Calculations with spatial data using pandas and numpy.

• Plot interactive map using Folium

• Create the website using Flask and Javascript.

Plot geochemical data using python: An example of analyzing Adakite rock data from GEOROC

We download the precompiled data of adakite from GEOROC database.
For a simple impression of adakite, the wikipedia page gives some clue: Adakites are volcanic rocks of intermediate to felsic composition that have geochemical characteristics of magma thought to have formed by partial melting of altered basalt that is subducted below volcanic arcs.

In this example, we demonstrate how to use python to simplify the data, discard the null data, classify and plot the geochemical properties.

First, let’s look at the data, it is quite large and differs in the available data: some elements are there, some are not.

Now we can start by importing some useful packages:

import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import warnings
warnings.filterwarnings("ignore")
plt.rcParams['figure.figsize'] = (20, 10)
plt.style.use('default')

We use pandas to import the data, the encoding is used to solve some font conflicts.

df = pd.read_csv("ADAKITE.csv", encoding="iso-8859-1")

We use dropna() method to discard columns with less than 50% data available.

df = df.dropna(axis=1, thresh=int(0.5*len(df)))

To have a look at the sample occurrence, let’s plot the data. We select the lon, lat column, drop null value and change all data to numeric format.

plt.figure(figsize=(10,10))
m = Basemap(lon_0=180,projection='hammer')
lon = df["LONGITUDE MIN"].dropna()
lat = df["LATITUDE MIN"].dropna()
lon = pd.to_numeric(lon, errors='ignore');
lat = pd.to_numeric(lat, errors='ignore');
lon_ = [];lat_ = []
for x, y in zip(lon,lat):
    try:
        xx, yy = m(float(x),float(y))
        lon_.append(xx);lat_.append(yy)
    except:
        pass
m.scatter(lon_, lat_, marker = "o" ,s=15, c="r" , edgecolors = "k", alpha = 1)
m.drawcoastlines()
plt.title('Adakite rocks sample')
plt.show()

We make a function called plot_harker() to plot Harker’s diagram:

def plot_harker(x,xlabel,y,ylabel,title=None,xlim=[40,80],ylim=None,color = "b",label=None):
    plt.scatter(x=x,y=y,marker="o", c=color,s=8,label = label)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.xlim(xlim)
    try:
        plt.ylim(ylim)
    except:
        pass
    if title != None:
        plt.title(title)

… and then plot some elements using that function:

plt.figure(figsize=(12,12))
plt.subplot(321)
plot_harker(x=df["SIO2(WT%)"],xlabel=r'$SiO_2$ (wt%)',
            y=df["AL2O3(WT%)"],ylabel=(r'$Al_2O_3$ (wt%)'),title=r'$SiO_2$ vs $Al_2O_3$')
plt.subplot(322)
plot_harker(x=df["SIO2(WT%)"],xlabel=r'$SiO_2$ (wt%)',
            y=df["MGO(WT%)"],ylabel=(r'$MgO$ (wt%)'),title=r'$SiO_2$ vs $MgO$')
plt.subplot(323)
plot_harker(x=df["SIO2(WT%)"],xlabel=r'$SiO_2$ (wt%)',
            y=df["FEOT(WT%)"],ylabel=(r'$FeOt$ (wt%)'),title=r'$SiO_2$ vs $FeOt$')
plt.subplot(324)
plot_harker(x=df["SIO2(WT%)"],xlabel=r'$SiO_2$ (wt%)',
            y=df["TIO2(WT%)"],ylabel=(r'$TiO_2$ (wt%)'),title=r'$SiO_2$ vs $TiO_2$')
plt.subplot(325)
plot_harker(x=df["SIO2(WT%)"],xlabel=r'$SiO_2$ (wt%)',
            y=df["NA2O(WT%)"],ylabel=(r'$Na_2O$ (wt%)'),title=r'$SiO_2$ vs $Na_2O$')
plt.subplot(326)
plot_harker(x=df["SIO2(WT%)"],xlabel=r'$SiO_2$ (wt%)',
            y=df["K2O(WT%)"],ylabel=(r'$MgO$ (wt%)'),title=r'$SiO_2$ vs $K_2O$')
plt.suptitle(r'Harker diagram of Adakite vs $SiO_2$',y=0.92,fontsize=15)
plt.subplots_adjust(hspace=0.3)
plt.show()

We can try to see the tectonic settings of the rock:

plt.figure(figsize=(8,8))
tec = df['TECTONIC SETTING'].dropna()
tec = tec.replace('ARCHEAN CRATON (INCLUDING GREENSTONE BELTS)','ARCHEAN CRATON')
tec_counts = tec.value_counts()
tec_counts.plot(kind="bar",fontsize=10)
plt.title('Tectonic settings of Adakite')
plt.ylim([0,500])
plt.show()

The following code demonstrates how to create new columns and divide the data. We divide the data in High Silica Adakite (SiO₂ > 60%) and Low Silica Adakite (SiO₂ < 60%)

df['SR/Y'] = (df["SR(PPM)"]/df["Y(PPM)"])
df['CAO+NA2O'] = df['CAO(WT%)'] + df['NA2O(WT%)']
df['CR/NI'] = df['CR(PPM)'] + df['NI(PPM)']
df_hsa = df[df["SIO2(WT%)"] > 60]
df_lsa = df[df["SIO2(WT%)"] < 60]

Bonus: Let’s see the publish about adakite in the GEOROC database by year:

cite = [df.CITATIONS[x] for x in range(0,len(df)) if len(df.CITATIONS[x]) > 20 and df.CITATIONS[x].count('[') < 3]
year = []
for i in range(0,len(cite)):
    year.append(int(cite[i].split('[')[2].split(']')[0]))

The python code below is the full program of this example:

Loading

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw Adakite.ipynb hosted with ❤ by GitHub

Nguyen Cong Nghia
IESAS

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
warnings.filterwarnings("ignore")

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)

plt.figure(figsize=(10,5))
# 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 
    m.drawgreatcircle(stlon,stlat,evlon,evlat,linewidth=1,color='b')
    xx,yy = m(stlon,stlat)
    m.scatter(xx, yy, marker = "^" ,s=150, c="g" , edgecolors = "k", alpha = 1)
    
m.drawcoastlines()
#m.fillcontinents()
m.drawparallels(np.arange(-90,90,20),labels=[1,1,0,1])
plt.title("Event-station map")
plt.show()

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

plt.figure(figsize=(10,5))
for tr in stream:
    tr.normalize()
    dist = tr.stats.sac.dist
    plt.plot(tr.times(),tr.data+dist*0.01,c="k",linewidth=0.5)
    plt.scatter(tr.stats.sac.t3,dist*0.01,marker="|",color="r")
plt.ylabel("x100 km")    
plt.ylim(84,77)
plt.xlim(650,800)
plt.show()

Plot the same seismogram but using filled colors (which is more suitable to plot other kind of seismic traces)

plt.figure(figsize=(10,5))
for tr in stream:
    tr.normalize()
    dist = tr.stats.sac.dist*0.01
    x = tr.times()
    y = tr.data+dist
    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")    
plt.ylim(84,77)
plt.xlim(650,800)
plt.show()

Nguyen Cong Nghia

IESAS

Time-series Analysis using Python I

Introduction

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.

Prerequisites

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.

Data

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"])
df.head()

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.

df['Year'].head()
df[['Year','Lat']].head()

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

df.loc[:,"Year"].head()
df.iloc[:,3].head()

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”).

Analysis

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)
df.head()

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

df.plot()

We can also customize the plot easily.

df.plot(y=['dN',"dE","dU"],grid=True)
plt.ylabel("Amplitude")
plt.suptitle("GPS Data Visualization")
plt.title("0498")

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

plt.savefig('0498Data.pdf',dpi=300,bbox_inches='tight')

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 https://www.python.org/downloads/.

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.

Simple 1D velocity model inversion from P arrival time

Refer to Chapter 5, Introduction to Seismology, Shearer 2009.

Problem:

From the P-wave travel time data below (note that the reduction velocity of 8km/s), inverse for the 1D velocity model using T(X) curve fitting (fit the T(X) curve with lines, then invert for the ray parameter p and delay time τ(p), then solve for the velocity and depth).

 

Fig 1. P-wave travel time from a seismic experiment.

 

Solution:

We construct the T(X) plot from the data: y = T  – X/8 => T = y + X/8 and then divide the T(X) curve into fitting lines (Fig. 2).

Fig 2. T(X) plot with lines fitting to data points.

The ray parameter and delay is calculated for each line by using:

p = dT/dX

τ = T – pX

These parameters p and τ are slope and intercept of a line, respectively.

The τ(p) can be written as:

as in the matrix form of:

or, in short,

We can invert for the thickness of each layer using least square inversion:

If we assume ray slowness u(i) = p(i) for each line, we can solve the h and v for the 1D velocity model.

Here is the result of the inversion:

The sharp increase in Vp from ~6.6 km/s to 7.7 km/s at ~ 30km depth suggest for the Moho discontinuity.

The codes for the problem is written in python and can be downloaded: source code, τ inversion code, data points.

Nguyen Cong Nghia, IESAS

 

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:

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.

The vertical slowness is defined as:

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

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?

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 main program is as follow:

from ttcal import layerxt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

v_model = pd.read_csv('v_model.dat')

p1 = 0.1236
p2 = 0.2217
n = 100

result = [['p','X','T']]
for p in np.linspace(p1, p2, n):
    X = 0; T = 0;
    for i in range(0, len(v_model)-1):
            p = p
            h = v_model.h[i+1] - v_model.h[i]
            utop = 1/v_model.vp[i]
            ubot = 1/v_model.vp[i+1]
            dx, dt, irtr = layerxt(p,h,utop,ubot)
            X += dx; T += dt
        #If the ray reflected at layer the dx, dt calculation stops
            if irtr == 2:
                    break
        result.append([p,X,T])

result_df = pd.DataFrame(result, columns=result.pop(0))
result_df = result_df[::-1]

#Multiply by 2 to get total surface-to-surface value ot X(p) and T(p)
result_df['T'] = 2*result_df['T']
result_df['X'] = 2*result_df['X']

plt.figure(1)
result_df['RT'] = result_df['T']-result_df['X']/8
result_df.plot(x='X', y='RT', legend=False)
plt.title('Reduction time')
plt.ylabel('Reduce time T - X/8 (s)')
plt.xlabel('Distance (km)')
plt.text(20,0.80,'Prograde', rotation =40)
plt.text(40,1.6,'Retrograde', rotation =35)
plt.text(80,1.00,'Prograde', rotation =0)
plt.savefig('ex4.8.1.png')

plt.figure(2)
result_df.plot(x='p', y='X', legend=False)
plt.title('X(p)')
plt.ylabel('Distance X (km)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.2.png')

plt.figure(3)
result_df['Taup'] = result_df['T'] - result_df['p']*result_df['X']
result_df.plot(x='p', y='Taup', legend=False)
plt.title('Tau(p)')
plt.ylabel('Tau(p) (s)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.3.png')

plt.figure(4)
v_model.h = v_model.h * (-1)
v_model.plot(x='vp', y='h', legend=False)
plt.ylabel('Depth (km)')
plt.xlabel('Velocity (km/s)')
plt.xlim([4,9])
plt.ylim([-10,0])
plt.savefig('ex4.8.4.png')

#plt.show()

LAYERXT subroutine written in python:

from math import sqrt, log


def layerxt(p,h,utop,ubot):
    ''' LAYERXT calculates dx and dt for ray in layer with linear velocity gradient
 
        Inputs:   p     =  horizontal slowness
             h     =  layer thickness
             utop  =  slowness at top of layer
             ubot  =  slowness at bottom of layer
        Returns:  dx    =  range offset
             dt    =  travel time
             irtr  =  return code
                   = -1,  zero thickness layer
                   =  0,  ray turned above layer
                   =  1,  ray passed through layer
                   =  2,  ray turned within layer, 1 segment counted in dx,dt'''
# ray turned above layer
    if p >= utop:
        dx = 0
        dt = 0
        irtr = 0
        return dx,dt,irtr

#Zero thickness layer
    elif h == 0:
        dx = 0
        dt = 0
        irtr = -1
        return dx,dt,irtr

#Calculate some parameters
    u1 = utop
    u2 = ubot
    v1 = 1.0/u1
    v2 = 1.0/u2
#Slope of velocity gradient
    b = (v2 - v1)/h
    eta1 = sqrt(u1**2 - p**2)

#Constant velocity layer
    if b == 0:
        dx = h*p/eta1
        dt = h*u1**2/eta1
        irtr = 1
        return dx,dt,irtr

    x1 = eta1/(u1*b*p)
    tau1=(log((u1+eta1)/p)-eta1/u1)/b

#Ray turned within layer, no contribution to integral from bottom point
    if p >= ubot:
        dx = x1
        dtau = tau1
        dt = dtau + p*dx
        irtr = 2    
        return dx,dt,irtr

        irtr=1
    eta2=sqrt(u2**2-p**2)
    x2=eta2/(u2*b*p)
    tau2=(log((u2+eta2)/p)-eta2/u2)/b

    dx=x1-x2
    dtau=tau1-tau2

    dt=dtau+p*dx
    return dx,dt,irtr

def flat(zsph,vsph):
    ''' Calculate the flatten transformation from a spherical Earth.
        zf, vf = flat(zsph,vsph)'''
    
# Radius of the Earth
    a = 6371.0

    rsph = a - zsph
    zf = -a*log(rsph/float(a))
    vf = (a/float(rsph))*vsph
    return zf, vf

Input velocity model:

h,vp,vs,den
0.0,4.5,2.4,2
1.5,6.8,3.75,2.8
6.0,7.0,3.5,2.9
6.5,8.0,4.6,3.1
10.0,8.1,4.7,3.1

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, https://jeffknupp.com/blog/2014/06/18/improve-your-python-python-classes-and-object-oriented-programming/. 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.

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).