From 3fb3d33bad66bbbba5c1122d21e51ff2862d4b93 Mon Sep 17 00:00:00 2001
From: jonnymaserati <jonnycorcutt@gmail.com>
Date: Mon, 6 Jun 2022 13:42:07 +0200
Subject: [PATCH] Add modified_tortuosity_index

---
 welleng/survey.py  | 290 +++++++++++++++++++++++++++++++++++++++++++--
 welleng/utils.py   |  68 +++++++++++
 welleng/version.py |   2 +-
 3 files changed, 352 insertions(+), 8 deletions(-)

diff --git a/welleng/survey.py b/welleng/survey.py
index 33a8de7..07e2b2e 100644
--- a/welleng/survey.py
+++ b/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(
diff --git a/welleng/utils.py b/welleng/utils.py
index 244f154..cfda41c 100644
--- a/welleng/utils.py
+++ b/welleng/utils.py
@@ -1,4 +1,6 @@
 import numpy as np
+from scipy.spatial.transform import Rotation as R
+
 try:
     from numba import njit
     NUMBA = True
@@ -457,3 +459,69 @@ def errors_from_cov(cov, data=False):
         }
 
     return np.array([nn, ne, nv, ee, ev, vv]).T
+
+
+class Arc:
+    def __init__(self, dogleg, radius):
+        """
+        Generates a generic arc that can be transformed with a specific pos
+        and vec via a transform method. The arc is initialized at a local
+        origin and kicks off down and to the north (assuming an NEV coordinate
+        system).
+
+        Parameters
+        ----------
+        dogleg: float
+            The sweep angle of the arc in radians.
+        radius: float
+            The radius of the arc in meters.
+
+        Returns
+        -------
+        arc: Arc object
+        """
+        self.dogleg = dogleg
+        self.radius = radius
+        self.delta_md = dogleg * radius
+
+        self.pos = np.array([
+            np.cos(dogleg),
+            0.,
+            np.sin(dogleg)
+        ]) * radius
+        self.pos[0] = radius - self.pos[0]
+
+        self.vec = np.array([
+            np.sin(dogleg),
+            0.,
+            np.cos(dogleg)
+        ])
+
+    def transform(self, toolface, pos=None, vec=None, target=False):
+        if vec is None:
+            vec = np.array([0., 0., 1.])
+        if target:
+            vec *= -1
+        inc, azi = get_angles(vec, nev=True).reshape(2)
+        angles = [
+            toolface,
+            inc,
+            azi
+        ]
+        r = R.from_euler('zyz', angles, degrees=False)
+
+        pos_new, vec_new = r.apply(np.vstack((self.pos, self.vec)))
+
+        if pos is not None:
+            pos_new += pos
+        if target:
+            vec_new *= -1
+
+        return (pos_new, vec_new)
+
+
+def get_arc(dogleg, radius, toolface, pos=None, vec=None, target=False):
+    arc = Arc(dogleg, radius)
+    pos_new, vec_new = arc.transform(toolface, pos, vec, target)
+
+    return (pos_new, vec_new, arc.delta_md)
diff --git a/welleng/version.py b/welleng/version.py
index 2d80271..b94cbb0 100644
--- a/welleng/version.py
+++ b/welleng/version.py
@@ -1 +1 @@
-__version__ = '0.4.9'
+__version__ = '0.4.10'