feat: derive basal contacts from geology map in thickness calculators

Basal contacts are now extracted automatically from the geology layer
and stratigraphic column when no basal-contacts layer is supplied, so
the separate "Extract Basal Contacts" step is no longer required
before running the thickness calculator (GUI widget and the
Processing Toolbox algorithm).

Author: lachlangrose

loopstructural/gui/map2loop_tools/thickness_calculator_widget.py

@@ -338,15 +338,30 @@ def _run_calculator(self):
         # Validate inputs based on calculator type
         calculator_type = self.calculatorTypeComboBox.currentText()
 
-        if calculator_type == "InterpolatedStructure":
-            if not self.geologyLayerComboBox.currentLayer():
-                QMessageBox.warning(self, "Missing Input", "Please select a geology layer.")
-                return False
-            if not self.basalContactsComboBox.currentLayer():
-                QMessageBox.warning(self, "Missing Input", "Please select a basal contacts layer.")
+        if not self.geologyLayerComboBox.currentLayer():
+            QMessageBox.warning(self, "Missing Input", "Please select a geology layer.")
+            return False
+
+        # A basal contacts layer is optional: when omitted, basal contacts are
+        # extracted automatically from the geology layer using the
+        # stratigraphic column, so a stratigraphic order is required instead.
+        if not self.basalContactsComboBox.currentLayer():
+            stratigraphic_order = (
+                self.data_manager.get_stratigraphic_unit_names()
+                if self.data_manager and hasattr(self.data_manager, 'get_stratigraphic_unit_names')
+                else None
+            )
+            if not stratigraphic_order:
+                QMessageBox.warning(
+                    self,
+                    "Missing Input",
+                    "Please select a basal contacts layer, or define a stratigraphic "
+                    "column so basal contacts can be extracted from the geology layer "
+                    "automatically.",
+                )
                 return False
 
-        elif calculator_type == "StructuralPoint":
+        if calculator_type == "StructuralPoint":
             if not self.structureLayerComboBox.currentLayer():
                 QMessageBox.warning(self, "Missing Input", "Please select a structure layer.")
                 return False

loopstructural/gui/map2loop_tools/thickness_calculator_widget.ui

@@ -63,12 +63,16 @@
      <item row="4" column="0">
       <widget class="QLabel" name="basalContactsLabel">
        <property name="text">
-        <string>Basal Contacts Layer:</string>
+        <string>Basal Contacts Layer (optional):</string>
        </property>
       </widget>
      </item>
      <item row="4" column="1">
-      <widget class="QgsMapLayerComboBox" name="basalContactsComboBox"/>
+      <widget class="QgsMapLayerComboBox" name="basalContactsComboBox">
+       <property name="allowEmptyLayer">
+        <bool>true</bool>
+       </property>
+      </widget>
      </item>
      <item row="6" column="0">
       <widget class="QLabel" name="sampledContactsLabel">

loopstructural/main/m2l_api.py

