Merge pull request #115 from jonnymaserati/feature/modified_tortuosit…

…y_index

Add modified_tortuosity_index

welleng/survey.py

@@ -25,6 +25,7 @@
 from .node import Node
 from .connector import Connector, interpolate_well
 from .visual import figure
+from .units import ureg
 
 
 AZI_REF = ["true", "magnetic", "grid"]
@@ -357,6 +358,15 @@ def __init__(
  self.n = np.array(n) if n is not None else n
  self.e = np.array(e) if e is not None else e
  self.tvd = np.array(tvd) if tvd is not None else tvd
+
+ # start_nev will be overwritten if n, e, tvd data provided
+ if not all((self.n is None, self.e is None, self.tvd is None)):
+ self.start_nev = np.array(
+ [self.n[0], self.e[0], self.tvd[0]]
+ )
+ else:
+ self.start_nev = np.array(start_nev)
+
  self.x = np.array(x) if x is not None else x
  self.y = np.array(y) if y is not None else y
  self.z = np.array(z) if z is not None else z
@@ -479,7 +489,7 @@ def _min_curve(self, vec):
  self.delta_md = mc.delta_md
  self.dls = mc.dls
  self.pos_xyz = mc.poss
- self.pos_nev = get_nev(self.pos_xyz)
+ self.pos_nev = get_nev(self.pos_xyz) + self.start_nev
 
  if self.x is None:
  # self.x, self.y, self.z = (mc.poss + self.start_xyz).T
@@ -716,7 +726,7 @@ def interpolate_survey(self, step=30, dls=1e-8):
  """
  Convenience method for interpolating a Survey object's MD.
  """
- survey_interpolated = interpolate_survey(self, step=30, dls=1e-8)
+ survey_interpolated = interpolate_survey(self, step=step, dls=dls)
  return survey_interpolated
 
  def figure(self, type='scatter3d', **kwargs):
@@ -858,6 +868,270 @@ def set_vertical_section(self, vertical_section_azimuth, deg=True):
  vertical_section_azimuth, deg
  )
 
+ def modified_tortuosity_index(
+ self, rtol=1.0, data=False, **kwargs
+ ):
+ """
+ Convenient method for calculating the Tortuosity Index (TI) using a
+ modified version of the method described in the International
+ Association of Directional Drilling presentation
+ (https://www.iadd-intl.org/media/files/files/47d68cb4/iadd-luncheon-february-22-2018-v2.pdf)
+ by Pradeep Ashok et al.
+
+ Parameters
+ ----------
+ rtol: float
+ Relative tolerance when determining closeness of normal vectors.
+ atol: float
+ Absolute tolerance when determining closeness of normal vectors.
+ data: boolean
+ If true returns a dictionary of properties.
+
+ Returns
+ ti: (n,1) array or dict
+ Array of tortuosity index or a dict of results where:
+ """
+ return modified_tortuosity_index(
+ self, rtol=rtol, data=data, **kwargs
+ )
+
+ def tortuosity_index(self, rtol=1.0, data=False, **kwargs):
+ """
+ A modified version of the Tortuosity Index function originally
+ referenced in an IADD presentation on "Measuring Wellbore
+ Tortuosity" by Pradeep Ashok - https://www.iadd-intl.org/media/
+ files/files/47d68cb4/iadd-luncheon-february-22-2018-v2.pdf with
+ reference to the original paper "A Novel Method for the Automatic
+ Grading of Retinal Vessel Tortuosity" by Enrico Grisan et al.
+ In SPE/IADC-194099-MS there's mention that a factor of 1e7 is
+ applied to the TI result since the results are otherwise very small
+ numbers.
+ Unlike the documented version that uses delta inc and azi for
+ determining continuity and directional intervals (which are not
+ independent and so values are double dipped in 3D), this method
+ determines the point around which the well bore is turning and tests
+ for continuity of these points. As such, this function takes account
+ of torsion of the well bore and demonstrates that actual/drilled
+ trajectories are significantly more tortuous than planned
+ trajectories (intuitively).
+
+ Parameters
+ ----------
+ rtol: float
+ Relative tolerance when determining closeness of normal vectors.
+ atol: float
+ Absolute tolerance when determining closeness of normal vectors.
+ data: boolean
+ If true returns a dictionary of properties.
+
+ Returns
+ ti: (n,1) array
+ Array of tortuosity index or a dict or results.
+ """
+
+ return tortuosity_index(
+ self, rtol=rtol, data=data, **kwargs
+ )
+
+ def directional_difficulty_index(self, **kwargs):
+ """
+ Taken from IADC/SPE 59196 The Directional Difficulty Index - A
+ New Approach to Performance Benchmarking by Alistair W. Oag et al.
+
+ Returns
+ -------
+ data: (n) array of floats
+ The ddi for each survey station.
+ """
+
+ return directional_difficulty_index(self, **kwargs)
+
+
+def modified_tortuosity_index(
+ survey, rtol=1.0, data=False, **kwargs
+):
+ """
+ Method for calculating the Tortuosity Index (TI) using a modified
+ version of the method described in the International Association of
+ Directional Drilling presentation (https://www.iadd-intl.org/media/files/files/47d68cb4/iadd-luncheon-february-22-2018-v2.pdf)
+ by Pradeep Ashok et al.
+ """
+ # set default params
+ coeff = kwargs.get('coeff', 1.0) # for testing dimensionlessness
+ kappa = kwargs.get('kapa', 1)
+
+ continuous, starts, mds, locs, n_sections, n_sections_arr = _get_ti_data(
+ survey, rtol
+ )
+
+ l_cs = (
+ survey.md[1:] - mds[n_sections_arr - 1]
+ ) / coeff
+ l_xs = np.linalg.norm(
+ survey.pos_nev[1:]
+ - np.array(locs)[n_sections_arr - 1],
+ axis=1
+ ) / coeff
+ b = (
+ (l_cs / l_xs) - 1
+ ) / l_cs
+ # )
+
+ cumsum = 0
+ a = []
+ for n in n_sections:
+ a.extend(
+ b[n_sections_arr == n]
+ + cumsum
+ )
+ cumsum = a[-1]
+ a = np.array(a)
+
+ mti = np.hstack((
+ np.array([0.0]),
+ (
+ # 1
+ (n_sections_arr / (n_sections_arr + 1))
+ * (kappa * ((survey.md[1:] - survey.md[0]) / coeff))
+ # * (kappa / (np.cumsum(survey.dogleg)[1:] + 1))
+ * a
+ )
+ ))
+
+ if data:
+ return {
+ 'continuous': continuous, 'starts': starts, 'mds': mds,
+ 'locs': locs, 'n_sections': n_sections,
+ 'n_sections_arr': n_sections_arr, 'l_cs': l_cs, 'l_xs': l_xs,
+ 'mti': mti
+ }
+
+ return mti
+
+
+def tortuosity_index(survey, rtol=1.0, data=False, **kwargs):
+ """
+ Method for calculating the Tortuosity Index (TI) as described in the 
+ International Association of Directional Drilling presentation
+ (https://www.iadd-intl.org/media/files/files/47d68cb4/iadd-luncheon-february-22-2018-v2.pdf)
+ by Pradeep Ashok et al.
+ """
+ # set default params
+ coeff = kwargs.get('coeff', 0.3048) # for testing dimensionlessness
+ kappa = kwargs.get('kapa', 1e7)
+
+ continuous, starts, mds, locs, n_sections, n_sections_arr = _get_ti_data(
+ survey, rtol
+ )
+
+ l_cs = (survey.md[1:] - mds[n_sections_arr - 1]) / coeff
+ l_xs = np.linalg.norm(
+ survey.pos_nev[1:]
+ - np.array(locs)[n_sections_arr - 1],
+ axis=1
+ ) / coeff
+ b = (
+ (l_cs / l_xs) - 1
+ )
+
+ cumsum = 0
+ a = []
+ for n in n_sections:
+ a.extend(
+ b[n_sections_arr == n]
+ + cumsum
+ )
+ cumsum = a[-1]
+ a = np.array(a)
+
+ ti = np.hstack((
+ np.array([0.0]),
+ (
+ (n_sections_arr / (n_sections_arr + 1))
+ * (kappa / (survey.md[1:] - survey.md[0] / coeff))
+ * a
+ )
+ ))
+
+ if data:
+ return {
+ 'continuous': continuous, 'starts': starts, 'mds': mds,
+ 'locs': locs, 'n_sections': n_sections,
+ 'n_sections_arr': n_sections_arr, 'l_cs': l_cs, 'l_xs': l_xs,
+ 'ti': ti
+ }
+
+ return ti
+
+
+def directional_difficulty_index(survey, **kwargs):
+ """
+ Taken from IADC/SPE 59196 The Directional Difficulty Index - A
+ New Approach to Performance Benchmarking by Alistair W. Oag et al.
+ Parameters
+ ----------
+ survey: welleng.survey.Survey object
+ data: bool
+ If True, returns the ddi at each survey station.
+ Returns
+ -------
+ ddi: float
+ The ddi for the well at well (at TD).
+ data: (n) array of floats
+ The ddi for each survey station.
+ """
+ with np.errstate(divide='ignore', invalid='ignore'):
+ ddi = np.nan_to_num(np.log10(
+ (
+ (survey.md * ureg.meters).to('ft').m
+ * (
+ np.linalg.norm(
+ (survey.n, survey.e), axis=0
+ ) * ureg.meters
+ ).to('ft').m
+ * np.cumsum(np.degrees(survey.dogleg))
+ )
+ / (survey.tvd * ureg.meters).to('ft').m
+ ), nan=0.0, posinf=0.0, neginf=0.0)
+
+ return ddi
+
+
+def _get_ti_data(survey, rtol):
+ # tol_0 = np.full((len(survey.md) - 2, 3), rtol)
+ # delta_md = (survey.md[2:] - survey.md[1:-1]).reshape(-1, 1)
+ # delta_norm = survey.normals[1:] - survey.normals[:-1]
+ # tol_norm = tol_0 * np.sin(delta_md / 30)
+
+ # continuous = (
+ # (delta_norm / 30 * delta_md <= tol_norm).all(axis=-1) | (
+ # np.isnan(survey.normals[1:]).all(axis=-1)
+ # & np.isnan(survey.normals[:-1]).all(axis=-1)
+ # )
+ # )
+ continuous = np.all(
+ np.isclose(
+ survey.normals[:-1],
+ survey.normals[1:],
+ equal_nan=True,
+ rtol=rtol, atol=rtol
+ ), axis=-1
+ )
+
+ starts = np.concatenate((
+ np.array([0]),
+ np.where(continuous == False)[0] + 1,
+ ))
+
+ mds = survey.md[starts]
+ locs = survey.pos_nev[starts]
+ n_sections = np.arange(1, len(starts) + 1, 1)
+ n_sections_arr = np.searchsorted(mds, survey.md[1:])
+
+ return (
+ continuous, starts, mds, locs, n_sections, n_sections_arr
+ )
+
 
 class TurnPoint:
  def __init__(
@@ -1280,14 +1554,16 @@ def get_sections(survey, rtol=1e-1, atol=1e-1, dls_cont=False, **targets):
 
  # check for DLS continuity
  if not dls_cont:
- dls_cont = [True] * (len(survey.dls) - 2)
+ # dls_cont = [True] * (len(survey.dls) - 2)
+ dls_cont = np.full(len(survey.dls) - 2, True)
  else:
  upper = np.around(survey.dls[1:-1], decimals=2)
  lower = np.around(survey.dls[2:], decimals=2)
- dls_cont = [
- True if u == l else False
- for u, l in zip(upper, lower)
- ]
+ # dls_cont = [
+ # True if u == l else False
+ # for u, l in zip(upper, lower)
+ # ]
+ dls_cont = np.equal(upper, lower)
 
  continuous = np.all((
  np.all(