#import paleomagnetic specific libraries import pmag, pmagplotlib, IPmag import pandas pandas.set_option('display.max_columns', 500) from IPython.core.display import HTML from mpl_toolkits.basemap import Basemap %pylab inline data = pandas.read_csv('../2014_Osler_Data/SimpsonIsland_OslerData.csv') #write these data to a latex file for the supplemental PDF with open('OslerData.txt','w') as f: f.write(data.to_latex(index=False,columns=['Site','Strat_Height','Site_Lat','Site_Long','Dec_Geo','Inc_Geo','Dec_TC','Inc_TC','n','VGP_lat','VGP_long'])) display(HTML(data.to_html())) data_file='../2014_Osler_Data/SimpsonIsland_OslerData.csv' SimpsonIsland_Strat = np.genfromtxt(data_file,delimiter=",", skip_header=1) strat_height=SimpsonIsland_Strat[:,1] Dec_TC=SimpsonIsland_Strat[:,6] Inc_TC=SimpsonIsland_Strat[:,7] A95=SimpsonIsland_Strat[:,8] plat=SimpsonIsland_Strat[:,10] abs_plat=SimpsonIsland_Strat[:,11] VGP_lat=SimpsonIsland_Strat[:,12] VGP_long=SimpsonIsland_Strat[:,13] SI_Directions=[] SI_Poles=[] for n in range(0,84): Dec,Inc=Dec_TC[n],Inc_TC[n] SI_Directions.append([Dec,Inc,1.]) Plong,Plat=VGP_long[n],VGP_lat[n] SI_Poles.append([Plong,Plat,1.]) SI_LowerThird_Directions=SI_Directions[0:30] SI_LowerThird_Poles=SI_Poles[0:30] SI_MiddleThird_Directions=SI_Directions[30:50] SI_MiddleThird_Poles=SI_Poles[30:50] SI_UpperThird_Directions=SI_Directions[50:84] SI_UpperThird_Poles=SI_Poles[50:84] #calculate and display the Fisher means for each subset of the directions lower_third_mean=pmag.fisher_mean(SI_LowerThird_Directions) middle_third_mean=pmag.fisher_mean(SI_MiddleThird_Directions) upper_third_mean=pmag.fisher_mean(SI_UpperThird_Directions) print 'The Fisher mean parameters for SI_LowerThird_Directions are: ' print 'Dec = ' + str(lower_third_mean['dec']) + ' Inc = ' + str(lower_third_mean['inc']) print 'alpha95 = ' + str(lower_third_mean['alpha95']) + ' k= ' + str(lower_third_mean['k']) print '' print 'The Fisher mean parameters for SI_MiddleThird_Directions are: ' print 'Dec = ' + str(middle_third_mean['dec']) + ' Inc = ' + str(middle_third_mean['inc']) print 'alpha95 = ' + str(middle_third_mean['alpha95']) + ' k= ' + str(middle_third_mean['k']) print '' print 'The Fisher mean parameters for SI_UpperThird_Directions are: ' print 'Dec = ' + str(upper_third_mean['dec']) + ' Inc = ' + str(upper_third_mean['inc']) print 'alpha95 = ' + str(upper_third_mean['alpha95']) + ' k= ' + str(upper_third_mean['k']) #plot the direction of each flow mean on an equal area plot fignum = 1 pylab.figure(num=fignum,figsize=(10,10),dpi=160) pmagplotlib.plotNET(fignum) IPmag.iplotDI(SI_LowerThird_Directions,color='r') IPmag.iplotDI(SI_MiddleThird_Directions,color='y') IPmag.iplotDI(SI_UpperThird_Directions,color='b') title('Directional data from Simpson Island lava flows:') #plot the group means and a95 ellipses on same equal area plot IPmag.iplotDImean(lower_third_mean['dec'],lower_third_mean['inc'],lower_third_mean["alpha95"], color='r',marker='s',label='lower third') IPmag.iplotDImean(middle_third_mean['dec'],middle_third_mean['inc'],middle_third_mean["alpha95"], color='y',marker='s',label='middle third') IPmag.iplotDImean(upper_third_mean['dec'],upper_third_mean['inc'],upper_third_mean["alpha95"], color='b',marker='s',label='upper third') IPmag.iWatsonV(SI_LowerThird_Directions,SI_MiddleThird_Directions) IPmag.iBootstrap(SI_LowerThird_Directions,SI_MiddleThird_Directions) IPmag.iWatsonV(SI_MiddleThird_Directions,SI_UpperThird_Directions) IPmag.iBootstrap(SI_MiddleThird_Directions,SI_UpperThird_Directions) IPmag.iWatsonV(SI_LowerThird_Directions,SI_UpperThird_Directions) IPmag.iBootstrap(SI_LowerThird_Directions,SI_UpperThird_Directions) IPmag.iWatsonV(SI_LowerThird_Poles,SI_UpperThird_Poles) IPmag.iBootstrap(SI_LowerThird_Poles,SI_UpperThird_Poles) SI_LowerThird_plat=numpy.abs(IPmag.lat_from_i(lower_third_mean['inc'])) SI_LowerThird_plat_max=numpy.abs(IPmag.lat_from_i(lower_third_mean['inc']- lower_third_mean['alpha95'])) SI_LowerThird_plat_min=numpy.abs(IPmag.lat_from_i(lower_third_mean['inc']+ lower_third_mean['alpha95'])) SI_MiddleThird_plat=numpy.abs(IPmag.lat_from_i(middle_third_mean['inc'])) SI_MiddleThird_plat_max=numpy.abs(IPmag.lat_from_i(middle_third_mean['inc']- middle_third_mean['alpha95'])) SI_MiddleThird_plat_min=numpy.abs(IPmag.lat_from_i(middle_third_mean['inc']+ middle_third_mean['alpha95'])) SI_UpperThird_plat=numpy.abs(IPmag.lat_from_i(upper_third_mean['inc'])) SI_UpperThird_plat_max=numpy.abs(IPmag.lat_from_i(upper_third_mean['inc']- upper_third_mean['alpha95'])) SI_UpperThird_plat_min=numpy.abs(IPmag.lat_from_i(upper_third_mean['inc']+ upper_third_mean['alpha95'])) SI_LowerThird_meanstrat=521 SI_MiddleThird_meanstrat=521+1041 SI_UpperThird_meanstrat=521+1041+1041 figure() errorbar(SI_LowerThird_plat, SI_LowerThird_meanstrat, yerr=[[521],[521]], xerr=[[SI_LowerThird_plat-SI_LowerThird_plat_min], [SI_LowerThird_plat_max-SI_LowerThird_plat]], fmt='ro',label="lower third (N=30)") errorbar(SI_MiddleThird_plat, SI_MiddleThird_meanstrat, yerr=[[521],[521]], xerr=[[SI_MiddleThird_plat-SI_MiddleThird_plat_min], [SI_MiddleThird_plat_max-SI_MiddleThird_plat]], fmt='yo',label="middle third (N=20)") errorbar(SI_UpperThird_plat, SI_UpperThird_meanstrat, yerr=[[521],[521]], xerr=[[SI_UpperThird_plat-SI_UpperThird_plat_min], [SI_UpperThird_plat_max-SI_UpperThird_plat]], fmt='bo',label="upper third (N=34)") xlabel('paleolatitude (degrees)') ylabel('stratigraphic thickness (meters)') legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) fig = matplotlib.pyplot.gcf() fig.set_size_inches(5,10) pylab.xlim([20,70]) pylab.ylim([0,3200]) plt.show() SI_Lower500m_Directions=SI_Directions[0:17] SI_Lower500m_Poles=SI_Poles[0:17] SI_Upper500m_Directions=SI_Directions[67:84] SI_Upper500m_Poles=SI_Poles[67:84] Lower500m_mean=pmag.fisher_mean(SI_Lower500m_Directions) Upper500m_mean=pmag.fisher_mean(SI_Upper500m_Directions) print 'The Fisher mean parameters for SI_Lower500m_Directions are: ' print 'Dec = ' + str(Lower500m_mean['dec']) + ' Inc = ' + str(Lower500m_mean['inc']) print 'alpha95 = ' + str(Lower500m_mean['alpha95']) + ' k= ' + str(Lower500m_mean['k']) print '' print 'The Fisher mean parameters for SI_Upper500m_Directions are: ' print 'Dec = ' + str(Upper500m_mean['dec']) + ' Inc = ' + str(Upper500m_mean['inc']) print 'alpha95 = ' + str(Upper500m_mean['alpha95']) + ' k= ' + str(Upper500m_mean['k']) IPmag.iWatsonV(SI_Lower500m_Directions,SI_Upper500m_Directions) IPmag.iBootstrap(SI_Lower500m_Directions,SI_Upper500m_Directions) SI_Directions=[] SI_Poles=[] for n in range(0,84): Dec,Inc=Dec_TC[n],Inc_TC[n] SI_Directions.append([Dec,Inc,1.]) Plong,Plat=VGP_long[n],VGP_lat[n] SI_Poles.append([Plong,Plat,1.]) SI_LowerHalf_Directions=SI_Directions[0:41] SI_LowerHalf_Poles=SI_Poles[0:41] SI_UpperHalf_Directions=SI_Directions[41:84] SI_UpperHalf_Poles=SI_Poles[41:84] LowerHalf_mean=pmag.fisher_mean(SI_LowerHalf_Directions) UpperHalf_mean=pmag.fisher_mean(SI_UpperHalf_Directions) print 'The fisher mean parameters for SI_LowerHalf_Directions are: ' print 'Dec = ' + str(LowerHalf_mean['dec']) + ' Inc = ' + str(LowerHalf_mean['inc']) print 'alpha95 = ' + str(LowerHalf_mean['alpha95']) + ' k= ' + str(LowerHalf_mean['k']) print '' print 'The fisher mean parameters for SI_UpperHalf_Directions are: ' print 'Dec = ' + str(UpperHalf_mean['dec']) + ' Inc = ' + str(UpperHalf_mean['inc']) print 'alpha95 = ' + str(UpperHalf_mean['alpha95']) + ' k= ' + str(UpperHalf_mean['k']) IPmag.iWatsonV(SI_LowerHalf_Directions,SI_UpperHalf_Directions) IPmag.iBootstrap(SI_LowerHalf_Directions,SI_UpperHalf_Directions) LowerThird_MeanPole=pmag.fisher_mean(SI_LowerThird_Poles) MiddleThird_MeanPole=pmag.fisher_mean(SI_MiddleThird_Poles) UpperThird_MeanPole=pmag.fisher_mean(SI_UpperThird_Poles) print 'The Fisher mean parameters for the LowerThird mean pole are: ' print 'Plong = ' + str(LowerThird_MeanPole['dec']) + ' Plat = ' + str(LowerThird_MeanPole['inc']) print 'A95 = ' + str(LowerThird_MeanPole['alpha95']) + ' k= ' + str(LowerThird_MeanPole['k']) print '' print 'The Fisher mean parameters for the MiddleThird mean pole are: ' print 'Plong = ' + str(MiddleThird_MeanPole['dec']) + ' Plat = ' + str(MiddleThird_MeanPole['inc']) print 'A95 = ' + str(MiddleThird_MeanPole['alpha95']) + ' k= ' + str(MiddleThird_MeanPole['k']) print '' print 'The Fisher mean parameters for the UpperThird mean pole are: ' print 'Plong = ' + str(UpperThird_MeanPole['dec']) + ' Plat = ' + str(UpperThird_MeanPole['inc']) print 'A95 = ' + str(UpperThird_MeanPole['alpha95']) + ' k= ' + str(UpperThird_MeanPole['k']) #setup the figure and map figure(figsize=(8, 8)) m = Basemap(projection='ortho',lat_0=35,lon_0=200,resolution='l') # draw coastlines, country boundaries, fill continents. m.drawcoastlines(linewidth=0.25) #map.drawcountries(linewidth=0.25) m.fillcontinents(color='coral',lake_color='white') m.drawmapboundary(fill_color='white') m.drawmeridians(np.arange(0,360,30)) m.drawparallels(np.arange(-90,90,30)) #plot the mean poles IPmag.poleplot(m,LowerThird_MeanPole['dec'],LowerThird_MeanPole['inc'], LowerThird_MeanPole['alpha95'],color='r',label='Osler Group Lower Reversed Pole') IPmag.poleplot(m,MiddleThird_MeanPole['dec'],MiddleThird_MeanPole['inc'], MiddleThird_MeanPole['alpha95'],color='y',label='Osler Group Middle Reversed Pole') IPmag.poleplot(m,UpperThird_MeanPole['dec'],UpperThird_MeanPole['inc'], UpperThird_MeanPole['alpha95'],color='b',label='Osler Group Upper Reversed Pole') #show the plot legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) plt.show() Bibtex citation info for Halls (1974) @article{Halls1974a, Author = {Halls, H.C.}, Journal = {Canadian Journal of Earth Science}, Pages = {1200-1207}, Title = {A paleomagnetic reversal in the {O}sler {V}olcanic {G}roup, northern {L}ake {S}uperior}, Volume = {11}, Year = {1974}} Halls1974_Osler_Data=pandas.read_csv('../2014_Osler_Data/Halls1974_data.csv',sep=',') Halls1974_Osler_Data IPmag.VGP_calc(Halls1974_Osler_Data) Halls1974_Osler_Data.head() Halls1974_Osler_Data_N = Halls1974_Osler_Data[:5] Halls1974_Osler_Data_R = Halls1974_Osler_Data[5:] Halls1974_Osler_Data_R=Halls1974_Osler_Data_R.reset_index() Halls1974_Osler_R_VGPs=[] Halls1974_Osler_R_Plong=[] Halls1974_Osler_R_Plat=[] for n in range(0,len(Halls1974_Osler_Data_R)): Plong,Plat=Halls1974_Osler_Data_R['pole_long_rev'][n],Halls1974_Osler_Data_R['pole_lat_rev'][n] Halls1974_Osler_R_Plong.append(Plong) Halls1974_Osler_R_Plat.append(Plat) Halls1974_Osler_R_VGPs.append([Plong,Plat,1.]) IPmag.iWatsonV(Halls1974_Osler_R_VGPs,SI_LowerThird_Poles) IPmag.iBootstrap(Halls1974_Osler_R_VGPs,SI_LowerThird_Poles) IPmag.iWatsonV(Halls1974_Osler_R_VGPs,SI_UpperThird_Poles) IPmag.iBootstrap(Halls1974_Osler_R_VGPs,SI_UpperThird_Poles) upperR_Osler_VGPs=Halls1974_Osler_R_VGPs+SI_UpperThird_Poles upperR_Osler_meanpole=pmag.fisher_mean(upperR_Osler_VGPs) print 'The fisher mean parameters for the pole calculated from upper reversed' print 'Osler Group VGPs (SI_UpperThird_Poles+Halls, 1974 reversed data) are: ' print 'Plong = ' + str(upperR_Osler_meanpole['dec']) + ' Plat = ' + str(upperR_Osler_meanpole['inc']) print 'A95 = ' + str(upperR_Osler_meanpole['alpha95']) + ' k= ' + str(upperR_Osler_meanpole['k']) + ' N= ' + str(upperR_Osler_meanpole['n']) 001 0.0 0.0 0.0 0.0 000 ! \\ 001 3900.0 0.0 0.0 0.0 000 ! \\ 199 0.0 0.0 0.0 0.0 001 ! \\ 199 1105.0 0.0 115.4 48.4 001 !Put Osler-upperthird at N for Laurentia \\ 199 1107.5 0.0 121.3 47.3 001 !Put Osler-middlethird at N for Laurentia \\ 199 1110.0 0.0 128.7 49.0 001 !Put Osler-lowerthird at N for Laurentia \\ from IPython.display import SVG SVG("Paleogeo.svg") #import a couple special functions from scipy from scipy import special from scipy import stats Tauxe_NSVG_Data=pandas.read_csv('../2014_Osler_Data/Tauxe2009a_data.csv',sep=',') #show first 5 rows Tauxe_NSVG_Data.head() NSVG_nswu=Tauxe_NSVG_Data.ix[Tauxe_NSVG_Data['sequence'] == 'nswu'] NSVG_nswu.reset_index(inplace=True) NSVG_nswu_VGPs=[] NSVG_nswu_Plong=[] NSVG_nswu_Plat=[] for n in range(0,len(NSVG_nswu)): Plong,Plat=NSVG_nswu['pole_long'][n],NSVG_nswu['pole_lat'][n] NSVG_nswu_Plong.append(Plong) NSVG_nswu_Plat.append(Plat) NSVG_nswu_VGPs.append([Plong,Plat,1.]) NSVG_nswu_mean=pmag.fisher_mean(NSVG_nswu_VGPs) print 'Here are the details for the calculated pole from the SW limb of the NSVG' print 'Pole longitude is: ' + str(round(NSVG_nswu_mean['dec'],1)) print 'Pole latitude is: ' + str(round(NSVG_nswu_mean['inc'],1)) print 'Pole kappa is: ' + str(round(NSVG_nswu_mean['k'],1)) print 'Pole A95 is: ' + str(round(NSVG_nswu_mean['alpha95'],1)) print 'Pole N is: ' + str(int(NSVG_nswu_mean['n'])) samplesize=100000 #set at 100,000 this code will take a long time to run #parameters for pole 1 (OVG pole calculated above) pole1_plong=201.6 pole1_plat=42.5 pole1_A95=3.7 pole1_kappa=25.7 pole1_N=59 pole1_age=1105 #1 sigma age uncertainty pole1_age_error=1 #parameters for pole 2 (NSVG pole calculated above) pole2_plong=182.1 pole2_plat=35.8 pole2_A95=3.1 pole2_kappa=45.7 pole2_N=47 #taking the average of the 40th Ave and Palisade pole2_age=1097.5 #2sigma uncertainty min age on Palisade is 1094.9 #2sigma uncertainty min age on 40th Ave is 1100.3 #a 2sigma of 2.7 on the pole2_age approx spans this range #giving a 1 sigma age uncertainty of pole2_age_error=1.35 #the longitude, latitude of Duluth, MN as a reference location Duluth=[267.9, 46.8] pole1=(pole1_plong, pole1_plat) pole1_paleolat=90-pmag.angle(pole1,Duluth) pole2=(pole2_plong, pole2_plat) pole2_paleolat=90-pmag.angle(pole2,Duluth) print "The paleolatitude for Duluth resulting from pole 1 is:" + str(pole1_paleolat) print "The paleolatitude for Duluth resulting from pole 2 is:" + str(pole2_paleolat) rate=((pole1_paleolat-pole2_paleolat)*111*100000)/((pole1_age-pole2_age)*1000000) print "The rate of paleolatitudinal change implied by the poles pairs in cm/yr is:" + str(rate) pole1_MCages=np.random.normal(pole1_age,pole1_age_error,samplesize) pole2_MCages=np.random.normal(pole2_age,pole2_age_error,samplesize) plt.hist(pole1_MCages,100,histtype='stepfilled',color='darkred',label='Pole 1 ages') plt.hist(pole2_MCages,100,histtype='stepfilled',color='darkblue',label='Pole 2 ages') plt.xlabel('Age (Ma)') plt.ylabel('n') plt.legend(loc=3) plt.show() pole1_MCpoles=[] pole1_MCpole_lat=[] pole1_MCpole_long=[] pole1_MCpaleolat=[] for n in range(samplesize): vgp_samples=[] for vgp in range(pole1_N): #pmag.dev returns a direction from a fisher distribution with specified kappa direction_atN=pmag.fshdev(pole1_kappa) #this direction is centered at latitude of 90º and needs to be rotated #to be centered on the mean pole position tilt_direction=pole1_plong tilt_amount=90-pole1_plat direction=pmag.dotilt(direction_atN[0],direction_atN[1],tilt_direction,tilt_amount) vgp_samples.append([direction[0],direction[1],1.]) mean=pmag.fisher_mean(vgp_samples) mean_pole_position=(mean['dec'],mean['inc']) pole1_MCpoles.append([mean['dec'],mean['inc'],1.]) pole1_MCpole_lat.append(mean['inc']) pole1_MCpole_long.append(mean['dec']) paleolat=90-pmag.angle(mean_pole_position,Duluth) pole1_MCpaleolat.append(paleolat[0]) pole2_MCpoles=[] pole2_MCpole_lat=[] pole2_MCpole_long=[] pole2_MCpaleolat=[] for n in range(samplesize): vgp_samples=[] for vgp in range(pole2_N): #pmag.dev returns a direction from a fisher distribution with specified kappa direction_atN=pmag.fshdev(pole2_kappa) #this direction is centered at latitude of 90º and needs to be rotated #to be centered on the mean pole position tilt_direction=pole2_plong tilt_amount=90-pole2_plat direction=pmag.dotilt(direction_atN[0],direction_atN[1],tilt_direction,tilt_amount) vgp_samples.append([direction[0],direction[1],1.]) mean=pmag.fisher_mean(vgp_samples) mean_pole_position=(mean['dec'],mean['inc']) pole2_MCpoles.append([mean['dec'],mean['inc'],1.]) pole2_MCpole_lat.append(mean['inc']) pole2_MCpole_long.append(mean['dec']) paleolat=90-pmag.angle(mean_pole_position,Duluth) pole2_MCpaleolat.append(paleolat[0]) plt.figure(figsize=(8, 8)) m = Basemap(projection='ortho',lat_0=35,lon_0=200,resolution='c',area_thresh=50000) m.drawcoastlines(linewidth=0.25) m.fillcontinents(color='bisque',lake_color='white',zorder=1) m.drawmapboundary(fill_color='white') m.drawmeridians(np.arange(0,360,30)) m.drawparallels(np.arange(-90,90,30)) IPmag.vgpplot(m,pole1_MCpole_long,pole1_MCpole_lat,color='b') IPmag.vgpplot(m,pole2_MCpole_long,pole2_MCpole_lat,color='g') #calculating the change in paleolatitude between the Monte Carlo pairs pole1_pole2_Delta_degrees=[] pole1_pole2_Delta_kilometers=[] pole1_pole2_Delta_myr=[] pole1_pole2_degrees_per_myr=[] pole1_pole2_cm_per_yr=[] for n in range(samplesize): Delta_degrees=pole1_MCpaleolat[n]-pole2_MCpaleolat[n] Delta_Myr=pole1_MCages[n]-pole2_MCages[n] pole1_pole2_Delta_degrees.append(Delta_degrees) degrees_per_myr=Delta_degrees/Delta_Myr cm_per_yr=((Delta_degrees*111)*100000)/(Delta_Myr*1000000) pole1_pole2_degrees_per_myr.append(degrees_per_myr) pole1_pole2_cm_per_yr.append(cm_per_yr) twopointfive_percentile=stats.scoreatpercentile(pole1_pole2_cm_per_yr,2.5) fifty_percentile=stats.scoreatpercentile(pole1_pole2_cm_per_yr,50) ninetysevenpointfive_percentile=stats.scoreatpercentile(pole1_pole2_cm_per_yr,97.5) print "2.5th percentile is: " + str(twopointfive_percentile) print "50th percentile is: " + str(fifty_percentile) print "97.5th percentile is: " + str(ninetysevenpointfive_percentile) plotnumber=5000 figure(num=None, figsize=(14, 4)) plt.subplot(1, 2, 1) for n in range(plotnumber): plt.plot([pole1_MCpaleolat[n],pole2_MCpaleolat[n]], [pole1_MCages[n],pole2_MCages[n]],'k-',linewidth=0.05,alpha=0.3) plt.scatter(pole1_MCpaleolat[:plotnumber],pole1_MCages[:plotnumber],color='b',s=3) plt.scatter(pole1_paleolat,pole1_age,color='lightblue',s=100, edgecolor='w', zorder=10000) plt.scatter(pole2_MCpaleolat[:plotnumber],pole2_MCages[:plotnumber],color='g',s=3) plt.scatter(pole2_paleolat,pole2_age,color='lightgreen',s=100, edgecolor='w', zorder=10000) plt.plot([pole1_paleolat,pole2_paleolat],[pole1_age,pole2_age],'w-',linewidth=2) plt.gca().invert_yaxis() plt.xlabel('paleolatitude at Duluth, MN (degrees)',size=14) plt.ylabel('time (Ma)',size=14) plt.text(pole1_paleolat,1099,'OVG', color='b',ha ='center',size=14) plt.text(pole2_paleolat,1104.5,'NSVG', color='g',ha ='center',size=14) plt.text(20.5,1092.5,'A', color='k',ha ='left',size=28) plt.subplot(1, 2, 2) plt.hist(pole1_pole2_cm_per_yr,bins=600) plt.vlines(twopointfive_percentile,0,4500,'r', linestyles='dashed') plt.text(twopointfive_percentile-0.5,4000,'15 cm/yr', color='r',ha ='right',size=12) plt.vlines(fifty_percentile,0,4500,'r', linestyles='solid') plt.text(fifty_percentile+0.5,4000,'24 cm/yr', color='r',ha ='left',size=12) plt.vlines(ninetysevenpointfive_percentile,0,4500,'r', linestyles='dashed') plt.text(ninetysevenpointfive_percentile+0.5,4000,'44 cm/yr', color='r',ha ='left',size=12) plt.text(1.5,3930,'B', color='k',ha ='left',size=28) plt.ylabel('n',size=14) plt.xlabel('latitudinal drift rate (cm/yr)',size=14) plt.xlim([0,60]) plt.ylim([0,4500]) plt.savefig('../2014_Osler_Manuscript_Files/2014_Osler_Figures/MonteCarlo.pdf', bbox_inches='tight') plt.show() lower_Osler_pole=(218.6,40.9) NSVG_nswu_pole=(182.1,35.8) lower_Osler_pole_paleolat=90-pmag.angle(lower_Osler_pole,Duluth) NSVG_nswu_paleolat=90-pmag.angle(NSVG_nswu_pole,Duluth) difference=lower_Osler_pole_paleolat-NSVG_nswu_paleolat print 'The difference in paleolatitude resulting from these poles is: ' + str(difference) print 'from ' + str(lower_Osler_pole_paleolat) + ' to ' + str(NSVG_nswu_paleolat) + ' using Duluth, MN as a reference location' #flow thicknesses where bottom and top are exposed without cover in between #or with minimal cover and assumption of continuity SI1_flowthickness = array([1.2, 1.1+0.7,0.2+2.9+0.3,5.1+5.4+0.6+3.4,0.6+1+.3,0.9,0.5,2.2,0.8,0.5,0.3+1.3+0.6+1.6,0.4+2.3,0.6,0.8+1.0,1.6+0.3,1.0+0.6+2.4,0.2+4.4+0.8+0.6,3.3+1.6,0.7,4+4.7+5.7+1.7+1.8,0.9,0.2+0.5,0.2+0.8+1.3,0.2+0.9,1.7,0.6,0.2+1.5,1.2+0.4+1.3,1.3+1.5,1.3,2.4,0.3+1.7,0.6+0.3,0.3+0.5+.3+1.3,0.5+1.1,0.1+0.3+1.3,1.5,0.1+1.9+1.7+1.2]) SI2_flowthickness = array([0.6+2.8,5.7+2+3.6,1.3+0.7+1.1+1.1,6.2+1,1.8+3.3+11.4]) SI3_flowthickness = array([0.6+0.8+1+.8,1+4.9+1.1,0.3+1.6+4.5+1.2,0.2+1.5+7.2+0.9,0.5+1.7+9.3,1.4+2.9+0.6,3.1+3+0.9,3.9+1.4+.3]) SI4_flowthickness = array([0.2+1.2+0.9,0.3+1.5+0.4+3.2,14+6.5,0.7+1+3.6,8.2+2.4+4.6,3.1+1.4+1.4,1.6,2.8+2.1,0.7+5.6+2.4+0.6,5.6+0.7,2.2+0.2+0.3,0.4+2.4+1.9+1,4.2+0.7]) SI5_flowthickness = array([6.6+9.8+7.4+1.4+2.3,10.2+0.8,5.6+2.8+1.4+2.2,0.2+2.8+2.4+0.7,0.6,0.2+0.5,0.2+2.2,1.7+1.7,1.6+0.6,0.2+3.8+1.4+1,7.6+20.4+6.3+4.2+1,2.8+1.3+1.2,10.2+9]) SI6_flowthickness = array([13.1+3.3,13.5+3.7,6.5+0.6,6.1+2.1+0.9,0.9+6.6+4.0,7.3+5.6+1.3+7.4,0.5+9.6+3.4+0.5,2.8+4.3,0.3+7.4+3.9]) SI7_flowthickness = array([0.4+0.9+1.2+1.4,0.3+1.2,0.5+2.5+5+1.9]) SI8_flowthickness = array([1.4+1.9+0.6,0.5+1+8.8+4+1.1,5+15.8+2.8+2.1,5+4.2,8.4+8.1+5,1.8+11.8+1+2.2+1.4]) SI9_flowthickness = array([5.0+4.2,14.0+1+2.1,2.8+2.1,3.6+12.3+3.6+0.7,4.2+2.3,1.3+1.4+2.2,0.8+1.8+9.8+4.2+0.8,0.2+2.0+0.6+2.3,0.3+0.4+3.3,7.5]) SIall_flowthickness = concatenate((SI1_flowthickness, SI2_flowthickness, SI3_flowthickness, SI4_flowthickness, SI5_flowthickness, SI6_flowthickness, SI7_flowthickness, SI8_flowthickness, SI9_flowthickness)) print "The number of flows for which thickness estimates could be obtained is:" print len(SIall_flowthickness) print "The mean flow thickness is:" print np.mean(SIall_flowthickness) print "The median flow thickness is:" print np.median(SIall_flowthickness) print "The first quartile thickness is:" print np.percentile(SIall_flowthickness,25) print "The third quartile thickness is:" print np.percentile(SIall_flowthickness,75) hist(SIall_flowthickness, 40, facecolor='green') title('Histogram of Osler flow thicknesses within measured sections') ylabel('number of flows') xlabel('flow thickness (meters)') plt.show() data=(SI1_flowthickness,SI2_flowthickness,SI3_flowthickness, SI4_flowthickness,SI5_flowthickness, SI6_flowthickness, SI7_flowthickness,SI8_flowthickness,SI9_flowthickness,SIall_flowthickness) figure(figsize=(12,4)) title('Boxplot of flow thicknesses grouped by stratigraphic section') boxplot(data) labels = ('SI1', 'SI2', 'SI3', 'SI4', 'SI5', 'SI6', 'SI7', 'SI8', 'SI9', 'all sections') xticks(range(1,11),labels) xlabel('stratigraphic section') ylabel('flow thickness (meters)') plt.show() %load IPmag.py import pmag, pmagplotlib import pylab import numpy as np import matplotlib.pyplot as plt def iflip(D): #function simplified from PmagPy pmag.flip function """ This function returns the antipode (flips) of the unit vectors in D (dec,inc,length). """ Dflip=[] for rec in D: d,i=(rec[0]-180.)%360.,-rec[1] Dflip.append([d,i,1.]) return Dflip def iBootstrap(Data1,Data2,NumSims=1000): """ Conduct a bootstrap test (Tauxe, 2010) for a common mean on two declination, inclination data sets This function modifies code from PmagPy for use calculating and plotting bootstrap statistics. Three plots are generated (one for x, one for y and one for z). If the 95 percent confidence bounds for each component overlap each other, the two directions are not significantly different. Parameters ---------- Data1 : a list of directional data [dec,inc] Data2 : a list of directional data [dec,inc] NumSims : number of bootstrap samples (default is 1000) """ counter=0 BDI1=pmag.di_boot(Data1) BDI2=pmag.di_boot(Data2) print "" print "===============" print "" print "Here are the results of the bootstrap test for a common mean" CDF={'X':1,'Y':2,'Z':3} pylab.figure(CDF['X'],figsize=(3,3),dpi=160) pylab.figure(CDF['Y'],figsize=(3,3),dpi=160) pylab.figure(CDF['Z'],figsize=(3,3),dpi=160) pmagplotlib.plotCOM(CDF,BDI1,BDI2,["",""]) def iWatsonV(Data1,Data2,NumSims=5000): """ Conduct a Watson V test for a common mean on two declination, inclination data sets This function calculates Watson's V statistic from input files through Monte Carlo simulation in order to test whether two populations of directional data could have been drawn from a common mean. The critical angle between the two sample mean directions and the corresponding McFadden and McElhinny (1990) classification is printed. Parameters ---------- Data1 : a list of directional data [dec,inc] Data2 : a list of directional data [dec,inc] NumSims : number of Monte Carlo simulations (default is 5000) """ pars_1=pmag.fisher_mean(Data1) pars_2=pmag.fisher_mean(Data2) cart_1=pmag.dir2cart([pars_1["dec"],pars_1["inc"],pars_1["r"]]) cart_2=pmag.dir2cart([pars_2['dec'],pars_2['inc'],pars_2["r"]]) Sw=pars_1['k']*pars_1['r']+pars_2['k']*pars_2['r'] # k1*r1+k2*r2 xhat_1=pars_1['k']*cart_1[0]+pars_2['k']*cart_2[0] # k1*x1+k2*x2 xhat_2=pars_1['k']*cart_1[1]+pars_2['k']*cart_2[1] # k1*y1+k2*y2 xhat_3=pars_1['k']*cart_1[2]+pars_2['k']*cart_2[2] # k1*z1+k2*z2 Rw=np.sqrt(xhat_1**2+xhat_2**2+xhat_3**2) V=2*(Sw-Rw) # keep weighted sum for later when determining the "critical angle" # let's save it as Sr (notation of McFadden and McElhinny, 1990) Sr=Sw # do monte carlo simulation of datasets with same kappas as data, # but a common mean counter=0 Vp=[] # set of Vs from simulations for k in range(NumSims): # get a set of N1 fisher distributed vectors with k1, # calculate fisher stats Dirp=[] for i in range(pars_1["n"]): Dirp.append(pmag.fshdev(pars_1["k"])) pars_p1=pmag.fisher_mean(Dirp) # get a set of N2 fisher distributed vectors with k2, # calculate fisher stats Dirp=[] for i in range(pars_2["n"]): Dirp.append(pmag.fshdev(pars_2["k"])) pars_p2=pmag.fisher_mean(Dirp) # get the V for these Vk=pmag.vfunc(pars_p1,pars_p2) Vp.append(Vk) # sort the Vs, get Vcrit (95th percentile one) Vp.sort() k=int(.95*NumSims) Vcrit=Vp[k] # equation 18 of McFadden and McElhinny, 1990 calculates the critical # value of R (Rwc) Rwc=Sr-(Vcrit/2) # following equation 19 of McFadden and McElhinny (1990) the critical # angle is calculated. If the observed angle (also calculated below) # between the data set means exceeds the critical angle the hypothesis # of a common mean direction may be rejected at the 95% confidence # level. The critical angle is simply a different way to present # Watson's V parameter so it makes sense to use the Watson V parameter # in comparison with the critical value of V for considering the test # results. What calculating the critical angle allows for is the # classification of McFadden and McElhinny (1990) to be made # for data sets that are consistent with sharing a common mean. k1=pars_1['k'] k2=pars_2['k'] R1=pars_1['r'] R2=pars_2['r'] critical_angle=np.degrees(np.arccos(((Rwc**2)-((k1*R1)**2) -((k2*R2)**2))/ (2*k1*R1*k2*R2))) D1=(pars_1['dec'],pars_1['inc']) D2=(pars_2['dec'],pars_2['inc']) angle=pmag.angle(D1,D2) print "Results of Watson V test: " print "" print "Watson's V: " '%.1f' %(V) print "Critical value of V: " '%.1f' %(Vcrit) if VVcrit: print '"Fail": Since V is greater than Vcrit, the two means can' print 'be distinguished at the 95% confidence level.' print "" print "M&M1990 classification:" print "" print "Angle between data set means: " '%.1f'%(angle) print "Critical angle for M&M1990: " '%.1f'%(critical_angle) if V>Vcrit: print "" elif V= 0: X_down.append(XY[0]) Y_down.append(XY[1]) else: X_up.append(XY[0]) Y_up.append(XY[1]) if len(X_up)>0: pylab.scatter(X_up,Y_up,facecolors='none', edgecolors=color) if len(X_down)>0: pylab.scatter(X_down,Y_down,facecolors=color, edgecolors=color) def iplotDImean(Dec,Inc,a95,color='k',marker='o',label=''): """ Plot a mean declination, inclination with alpha_95 ellipse on an equal area plot. Before this function is called a plot needs to be initialized with code that looks something like: >fignum = 1 >pylab.figure(num=fignum,figsize=(10,10),dpi=160) >pmagplotlib.plotNET(fignum) Parameters ---------- Dec : declination of mean being plotted Inc : inclination of mean being plotted a95 : a95 confidence ellipse of mean being plotted color : the default color is black. Other colors can be chosen (e.g. 'r') marker : the default is a circle. Other symbols can be chose (e.g. 's') label : the default is no label. Labels can be assigned """ DI_dimap=pmag.dimap(Dec,Inc) pylab.plot(DI_dimap[0],DI_dimap[1],color=color,marker=marker,label=label) Xcirc,Ycirc=[],[] Da95,Ia95=pmag.circ(Dec,Inc,a95) pylab.legend(loc=2) for k in range(len(Da95)): XY=pmag.dimap(Da95[k],Ia95[k]) Xcirc.append(XY[0]) Ycirc.append(XY[1]) pylab.plot(Xcirc,Ycirc,color) def shoot(lon, lat, azimuth, maxdist=None): """ This function enables A95 error ellipses to be drawn in basemap around paleomagnetic poles in conjunction with equi (from: http://www.geophysique.be/2011/02/20/matplotlib-basemap-tutorial-09-drawing-circles/) """ glat1 = lat * np.pi / 180. glon1 = lon * np.pi / 180. s = maxdist / 1.852 faz = azimuth * np.pi / 180. EPS= 0.00000000005 if ((np.abs(np.cos(glat1)) EPS): sy = np.sin(y) cy = np.cos(y) cz = np.cos(b + y) e = 2. * cz * cz - 1. c = y x = e * cy y = e + e - 1. y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) * d / 4. - cz) * sy * d + tu b = cu * cy * cf - su * sy c = r * np.sqrt(sa * sa + b * b) d = su * cy + cu * sy * cf glat2 = (np.arctan2(d, c) + np.pi) % (2*np.pi) - np.pi c = cu * cy - su * sy * cf x = np.arctan2(sy * sf, c) c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16. d = ((e * cy * c + cz) * sy * c + y) * sa glon2 = ((glon1 + x - (1. - c) * d * f + np.pi) % (2*np.pi)) - np.pi baz = (np.arctan2(sa, b) + np.pi) % (2 * np.pi) glon2 *= 180./np.pi glat2 *= 180./np.pi baz *= 180./np.pi return (glon2, glat2, baz) def equi(m, centerlon, centerlat, radius, color): """ This function enables A95 error ellipses to be drawn in basemap around paleomagnetic poles in conjunction with shoot (from: http://www.geophysique.be/2011/02/20/matplotlib-basemap-tutorial-09-drawing-circles/). """ glon1 = centerlon glat1 = centerlat X = [] Y = [] for azimuth in range(0, 360): glon2, glat2, baz = shoot(glon1, glat1, azimuth, radius) X.append(glon2) Y.append(glat2) X.append(X[0]) Y.append(Y[0]) X,Y = m(X,Y) plt.plot(X,Y,color) def poleplot(mapname,plong,plat,A95,label='',color='k',marker='o'): """ This function plots a paleomagnetic pole and A95 error ellipse on whatever current map projection has been set using the basemap plotting library. Parameters ----------- mapname : the name of the current map that has been developed using basemap plong : the longitude of the paleomagnetic pole being plotted (in degrees E) plat : the latitude of the paleomagnetic pole being plotted (in degrees) A95 : the A_95 confidence ellipse of the paleomagnetic pole (in degrees) label : a string that is the label for the paleomagnetic pole being plotted color : the color desired for the symbol and its A95 ellipse (default is 'k' aka black) marker : the marker shape desired for the pole mean symbol (default is 'o' aka a circle) """ centerlon, centerlat = mapname(plong,plat) A95_km=A95*111.32 mapname.scatter(centerlon,centerlat,20,marker=marker,color=color,label=label,zorder=101) equi(mapname, plong, plat, A95_km,color) def vgpplot(mapname,plong,plat,label='',color='k',marker='o'): """ This function plots a paleomagnetic pole on whatever current map projection has been set using the basemap plotting library. Parameters ----------- mapname : the name of the current map that has been developed using basemap plong : the longitude of the paleomagnetic pole being plotted (in degrees E) plat : the latitude of the paleomagnetic pole being plotted (in degrees) color : the color desired for the symbol and its A95 ellipse (default is 'k' aka black) marker : the marker shape desired for the pole mean symbol (default is 'o' aka a circle) """ centerlon, centerlat = mapname(plong,plat) mapname.scatter(centerlon,centerlat,20,marker=marker,color=color,label=label,zorder=100) def VGP_calc(dataframe): """ This function calculates paleomagnetic poles from directional data within a pandas.DataFrame Parameters ----------- dataframe : the name of the pandas.DataFrame containing the data dataframe['site_lat'] : the latitude of the site dataframe['site_long'] : the longitude of the site dataframe['inc_tc'] : the tilt-corrected inclination dataframe['dec_tc'] : the tilt-corrected declination """ #calculate the paleolatitude/colatitude dataframe['paleolatitude']=np.degrees(np.arctan(0.5*np.tan(np.radians(dataframe['inc_tc'])))) dataframe['colatitude']=90-dataframe['paleolatitude'] #calculate the latitude of the pole dataframe['pole_lat']=np.degrees(np.arcsin(np.sin(np.radians(dataframe['site_lat']))* np.cos(np.radians(dataframe['colatitude']))+ np.cos(np.radians(dataframe['site_lat']))* np.sin(np.radians(dataframe['colatitude']))* np.cos(np.radians(dataframe['dec_tc'])))) #calculate the longitudinal difference between the pole and the site (beta) dataframe['beta']=np.degrees(np.arcsin((np.sin(np.radians(dataframe['colatitude']))* np.sin(np.radians(dataframe['dec_tc'])))/ (np.cos(np.radians(dataframe['pole_lat']))))) #generate a boolean array (mask) to use to distinguish between the two possibilities for pole longitude #and then calculate pole longitude using the site location and calculated beta mask = np.cos(np.radians(dataframe['colatitude']))>np.sin(np.radians(dataframe['site_lat']))*np.sin(np.radians(dataframe['pole_lat'])) dataframe['pole_long']=np.where(mask,(dataframe['site_long']+dataframe['beta'])%360.,(dataframe['site_long']+180-dataframe['beta'])%360.) #calculate the antipode of the poles dataframe['pole_lat_rev']=-dataframe['pole_lat'] dataframe['pole_long_rev']=(dataframe['pole_long']-180.)%360. #the 'colatitude' and 'beta' columns were created for the purposes of the pole calculations #but aren't of further use and are deleted del dataframe['colatitude'] del dataframe['beta']