bugfix : crs for vertical wellbore

Author: valentin-gauthier-geosiris

energyml-utils/example/attic/arrays_test.py

@@ -378,6 +378,40 @@ def read_pointset(
     return meshes
 
 
+def read_wellbore_frame_repr_demo_jfr_02_26(
+    epc_path: str = r"rc/epc/out-galaxy-12-pts.epc",
+    well_uuid: str = "cfad9cb6-99fe-4172-b560-d2feca75dd9f",
+) -> List[AbstractMesh]:
+    epc = Epc.read_file(f"{epc_path}", read_rels_from_files=False, recompute_rels=False)
+
+    frame_repr = epc.get_object_by_uuid(well_uuid)[0]
+    # print(frame_repr)
+    # print(epc.get_h5_file_paths(frame_repr))
+
+    print(epc.get_h5_file_paths())
+
+    print(epc.get_h5_file_paths(frame_repr))
+
+    print("Object type: ", type(frame_repr))
+
+    meshes = read_mesh_object(energyml_object=frame_repr, workspace=epc)
+
+    # Previous result :
+    # points:
+    #     [[   0.    0.    0.]
+    #     [   0.    0.  250.]
+    #     [   0.    0.  500.]
+    #     [   0.    0.  750.]
+    #     [   0.    0. 1000.]]
+    # line indices:
+    #     [[0 1]
+    #     [1 2]
+    #     [2 3]
+    #     [3 4]]
+
+    return meshes
+
+
 if __name__ == "__main__":
     logging.basicConfig(level=logging.DEBUG)
 
@@ -386,7 +420,8 @@ def read_pointset(
     # meshes = read_wellbore_frame_repr()
     # meshes = read_representation_set_representation()
     # meshes = read_trset()
-    meshes = read_pointset()
+    # meshes = read_pointset()
+    meshes = read_wellbore_frame_repr_demo_jfr_02_26()
 
     print(f"Number of meshes read: {len(meshes)}")
 

energyml-utils/src/energyml/utils/data/helper.py

@@ -311,10 +311,12 @@ def get_crs_origin_offset(crs_obj: Any) -> List[float | int]:
     return crs_point_offset
 
 
-def get_datum_information(datum_obj: Any, workspace: Optional[EnergymlStorageInterface] = None):
-    "From a ObjMdDatum or a ReferencePointInACrs, return x, y, z, z_increas_downward, projected_epsg_code, vertical_epsg_code"
+def get_datum_information(
+    datum_obj: Any, workspace: Optional[EnergymlStorageInterface] = None
+) -> Tuple[float, float, float, bool, Optional[str], Optional[str], Optional[Any]]:
+    "From a ObjMdDatum or a ReferencePointInACrs, return x, y, z, z_increas_downward, projected_epsg_code, vertical_epsg_code, crs object"
     if datum_obj is None:
-        return 0.0, 0.0, 0.0, False, None, None
+        return 0.0, 0.0, 0.0, False, None, None, None
 
     t_lw = type(datum_obj).__name__.lower()
 
@@ -333,12 +335,18 @@ def get_datum_information(datum_obj: Any, workspace: Optional[EnergymlStorageInt
             z_increasing_downward,
             projected_epsg_code,
             vertical_epsg_code,
+            datum_obj,
         )
     elif "referencepointinacrs" in t_lw:
         x = get_object_attribute_rgx(datum_obj, "horizontal_coordinates.coordinate1")
         y = get_object_attribute_rgx(datum_obj, "horizontal_coordinates.coordinate2")
         z = get_object_attribute_rgx(datum_obj, "vertical_coordinate")
-        z_increasing_downward = get_object_attribute(datum_obj, "ZIncreasingDownward") or False
+        z_increasing_downward = False
+        v_crs_dor = get_object_attribute_rgx(datum_obj, "vertical_crs")
+        if v_crs_dor is not None and workspace is not None:
+            v_crs = workspace.get_object(get_obj_uri(v_crs_dor))
+            if v_crs is not None:
+                z_increasing_downward = is_z_reversed(v_crs)
         p_crs = get_object_attribute(datum_obj, "horizontal_coordinates.crs")
         projected_epsg_code = (
             get_projected_epsg_code(workspace.get_object(get_obj_uri(p_crs)), workspace)
@@ -354,22 +362,26 @@ def get_datum_information(datum_obj: Any, workspace: Optional[EnergymlStorageInt
             z_increasing_downward,
             projected_epsg_code,
             vertical_epsg_code,
+            p_crs,
         )
     elif "mddatum" in t_lw:
         x = get_object_attribute_rgx(datum_obj, "location.coordinate1")
         y = get_object_attribute_rgx(datum_obj, "location.coordinate2")
         z = get_object_attribute_rgx(datum_obj, "location.coordinate3")
         crs = get_object_attribute(datum_obj, "LocalCrs")
-        _, _, _, z_increasing_downward, projected_epsg_code, vertical_epsg_code = get_datum_information(crs, workspace)
+        _, _, _, z_increasing_downward, projected_epsg_code, vertical_epsg_code, _ = get_datum_information(
+            crs, workspace
+        )
         return (
             float(x) if x is not None else 0.0,
             float(y) if y is not None else 0.0,
             float(z) if z is not None else 0.0,
             z_increasing_downward,
             projected_epsg_code,
             vertical_epsg_code,
+            crs,
         )
-    return 0.0, 0.0, 0.0, False, None, None
+    return 0.0, 0.0, 0.0, False, None, None, None
 
 
 # ==================================================
@@ -659,7 +671,9 @@ def generate_smooth_trajectory(
     return np.array(smooth_points)
 
 
-def generate_vertical_well_points(wellbore_mds: np.ndarray, head_x: float, head_y: float, head_z: float) -> np.ndarray:
+def generate_vertical_well_points(
+    wellbore_mds: np.ndarray, head_x: float, head_y: float, head_z: float, z_increasing_downward: bool = False
+) -> np.ndarray:
     """
     Generates local 3D coordinates for a perfectly vertical wellbore.
 
@@ -684,8 +698,12 @@ def generate_vertical_well_points(wellbore_mds: np.ndarray, head_x: float, head_
     # In a vertical well, Z_point = Z_datum + (MD_point - MD_datum_at_0)
     # Most of the time, MD at head is 0.
     # If wellbore_mds start at 0, Z starts at head_z.
+    # if z_increasing_downward is False, we add the MD to head_z, otherwise we subtract it.
     md_start = wellbore_mds[0]
-    local_points[:, 2] = head_z + (wellbore_mds - md_start)
+    if z_increasing_downward:
+        local_points[:, 2] = head_z - (wellbore_mds - md_start)
+    else:
+        local_points[:, 2] = head_z + (wellbore_mds - md_start)
 
     return local_points
 

energyml-utils/src/energyml/utils/data/mesh.py

@@ -6,6 +6,7 @@
 import os
 import re
 import sys
+import traceback
 import numpy as np
 from dataclasses import dataclass, field
 from enum import Enum
@@ -780,6 +781,8 @@ def read_wellbore_frame_representation(
         trajectory_dor = search_attribute_matching_name(obj=energyml_object, name_rgx="Trajectory")[0]
         trajectory_obj = workspace.get_object(get_obj_uri(trajectory_dor))
 
+        print(f"Mds {wellbore_frame_mds}")
+
         meshes = read_wellbore_trajectory_representation(
             energyml_object=trajectory_obj,
             workspace=workspace,
@@ -831,8 +834,8 @@ def read_wellbore_trajectory_representation(
     # Get CRS from trajectory geometry if available
     try:
         crs = workspace.get_object(get_obj_uri(get_object_attribute(energyml_object, "geometry.LocalCrs")))
-    except Exception as e:
-        logging.debug(f"Could not get CRS from trajectory geometry")
+    except Exception:
+        logging.debug("Could not get CRS from trajectory geometry")
 
     # ==========
     # MD Datum
@@ -853,14 +856,24 @@ def read_wellbore_trajectory_representation(
             md_datum_obj = workspace.get_object(md_datum_identifier)
 
             if md_datum_obj is not None:
-                head_x, head_y, head_z, z_increasing_downward, projected_epsg_code, vertical_epsg_code = (
+                head_x, head_y, head_z, z_increasing_downward, projected_epsg_code, vertical_epsg_code, crs = (
                     get_datum_information(md_datum_obj, workspace)
                 )
+                # if crs is None:
+                #     crs = get_crs_obj(
+                #         context_obj=md_datum_obj,
+                #         path_in_root=".",
+                #         root_obj=energyml_object,
+                #         workspace=workspace,
+                #     )
     except Exception as e:
         logging.debug(f"Could not get reference point / Datum from trajectory: {e}")
 
     # ==========
     well_points = None
+    logging.debug(
+        f"wellbore mds : {wellbore_frame_mds}\n\tCRs : {crs}\n\thead x,y,z : {head_x}, {head_y}, {head_z}\n\tz increasing downward : {z_increasing_downward}"
+    )
     try:
         x_offset, y_offset, z_offset, (azimuth, azimuth_uom) = get_crs_offsets_and_angle(crs, workspace)
         # Try to read parametric Geometry from the trajectory.
@@ -886,11 +899,13 @@ def read_wellbore_trajectory_representation(
                 head_y=head_y,
                 head_z=head_z,
                 wellbore_mds=wellbore_frame_mds,
+                z_increasing_downward=z_increasing_downward,
             )
         else:
+            traceback.print_exc()
             raise ValueError(
                 "Cannot read wellbore trajectory representation: no parametric geometry and no measured depth information available to generate points"
-            ) from e
+            )
 
     meshes = []
     if well_points is not None and len(well_points) > 0:

energyml-utils/src/energyml/utils/epc.py

@@ -1178,11 +1178,10 @@ def write_array(
             obj = self.get_object_by_identifier(proxy)
 
         # Determine which external files to use
-        file_paths = (
-            [make_path_relative_to_other_file(external_uri, self.epc_file_path)]
-            if external_uri
-            else self.get_h5_file_paths(obj)
-        )
+        file_paths = self.get_h5_file_paths(obj)
+        if external_uri:
+            file_paths.insert(0, make_path_relative_to_other_file(external_uri, self.epc_file_path))
+
         if not file_paths or len(file_paths) == 0:
             file_paths = self.external_files_path
 
@@ -1284,39 +1283,37 @@ def get_array_metadata(
 
         return None
 
-    def get_h5_file_paths(self, obj: Any) -> List[str]:
+    def get_h5_file_paths(self, obj_or_id: Optional[Any] = None) -> List[str]:
         """
         Get all HDF5 file paths referenced in the EPC file (from rels to external resources)
         :return: list of HDF5 file paths
         """
 
         if self.force_h5_path is not None:
             return [self.force_h5_path]
+        h5_paths = set()
 
-        is_uri = (isinstance(obj, str) and parse_uri(obj) is not None) or isinstance(obj, Uri)
-        if is_uri:
-            obj = self.get_object_by_identifier(obj)
+        if obj_or_id is None:
+            return [self.epc_file_path.replace(".epc", ".h5")] if self.epc_file_path else []
 
-        h5_paths = set()
+        obj = self.get_object(obj_or_id) if isinstance(obj_or_id, (str, Uri)) else obj_or_id
 
-        if isinstance(obj, str):
-            obj = self.get_object_by_identifier(obj)
         # for rels in self.additional_rels.get(get_obj_identifier(obj), []):
         for rels in self._rels_cache.get_supplemental_rels(obj):
             if rels.type_value == EPCRelsRelationshipType.EXTERNAL_RESOURCE.get_type():
                 h5_paths.add(rels.target)
 
+        h5_paths = set(make_path_relative_to_filepath_list(list(h5_paths), self.epc_file_path))
+
         if len(h5_paths) == 0:
-            # search if an h5 file has the same name than the epc file
+            # Collect all .h5 files in the EPC file's folder
             epc_folder = self.get_epc_file_folder()
-            if epc_folder is not None and self.epc_file_path is not None:
-                epc_file_name = os.path.basename(self.epc_file_path)
-                epc_file_base, _ = os.path.splitext(epc_file_name)
-                possible_h5_path = os.path.join(epc_folder, epc_file_base + ".h5")
-                if os.path.exists(possible_h5_path):
-                    h5_paths.add(possible_h5_path)
-
-        return make_path_relative_to_filepath_list(list(h5_paths), self.epc_file_path)
+            if epc_folder is not None and os.path.isdir(epc_folder):
+                for fname in os.listdir(epc_folder):
+                    if fname.lower().endswith(".h5"):
+                        h5_paths.add(os.path.join(epc_folder, fname))
+
+        return list(h5_paths)
 
     def get_object_as_dor(self, identifier: str, dor_qualified_type) -> Optional[Any]:
         """

energyml-utils/src/energyml/utils/epc_stream.py

@@ -49,6 +49,7 @@
     in_epc_file_path_to_mime_type,
     is_core_prop_or_extension_path,
     repair_epc_structure_if_not_valid,
+    get_file_folder,
 )
 from energyml.utils.storage_interface import (
     DataArrayMetadata,
@@ -1391,20 +1392,18 @@ def get_h5_file_paths(self, obj: Union[str, Uri, Any], make_path_absolute_from_e
                 except KeyError:
                     pass
 
-        if len(h5_paths) == 0:
-            # search if an h5 file has the same name than the epc file
-            epc_folder = os.path.dirname(self.epc_file_path)
-            if epc_folder is not None and self.epc_file_path is not None:
-                epc_file_name = os.path.basename(self.epc_file_path)
-                epc_file_base, _ = os.path.splitext(epc_file_name)
-                possible_h5_path = os.path.join(epc_folder, epc_file_base + ".h5")
-                if os.path.exists(possible_h5_path):
-                    h5_paths.add(possible_h5_path)
-
         if make_path_absolute_from_epc_path:
-            return make_path_relative_to_filepath_list(list(h5_paths), self.epc_file_path)
-        else:
-            return list(h5_paths)
+            h5_paths = set(make_path_relative_to_filepath_list(list(h5_paths), self.epc_file_path))
+
+        if len(h5_paths) == 0:
+            # Collect all .h5 files in the EPC file's folder
+            epc_folder = get_file_folder(self.epc_file_path)
+            if epc_folder is not None and os.path.isdir(epc_folder):
+                for fname in os.listdir(epc_folder):
+                    if fname.lower().endswith(".h5"):
+                        h5_paths.add(os.path.join(epc_folder, fname))
+
+        return list(h5_paths)
 
     #    ________    ___   __________    __  _________________  ______  ____  _____
     #   / ____/ /   /   | / ___/ ___/   /  |/  / ____/_  __/ / / / __ \/ __ \/ ___/