Almost finished Resolution stuff. misfit for functions finished and debugged, plots made, etc.

Author: Murray2015

Resolution_forward_model.pdf

@@ -11,48 +11,6 @@
 import matplotlib.patches as patches
 import numpy as np
 
-## Read in wireline log data from .las files
-log_data = pd.read_csv("Wireline/resolution-1_final.las", delim_whitespace=True, 
-                       header=None, names=['DEPTH', 'BS', 'CALI', 'DTC', 'GR', 'GR_CORR', 
-                       'RESD', 'RESS', 'SP'], na_values=-999.2500, skiprows=56)
-                       
-## Keep a copy incase of need to check unaltered params. 
-log_data_original = log_data.copy()
-
-## Fill log data with dummy values for water and sed. Create a dataframe of the filled values, 
-## then concat to top of log_data.
-# Water depth is looked up from paper logs
-water_depth = 64.0 # meters
-# Get the first recorded depth in the log (note shallow unconsolidate section rarely logged) 
-first_depth_in_well = log_data['DEPTH'][0]
-# Find the depth sampling interval 
-depth_sample_interval = np.round(log_data['DEPTH'][1] - log_data['DEPTH'][0], decimals=4)
-# Create a vector of depths to fill 
-depths = np.arange(start=0,stop=first_depth_in_well, step=depth_sample_interval)
-# Create a vector of bit sizes
-BS = np.repeat(log_data['BS'][0],len(depths))
-# Create a vector of NaNs for logs not used in analysis.
-CALI= np.repeat(np.nan,len(depths))
-GR= np.repeat(np.nan,len(depths))
-GR_CORR= np.repeat(np.nan,len(depths))
-RESD= np.repeat(np.nan,len(depths))
-RESS= np.repeat(np.nan,len(depths))
-SP= np.repeat(np.nan,len(depths))
-# Create two vectors of velocity in the water and shalow sed, and concat.
-DTC_water = np.repeat(1/(1500*3.28084/1e6),water_depth)
-DTC_sed = np.repeat(1/(1800*3.28084/1e6),len(depths)-water_depth) 
-DTC=np.insert(DTC_sed, 0, DTC_water)
-# Make a temp dataframe of the newly generated data and NaNs, and then concat to make the full data.
-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)
-
-# Interpolate missing values in the logs (different tools turned on at different times) 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()
-log_data['DTC'] = log_data['DTC'].interpolate()
-
-# Make a Vp log, to plot and also to later find the sill velocity from. 
-log_data['VP_m_per_s'] = 1/(log_data['DTC'] * (1e-6) / 0.3048)
 
 #########################################
 ############### Load files ##############
@@ -206,23 +164,21 @@ def detrend_fold_profiles(fold):
 ##################################
 ######### Forward model ##########
 ##################################
+## First we make forward models. We use the Southland horizon just 
+## below the forced fold, as it shows no signs of erosion. 
 greatest_sill_vel = 7000
 least_sill_vel = 4500
 greatest_phi0 = 0.7
-least_phi0=0.25
+least_phi0=0.20
 greatest_lambda = 3.7
 least_lambda = 1.4
 
 ## Make upper limit of fold profile 
-greatest_fold_decompacted = decompact_fold(top_forced_fold_clip, tolerance=0.0001, phi0=greatest_phi0, Lambda=greatest_lambda)
 greatest_southland_decompacted = decompact_fold(southland_clip, tolerance=0.0001, phi0=greatest_phi0, Lambda=greatest_lambda)
-greatest_fold_detrended = detrend_fold_profiles(greatest_fold_decompacted)
 greatest_southland_detrended = detrend_fold_profiles(greatest_southland_decompacted)
 
 ## Make lower limit of fold profile 
-least_fold_decompacted = decompact_fold(top_forced_fold_clip, tolerance=0.0001, phi0=least_phi0, Lambda=least_lambda)
 least_southland_decompacted = decompact_fold(southland_clip, tolerance=0.0001, phi0=least_phi0, Lambda=least_lambda)
-least_fold_detrended = detrend_fold_profiles(least_fold_decompacted)
 least_southland_detrended = detrend_fold_profiles(least_southland_decompacted)
 
 ## Make upper limit of sill thickness 
@@ -232,7 +188,6 @@ def detrend_fold_profiles(fold):
 least_sill_amp_from_time = time_sill_amp_extraction_2d(top_sill_time, base_sill_time, sill_vel=least_sill_vel)
 
 ## Plot 
-plt.fill_between(greatest_fold_detrended['trace'], least_fold_detrended['decompacted_profile_detrended_km']*1000, greatest_fold_detrended['decompacted_profile_detrended_km']*1000, color='b', alpha=0.2, label='Decompacted fold')
 plt.fill_between(greatest_southland_detrended['trace'], least_southland_detrended['decompacted_profile_detrended_km']*1000, greatest_southland_detrended['decompacted_profile_detrended_km']*1000, color='g', alpha=0.2, label='Decompacted Southland')
 plt.fill_between(greatest_sill_amp_from_time['trace'], least_sill_amp_from_time['sill_amp_meters'], greatest_sill_amp_from_time['sill_amp_meters'], color='r', alpha=0.2, label='Sill thickness')
 plt.ylabel("Sill thickness / Fold amplitude (meters)")
@@ -246,22 +201,24 @@ def detrend_fold_profiles(fold):
 ##################################
 ######### Inverse model ##########
 ##################################
+## Secondly, we solve the inverse problem to find the best fitting decompaction
+## parameters to fit the sill. 
 
 ## Set vars 
 greatest_sill_vel = 7000
 least_sill_vel = 4500
 greatest_phi0 =0.7
-least_phi0=0.25
-phi0_decimation = 0.05
+least_phi0=0.20
+phi0_decimation = 0.01
 greatest_lambda = 3.7
 least_lambda = 1.4
-lambda_decimation = 0.1
+lambda_decimation = 0.05
 phi0_vec = np.arange(start=least_phi0, stop=greatest_phi0, step=phi0_decimation)
 lambda_vec = np.arange(start=least_lambda, stop=greatest_lambda, step=lambda_decimation)
 grid_search = np.zeros((len(phi0_vec), len(lambda_vec)))
 
 
-## Misfit function 
+## Misfit functions 
 def misfit(sill, forced_fold):
     ''' 
     Function which returns the misfit between a sill curve and a 
@@ -272,54 +229,172 @@ def misfit(sill, forced_fold):
         if sill['trace'][i] in forced_fold['trace'].values:
             total_misfit += (sill['sill_amp_meters'][i]/1000.0) - forced_fold['decompacted_profile_detrended_km'].loc[forced_fold['trace']==sill['trace'][i]].values
     return total_misfit
-    
 #misfit(sill_amp_from_time, fold_detrended)
 
+## Misfit function 
+def RMS_misfit(sill, forced_fold):
+    ''' 
+    Function which returns the misfit between a sill curve and a 
+    forced fold curve.
+    '''
+    total_misfit = 0 
+    counter = 0
+    for i in range(len(sill)):
+        if sill['trace'][i] in forced_fold['trace'].values:
+            counter += 1
+            total_misfit += np.power((sill['sill_amp_meters'][i]/1000.0) - forced_fold['decompacted_profile_detrended_km'].loc[forced_fold['trace']==sill['trace'][i]].values, 2)
+    rms_misfit = np.power((total_misfit / counter), 0.5)
+    return rms_misfit
+RMS_misfit(sill_amp_from_time, fold_detrended)
+
 ## Grid search
 ## Outer loop - phi0
+sill_amp_gs = time_sill_amp_extraction_2d(top_sill_time, base_sill_time, sill_vel=6500)
 for i_loc, i in enumerate(phi0_vec):
     ## Inner loop - lambda
     for j_loc, j in enumerate(lambda_vec):
         ## Note, could do a third loop for sill vel
-        sill_amp_gs = time_sill_amp_extraction_2d(top_sill_time, base_sill_time, sill_vel=5500)
-        fold_gs = detrend_fold_profiles(decompact_fold(top_forced_fold_clip, tolerance=0.0001, phi0=i, Lambda=j))
-        grid_search[i_loc,j_loc] = misfit(sill_amp_gs, fold_gs) ## fill this with the misfit function. 
+        fold_gs = detrend_fold_profiles(decompact_fold(southland_clip, tolerance=0.0001, phi0=i, Lambda=j))
+        grid_search[len(phi0_vec)-(i_loc+1),j_loc] = RMS_misfit(sill_amp_gs, fold_gs) ## fill this with the misfit function. 
         print('i = ', i_loc, ', j = ', j_loc, ', total = ', j_loc+(len(lambda_vec)*i_loc), 'out of ', len(phi0_vec)*len(lambda_vec))
 
-
-
-## Plot the grid searc
-g2 = abs(grid_search)
+## Find the optimum values 
+#g2 = abs(grid_search)
+g2 = grid_search
 max_loc = np.where(g2 == g2.min())
 opt_phi0 = phi0_vec[::-1][max_loc[0]]
 opt_lambda = lambda_vec[max_loc[1]]
-plt.figure()
-plt.imshow(g2, aspect='equal', cmap='jet_r', extent=[least_lambda, greatest_lambda, least_phi0, greatest_phi0])
+
+## Plot the grid searc
+plt.figure(figsize=(4,3))
+plt.imshow(g2, aspect='auto', cmap='jet_r', extent=[least_lambda, greatest_lambda, least_phi0, greatest_phi0])
 #plt.scatter(y=[opt_phi0], x=[opt_lambda], c='r', s=40)
-plt.colorbar(cmap='jet_r')
+cbar = plt.colorbar(cmap='jet_r')
+cbar.ax.set_ylabel('RMS Misfit')
 plt.ylabel(r'$\Phi_0$')
 plt.xlabel(r'$\lambda$')
-plt.savefig("/home/murray/Documents/thesis_python_code/Resolution_inverse_grid.pdf",bbox_inches='tight')
+#plt.savefig("/home/murray/Documents/thesis_python_code/Resolution_inverse_grid.pdf",bbox_inches='tight')
 plt.show()
 
-## Next to do - change code so the white padding at the base of the image goes. 
-# add contours
-# replot with more squares 
-# find why lambda axes doesn't go to lambda max (def a bug somewhere)
+plt.figure(figsize=(5,4))
+plt.contour(g2, 20, linewidths=0.5, colors='k', aspect='auto', extent=[least_lambda, greatest_lambda, greatest_phi0, least_phi0])
+plt.contourf(g2, 20, aspect='auto', cmap='jet_r', extent=[least_lambda, greatest_lambda, greatest_phi0, least_phi0])
+cbar = plt.colorbar(cmap='jet_r')
+cbar.ax.set_ylabel('RMS Misfit')
+plt.scatter(opt_lambda, opt_phi0, marker='o', s=100, c='k')
+plt.ylabel(r'$\Phi_0$')
+plt.xlabel(r'$\lambda$')
+#plt.savefig("/home/murray/Documents/thesis_python_code/Resolution_inverse_grid.pdf",bbox_inches='tight')
+plt.show()
 
 ## Plot fold profile with minimum misfit 
-fold_decompacted = decompact_fold(top_forced_fold_clip, tolerance=0.0001, phi0=opt_phi0, Lambda=opt_lambda)
 southland_decompacted = decompact_fold(southland_clip, tolerance=0.0001, phi0=opt_phi0, Lambda=opt_lambda)
-
-fold_detrended = detrend_fold_profiles(fold_decompacted)
 southland_detrended = detrend_fold_profiles(southland_decompacted)
-
-plt.plot(fold_detrended['trace'], fold_detrended['decompacted_profile_detrended_km']*1000, label='Decompacted fold')
 plt.plot(southland_detrended['trace'], southland_detrended['decompacted_profile_detrended_km']*1000, label='Decompacted Southland')
 plt.fill_between(greatest_sill_amp_from_time['trace'], least_sill_amp_from_time['sill_amp_meters'], greatest_sill_amp_from_time['sill_amp_meters'], color='r', alpha=0.2, label='Sill thickness')
 plt.ylabel("Sill thickness / Fold amplitude (meters)")
 plt.xlabel("Trace")
 plt.legend(loc=8)
-plt.savefig("/home/murray/Documents/thesis_python_code/Resolution_inverse_model.pdf",bbox_inches='tight')
+plt.savefig("/home/murray/Documents/thesis_python_code/Resolution_optimum_inverse_model.pdf",bbox_inches='tight')
 plt.show()
 
+### Finally, make a theoretical porosity curve from the above results, fit to the 
+### porosity curve with a 1D gridsearch, and plot the misfit. 
+
+### Read in well logs and apply corrections 
+## Read in wireline log data from .las files
+log_data = pd.read_csv("Wireline/resolution-1_final.las", delim_whitespace=True, 
+                       header=None, names=['DEPTH', 'BS', 'CALI', 'DTC', 'GR', 'GR_CORR', 
+                       'RESD', 'RESS', 'SP'], na_values=-999.2500, skiprows=56)
+                       
+## Keep a copy incase of need to check unaltered params. 
+log_data_original = log_data.copy()
+
+## Fill log data with dummy values for water and sed. Create a dataframe of the filled values, 
+## then concat to top of log_data.
+# Water depth is looked up from paper logs
+water_depth = 64.0 # meters
+# Get the first recorded depth in the log (note shallow unconsolidate section rarely logged) 
+first_depth_in_well = log_data['DEPTH'][0]
+# Find the depth sampling interval 
+depth_sample_interval = np.round(log_data['DEPTH'][1] - log_data['DEPTH'][0], decimals=4)
+# Create a vector of depths to fill 
+depths = np.arange(start=0,stop=first_depth_in_well, step=depth_sample_interval)
+# Create a vector of bit sizes
+BS = np.repeat(log_data['BS'][0],len(depths))
+# Create a vector of NaNs for logs not used in analysis.
+CALI= np.repeat(np.nan,len(depths))
+GR= np.repeat(np.nan,len(depths))
+GR_CORR= np.repeat(np.nan,len(depths))
+RESD= np.repeat(np.nan,len(depths))
+RESS= np.repeat(np.nan,len(depths))
+SP= np.repeat(np.nan,len(depths))
+# Create two vectors of velocity in the water and shalow sed, and concat.
+DTC_water = np.repeat(1/(1500*3.28084/1e6),water_depth)
+DTC_sed = np.repeat(1/(1800*3.28084/1e6),len(depths)-water_depth) 
+DTC=np.insert(DTC_sed, 0, DTC_water)
+# Make a temp dataframe of the newly generated data and NaNs, and then concat to make the full data.
+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)
+
+# Interpolate missing values in the logs (different tools turned on at different times) 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()
+log_data['DTC'] = log_data['DTC'].interpolate()
+
+# Get array of depths
+log_depths = log_data['DEPTH'].values/1000.00
+
+# Apply formula to find the poro at those depths 
+theoretical_poro = opt_phi0 * np.exp(-(log_depths) / opt_lambda)
+#plt.figure(figsize=(2,4))
+#plt.plot(theoretical_poro, log_depths)
+#plt.gca().invert_yaxis()
+#plt.xlim(0,0.6)
+#plt.show()
+
+# Define RHG function to change sonic to poro 
+def RHG(well_log, matrix_vel):
+    '''Function to calculate Raymer-Hunt-Gardner porosity from Sonic log'''
+    return (5.0/8)*((log_data['DTC'] - matrix_vel)/(log_data['DTC']))
+    
+## Define misfit function for RHG vs. theoretical poro
+def RHG_misfit(log_poro, theoretical_poro, range_of_interest=(2145, 12100)):
+    '''Function to calculate the misfit between a theoretical porosity log and a calculated
+    porosity log ''' 
+    misfit_vec = np.power((log_poro[range_of_interest[0]:range_of_interest[1]] - theoretical_poro[range_of_interest[0]:range_of_interest[1]]), 2)
+    return np.sqrt(np.divide(np.sum(misfit_vec), len(misfit_vec)))
+#RHG_misfit(rhg, theoretical_poro)
+    
+## Plot test curves of theoretical poro and cal poro, to check everything
+rhg = RHG(log_data, 55)
+plt.figure(figsize=(4,6))
+plt.plot(theoretical_poro, log_depths, label='Theoretical poro')
+plt.plot(rhg, log_depths, label='Calc poro')
+plt.plot(rhg[2145:12100], log_depths[2145:12100], label='Fitted region')
+plt.legend()
+plt.xlabel(r"$\phi$")
+plt.ylabel('Depth, km')
+plt.gca().invert_yaxis()
+plt.xlim(0,0.6)
+plt.show()
+
+## Define range of matrix velocities 
+matrix_vels = np.arange(start=40, stop=80, step=1)
+
+## Define vector of zeroes for misfit results 
+misfits = np.zeros(len(matrix_vels))
+
+## Loop over matrix vels to find the misfit for the matrix velocities 
+for i_loc, i in enumerate(matrix_vels):
+    misfits[i_loc] = RHG_misfit(RHG(log_data, i), theoretical_poro, range_of_interest=(0, 12100))
+
+plt.plot(matrix_vels, abs(misfits))
+plt.axvline(x=matrix_vels[np.argmin(misfits)], c='r')
+plt.text(54, 0.1, r'$Minimum = 52 \mu s / ft $', color='r')
+plt.ylabel("RMSE misfit")
+plt.xlabel(r"$Matrix \ velocity \ (\mu s / ft)$")
+plt.show()
+
+
+# Plot curve of misfit vs matrix sonic velocity, and add static ranges for sensible velcities