One of the primary purposes of programming is automatization. Whenever you work on data, you should try to avoid any manual intervention such as copying and pasting from one file to another, for example into an excel spreadsheet. This introduces a source of error and severely limits the amount of data you can process. Further, whenever you change parts of your program, it is likely that the manual steps have to be repeated. For all your work on data, good practice is to work on the rawest form of files you get, i.e., the ones supplied by a contractor or produced by the measurement apparatus.
One of the most important data sources in subsurface geoscience are well logs. These are cylindrical logging tools of a length between one and several dozens of meters that are being run up and down wellbores along a wireline and measure, collect and transmit rock formation data such as electrical resistivity, radiation, temperature and more, as a function of depth. This allows for accurate chracterization of lithologic layers, rock properties, oil saturation and more, to a high resolution.
In [1]:
%pylab inline
import os
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
A common well log file in the raw LAS format looks as follows (we print the first 55 lines and skip the rest of the data section):
In [2]:
fname = '1044944685.las'
lines_to_display = 55
with open(fname, 'r') as f:
for i in range(lines_to_display):
print f.readline().rstrip()
We'll try to keep things simpel here. Say we know that the data array always starts at line 49, then we can load the entire dataset by using loadtxt, telling the method to skip everything before that line/row.
In [3]:
las_data = np.loadtxt(fname, skiprows=49)
print 'readings/rows, logs/columns:'
print las_data.shape
print 'values:'
print las_data
Which creates a two-dimensional array containing all logs (columns) and their readings (rows). If we further know that Depth is in column 0, and our log of interest RHOB is in column 4, we can extract the corresponding data:
In [4]:
las_depth = las_data[:,0] # all readings of log 0
las_rhob = las_data[:,4] # all readings of log 4
print las_depth
print las_rhob
We have now extracted our arrays of interest, and can plot the measured rock density as a function of depth:
In [5]:
plt.plot(las_rhob, las_depth)
# so that depth, albeit positive, increases downward
plt.gca().invert_yaxis()
# cosmetics
plt.xlim(0.,3.)
plt.xlabel('Specific Bulk Density [-]')
plt.ylabel('Depth [ft]')
plt.show()
We will operate on the extracted data now, and plot the result. The simple operation is
$\sigma(z) = g \int \rho dz$
to compute the overburden stress resulting from the lithostatic (rock) weight. We will therefore first convert the density from specific to [kg m-3] and the depth from [ft] to [m]
In [6]:
rhob = las_rhob*1.0E3
depth = las_depth*0.3048
and then use a scipy function to perform the integration
In [7]:
# sigma(z) = g int rho dz
overburden_stress_pa = 9.81 * integrate.cumtrapz(rhob, depth, initial=0)
overburden_stress_mpa = overburden_stress_pa/1.0E6
In [8]:
plt.plot(overburden_stress_mpa, depth)
# so that depth, albeit positive, increases downward
plt.gca().invert_yaxis()
# cosmetics
plt.xlabel('$\sigma_v$ [MPa]')
plt.ylabel('Depth [m]')
plt.show()
DCAL measures the diamater of the wellbore in [in]. Plot its readings vs depth, and integrate the cross-sectional area (use np.pi as $\pi$) to plot the cumulative borehole volume in [m3] vs [m]. (Download the notebook and las file here)
In [11]:
las_dcal = las_data[:,3] # all readings of log 3
In [12]:
plt.plot(las_dcal, las_depth)
# so that depth, albeit positive, increases downward
plt.gca().invert_yaxis()
# cosmetics
plt.xlabel('Wellbor Diameter [in]')
plt.ylabel('Depth [ft]')
plt.show()
In [14]:
dcal = las_dcal*0.0254 # from [in] to [m]
xarea = np.pi * dcal**2/4. # cross-sectional area [m2], dcal=diameter
In [15]:
# V(z) = int area dz
volume_m3 = integrate.cumtrapz(xarea, depth, initial=0)
In [18]:
plt.plot(volume_m3, depth)
# so that depth, albeit positive, increases downward
plt.gca().invert_yaxis()
# cosmetics
plt.xlabel('Wellbore Volume [m3]')
plt.ylabel('Depth [m]')
plt.show()
In [ ]: