feature: add multi cone cpt handling (#5)

* add simplify from simplify py

* add url from original implementation

* add support for multi cone

* improve unique cpt identifier

---------

Co-authored-by: Julien Hericher <Julien.Hericher@ngi.no>

Author: Julien76

src/ngi_calculations/cpt_correlations/methods/cpt_process/calculations.py

@@ -8,6 +8,7 @@
 from ngi_calculations.cpt_correlations.definitions.geo import GEO
 from ngi_calculations.cpt_correlations.definitions.physics import PhysicParameters as PHY
 from ngi_calculations.cpt_correlations.methods.cpt_process.options import CptProcessOptions
+from ngi_calculations.cpt_correlations.models.cpt_cone import CptCone
 from ngi_calculations.cpt_correlations.models.cpt_raw import RawCPT
 from ngi_calculations.cpt_correlations.models.lab_data import LabData
 from ngi_calculations.cpt_correlations.utils.interpolation import interpolate_missing_values
@@ -30,7 +31,6 @@ def __init__(self, raw_cpt: RawCPT, lab_data: LabData, options: CptProcessOption
         self.raw_cpt = raw_cpt
         self.lab_data = lab_data
         self.options = options
-        # self.calculate()
 
     @measure("ProcessCPT_Calculation", log_time=False, log_child=False)
     @track_execution
@@ -50,6 +50,7 @@ def calculate(self):
         self._differential_pressure()
         self._normalized_differential_pressure()
         self._elevation()
+        self._set_cone_info()
         self._total_cone_resistance()
         self._net_cone_resistance()
         self._normalized_cone_resistance()
@@ -63,7 +64,7 @@ def calculate(self):
         self._filter_data()
 
     def _set_initial_data(self):
-        self.data = self.raw_cpt.data[self.raw_cpt.columns.all].copy()
+        self.data = self.raw_cpt.data.copy()
         self.data[GEO.u2_raw.key] = self.data[GEO.u2.key]
         self.data[GEO.depth_raw.key] = self.data[GEO.depth.key]
         self.depth_raw = self.data[GEO.depth.key].values
@@ -164,7 +165,14 @@ def _integrate_lab_profile(self):
         _df = _df.drop([x for x in lab_cols if x in _df.columns], axis=1)
 
         # Merge the two Dataframes
-        _df.astype(np.float64, copy=False)
+        # Create a type dictionary for all columns
+        types_dict = {col: np.float64 for col in _df.columns if col != self.options.cpt_identifier and col != depth}
+        # Add depth type if needed
+        types_dict[depth] = np.float64
+        # Apply types
+        _df = _df.astype(types_dict)
+
+        # _df.astype(np.float64, copy=False)
         _df = pd.merge(_df, lab_df, on=depth, how="outer")
         _df.sort_values(by=depth, inplace=True)
         _df.set_index(depth, drop=False, inplace=True)
@@ -231,16 +239,36 @@ def _elevation(self):
         # TODO: Where to make available the self.options.elevation?
         self.data[GEO.elevation.key] = self.options.elevation - self.data[GEO.depth.key]
 
+    def _set_cone_info(self):
+        if self.options.cpt_identifier in self.data.columns:
+            # Create a function to map each identifier to the corresponding cone area ratio value
+            def get_cone_area_ratio(identifier):
+                cone = self.raw_cpt.cone.get(identifier)
+                # Return the numeric cone_area_ratio attribute, not the entire CptCone object
+                return cone.cone_area_ratio if cone else CptCone().cone_area_ratio
+
+            # Apply the function to each row's identifier
+            self.data[GEO.cone_area_ratio.key] = self.data[self.options.cpt_identifier].apply(get_cone_area_ratio)
+        else:
+            # If the identifier column doesn't exist, use the default value
+            self.data[GEO.cone_area_ratio.key] = CptCone().cone_area_ratio
+
+        # Do the same for sleeve_area_ratio
+        if self.options.cpt_identifier in self.data.columns:
+
+            def get_sleeve_area_ratio(identifier):
+                cone = self.raw_cpt.cone.get(identifier)
+                # Return the numeric sleeve_area_ratio attribute, not the entire CptCone object
+                return cone.sleeve_area_ratio if cone else CptCone().sleeve_area_ratio
+
+            self.data[GEO.sleeve_area_ratio.key] = self.data[self.options.cpt_identifier].apply(get_sleeve_area_ratio)
+        else:
+            self.data[GEO.sleeve_area_ratio.key] = CptCone().sleeve_area_ratio
+
     def _total_cone_resistance(self):
-        self.data[GEO.cone_area_ratio.key] = self.raw_cpt.cone.cone_area_ratio
-        self.data[GEO.sleeve_area_ratio.key] = self.raw_cpt.cone.sleeve_area_ratio
         self.data[GEO.qt.key] = (
             1000 * self.data[GEO.qc.key] + self.data[GEO.u2.key] * (1 - self.data[GEO.cone_area_ratio.key])
         ) / 1000
-        # self.df[GEO.qt.key] = (
-        #     1000 * self.df[GEO.qc.key]
-        #     + self.df[GEO.u2.key] * (1 - self.df[GEO.cone_area_ratio.key]) * self.df[GEO.sleeve_area_ratio.key]
-        # ) / 1000
 
     def _net_cone_resistance(self):
         self.data[GEO.qn.key] = self.data[GEO.qt.key] - self.data[GEO.sigVtTotal.key] / 1000

src/ngi_calculations/cpt_correlations/methods/cpt_process/options.py

@@ -9,3 +9,4 @@ class CptProcessOptions(BaseModel):
     adjust_depth_tilt: bool = False  # Adjust depth due to tilt
     compensate_atm_pressure: bool = False
     interpolation_mode: Literal["linear", "padding"] = "linear"
+    cpt_identifier: str = "method_id"

src/ngi_calculations/cpt_correlations/models/cpt_raw.py

@@ -34,7 +34,7 @@ def all(self):
 class RawCPT(CustomBaseModel):
     data: pd.DataFrame = Field(..., description="Pandas DataFrame containing raw CPT data")
     columns: RawCPTColumns = RawCPTColumns()
-    cone: CptCone = CptCone()
+    cone: dict[str, CptCone] | None = None
 
     @field_validator("data")
     def validate_data(cls, value):

src/ngi_calculations/cpt_correlations/utils/simplify.py

@@ -0,0 +1,121 @@
+# Description: inspired from https://github.com/omarestrella/simplify.py
+
+
+def getSquareDistance(p1, p2, xKey, yKey):
+    """
+    Square distance between two points
+    """
+    dx = p1[xKey] - p2[xKey]
+    dy = p1[yKey] - p2[yKey]
+
+    return dx * dx + dy * dy
+
+
+def getSquareSegmentDistance(p, p1, p2, xKey, yKey):
+    """
+    Square distance between point and a segment
+    """
+    x = p1[xKey]
+    y = p1[yKey]
+
+    dx = p2[xKey] - x
+    dy = p2[yKey] - y
+
+    if dx != 0 or dy != 0:
+        t = ((p[xKey] - x) * dx + (p[yKey] - y) * dy) / (dx * dx + dy * dy)
+
+        if t > 1:
+            x = p2[xKey]
+            y = p2[yKey]
+        elif t > 0:
+            x += dx * t
+            y += dy * t
+
+    dx = p[xKey] - x
+    dy = p[yKey] - y
+
+    return dx * dx + dy * dy
+
+
+def simplifyRadialDistance(points, tolerance, xKey, yKey):
+    length = len(points)
+    prev_point = points[0]
+    new_points = [prev_point]
+
+    for i in range(length):
+        point = points[i]
+
+        if getSquareDistance(point, prev_point, xKey, yKey) > tolerance:
+            new_points.append(point)
+            prev_point = point
+
+    if prev_point != point:
+        new_points.append(point)
+
+    return new_points
+
+
+def simplifyDouglasPeucker(points, tolerance, xKey="x", yKey="y"):
+    length = len(points)
+    markers = [0] * length  # Maybe not the most efficent way?
+
+    first = 0
+    last = length - 1
+
+    first_stack = []
+    last_stack = []
+
+    new_points = []
+
+    markers[first] = 1
+    markers[last] = 1
+
+    while last:
+        max_sqdist = 0
+
+        for i in range(first, last):
+            sqdist = getSquareSegmentDistance(points[i], points[first], points[last], xKey, yKey)
+
+            if sqdist > max_sqdist:
+                index = i
+                max_sqdist = sqdist
+
+        if max_sqdist > tolerance:
+            markers[index] = 1
+
+            first_stack.append(first)
+            last_stack.append(index)
+
+            first_stack.append(index)
+            last_stack.append(last)
+
+        # Can pop an empty array in Javascript, but not Python, so check
+        # the length of the list first
+        if len(first_stack) == 0:
+            first = None
+        else:
+            first = first_stack.pop()
+
+        if len(last_stack) == 0:
+            last = None
+        else:
+            last = last_stack.pop()
+
+    for i in range(length):
+        if markers[i]:
+            new_points.append(points[i])
+
+    return new_points
+
+
+def simplify(
+    points: list[dict], tolerance: float = 0.1, highestQuality: bool = True, xKey: str = "x", yKey: str = "y"
+) -> list[dict]:
+    sqtolerance = tolerance * tolerance
+
+    if not highestQuality:
+        points = simplifyRadialDistance(points, sqtolerance, xKey, yKey)
+
+    points = simplifyDouglasPeucker(points, sqtolerance, xKey, yKey)
+
+    return points

tests/conftest.py

@@ -3,12 +3,11 @@
 
 import pandas as pd
 import pytest
-from tests.data.calculation_data import Files2, get_data_from_excel_calculation_file
 
 from ngi_calculations.cpt_correlations.methods.cpt_process.calculations import CPTProcessCalculation
-from ngi_calculations.cpt_correlations.models.cpt_cone import CptCone
 from ngi_calculations.cpt_correlations.models.cpt_raw import RawCPT
 from ngi_calculations.cpt_correlations.models.lab_data import LabData
+from tests.data.calculation_data import Files2, get_data_from_excel_calculation_file
 
 
 def dict_parametrize(data, **kwargs) -> Callable:
@@ -29,7 +28,7 @@ class FullCase:
 @pytest.fixture(params=[Files2.testA], scope="session")
 def analysis_from_excel(request):
     excel_data = get_data_from_excel_calculation_file(request.param)
-    raw_cpt = RawCPT(data=excel_data.raw_cpt, cone=CptCone(cone_area_ratio=excel_data.cone_area_ratio))
+    raw_cpt = RawCPT(data=excel_data.raw_cpt, cone=excel_data.cone)
     lab_data = LabData(data=excel_data.lab_data)
     processor = CPTProcessCalculation(raw_cpt=raw_cpt, lab_data=lab_data)
     return FullCase(raw_cpt=raw_cpt, processor=processor, expected=excel_data.processed_cpt)

tests/data/calculation_data.py

@@ -4,6 +4,7 @@
 
 import pandas as pd
 
+from ngi_calculations.cpt_correlations.models.cpt_cone import CptCone
 from tests.data.excel_data_column_mapping import (
     CptProcessedColumns,
     CptRawColumns,
@@ -83,8 +84,7 @@ class ExcelFileConfig:
     lab_data: ExcelDataConfig
     raw_cpt: ExcelDataConfig
     processed_cpt: ExcelDataConfig
-    # interpretation: ExcelInterpretation
-    cone_area_ratio: float
+    cone: dict[str, CptCone]
 
 
 @dataclass
@@ -104,7 +104,7 @@ class Files2(Enum):
         raw_cpt=ExcelDataConfig(
             sheetname="cpt_raw",
             start_column="A",
-            end_column="H",
+            end_column="I",
             header_row=2,
             data_start_row=4,
             data_end_row=439,
@@ -122,13 +122,17 @@ class Files2(Enum):
         processed_cpt=ExcelDataConfig(
             sheetname="cpt_processed",
             start_column="A",
-            end_column="AC",
+            end_column="AF",
             header_row=2,
             data_start_row=4,
             data_end_row=439,
             column_mapping=CptProcessedColumns().columns,
         ),
-        cone_area_ratio=0.846,
+        cone={
+            "aaa": CptCone(cone_area_ratio=0.85, sleeve_area_ratio=1.0),
+            "bbb": CptCone(cone_area_ratio=0.92, sleeve_area_ratio=0.9),
+            "ccc": CptCone(cone_area_ratio=0.75, sleeve_area_ratio=0.8),
+        },
     )
 
 
@@ -137,8 +141,7 @@ class ExcelData2:
     lab_data: pd.DataFrame
     raw_cpt: pd.DataFrame
     processed_cpt: pd.DataFrame
-    # interpretation: ExcelInterpretation
-    cone_area_ratio: float
+    cone: dict[str, CptCone]
 
 
 CURRENT_DIR = os.getcwd()
@@ -156,5 +159,5 @@ def get_data_from_excel_calculation_file(file: Files2):
         raw_cpt=cpt_raw_data,
         lab_data=lab_data,
         processed_cpt=cpt_processed_data,
-        cone_area_ratio=settings.cone_area_ratio,
+        cone=settings.cone,
     )

tests/data/calculations/testA.xlsx

@@ -7,14 +7,15 @@
 class CptRawColumns:
     columns: dict = field(
         default_factory=lambda: {
-            "A": GEO.depth.key,
-            "B": GEO.qc.key,
-            "C": GEO.fs.key,
-            "D": GEO.u2.key,
-            "E": GEO.temperature.key,
-            "F": GEO.penetration_rate.key,
-            "G": GEO.penetration_force.key,
-            "H": GEO.tilt.key,
+            "A": "method_id",
+            "B": GEO.depth.key,
+            "C": GEO.qc.key,
+            "D": GEO.fs.key,
+            "E": GEO.u2.key,
+            "F": GEO.temperature.key,
+            "G": GEO.penetration_rate.key,
+            "H": GEO.penetration_force.key,
+            "I": GEO.tilt.key,
         }
     )
 
@@ -40,26 +41,28 @@ class CptProcessedColumns:
     columns: dict = field(
         default_factory=lambda: {
             **CptRawColumns().columns,
-            "I": GEO.wc.key,
-            "J": GEO.WP.key,
-            "K": GEO.LL.key,
-            "L": GEO.Ip.key,
-            "M": GEO.St.key,
-            "N": GEO.uw.key,
-            "O": GEO.u0.key,
-            "P": GEO.sigVtTotal.key,
-            "Q": GEO.sigVtEff.key,
-            "R": GEO.qt.key,
-            "S": GEO.qn.key,
-            "T": GEO.u_delta.key,
-            "U": GEO.Fr.key,
-            "V": GEO.Qt.key,
-            "W": GEO.Bq.key,
-            "X": GEO.u_delta_norm.key,
-            "Y": GEO.Ic.key,
-            "Z": GEO.n.key,
-            "AA": GEO.Qtn.key,
-            "AB": GEO.Icn.key,
-            "AC": GEO.Rf.key,
+            "J": GEO.wc.key,
+            "K": GEO.WP.key,
+            "L": GEO.LL.key,
+            "M": GEO.Ip.key,
+            "N": GEO.St.key,
+            "O": GEO.uw.key,
+            "P": GEO.u0.key,
+            "Q": GEO.cone_area_ratio.key,
+            "R": GEO.sleeve_area_ratio.key,
+            "S": GEO.sigVtTotal.key,
+            "T": GEO.sigVtEff.key,
+            "U": GEO.qt.key,
+            "V": GEO.qn.key,
+            "W": GEO.u_delta.key,
+            "X": GEO.Fr.key,
+            "Y": GEO.Qt.key,
+            "Z": GEO.Bq.key,
+            "AA": GEO.u_delta_norm.key,
+            "AB": GEO.Ic.key,
+            "AC": GEO.n.key,
+            "AD": GEO.Qtn.key,
+            "AE": GEO.Icn.key,
+            "AF": GEO.Rf.key,
         }
     )