In [1]:
import numpy as np
import pandas as pd
Let's write a little function to do the conversion.
In [2]:
def dir2cart(data):
""" Converts data array with [Declination, Inclination, Intensity]
to cartesian coordinates, X=North, Y=East, Z=Down
Returns array with [X,Y,Z]
"""
# convert degrees to radians for declination and inclination
decs,incs,ints=np.radians(data[0]),np.radians(data[1]),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 [3]:
# 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))
To solve this problem, 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 [1]:
import pmagpy.pmag as pmag
help(pmag.get_unf)
Now we can use that function to generate a list of random points on the Earth's surface.
In [4]:
places=pmag.get_unf(10)
print (places)
Now let's find out about ipmag.igrf
In [5]:
import pmagpy.ipmag as ipmag
help(ipmag.igrf)
And now we can ship the data in places to ipmag.igrf.
In [6]:
for place in places:
print (ipmag.igrf([2006,0,place[1],place[0]]))
In [8]:
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 [9]:
import pmagpy.pmagplotlib as pmagplotlib
import matplotlib.pyplot as plt
# this 'magic command' (starts with %) let's us plot things in the notebook
%matplotlib inline
In [10]:
ipmag.plot_net(1) # make an equal angle net
ipmag.plot_di(data[0],data[1]) # put on the dots
In [11]:
lat = np.radians(36.) # remember to convert to radians!
inc = lambda lat: np.degrees(np.arctan(2.*np.tan(lat))) # and back!
print ('%7.1f'%(inc(lat))) # and print it out
Let's use the pmag function dia_vgp. First let's figure out what it does:
In [7]:
help(pmag.dia_vgp)
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 [13]:
vgp_lat,vgp_lon,dp,dp= pmag.dia_vgp(345,47,0.,36,-112)
print ('%7.1f %7.1f'%(vgp_lat,vgp_lon))