import welly
welly.__version__

"""
Explanation: Location and the deviation survey
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.
End of explanation
"""


from welly import Well

w = Well.from_las("data/P-130_out.LAS")
w

w.plot()

"""
Explanation: Adding deviation to an existing well
First we'll read a LAS and instantiate a well w
End of explanation
"""


with open('data/P-130_deviation_survey.csv') as f:
    lines = f.readlines()

for line in lines[:6]:
    print(line, end='')

"""
Explanation: 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:
End of explanation
"""


import numpy as np

dev = np.loadtxt('data/P-130_deviation_survey.csv', delimiter=',', skiprows=1, usecols=[0,1,2])
dev[:5]

"""
Explanation: Then we can turn that into an ndarray:
End of explanation
"""


w.location.add_deviation(dev, td=w.location.tdd)

"""
Explanation: 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.
End of explanation
"""


w.location.position[:5]

"""
Explanation: Now you have the position log:
End of explanation
"""


w.location.md2tvd(1000)

w.location.tvd2md(998.78525)

"""
Explanation: Note that it is irregularly sampled &mdash; 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. 
MD to TVD and vice versa
We now have the methods md2tvd and tvd2md available to us:
End of explanation
"""


md = np.linspace(0, 300, 31)

w.location.md2tvd(md)

"""
Explanation: These can also accept an array:
End of explanation
"""


w.location.md2tvd([0, 10, 20, 30])

"""
Explanation: Note that these are linear in MD, but not in TVD.
End of explanation
"""


w.location.dogleg[:10]

"""
Explanation: If you have the position log, but no deviation survey
In general, deviation surveys are considered 'canonical'. That is, they are data recorded in the well. The position log &mdash; a set of (x, y, z) points in a linear Euclidean space like (X_UTM, Y_UTM, TVDSS) &mdash; 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.
Dogleg severity
The dogleg severity array is captured in the dogleg attribute:
End of explanation
"""


import pandas as pd

dev = pd.read_csv('data/deviation.csv')
dev.head(10)

dev.tail()

"""
Explanation: Starting from new well
Data from Rob:
End of explanation
"""


x = Well(params={'header': {'name': 'foo'}})

"""
Explanation: First we'll create an 'empty' well.
End of explanation
"""


from welly import Location

x.location = Location(params={'kb': 100})

x.location.add_deviation(dev[['MD[m]', 'Inc[deg]', 'Azi[deg]']].values)

"""
Explanation: Now add the Location object to the well's location attribute, finally calling its add_deviation() method on the deviation data:
End of explanation
"""


np.set_printoptions(suppress=True)

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()

"""
Explanation: Let's see how our new position data compares to what was in the deviation.csv data file:
Compare x, y, and dogleg
End of explanation
"""


x.location.position[:7]

"""
Explanation: 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:
End of explanation
"""


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)

"""
Explanation: In plan view, the wells match:
End of explanation
"""


np.set_printoptions(suppress=True, precision=2)

x.location.trajectory(datum=[111000, 2222000, 100], elev=True)

"""
Explanation: Fit a spline to the position log
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).
End of explanation
"""


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()

"""
Explanation: We can make a 3D plot with this trajectory:
End of explanation
"""


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)

"""
Explanation: Compare doglegs
The deviation.csv file also contains a measure of dogleg severity, which welly also generates (since v0.4.2).
Note that in the current version dogleg severity is in radians, whereas the usual units are degrees per 100 ft or degrees per 30 m. The next release of welly, v0.5, will start using degrees per 30 m by default.
End of explanation
"""


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],
      ]

z = welly.Well()
z.location = welly.Location(params={'kb': 10})
z.location.add_deviation(dev, td=1000, azimuth_datum=20)

z.location.plot_plan()

z.location.plot_3d()

"""
Explanation: Apart from the scaling, they agree. 
Implementation details
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.

Azimuth datum
You can adjust the angle of the azimuth datum with the azimuth_datum keyword argument. The default is zero, which means the azimuths in your survey are in degrees relative to grid north (of your UTM grid, say).
Let's make some fake data like
MD, INCL, AZI
End of explanation
"""


z.location.trajectory(points=20)

"""
Explanation: Trajectory
Get regularly sampled well trajectory with a specified number of points. Assumes there is a position log already, e.g. resulting from calling add_deviation() on a deviation survey. Computed from the position log by scipy.interpolate.splprep().
End of explanation
"""


