Add seismic intersection module - seismic fence (#449)

Co-authored-by: Hans Kallekleiv <[email protected]>

Author: jorgenherje

.github/workflows/webviz.yml

@@ -71,7 +71,7 @@ jobs:
             - name: 🕵️ Check auto-generated frontend code is in sync with backend
               run: |
                   docker build -f backend.Dockerfile -t backend:latest .
-                  CONTAINER_ID=$(docker run --detach -p 5000:5000 --env UVICORN_PORT=5000 --env UVICORN_ENTRYPOINT=src.backend.primary.main:app --env WEBVIZ_CLIENT_SECRET=0 --env WEBVIZ_SMDA_SUBSCRIPTION_KEY=0 --env WEBVIZ_SMDA_RESOURCE_SCOPE=0 backend:latest)
+                  CONTAINER_ID=$(docker run --detach -p 5000:5000 --env UVICORN_PORT=5000 --env UVICORN_ENTRYPOINT=src.backend.primary.main:app --env WEBVIZ_CLIENT_SECRET=0 --env WEBVIZ_SMDA_SUBSCRIPTION_KEY=0 --env WEBVIZ_SMDA_RESOURCE_SCOPE=0 --env WEBVIZ_VDS_HOST_ADDRESS=0 backend:latest)
                   sleep 5  # Ensure the backend server is up and running exposing /openapi.json
                   npm run generate-api --prefix ./frontend
                   docker stop $CONTAINER_ID
@@ -113,6 +113,7 @@ jobs:
                   WEBVIZ_CLIENT_SECRET: 0
                   WEBVIZ_SMDA_SUBSCRIPTION_KEY: 0
                   WEBVIZ_SMDA_RESOURCE_SCOPE: 0
+                  WEBVIZ_VDS_HOST_ADDRESS: 0
               run: |
                   pytest ./tests/unit
 

backend/poetry.lock

@@ -26,6 +26,7 @@ vtk = "^9.2.6"
 fmu-sumo = "1.0.3"
 sumo-wrapper-python = "1.0.6"
 azure-monitor-opentelemetry = "^1.1.0"
+requests-toolbelt = "^1.0.0"
 
 
 [tool.poetry.group.dev.dependencies]

backend/pyproject.toml

@@ -1,31 +1,100 @@
 import logging
-from typing import List
+from typing import List, Optional
 
-from fastapi import APIRouter, Depends, HTTPException, Query
+from fastapi import APIRouter, Depends, HTTPException, Query, Body
 
-from src.services.sumo_access.seismic_access import SeismicAccess
+from src.services.sumo_access.seismic_access import SeismicAccess, VdsHandle
+from src.services.vds_access.vds_access import VdsAccess
 from src.services.utils.authenticated_user import AuthenticatedUser
 from src.backend.auth.auth_helper import AuthHelper
+from src.services.utils.b64 import b64_encode_float_array_as_float32
+from src.services.vds_access.response_types import VdsMetadata
+from src.services.vds_access.request_types import VdsCoordinateSystem, VdsCoordinates
 
 from . import schemas
 
+
 LOGGER = logging.getLogger(__name__)
 
 router = APIRouter()
 
 
-@router.get("/seismic_directory/")
-async def get_seismic_directory(
+@router.get("/seismic_cube_meta_list/")
+async def get_seismic_cube_meta_list(
     authenticated_user: AuthenticatedUser = Depends(AuthHelper.get_authenticated_user),
     case_uuid: str = Query(description="Sumo case uuid"),
     ensemble_name: str = Query(description="Ensemble name"),
 ) -> List[schemas.SeismicCubeMeta]:
     """
-    Get a directory of seismic cubes.
+    Get a list of seismic cube meta.
     """
     access = await SeismicAccess.from_case_uuid(authenticated_user.get_sumo_access_token(), case_uuid, ensemble_name)
-    seismic_cube_metas = await access.get_seismic_directory()
+    seismic_cube_meta_list = await access.get_seismic_cube_meta_list_async()
     try:
-        return [schemas.SeismicCubeMeta(**meta.__dict__) for meta in seismic_cube_metas]
+        return [schemas.SeismicCubeMeta(**meta.__dict__) for meta in seismic_cube_meta_list]
     except ValueError as exc:
         raise HTTPException(status_code=400, detail=str(exc)) from exc
+
+
+@router.post("/get_seismic_fence/")
+async def post_get_seismic_fence(
+    authenticated_user: AuthenticatedUser = Depends(AuthHelper.get_authenticated_user),
+    case_uuid: str = Query(description="Sumo case uuid"),
+    ensemble_name: str = Query(description="Ensemble name"),
+    realization_num: int = Query(description="Realization number"),
+    seismic_attribute: str = Query(description="Seismic cube attribute"),
+    time_or_interval_str: str = Query(description="Timestamp or timestep"),
+    observed: bool = Query(description="Observed or simulated"),
+    polyline: schemas.SeismicFencePolyline = Body(embed=True),
+) -> schemas.SeismicFenceData:
+    """Get a fence of seismic data from a polyline defined by a set of (x, y) coordinates in domain coordinate system.
+
+    The fence data contains a set of traces perpendicular to the polyline, with one trace per (x, y)-point in polyline.
+    Each trace has equal number of samples, and is a set of sample values along the depth direction of the seismic cube.
+
+    Returns:
+    A SeismicFenceData object with fence traces in encoded 1D array, metadata for trace array decoding and fence min/max depth.
+    """
+    seismic_access = await SeismicAccess.from_case_uuid(
+        authenticated_user.get_sumo_access_token(), case_uuid, ensemble_name
+    )
+
+    vds_handle: Optional[VdsHandle] = None
+    try:
+        vds_handle = await seismic_access.get_vds_handle_async(
+            realization=realization_num,
+            seismic_attribute=seismic_attribute,
+            time_or_interval_str=time_or_interval_str,
+            observed=observed,
+        )
+    except ValueError as err:
+        raise HTTPException(status_code=404, detail=str(err)) from err
+
+    if vds_handle is None:
+        raise HTTPException(status_code=404, detail="Vds handle not found")
+
+    vds_access = VdsAccess(sas_token=vds_handle.sas_token, vds_url=vds_handle.vds_url)
+
+    # Retrieve fence and post as seismic intersection using cdp coordinates for vds-slice
+    # NOTE: Correct coordinate format and scaling - see VdsCoordinateSystem?
+    [
+        flattened_fence_traces_array,
+        num_traces,
+        num_samples_per_trace,
+    ] = await vds_access.get_flattened_fence_traces_array_and_metadata_async(
+        coordinates=VdsCoordinates(polyline.x_points, polyline.y_points),
+        coordinate_system=VdsCoordinateSystem.CDP,
+    )
+
+    meta: VdsMetadata = await vds_access.get_metadata_async()
+    if len(meta.axis) != 3:
+        raise ValueError(f"Expected 3 axes, got {len(meta.axis)}")
+    depth_axis_meta = meta.axis[2]
+
+    return schemas.SeismicFenceData(
+        fence_traces_b64arr=b64_encode_float_array_as_float32(flattened_fence_traces_array),
+        num_traces=num_traces,
+        num_samples_per_trace=num_samples_per_trace,
+        min_fence_depth=depth_axis_meta.min,
+        max_fence_depth=depth_axis_meta.max,
+    )

backend/src/backend/primary/routers/seismic/router.py

@@ -1,5 +1,9 @@
+from typing import List
+
 from pydantic import BaseModel
 
+from src.services.utils.b64 import B64FloatArray
+
 
 class SeismicCubeMeta(BaseModel):
     seismic_attribute: str
@@ -8,6 +12,50 @@ class SeismicCubeMeta(BaseModel):
     is_depth: bool
 
 
-class VdsHandle(BaseModel):
-    sas_token: str
-    vds_url: str
+class SeismicFencePolyline(BaseModel):
+    """
+    (x, y) points defining a polyline in domain coordinate system, to retrieve fence of seismic data.
+
+    Expect equal number of x- and y-points.
+
+    Note: Coordinates are in domain coordinate system (UTM).
+
+    NOTE:
+    - Verify coordinates are in domain coordinate system (UTM)?
+    - Consider points_xy: List[float] - i.e. array with [x1, y1, x2, y2, ..., xn, yn] instead of x_points and y_points arrays?
+    - Ensure equal length of x_points and y_points arrays?
+    """
+
+    x_points: List[float]
+    y_points: List[float]
+
+
+class SeismicFenceData(BaseModel):
+    """
+    Definition of a fence of seismic data from a set of (x, y) coordinates in domain coordinate system.
+    Each (x, y) point provides a trace perpendicular to the x-y plane, with number of samples equal to the depth of the seismic cube.
+
+    Each trace is defined to be a set of depth value samples along the length direction of the fence.
+
+    `Properties:`
+    - `fence_traces_b64arr`: The fence trace array is base64 encoded 1D float array - where data is stored trace by trace.
+    - `num_traces`: The number of traces in the fence trace array. Equals the number of (x, y) coordinates in requested polyline.
+    - `num_samples_per_trace`: The number of samples in each trace.
+    - `min_fence_depth`: The minimum depth value of the fence.
+    - `max_fence_depth`: The maximum depth value of the fence.
+
+    `Description - fence_traces_b64arr:`\n
+    The encoded fence trace array is a flattened array of traces, where data is stored trace by trace.
+    With `m = num_traces`, and `n = num_samples_per_trace`, the flattened array has length `mxn`.
+
+    Fence traces 1D array: [trace_1_sample_1, trace_1_sample_2, ..., trace_1_sample_n, ..., trace_m_sample_1, trace_m_sample_2, ..., trace_m_sample_n] \n
+
+    See:
+    - VdsAxis: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/internal/core/core.go#L37-L55
+    """
+
+    fence_traces_b64arr: B64FloatArray
+    num_traces: int
+    num_samples_per_trace: int
+    min_fence_depth: float
+    max_fence_depth: float

backend/src/backend/primary/routers/seismic/schemas.py

@@ -17,6 +17,7 @@
 SMDA_RESOURCE_SCOPE = os.environ["WEBVIZ_SMDA_RESOURCE_SCOPE"]
 SUMO_ENV = os.getenv("WEBVIZ_SUMO_ENV", "prod")
 GRAPH_SCOPES = ["User.Read", "User.ReadBasic.All"]
+VDS_HOST_ADDRESS = os.environ["WEBVIZ_VDS_HOST_ADDRESS"]
 
 RESOURCE_SCOPES_DICT = {
     "sumo": [f"api://{sumo_app_reg[SUMO_ENV]['RESOURCE_ID']}/access_as_user"],

backend/src/config.py

@@ -12,9 +12,9 @@
 
 
 class SeismicAccess(SumoEnsemble):
-    async def get_seismic_directory(self) -> List[SeismicCubeMeta]:
+    async def get_seismic_cube_meta_list_async(self) -> List[SeismicCubeMeta]:
         seismic_cube_collection: CubeCollection = self._case.cubes.filter(iteration=self._iteration_name, realization=0)
-        seismic_cube_metas: List[SeismicCubeMeta] = []
+        seismic_cube_meta_list: List[SeismicCubeMeta] = []
         async for cube in seismic_cube_collection:
             t_start = cube["data"].get("time", {}).get("t0", {}).get("value", None)
             t_end = cube["data"].get("time", {}).get("t1", {}).get("value", None)
@@ -34,10 +34,10 @@ async def get_seismic_directory(self) -> List[SeismicCubeMeta]:
                 is_observation=cube["data"]["is_observation"],
                 is_depth=cube["data"]["vertical_domain"] == "depth",
             )
-            seismic_cube_metas.append(seismic_meta)
-        return seismic_cube_metas
+            seismic_cube_meta_list.append(seismic_meta)
+        return seismic_cube_meta_list
 
-    async def get_vds_handle(
+    async def get_vds_handle_async(
         self,
         seismic_attribute: str,
         realization: int,

backend/src/services/sumo_access/seismic_access.py

@@ -0,0 +1,121 @@
+from dataclasses import dataclass
+from enum import StrEnum
+from typing import List
+
+######################################################################################################
+#
+# This file contains the request types for the vds-slice service found in the following file:
+#
+# https://github.com/equinor/vds-slice/blob/master/api/request.go
+#
+# Master commit hash: ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3
+#
+# https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go
+#
+######################################################################################################
+
+
+class VdsInterpolation(StrEnum):
+    """
+    Interpolation options for vds fence
+
+    Source: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go#L98
+    """
+
+    NEAREST = "nearest"
+    LINEAR = "linear"
+    CUBIC = "cubic"
+    ANGULAR = "angular"
+    TRIANGULAR = "triangular"
+
+
+class VdsCoordinateSystem(StrEnum):
+    """
+    Coordinate system options for vds fence
+
+    * ilxl: inline, crossline pairs
+    * ij: Coordinates are given as in 0-indexed system, where the first
+          line in each direction is 0 and the last is number-of-lines - 1.
+    * cdp: Coordinates are given as cdpx/cdpy pairs. In the original SEGY
+           this would correspond to the cdpx and cdpy fields in the
+           trace-headers after applying the scaling factor.
+
+    Source: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go#L86C3-L86C3
+    """
+
+    CDP = "cdp"
+    IJ = "ij"
+    ILXL = "ilxl"
+
+
+@dataclass
+class VdsCoordinates:
+    """
+    A list of coordinates in the selected VdsCoordinateSystem, as (x, y) points.
+
+    Convert coordinates to format for query request parameter - [[x1,y1], [x2,y2], ..., [xn,yn]]
+
+    Source: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go#L90
+    """
+
+    x_points: List[float]
+    y_points: List[float]
+
+    def __init__(self, x_points: List[float], y_points: List[float]) -> None:
+        if len(x_points) != len(y_points):
+            raise ValueError("x_points and y_points must be of equal length")
+
+        self.x_points = x_points
+        self.y_points = y_points
+
+    def to_list(self) -> List[List[float]]:
+        return [[x, y] for x, y in zip(self.x_points, self.y_points)]
+
+
+@dataclass
+class VdsRequestedResource:
+    """
+    Definition of requested vds resource for vds-slice
+    This is a base class for request types for vds-slice requests
+
+    See: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go#L13-L35
+    """
+
+    vds: str  # blob url
+    sas: str  # sas-token
+
+    def request_parameters(self) -> dict:
+        return {"vds": self.vds, "sas": self.sas}
+
+
+@dataclass
+class VdsMetadataRequest(VdsRequestedResource):
+    """
+    Definition of metadata request for vds-slice
+
+    See: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go#L62-L64
+    """
+
+
+@dataclass
+class VdsFenceRequest(VdsRequestedResource):
+    """
+    Definition of a fence request struct for vds-slice
+
+    See: https://github.com/equinor/vds-slice/blob/ab6f39789bf3d3b59a8df14f1c4682d340dc0bf3/api/request.go#L76-L105
+    """
+
+    coordinate_system: VdsCoordinateSystem
+    coordinates: VdsCoordinates
+    interpolation: VdsInterpolation
+    fill_value: float
+
+    def request_parameters(self) -> dict:
+        return {
+            "vds": self.vds,
+            "sas": self.sas,
+            "coordinateSystem": self.coordinate_system.value,
+            "coordinates": self.coordinates.to_list(),
+            "interpolation": self.interpolation.value,
+            "fillValue": self.fill_value,
+        }