Made resolution 1 base map, and changed well plot

Author: Murray2015

Resolution-1_well_plot.jpg

@@ -6,7 +6,9 @@
 """
 
 import pandas as pd 
+import numpy as np
 import matplotlib.pyplot as plt
+import statsmodels.formula.api as sm
 
 # Read in data 
 data = pd.read_csv("opal_ct.csv")
@@ -18,6 +20,7 @@
 ax1.set_ylabel('Depth below seafloor')
 # Make the y-axis label, ticks and tick labels match the line color.
 ax1.set_xlabel('Count', color='b')
+plt.xlim(0,20)
 ax1.tick_params('x', colors='b')
 ax2 = ax1.twiny()
 ax2.hist(data2['Depth_below_seafloor'], normed=1, histtype='step',cumulative=True, color='r', orientation='horizontal')
@@ -31,15 +34,110 @@
 data2['burial_rate'] = data2['Depth_below_seafloor'] / data2['Age_host_seds']
 
 # Plot data on a semilog plot 
-plt.scatter(data2['burial_rate'][data2['burial_rate'] < 3], data2['Depth_below_seafloor'][data2['burial_rate'] < 3], color='b') 
-plt.scatter(data2['burial_rate'][data2['burial_rate'] > 3], data2['Depth_below_seafloor'][data2['burial_rate'] > 3], color='c') 
-plt.scatter(data2['burial_rate'][data2['burial_rate']> 30], data2['Depth_below_seafloor'][data2['burial_rate'] > 30], color='r') 
+plt.scatter(data2['burial_rate'][data2['burial_rate'] < 4], data2['Depth_below_seafloor'][data2['burial_rate'] < 4], color='b') 
+plt.scatter(data2['burial_rate'][data2['burial_rate'] > 4], data2['Depth_below_seafloor'][data2['burial_rate'] > 4], color='c') 
+plt.scatter(data2['burial_rate'][data2['burial_rate']> 40], data2['Depth_below_seafloor'][data2['burial_rate'] > 40], color='r') 
 plt.xlabel('Burial rate (m/Ma)')
 plt.ylabel('Depth (m)')
+plt.ylim(0,1000)
 ax = plt.gca()
 ax.invert_yaxis()
 ax.set_xscale('log')
 plt.show()
 
+# Make cumulative hist for subsets of data
+fig, ax1 = plt.subplots(figsize=(4,5))
+ax1.hist(data2['Depth_below_seafloor'][data2['burial_rate'] < 4], color='b', orientation='horizontal', bins=range(0,1000,100))
+ax1.set_ylabel('Depth below seafloor')
+# Make the y-axis label, ticks and tick labels match the line color.
+ax1.set_xlabel('Count', color='b')
+ax1.set_xlim(0,15)
+ax1.set_ylim([0,1000])
+ax1.tick_params('x', colors='b')
+ax2 = ax1.twiny()
+ax2.hist(data2['Depth_below_seafloor'][data2['burial_rate'] < 4], normed=1, histtype='step',cumulative=True, color='r', orientation='horizontal')
+ax2.set_xlabel('Cumulative frequency', color='r')
+ax2.tick_params('x', colors='r')
+ax2.set_ylim([0,1000])
+plt.gca().invert_yaxis()
+fig.tight_layout()
+plt.show()
+
+fig, ax1 = plt.subplots(figsize=(4,5))
+ax1.hist(data2.loc[(data2["burial_rate"] > 4) & (data2["burial_rate"] < 40), "Depth_below_seafloor"], color='c', orientation='horizontal', bins=range(0,1000,100))
+ax1.set_ylabel('Depth below seafloor')
+# Make the y-axis label, ticks and tick labels match the line color.
+ax1.set_xlabel('Count', color='c')
+ax1.set_xlim(0,15)
+ax1.set_ylim([0,1000])
+ax1.tick_params('x', colors='c')
+ax2 = ax1.twiny()
+ax2.hist(data2.loc[(data2["burial_rate"] > 4) & (data2["burial_rate"] < 40), "Depth_below_seafloor"], normed=1, histtype='step',cumulative=True, color='r', orientation='horizontal')
+ax2.set_xlabel('Cumulative frequency', color='r')
+ax2.tick_params('x', colors='r')
+ax2.set_ylim([0,1000])
+plt.gca().invert_yaxis()
+fig.tight_layout()
+plt.show()
+
+fig, ax1 = plt.subplots(figsize=(4,5))
+ax1.hist(data2['Depth_below_seafloor'][data2['burial_rate'] >40], color='r', orientation='horizontal', bins=range(0,1000,100))
+ax1.set_ylabel('Depth below seafloor')
+# Make the y-axis label, ticks and tick labels match the line color.
+ax1.set_xlabel('Count', color='r')
+ax1.set_xlim(0,15)
+ax1.set_ylim([0,1000])
+ax1.tick_params('x', colors='r')
+ax2 = ax1.twiny()
+ax2.hist(data2['Depth_below_seafloor'][data2['burial_rate'] >40], normed=1, histtype='step',cumulative=True, color='r', orientation='horizontal')
+ax2.set_xlabel('Cumulative frequency', color='r')
+ax2.tick_params('x', colors='r')
+ax2.set_ylim([0,1000])
+plt.gca().invert_yaxis()
+fig.tight_layout()
+plt.show()
+
+# 0th order 
+## plot cumulative freq against time
+# time = depth / burial rate 
+# fit regression line, slope is reaction rate k. 
+
+# For inversion, vector of A and vector of E (to make grid search). 
+A = np.linspace(start=10e-2, stop=10e4, num=10)
+E = np.linspace(start=70, stop=110, num=10)
+misfit_grid = np.zeros((10,10))
+# find temp at each point with a geothermal gradient of 30 degrees c and a surface temp of 0. 
+data2['temp'] = data2['Depth_below_seafloor']/1000 * 30
+# Find misfit as abs(predicted - real). Fill it into grid search. 
+# Sort by depth 
+data3 = data2.sort_values(by='Depth_below_seafloor')
+# Make column of cumulative percentage
+data3['cum_perc'] = np.cumsum(data3['Depth_below_seafloor'].values)/np.max(np.cumsum(data3['Depth_below_seafloor'].values)) * 100
+result = sm.ols(formula="cum_perc ~ Depth_below_seafloor", data=data3).fit()
+k = result.params['Depth_below_seafloor']
+
+
+
+
+
+# 1st order 
+## plot ln(cumulative freq) against time
+# time = depth / burial rate 
+# fit regression line, slope is -ve reaction rate k. 
+
+# For inversion, vector of A and vector of E (to make grid search). 
+# find temp at each point with a geothermal gradient of 30 degrees c and a surface temp of 0. 
+# Find misfit as abs(predicted - real). Fill it into grid search. 
+
+
+# Nucleation and growth reaction 
+## plot ln(cumulative freq) against time
+# time = depth / burial rate 
+# fit regression line, slope is -ve reaction rate k. 
+
+# For inversion, vector of A and vector of E (to make grid search). 
+# find temp at each point with a geothermal gradient of 30 degrees c and a surface temp of 0. 
+# Find misfit as abs(predicted - real). Fill it into grid search. 
+
 
 

opal_ct/opal_ct_inversion.py

@@ -0,0 +1 @@
+,murray,murray-Latitude-E5430-non-vPro,11.12.2017 21:55,file:///home/murray/.config/libreoffice/4;

resolution_1/Resolution-1 well/.~lock.Sonic log data.xlsx#

@@ -49,7 +49,7 @@
 ax5.axvline(x=1900, color='red')
 ax5.axvline(x=1220, color='blue')
 fig.subplots_adjust(wspace=0, hspace=0)
-plt.savefig("/home/murray/Documents/thesis_python_code/Resolution-1_well_plot.pdf")
+#plt.savefig("/home/murray/Documents/thesis_python_code/Resolution-1_well_plot.pdf")
 plt.show()
 
 ## Fill log data with dummy values for water and sed. Create a dataframe of the filled values, 
@@ -70,6 +70,8 @@
 temp_dataframe = pd.DataFrame({'DEPTH':depths, 'BS':BS, 'CALI':CALI, 'DTC':DTC, 'GR':GR, 'GR_CORR':GR_CORR, 'RESD':RESD, 'RESS':RESS, 'SP':SP})
 log_data = pd.concat((temp_dataframe,log_data), ignore_index=True)
 
+# Add the water column 
+
 # Interpolate missing values to enable the cumulative integration 
 log_data.loc[0,'DT_s_per_m'] = (1.0/1500)
 log_data['DT_s_per_m'] = log_data['DT_s_per_m'].interpolate()
@@ -89,7 +91,7 @@
 plt.xticks(np.arange(1500, 6000, 1000))
 plt.xlabel(r'$Vp \quad ms^{-1}$', fontsize=18)
 plt.ylabel("Depth (m)", fontsize=14)
-plt.savefig("/home/mxh909/Desktop/magee_resolution_images/Resolution-1_Vp_plot.pdf",bbox_inches='tight')
+#plt.savefig("/home/mxh909/Desktop/magee_resolution_images/Resolution-1_Vp_plot.pdf",bbox_inches='tight')
 plt.show()
 
 
@@ -99,14 +101,16 @@
 # Cumulativly integrate to form the time depth curve
 from scipy.integrate import cumtrapz 
 time_s = 2*(cumtrapz(y=log_data['DT_s_per_m'], x=log_data['DEPTH'])) # x2 as log is in one-way-time
-time_s += 0.1 # add in the value of the seafloor from the seismic. Entered here rather than calculated, as well is a little upslope.
+#time_s += 0.1 # add in the value of the seafloor from the seismic. Entered here rather than calculated, as well is a little upslope.
 plt.figure()
 plt.plot(time_s, log_data['DEPTH'][:-1])
-plt.gca().invert_yaxis()
 plt.title("Time-depth curve")
 plt.xlabel("Time (s)")
 plt.ylabel("Depth (m)")
-plt.savefig("/home/mxh909/Desktop/magee_resolution_images/Resolution-TD_curve.pdf",bbox_inches='tight')
+plt.xlim(0,1.8)
+plt.ylim(0,2100)
+plt.gca().invert_yaxis()
+plt.savefig("/home/murray/Documents/thesis_python_code/resolution_1/Resolution-TD_curve.pdf",bbox_inches='tight')
 plt.show()
 
 np.savetxt("/home/mxh909/Desktop/magee_resolution_images/Resolution_1_time_depth_curve.txt", np.vstack((log_data['DEPTH'].values[:-1], time_s)).T, header="Depth_m Time_s")