In [2]:
import numpy as np
import pandas as pd
deg2rad=np.pi/180. # converts degrees to radians
Let's write a little function to do the conversion.
In [3]:
def dir2cart(data):
decs,incs,ints=data[0]*deg2rad,data[1]*deg2rad,data[2]
X=ints*np.cos(decs)*np.cos(incs)
Y=ints*np.sin(decs)*np.cos(incs)
Z=ints*np.sin(incs)
cart=np.array([X,Y,Z]).transpose()
return cart
Now let's read in a data file with some geomagnetic field vectors in it.
In [4]:
# read in the data and transpose it to rows of dec, inc, int
data=np.loadtxt('../Chapter_2/ps2_prob1_data.txt').transpose()
print dir2cart(data)
First we have to understand how the function pmag.get_unf() works. To do this, we need to tell the notebook where the pmag module lives, import it and print out the doc string for get_unf():
In [5]:
import sys
sys.path.insert(0,'/Users/ltauxe/PmagPy')
import pmag
print pmag.get_unf.__doc__
Now we can use that function to generate a list of random points on the Earth's surface.
In [6]:
places=pmag.get_unf(10)
print places
Now let's find out about ipmag.igrf
In [7]:
import ipmag
print ipmag.igrf.__doc__
And now we can ship the data in places to doigrf.
In [8]:
for place in places:
print ipmag.igrf([2006,0,place[1],place[0]])
In [9]:
data=[] # make a blank list
for place in places:
Dir=ipmag.igrf([2006,0,place[1],place[0]])
data.append(Dir) # append to the data list
data=np.array(data).transpose() # dir2cart takes arrays of data
print dir2cart(data)
In [10]:
import pmagplotlib
import matplotlib.pyplot as plt
%matplotlib inline
In [11]:
ipmag.plot_net(1) # make an equal angle net
ipmag.plot_di(data[0],data[1]) # put on the dots
In [12]:
lat = 36.*deg2rad # remember to convert to radians!
inc = np.arctan(2.*np.tan(lat)) /deg2rad # and back!
print '%7.1f'%(inc) # and print it out
Let's use the pmag function dia_vgp. First let's figure out what it does:
In [13]:
print pmag.dia_vgp.__doc__
Now we can use it to convert our directions to VGPs. Note that alpha95 is require but is not given here so we supply a zero in its place. Note also that westward longitudes are indicated by minus signs
In [14]:
vgp_lat,vgp_lon,dp,dp= pmag.dia_vgp(345,47,0.,36,-112)
print '%7.1f %7.1f'%(vgp_lat,vgp_lon)