Most wells are vertical, but many are not. All modern wells have a deviation survey, which is converted into a position log, giving the 3D position of the well in space.
welly has a simple way to add a position log in a specific format, and computes a position log from it. You can use the position log to convert between MD and TVD.
First, version check.
In [1]:
import welly
welly.__version__
Out[1]:
In [2]:
from welly import Well
w = Well.from_las("P-130_out.LAS")
w
Out[2]:
In [3]:
%matplotlib inline
w.plot()
There aren't a lot of tricks for handling the input data, which is assumed to be a CSV-like file containing columns like:
MD, inclination, azimuth
For example:
In [4]:
with open('P-130_deviation_survey.csv') as f:
lines = f.readlines()
for line in lines[:6]:
print(line, end='')
Then we can turn that into an ndarray:
In [5]:
import numpy as np
dev = np.loadtxt('P-130_deviation_survey.csv', delimiter=',', skiprows=1, usecols=[0,1,2])
dev[:5]
Out[5]:
You can use any other method to get to an array or pandas.DataFrame like this one.
Then we can add the deviation survey to the well's location attribute. This will automatically convert it into a position log, which is an array containing the x-offset, y-offset, and TVD of the well, in that order.
In [6]:
w.location.add_deviation(dev, td=w.location.tdd)
Now you have the position log:
In [7]:
w.location.position[:5]
Out[7]:
Note that it is irregularly sampled — this is nothing more than the deviation survey (which is MD, INCL, AZI) converted into relative positions (i.e. deltaX, deltaY, deltaZ). These positions are relative to the tophole location.
In [8]:
w.location.md2tvd(1000)
Out[8]:
In [9]:
w.location.tvd2md(998.78525)
Out[9]:
These can also accept an array:
In [10]:
md = np.linspace(0, 300, 31)
w.location.md2tvd(md)
Out[10]:
Note that these are linear in MD, but not in TVD.
In [11]:
w.location.md2tvd([0, 10, 20, 30])
Out[11]:
In general, deviation surveys are considered 'canonical'. That is, they are data recorded in the well. The position log — a set of (x, y, z) points in a linear Euclidean space like (X_UTM, Y_UTM, TVDSS) — is then computed from the deviation survey.
If you have deviation and position log, I recommend loading the deviation survey as above.
If you only have position, in a 3-column array-like called position (say), then you can add it to the well like so:
w.location.position = np.array(position)You can still use the MD-to-TVD and TVD-to-MD converters above, and w.position.trajectory() will work as usual, but you won't have w.position.dogleg or w.position.deviation.
In [12]:
w.location.dogleg[:10]
Out[12]:
Data from Rob:
In [13]:
import pandas as pd
dev = pd.read_csv('deviation.csv')
dev.head(10)
Out[13]:
In [14]:
dev.tail()
Out[14]:
First we'll create an 'empty' well.
In [18]:
x = Well(params={'header': {'name': 'foo'}})
Now add the Location object to the well's location attribute, finally calling its add_deviation() method on the deviation data:
In [19]:
from welly import Location
x.location = Location(params={'kb': 100})
x.location.add_deviation(dev[['MD[m]', 'Inc[deg]', 'Azi[deg]']].values)
In [20]:
np.set_printoptions(suppress=True)
In [21]:
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(15,5))
md_welly = x.location.deviation[:, 0]
# Plot x vs depth
ax.plot(x.location.position[:, 0], md_welly, lw=6, label="welly")
ax.plot(dev['East[m]'], dev['MD[m]'], c='limegreen', label="file")
ax.invert_yaxis()
ax.legend()
Out[21]:
They seem to match well. There's a difference at the top because welly always adds a (0, 0, 0) point to both the deviation and position logs:
In [22]:
x.location.position[:7]
Out[22]:
In plan view, the wells match:
In [19]:
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(*x.location.position[:, :2].T, c='c', lw=5, label="welly")
ax.plot(dev['East[m]'], dev['North[m]'], c='yellow', ls='--', label="file")
#ax.set_xlim(-20, 800); ax.set_ylim(-820, 20)
ax.grid(color='black', alpha=0.2)
To make things a bit more realistic, we can shift to the correct spatial datum, i.e. the (x, y, z) of the top hole, where z is the KB elevation.
We can also adjust the z value to elevation (i.e. negative downwards).
In [20]:
np.set_printoptions(suppress=True, precision=2)
x.location.trajectory(datum=[111000, 2222000, 100], elev=True)
Out[20]:
We can make a 3D plot with this trajectory:
In [21]:
from mpl_toolkits.mplot3d import Axes3D
fig, ax = plt.subplots(figsize=(12, 7), subplot_kw={'projection': '3d'})
ax.plot(*x.location.trajectory().T, lw=3, alpha=0.75)
plt.show()
In [22]:
fig, ax = plt.subplots(figsize=(15,4))
ax.plot(x.location.dogleg, lw=5, label="welly")
ax = plt.twinx(ax=ax)
ax.plot(dev['Dogleg [deg/30m]'], c='limegreen', ls='--', label="file")
ax.text(80, 4, 'file', color='limegreen', ha='right', va='top', size=16)
ax.text(80, 3.5, 'welly', color='C0', ha='right', va='top', size=16)
Out[22]:
They mostly agree.
The position log is computed from the deviation survey with the minimum curvature algorithm, which is fairly standard in the industry. To use a different method, pass method='aa' (average angle) or method='bt' (balanced tangent) directly to Location.compute_position_log() yourself.
Once we have the position log, we still need a way to look up arbitrary depths. To do this, we use a cubic spline fitted to the position log. This should be OK for most 'natural' well paths, but it might break horribly. If you get weird results, you can pass method='linear' to the conversion functions — less accurate but more stable.
In [23]:
dev = [[100, 0, 0],
[200, 10, 45],
[300, 20, 45],
[400, 20, 45],
[500, 20, 60],
[600, 20, 75],
[700, 90, 90],
[800, 90, 90],
[900, 90, 90],
]
In [24]:
z = welly.Well(params={'name': 'foo'})
z.location = welly.Location(params={'kb': 10})
z.location.add_deviation(dev, td=1000, azimuth_datum=20)
In [25]:
z.location.plot_plan()
In [26]:
z.location.plot_3d()