First let's set things up for business.
In [1]:
import pmagpy.ipmag as ipmag
import pmagpy.pmag as pmag
import pmagpy.pmagplotlib as pmagplotlib
from IPython.display import Image
import pandas as pd
import numpy as np
import matplotlib
%matplotlib inline
import pylab as plt
Having downloaded the Tauxe and Hartl (1997) contribution from http://earthref.org/doi/10.1111/j.1365-246X.1997.tb04082.x we can 'unpack' it using the command line download_magic.py program. But first, let's find out how to use it...
In [2]:
help(ipmag.download_magic)
So we just specify the file name and input and output directories, making a destination directory (Chapter_15/Problem_1) first.
In [5]:
ipmag.download_magic('magic_contribution_14132.txt',dir_path='Chapter_15/Problem_1',\
input_dir_path='Chapter_15')
Out[5]:
Read in the two data files specimens.txt and sites.txt to see what's in there.
In [7]:
specimens=pd.read_csv('Chapter_15/Problem_1/specimens.txt',sep='\t',header=1)
specimens.head()
Out[7]:
In [9]:
sites=pd.read_csv('Chapter_15/Problem_1/sites.txt',sep='\t',header=1)
sites.head()
Out[9]:
Let's do things this way:
In [12]:
help(pmag.pinc)
In [13]:
GAD_inc=pmag.pinc(26)
In [14]:
specimens=specimens.dropna(subset=['analysts']) # kill unwanted lines with duplicate or irrelevent info
specimens['site']=specimens['specimen'] # make a column with site name
data=pd.merge(specimens,sites,on='site') # merge the two data frames on site
In [15]:
plt.figure(1,(10,4))
plt.plot(data['core_depth'],data['dir_inc'],'b-')
plt.plot(data['core_depth'],data['dir_inc'],'ro',markeredgecolor='black') # plot the data
plt.axhline(0,color='black') # decorate with zero line
plt.axhline(GAD_inc,color='green',linestyle='dashed',linewidth=2) #GAD inclinations
plt.axhline(-GAD_inc,color='green',linestyle='dashed',linewidth=2)
plt.xlabel('Depth (m)')
plt.ylabel('Inclination')
plt.ylim(-90,90); # set bounds for inclination
Hmm. The inlinations look a little steeper than expected from GAD. Oh, I looked up where it was and googled around about Site 522 and found this paper by Tauxe et al. (1983) [https://doi.org/10.1111/j.1365-246X.1983.tb03318.x] which says that it is the northward motion of the African plate. Makes sense.
We are asked to plot inclination versus age and plot the magnetic timescale for the Oligocene. I went to the stratigraphy.org website and downloaded this neat figure from which I learned that the Oligocene lasted from about 34 to 23 Ma.
In [16]:
Image(filename='ChronostratChart2017-02.jpg')
Out[16]:
So now I look in my data DataFrame to find inclination (dir_inc) and age and to plot inclination versus age and the time scale from 23 to 34 Ma.
In [22]:
fig=plt.figure(1,(6,12))
ax=fig.add_subplot(121)
agemin,agemax=0,10
pmagplotlib.plot_ts(ax,23,34,timescale='gts12')
ax1=fig.add_subplot(122)
plt.plot(data.dir_inc,data.age,'b-')
plt.plot(data.dir_inc,data.age,'ro',markeredgecolor='black')
plt.ylim(35,23)
plt.ylabel('Inclination');
One more little detail - looks like the ages aren't exactly in order, so let us sort the DataFrame by age and replot.
In [23]:
data.sort_values(by='age',inplace=True)
fig=plt.figure(1,(6,12))
ax=fig.add_subplot(121)
agemin,agemax=0,10
pmagplotlib.plot_ts(ax,23,34,timescale='gts12')
ax1=fig.add_subplot(122)
plt.plot(data.dir_inc,data.age,'b-')
plt.plot(data.dir_inc,data.age,'ro',markeredgecolor='black')
plt.ylim(35,23)
plt.ylabel('Inclination');
So... it does look like the ages are offset a bit between the author assigned ages and the GTS12 and in fact as this paper was published in 1997, it is impossible that they used GTS12 (published in 2012, duh). So let's replot this using the 'ck95' time scale.
In [24]:
fig=plt.figure(1,(6,12))
ax=fig.add_subplot(121)
agemin,agemax=0,10
pmagplotlib.plot_ts(ax,23,34,timescale='ck95')
ax1=fig.add_subplot(122)
plt.plot(data.dir_inc,data.age,'b-')
plt.plot(data.dir_inc,data.age,'ro',markeredgecolor='black')
plt.ylim(35,23)
plt.ylabel('Inclination');
And yes, the CK95 timescale works better. And as the paper was published in 1997, a 1995 time scale would have been the standard.
We already have the data in a Pandas Dataframe (data) and can just plot the 'age' column' against the 'core_depth' column. Because core_depth increases down in the core, we also have to reverse the sense of the y-axis to make a plot that is geologically sensible.
In [25]:
plt.plot(data['age'].values,data['core_depth'],'k-')
plt.ylabel('Depth below sea floor (m)')
plt.xlabel('Age (Ma)')
plt.ylim(data['core_depth'].max()+1,data['core_depth'].min()-1);
To get average sedimentation rate, I need to subtract the top and bottom meter levels and divide by the top and bottom age.
In [76]:
DepthRange=data.core_depth.max()-data.core_depth.min()
AgeRange=data.age.max()-data.age.min()
print ('Average sediment accumulation rate was: ','%7.1f'%(DepthRange/AgeRange), 'm/myr')
It gets higher down core.
When the ridge crest is shallow, it is much above the Carbonate Compensation Depth (CCD). As the crust ages, it sinks and carbonate dissolution increases. Also, the CCD was higher in the Miocene in general, so there are probably climatic reasons as well.
So this is a tricky question. If we use inclination only, we can use pmag.doincfish( ) (from Chapter 11, Problem 2) to calculate the average inclination.
In [80]:
help(pmag.doincfish)
So we need to send pmag.doincfish( ) an array of the inclination values.
In [81]:
pmag.doincfish(data.dir_inc.values)
Out[81]:
But from the description column we glimpsed Problem 1b it appears that someone oriented the declination data....
So we can use the declination,inclination data and run our old favorite, pmag.doprinc( ) from Chapter 12. (We can't use ipmag.fisher_mean( ) without separating data by polarity and I'm feeling lazy.)
In [82]:
help(pmag.doprinc)
In [84]:
DIs=np.array([data['dir_dec'].values,data['dir_inc'].values]).transpose()
pmag.doprinc(DIs)
Out[84]:
So the principle inclination is 53.5 and the incfish estimate was 52.7. Not bad. But now we need to calculate the paleolatitude from the inclination - the opposite of what we did in Problem 1b.
In [86]:
inc=np.radians(53) # remember to change to radians.
print ('Paleolatitude for inclination of 53: ','%7.1f'%(np.degrees(np.arctan(np.tan(inc)/2.))))
It is now at 23(S) and this was presumably 33.6 South, so it was heading North. I know that it was always in the southern hemisphere because if you look at the magstrat - it is the negative inclination intervals that correlate with the normal polarities!