@@ -429,9 +429,9 @@ def sample_contacts(
 
 def calculate_thickness(
     geology,
-    basal_contacts,
     sampled_contacts,
     structure,
+    basal_contacts=None,
     cross_sections=None,
     calculator_type="InterpolatedStructure",
     dtm=None,
@@ -451,12 +451,15 @@ def calculate_thickness(
     ----------
     geology : QgsVectorLayer or GeoDataFrame
         Geology polygon layer.
-    basal_contacts : QgsVectorLayer or GeoDataFrame
-        Basal contacts line layer.
     sampled_contacts : QgsVectorLayer or GeoDataFrame
         Sampled contacts point layer.
     structure : QgsVectorLayer or GeoDataFrame
         Structure point layer with orientation data.
+    basal_contacts : QgsVectorLayer or GeoDataFrame, optional
+        Basal contacts line layer, by default None. When omitted, basal
+        contacts are extracted automatically from ``geology`` using
+        ``stratigraphic_order``, so a separate basal-contacts extraction
+        step is not required.
     cross_sections : QgsVectorLayer or GeoDataFrame, optional
         Cross-sections line layer, by default None.
     calculator_type : str, optional
@@ -474,7 +477,9 @@ def calculate_thickness(
     max_line_length : float, optional
         Maximum line length for StructuralPoint calculator, by default None.
     stratigraphic_order : list, optional
-        List of unit names in stratigraphic order, by default None.
+        List of unit names in stratigraphic order, by default None. Required
+        when ``basal_contacts`` is omitted, so basal contacts can be derived
+        from the geology layer.
     updater : callable, optional
         Callback function for progress updates, by default None.
 
@@ -486,14 +491,35 @@ def calculate_thickness(
     if updater:
         updater(f"Calculating thickness using {calculator_type}...")
 
+    auto_extracted_basal_contacts = basal_contacts is None
+    if auto_extracted_basal_contacts:
+        if not stratigraphic_order:
+            raise ValueError(
+                "A basal contacts layer or a stratigraphic order is required to calculate "
+                "thickness: without a basal contacts layer, the stratigraphic order is "
+                "needed to extract basal contacts from the geology layer."
+            )
+        if updater:
+            updater("No basal contacts layer supplied; extracting basal contacts from geology...")
+        target_crs = geology.crs() if hasattr(geology, 'crs') else None
+        basal_contacts_gdf = extract_basal_contacts(
+            geology=geology,
+            stratigraphic_order=stratigraphic_order,
+            unit_name_field=unit_name_field,
+            updater=updater,
+            debug_manager=debug_manager,
+            target_crs=target_crs,
+        )['basal_contacts']
+
     # Convert layers to GeoDataFrames
     geology_gdf = qgsLayerToGeoDataFrame(geology)
-    basal_contacts_gdf = qgsLayerToGeoDataFrame(basal_contacts)
-    basal_contacts_gdf = (
-        basal_contacts_gdf.rename(columns={basal_contacts_unit_name: 'basal_unit'})
-        if basal_contacts_unit_name
-        else basal_contacts_gdf
-    )
+    if not auto_extracted_basal_contacts:
+        basal_contacts_gdf = qgsLayerToGeoDataFrame(basal_contacts)
+        basal_contacts_gdf = (
+            basal_contacts_gdf.rename(columns={basal_contacts_unit_name: 'basal_unit'})
+            if basal_contacts_unit_name
+            else basal_contacts_gdf
+        )
     sampled_contacts_gdf = qgsLayerToGeoDataFrame(sampled_contacts)
     structure_gdf = qgsLayerToGeoDataFrame(structure)
     cross_sections_gdf = qgsLayerToGeoDataFrame(cross_sections)

loopstructural/processing/algorithms/thickness_calculator.py

@@ -33,6 +33,7 @@
 )
 
 # Internal imports
+from ...main.m2l_api import extract_basal_contacts
 from ...main.vectorLayerWrapper import (
     dataframeToQgsTable,
     matrixToDict,
@@ -132,9 +133,9 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None:
         self.addParameter(
             QgsProcessingParameterFeatureSource(
                 self.INPUT_BASAL_CONTACTS,
-                "Basal Contacts",
+                "Basal Contacts (optional, extracted from geology if omitted)",
                 [QgsProcessing.TypeVectorLine],
-                defaultValue='Basal Contacts',
+                optional=True,
             )
         )
         self.addParameter(
@@ -261,9 +262,9 @@ def processAlgorithm(
 
         sampled_contacts = self.parameterAsSource(parameters, self.INPUT_SAMPLED_CONTACTS, context)
         unit_name_field = self.parameterAsString(parameters, self.INPUT_UNIT_NAME_FIELD, context)
+        geology_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context)
 
         if bounding_box_type == 0:
-            geology_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context)
             extent = geology_layer.extent()
             bounding_box = {
                 'minx': extent.xMinimum(),
@@ -315,7 +316,19 @@ def processAlgorithm(
         # convert layers to dataframe or geodataframe
         units = qgsLayerToDataFrame(geology_data)
         geology_data = qgsLayerToGeoDataFrame(geology_data)
-        basal_contacts = qgsLayerToGeoDataFrame(basal_contacts)
+        if basal_contacts is not None:
+            basal_contacts = qgsLayerToGeoDataFrame(basal_contacts)
+        else:
+            feedback.pushInfo(
+                "No basal contacts layer supplied; extracting basal contacts from geology layer..."
+            )
+            basal_contacts = extract_basal_contacts(
+                geology=geology_layer,
+                stratigraphic_order=stratigraphic_order,
+                unit_name_field=unit_name_field,
+                updater=feedback.pushInfo,
+                target_crs=geology_layer.crs() if geology_layer else None,
+            )['basal_contacts']
         structure_data = qgsLayerToDataFrame(structure_data)
         rename_map = {}
         missing_fields = []