Some of the common data conversion and analysis functionality of PmagPy is exposed within the GUI programs as we saw yesterday. However, the capabilities of PmagPy extend well beyond what is in the GUIs and can access within the command line programs and as functions that can imported as Python modules.
As an example, let's consider a situation where we want to draw random samples from a specified Fisher distribution. This can be done with the command line program fishrot.py. We can learn a little bit about fishrot.py in the PmagPy Cookbook. The command line programs take flags. Type this command at the terminal to learn about them:
fishrot.py -hFor example, if we want to sample 10 directions from a Fisher distribution with a mean declination of 45, a mean inclination of 30 and a kappa of 20, we would use this command in terminal:
fishrot.py -k 20 -n 10 -D 45 -I 30Let's look at the code of fishrot.py on Github: https://github.com/PmagPy/PmagPy/blob/master/programs/fishrot.py
We can see within the program that the program uses two functions within pmag.py (pmag.fshdev and pmag.dodirot). Let's look at the source of these functions which is taken from pmag.py into the code cell below:
In [1]:
def fshdev(k):
"""
Generate a random draw from a Fisher distribution with mean declination
of 0 and inclination of 90 with a specified kappa.
Parameters
----------
k : kappa (precision parameter) of the distribution
Returns
----------
dec, inc : declination and inclination of random Fisher distribution draw
"""
R1=random.random()
R2=random.random()
L=numpy.exp(-2*k)
a=R1*(1-L)+L
fac=numpy.sqrt((-numpy.log(a))/(2*k))
inc=90.-2*numpy.arcsin(fac)*180./numpy.pi
dec=2*numpy.pi*R2*180./numpy.pi
return dec,inc
def dodirot(D,I,Dbar,Ibar):
"""
Rotate a declination/inclination pair by the difference between dec=0 and
inc = 90 and the provided desired mean direction
Parameters
----------
D : declination to be rotated
I : inclination to be rotated
Dbar : declination of desired mean
Ibar : inclination of desired mean
Returns
----------
drot, irot : rotated declination and inclination
"""
d,irot=dogeo(D,I,Dbar,90.-Ibar)
drot=d-180.
if drot<360.:drot=drot+360.
if drot>360.:drot=drot-360.
return drot,irot
To use functions in the notebook we can import PmagPy function modules as in the code cell below. These function modules are: pmag, a module with functions for analyzing paleomagnetic and rock magnetic data that are used within the command lines programs, the GUIs and within ipmag, a module with functions that combine and extend pmag functions and exposes pmagplotlib functions in order to generate output that works well within the Jupyter notebook environment.
In [12]:
import pmagpy.pmag as pmag
import pmagpy.ipmag as ipmag
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Putting a question mark after a function results in the docstring being displayed.
In [13]:
ipmag.fishrot?
In [14]:
ipmag.fishrot??
Putting two questions marks displays the docstring and the source code.
In [15]:
ipmag.fishrot??
The docstring helps guide us to know what parameters to put into the function. Here let's take 10 samples from a distribution with mean of dec=45, inc= 30 and kappa precision parameter of 20. If simply call the function it will print to the console as in the code cell below. In the code cell below that, we assign the variable to be this this block of declination/inclination values which we call a di_block.
In [5]:
ipmag.fishrot(k=20, n=10, dec=45, inc=30)
Out[5]:
In [7]:
directions = ipmag.fishrot(k=20, n=10, dec=45, inc=30)
directions
Out[7]:
In [8]:
ipmag.plot_di()
Following the instructions in the documentation string, let's make a plot and plot the directions.
In [15]:
fignum = 1
plt.figure(num=fignum,figsize=(4,4),dpi=160)
ipmag.plot_net(fignum)
ipmag.plot_di(di_block=directions)
In [11]:
decs = [100,120,115]
incs = [24,26,10]
mean = ipmag.fisher_mean(dec=decs,inc=incs)
print(mean)
In [18]:
lots_of_directions = ipmag.fishrot(k=50, n=100, dec=45, inc=30)
mean = ipmag.fisher_mean(di_block=lots_of_directions)
print(mean)
In [19]:
ipmag.print_direction_mean(mean)
To access a single element in the dictionary, say the declination ('dec'), we would do it by putting the name of a dictionary and following it with the key related to the value we want to see. We will use this approach to extract and plot the mean dec, inc and $\alpha_{95}$ in the plot below.
In [20]:
mean['dec']
Out[20]:
Let's plot both the directions (using ipmag.plot_di) and the mean (using ipmag.plot_di_mean).
In [21]:
fignum = 1
plt.figure(num=fignum,figsize=(4,4),dpi=160)
ipmag.plot_net(fignum)
ipmag.plot_di(di_block=lots_of_directions)
ipmag.plot_di_mean(mean['dec'],mean['inc'],mean['alpha95'],color='r')
plt.show()
Let's use another ipmag function to determine the field predicted by the IGRF model here in La Jolla over the past 100 hundred years.
*Before we do so, let's all open up ipmag.py in a text editor. Find the function ipmag.igrf* and report back what functions it uses (hint they are from the pmag.py module).
To determine the IGRF field at a bunch of different years we will first create a list of years using the numpy function np.linspace. We can then take these years along with the igrf function and the location of La Jolla to calculate the local direction through time using a for loop.
In [22]:
pmag.doigrf??
In [25]:
years = np.arange(1900,2020,1)
In [ ]:
In [26]:
years
Out[26]:
In [36]:
field = ipmag.igrf([year,0,La_Jolla_lat,La_Jolla_lon])
field
Out[36]:
In [32]:
local_field = []
local_field_intensity = []
La_Jolla_lat = 32.8328
La_Jolla_lon = -117.2713
for year in years:
field = ipmag.igrf([year,0,La_Jolla_lat,La_Jolla_lon])
local_field.append(field)
local_field_intensity.append(field[2])
In [33]:
fignum = 1
plt.figure(num=fignum,figsize=(8,8),dpi=160)
ipmag.plot_net(fignum)
ipmag.plot_di(di_block=local_field,markersize=10,label='field since 1900 in La Jolla')
ipmag.plot_di(di_block=local_field[-2:],color='red',label='field today in La Jolla')
plt.legend()
plt.show()
In [34]:
plt.plot(years,local_field_intensity)
Out[34]:
In [ ]: