[4D]: New rich wellbore trajectory provider

backend_py/libs/services/src/webviz_services/smda_access/drogon/_drogon_well_data.py

@@ -1,4 +1,5 @@
 from typing import List
+import numpy as np
 
 from ..types import WellboreTrajectory, WellboreHeader, WellborePick, StratigraphicColumn
 
@@ -28,8 +29,12 @@ def get_drogon_well_headers() -> List[WellboreHeader]:
             wellbore_purpose="production",
             wellbore_status="active",
             current_track=1,
-            tvd_max=1774.5,
+            md_min=0.0,
             md_max=1799.5,
+            md_unit="m",
+            tvd_min=-25.0,
+            tvd_max=1774.5,
+            tvd_unit="m",
             kickoff_depth_md=None,
             kickoff_depth_tvd=None,
             parent_wellbore=None,
@@ -46,8 +51,12 @@ def get_drogon_well_headers() -> List[WellboreHeader]:
             wellbore_purpose="production",
             wellbore_status="active",
             current_track=1,
-            tvd_max=1656.9874,
+            md_min=0.0,
             md_max=3578.5,
+            md_unit="m",
+            tvd_min=-49.0,
+            tvd_max=1656.9874,
+            tvd_unit="m",
             kickoff_depth_md=None,
             kickoff_depth_tvd=None,
             parent_wellbore=None,
@@ -56,6 +65,22 @@ def get_drogon_well_headers() -> List[WellboreHeader]:
 
 
 def get_drogon_well_trajectories() -> List[WellboreTrajectory]:
+    # Original second wellbore data
+    original_tvd_msl = np.array([-49.0, 1293.4185, 1536.9384, 1616.4998, 1630.5153, 1656.9874])
+    original_easting = np.array([463256.911, 463564.402, 463637.925, 463690.658, 463910.452, 464465.876])
+    original_northing = np.array([5930542.294, 5931057.803, 5931184.235, 5931278.837, 5931688.122, 5932767.761])
+    original_md = np.array([0.0, 1477.0, 1761.5, 1899.2601, 2363.9988, 3578.5])
+
+    # Create 100x more sample points using linear interpolation
+    num_original_points = len(original_md)
+    num_interpolated_points = (num_original_points - 1) * 100 + 1
+
+    # Interpolate based on MD (measured depth) as the independent variable
+    md_interp = np.linspace(original_md[0], original_md[-1], num_interpolated_points)
+    tvd_msl_interp = np.interp(md_interp, original_md, original_tvd_msl)
+    easting_interp = np.interp(md_interp, original_md, original_easting)
+    northing_interp = np.interp(md_interp, original_md, original_northing)
+
     return [
         WellboreTrajectory(
             wellbore_uuid="drogon_vertical",
@@ -68,10 +93,10 @@ def get_drogon_well_trajectories() -> List[WellboreTrajectory]:
         WellboreTrajectory(
             wellbore_uuid="drogon_horizontal",
             unique_wellbore_identifier="55/33-A-4",
-            tvd_msl_arr=[-49.0, 1293.4185, 1536.9384, 1616.4998, 1630.5153, 1656.9874],
-            md_arr=[0.0, 1477.0, 1761.5, 1899.2601, 2363.9988, 3578.5],
-            easting_arr=[463256.911, 463564.402, 463637.925, 463690.658, 463910.452, 464465.876],
-            northing_arr=[5930542.294, 5931057.803, 5931184.235, 5931278.837, 5931688.122, 5932767.761],
+            tvd_msl_arr=tvd_msl_interp.tolist(),
+            md_arr=md_interp.tolist(),
+            easting_arr=easting_interp.tolist(),
+            northing_arr=northing_interp.tolist(),
         ),
     ]
 

backend_py/libs/services/src/webviz_services/smda_access/smda_access.py

@@ -175,8 +175,12 @@ async def get_wellbore_headers_async(self, field_identifier: str) -> List[Wellbo
             "well_northing",
             "depth_reference_point",
             "depth_reference_elevation",
+            "tvd_min",
             "tvd_max",
+            "tvd_unit",
+            "md_min",
             "md_max",
+            "md_unit",
         ]
         params = {
             "_projection": ",".join(projection),
@@ -200,6 +204,7 @@ async def get_wellbore_headers_async(self, field_identifier: str) -> List[Wellbo
             "kickoff_depth_tvd",
             "parent_wellbore",
         ]
+
         params = {
             "_projection": ",".join(projection),
             "_sort": "unique_wellbore_identifier",

backend_py/libs/services/src/webviz_services/smda_access/types.py

@@ -87,8 +87,12 @@ class WellboreHeader(BaseModel):
     wellbore_purpose: str | None
     wellbore_status: str | None
     current_track: int
-    tvd_max: float
-    md_max: float
+    md_min: float | None = None
+    md_max: float | None = None
+    md_unit: str | None = None
+    tvd_min: float | None = None
+    tvd_max: float | None = None
+    tvd_unit: str | None = None
     kickoff_depth_md: float | None
     kickoff_depth_tvd: float | None
     parent_wellbore: str | None

backend_py/libs/services/src/webviz_services/utils/surface_helpers.py

@@ -1,9 +1,15 @@
 from dataclasses import dataclass
+from enum import StrEnum
 
 import numpy as np
 import xtgeo
+from xtgeo import _cxtgeo
+from xtgeo.common.constants import UNDEF_LIMIT
+
 from numpy.typing import NDArray
 
+from webviz_services.service_exceptions import InvalidParameterError, Service
+
 
 def surface_to_float32_numpy_array(surface: xtgeo.RegularSurface) -> NDArray[np.float32]:
     masked_values = surface.values.astype(np.float32)
@@ -48,3 +54,120 @@ def get_min_max_surface_values(surface: xtgeo.RegularSurface) -> MinMax | None:
         return None
 
     return MinMax(min=masked_min_val, max=masked_max_val)
+
+
+class PickDirection(StrEnum):
+    """Direction of the pick relative to the surface"""
+
+    UPWARD = "UPWARD"
+    DOWNWARD = "DOWNWARD"
+
+
+@dataclass
+class SurfaceWellPick:
+    """Surface pick data along a well trajectory"""
+
+    unique_wellbore_identifier: str
+    x: float
+    y: float
+    z: float
+    md: float
+    direction: PickDirection
+
+
+@dataclass
+class WellTrajectory:
+    """
+    Well trajectory defined by a set of (x, y, z) coordinates and measured depths (md).
+
+    unique_wellbore_identifier: str
+    x_points: X-coordinates of well trajectory points.
+    y_points: Y-coordinates of well trajectory points.
+    z_points: Z-coordinates (depth values) of well trajectory points.
+    md_points: Measured depth values at each well trajectory point.
+    """
+
+    unique_wellbore_identifier: str
+    x_points: list[float]
+    y_points: list[float]
+    z_points: list[float]
+    md_points: list[float]
+
+
+def get_surface_picks_for_well_trajectory_from_xtgeo(
+    surf: xtgeo.RegularSurface,
+    well_trajectory: WellTrajectory,
+) -> list[SurfaceWellPick] | None:
+    """
+    Calculate intersections (wellpicks) between a surface and a well trajectory.
+
+    Uses the underlying xtgeo C-extension directly for performance.
+
+    Note that this function performs interpolation internally to find the intersections. Thereby
+    it can be inaccurate calculations of picks if the length between a well trajectory point and
+    the next is large and the surface intersects between these points.
+
+    """
+
+    # Ensure equal length arrays
+    if not (
+        well_trajectory.x_points
+        and well_trajectory.y_points
+        and well_trajectory.z_points
+        and well_trajectory.md_points
+        and len(well_trajectory.x_points)
+        == len(well_trajectory.y_points)
+        == len(well_trajectory.z_points)
+        == len(well_trajectory.md_points)
+    ):
+        raise InvalidParameterError(
+            "Well trajectory point arrays must be non-empty and of equal length", Service.GENERAL
+        )
+
+    xarray = np.array(well_trajectory.x_points, dtype=np.float32)
+    yarray = np.array(well_trajectory.y_points, dtype=np.float32)
+    zarray = np.array(well_trajectory.z_points, dtype=np.float32)
+    mdarray = np.array(well_trajectory.md_points, dtype=np.float32)
+
+    # nval = number of valid picks
+    # xres, yres, zres = arrays of x,y,z coordinates of picks
+    # mres = array of measured depth values of picks
+    # dres = array of direction indicators of picks (1=downward, 0=upward)
+    nval, xres, yres, zres, mres, dres = _cxtgeo.well_surf_picks(
+        xarray,
+        yarray,
+        zarray,
+        mdarray,
+        surf.ncol,
+        surf.nrow,
+        surf.xori,
+        surf.yori,
+        surf.xinc,
+        surf.yinc,
+        surf.yflip,
+        surf.rotation,
+        surf.npvalues1d,
+        xarray.size,
+        yarray.size,
+        zarray.size,
+        mdarray.size,
+        mdarray.size,
+    )
+    if nval < 1:
+        return None
+
+    mres[mres > UNDEF_LIMIT] = np.nan
+
+    res: list[SurfaceWellPick] = []
+    for i in range(nval):
+        res.append(
+            SurfaceWellPick(
+                unique_wellbore_identifier=well_trajectory.unique_wellbore_identifier,
+                x=xres[i],
+                y=yres[i],
+                z=zres[i],
+                md=mres[i],
+                direction=PickDirection.DOWNWARD if dres[i] == 1 else PickDirection.UPWARD,
+            )
+        )
+    return res