From 19c49ba0a05448566b159c3a86e5d456eae6c697 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Fri, 30 Oct 2026 19:21:24 +1030 Subject: [PATCH 1/4] feat: derive basal contacts from geology map in thickness algorithm The Processing thickness calculator now extracts basal contacts from the geology layer and the stratigraphic column when no basal contacts layer is given. The basal contacts input is optional. --- .../algorithms/thickness_calculator.py | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/loopstructural/processing/algorithms/thickness_calculator.py b/loopstructural/processing/algorithms/thickness_calculator.py index 299a6c5..433747f 100644 --- a/loopstructural/processing/algorithms/thickness_calculator.py +++ b/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.SourceType.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 = [] From 7ed4ab6ebaa03232d3887676affe5cd4fc0904fe Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Fri, 30 Oct 2026 19:24:40 +1030 Subject: [PATCH 2/4] fix: show the error when basal contact extraction fails Before, an error in map2loop gave an empty result and the user did not see the cause. Now extract_basal_contacts raises a RuntimeError with the cause. The GUI and the Processing algorithms show this message. --- loopstructural/main/m2l_api.py | 12 +++++++++--- tests/qgis/test_basal_contacts.py | 10 ++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/loopstructural/main/m2l_api.py b/loopstructural/main/m2l_api.py index 6bb22a4..734d089 100644 --- a/loopstructural/main/m2l_api.py +++ b/loopstructural/main/m2l_api.py @@ -295,6 +295,11 @@ def extract_basal_contacts( ------- dict Dictionary containing 'basal_contacts' GeoDataFrame and optionally 'all_contacts' GeoDataFrame. + + Raises + ------ + RuntimeError + If map2loop cannot extract the contacts. """ geology = qgsLayerToGeoDataFrame(geology, target_crs=target_crs) if unit_name_field and unit_name_field in geology.columns: @@ -400,10 +405,11 @@ def extract_basal_contacts( all_contacts_result.shape, basal_contacts.shape, ) - except Exception: + except Exception as err: + # Do not hide the error with an empty result. The caller shows the + # message to the user. logger.exception("Error during contact extraction") - basal_contacts = pd.DataFrame() - all_contacts_result = pd.DataFrame() + raise RuntimeError(f"Could not extract the basal contacts: {err}") from err if ignore_units and basal_contacts.empty is False: basal_contacts = basal_contacts[ diff --git a/tests/qgis/test_basal_contacts.py b/tests/qgis/test_basal_contacts.py index 385a8de..bbf2950 100644 --- a/tests/qgis/test_basal_contacts.py +++ b/tests/qgis/test_basal_contacts.py @@ -1,5 +1,6 @@ import unittest from pathlib import Path +from unittest.mock import patch from qgis.core import ( QgsFeature, QgsGeometry, @@ -212,6 +213,15 @@ def test_intrusion_in_column(self): self.assertAlmostEqual(self._length(contacts, 'D'), 300, delta=10) self.assertAlmostEqual(self._length(contacts, 'A'), 100, delta=5) + def test_extraction_error_is_raised(self): + # An error in map2loop must not give an empty result. + with patch( + 'loopstructural.main.m2l_api.ContactExtractor.extract_all_contacts', + side_effect=ValueError('bad geology'), + ): + with self.assertRaisesRegex(RuntimeError, 'bad geology'): + self._extract(['A', 'B', 'C']) + if __name__ == '__main__': unittest.main() From 5cd2a60dde00e3ea4bebab988430d3d933ae9729 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Fri, 30 Oct 2026 19:24:40 +1030 Subject: [PATCH 3/4] fix: use one geology source in the thickness algorithm The algorithm now uses the geology feature source for the units, the bounding box and the basal contact extraction. Thus "selected features only" applies to all of them, and an input that is not a map layer works. Add an optional faults input. The faults are used when the algorithm extracts the basal contacts from the geology. Show an error if the extraction finds no basal contacts. Add tests for basal contacts from the geology layer and from the selected features only. --- .../algorithms/thickness_calculator.py | 59 +++++-- tests/qgis/test_thickness_calculator.py | 148 ++++++++++++++++++ 2 files changed, 191 insertions(+), 16 deletions(-) create mode 100644 tests/qgis/test_thickness_calculator.py diff --git a/loopstructural/processing/algorithms/thickness_calculator.py b/loopstructural/processing/algorithms/thickness_calculator.py index 433747f..fe1cc56 100644 --- a/loopstructural/processing/algorithms/thickness_calculator.py +++ b/loopstructural/processing/algorithms/thickness_calculator.py @@ -53,6 +53,7 @@ class ThicknessCalculatorAlgorithm(QgsProcessingAlgorithm): INPUT_MAX_LINE_LENGTH = 'MAX_LINE_LENGTH' INPUT_STRATI_COLUMN = 'STRATIGRAPHIC_COLUMN' INPUT_BASAL_CONTACTS = 'BASAL_CONTACTS' + INPUT_FAULTS = 'FAULTS' INPUT_STRUCTURE_DATA = 'STRUCTURE_DATA' INPUT_DIPDIR_FIELD = 'DIPDIR_FIELD' INPUT_DIP_FIELD = 'DIP_FIELD' @@ -145,6 +146,16 @@ def initAlgorithm(self, config: Optional[dict[str, Any]] = None) -> None: [QgsProcessing.SourceType.TypeVectorPolygon], ) ) + # Faults are only used when the basal contacts are extracted from the + # geology layer. + self.addParameter( + QgsProcessingParameterFeatureSource( + self.INPUT_FAULTS, + "Faults (optional, used to extract basal contacts)", + [QgsProcessing.SourceType.TypeVectorLine], + optional=True, + ) + ) self.addParameter( QgsProcessingParameterField( @@ -251,7 +262,10 @@ def processAlgorithm( bounding_box_type = self.parameterAsEnum(parameters, self.INPUT_BOUNDING_BOX_TYPE, context) max_line_length = self.parameterAsSource(parameters, self.INPUT_MAX_LINE_LENGTH, context) basal_contacts = self.parameterAsSource(parameters, self.INPUT_BASAL_CONTACTS, context) - geology_data = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + geology_source = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) + if geology_source is None: + raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_GEOLOGY)) + faults_source = self.parameterAsSource(parameters, self.INPUT_FAULTS, context) structure_data = self.parameterAsSource(parameters, self.INPUT_STRUCTURE_DATA, context) orientation_type = self.parameterAsEnum(parameters, self.INPUT_ORIENTATION_TYPE, context) is_strike = orientation_type == 1 @@ -262,15 +276,19 @@ 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) + # Use the geology source for the units, the bounding box and the basal + # contacts. The source follows "selected features only", and it is + # also correct when the input is not a map layer. + units = qgsLayerToDataFrame(geology_source) + geology_data = qgsLayerToGeoDataFrame(geology_source) if bounding_box_type == 0: - extent = geology_layer.extent() + minx, miny, maxx, maxy = geology_data.total_bounds bounding_box = { - 'minx': extent.xMinimum(), - 'miny': extent.yMinimum(), - 'maxx': extent.xMaximum(), - 'maxy': extent.yMaximum(), + 'minx': minx, + 'miny': miny, + 'maxx': maxx, + 'maxy': maxy, } feedback.pushInfo("Using bounding box from geology layer") else: @@ -314,21 +332,30 @@ def processAlgorithm( strati_column_settings = QgsSettings() strati_column_settings.setValue('m2l/strati_column', matrix) # convert layers to dataframe or geodataframe - units = qgsLayerToDataFrame(geology_data) - geology_data = qgsLayerToGeoDataFrame(geology_data) 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'] + # Reproject the faults to the CRS of the geology, so that the two + # layers align. + try: + basal_contacts = extract_basal_contacts( + geology=geology_source, + stratigraphic_order=stratigraphic_order, + faults=faults_source, + unit_name_field=unit_name_field, + updater=feedback.pushInfo, + target_crs=geology_source.sourceCrs(), + )['basal_contacts'] + except Exception as err: + raise QgsProcessingException(str(err)) from err + if basal_contacts.empty: + raise QgsProcessingException( + "No basal contacts were found in the geology layer. Make sure that " + "the units in the stratigraphic column are in the geology layer." + ) structure_data = qgsLayerToDataFrame(structure_data) rename_map = {} missing_fields = [] diff --git a/tests/qgis/test_thickness_calculator.py b/tests/qgis/test_thickness_calculator.py new file mode 100644 index 0000000..176eba9 --- /dev/null +++ b/tests/qgis/test_thickness_calculator.py @@ -0,0 +1,148 @@ +"""The thickness calculator algorithm extracts basal contacts from geology. + +When no basal contacts layer is given, the algorithm gets the basal contacts +from the geology layer and the stratigraphic column. The map2loop calculator +is replaced with a fake, so these tests examine only the data that the +algorithm gives to it. +""" + +import unittest +from unittest.mock import patch + +import pandas as pd +from qgis.core import ( + QgsFeature, + QgsGeometry, + QgsProcessing, + QgsProcessingContext, + QgsProcessingFeatureSourceDefinition, + QgsProcessingFeedback, + QgsProject, + QgsVectorLayer, +) +from qgis.testing import start_app + +from loopstructural.processing.algorithms.thickness_calculator import ( + ThicknessCalculatorAlgorithm, +) + +MODULE = 'loopstructural.processing.algorithms.thickness_calculator' + + +class FakeCalculator: + """Keep the arguments of compute and return a thickness for each unit.""" + + calls = [] + + def __init__(self, **kwargs): + self.kwargs = kwargs + + def compute(self, units, stratigraphic_order, basal_contacts, *args): + FakeCalculator.calls.append({'kwargs': self.kwargs, 'basal_contacts': basal_contacts}) + return pd.DataFrame( + { + 'name': list(units['name']), + 'ThicknessMean': 1.0, + 'ThicknessMedian': 1.0, + 'ThicknessStdDev': 0.0, + } + ) + + +class TestThicknessFromGeology(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.qgs = start_app() + + def setUp(self): + FakeCalculator.calls = [] + self.geology = self._geology() + QgsProject.instance().addMapLayer(self.geology) + + def tearDown(self): + QgsProject.instance().removeMapLayer(self.geology.id()) + + def _geology(self): + # Two blocks of the layers A, B and C (top to bottom). The first block + # is at 0 <= x <= 100 and the second block is at 1000 <= x <= 1100. + layer = QgsVectorLayer("Polygon?crs=EPSG:28350&field=unitname:string", "geology", "memory") + features = [] + for x0 in (0, 1000): + x1 = x0 + 100 + for name, (y0, y1) in {'A': (200, 300), 'B': (100, 200), 'C': (0, 100)}.items(): + feature = QgsFeature(layer.fields()) + feature.setAttribute('unitname', name) + feature.setGeometry( + QgsGeometry.fromWkt( + f"POLYGON(({x0} {y0}, {x1} {y0}, {x1} {y1}, {x0} {y1}, {x0} {y0}))" + ) + ) + features.append(feature) + layer.dataProvider().addFeatures(features) + return layer + + def _points(self, fields, values): + layer = QgsVectorLayer(f"Point?crs=EPSG:28350&{fields}", "points", "memory") + features = [] + for x, y, attributes in values: + feature = QgsFeature(layer.fields()) + feature.setAttributes(attributes) + feature.setGeometry(QgsGeometry.fromWkt(f"POINT({x} {y})")) + features.append(feature) + layer.dataProvider().addFeatures(features) + return layer + + def _run(self, geology): + structure = self._points("field=DIPDIR:double&field=DIP:double", [(50, 150, [90.0, 10.0])]) + sampled_contacts = self._points("field=unit:string", [(50, 200, ['A'])]) + algorithm = ThicknessCalculatorAlgorithm() + algorithm.initAlgorithm() + parameters = { + 'THICKNESS_CALCULATOR_TYPE': 0, + 'BOUNDING_BOX_TYPE': 0, + 'MAX_LINE_LENGTH': 1000, + 'GEOLOGY': geology, + 'UNIT_NAME_FIELD': 'unitname', + 'STRATIGRAPHIC_COLUMN': ['A', 'B', 'C'], + 'SAMPLED_CONTACTS': sampled_contacts, + 'STRUCTURE_DATA': structure, + 'ORIENTATION_TYPE': 0, + 'DIPDIR_FIELD': 'DIPDIR', + 'DIP_FIELD': 'DIP', + 'THICKNESS': QgsProcessing.TEMPORARY_OUTPUT, + } + context = QgsProcessingContext() + context.setProject(QgsProject.instance()) + with patch(f'{MODULE}.InterpolatedStructure', FakeCalculator): + algorithm.processAlgorithm(parameters, context, QgsProcessingFeedback()) + self.assertEqual(len(FakeCalculator.calls), 1) + return FakeCalculator.calls[0] + + def _length(self, contacts, unit): + return contacts[contacts['basal_unit'] == unit].geometry.length.sum() + + def test_contacts_from_geology_layer(self): + call = self._run(self.geology) + contacts = call['basal_contacts'] + self.assertFalse(contacts.empty) + # Each block has one contact of 100 m for A (on B) and for B (on C). + self.assertAlmostEqual(self._length(contacts, 'A'), 200, delta=5) + self.assertAlmostEqual(self._length(contacts, 'B'), 200, delta=5) + self.assertAlmostEqual(call['kwargs']['bounding_box']['maxx'], 1100) + + def test_contacts_from_selected_features_only(self): + # Select only the first block. + self.geology.selectByExpression('x(centroid($geometry)) < 500') + self.assertEqual(self.geology.selectedFeatureCount(), 3) + call = self._run(QgsProcessingFeatureSourceDefinition(self.geology.id(), True)) + contacts = call['basal_contacts'] + self.assertFalse(contacts.empty) + self.assertLessEqual(contacts.total_bounds[2], 100 + 1e-6) + self.assertAlmostEqual(self._length(contacts, 'A'), 100, delta=5) + self.assertAlmostEqual(self._length(contacts, 'B'), 100, delta=5) + self.assertAlmostEqual(call['kwargs']['bounding_box']['maxx'], 100) + + +if __name__ == '__main__': + unittest.main() From d33c2dfb57b341b116c4a7bbd02dfe48001dd26d Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Sat, 31 Oct 2026 11:16:50 +1030 Subject: [PATCH 4/4] fix: two small bugs of the thickness and basal contacts tools - extract_basal_contacts had a mutable default (ignore_units=[]) and added None to it, so the list grew at each call. It also changed the list of the caller. Use None and a new list. - The thickness algorithm read MAX_LINE_LENGTH with parameterAsSource, so StructuralPoint always got None. Use parameterAsDouble. --- loopstructural/main/m2l_api.py | 5 +++-- .../processing/algorithms/thickness_calculator.py | 2 +- tests/qgis/test_basal_contacts.py | 15 +++++++++++++++ 3 files changed, 19 insertions(+), 3 deletions(-) diff --git a/loopstructural/main/m2l_api.py b/loopstructural/main/m2l_api.py index 734d089..884c87b 100644 --- a/loopstructural/main/m2l_api.py +++ b/loopstructural/main/m2l_api.py @@ -249,7 +249,7 @@ def extract_basal_contacts( geology, stratigraphic_order, faults=None, - ignore_units=[], + ignore_units=None, unit_name_field=None, all_contacts=False, updater=None, @@ -320,7 +320,8 @@ def extract_basal_contacts( geology = geology.drop(columns=['UNITNAME']) geology = geology.rename(columns={unit_name_field: 'UNITNAME'}) # Log parameters via DebugManager if provided - ignore_units += [None] + # a new list: do not change the list of the caller + ignore_units = list(ignore_units or []) + [None] # Units in the geology layer that are not in the stratigraphic column are # kept in the geology (so they still bound the other units and appear in diff --git a/loopstructural/processing/algorithms/thickness_calculator.py b/loopstructural/processing/algorithms/thickness_calculator.py index fe1cc56..f280ec8 100644 --- a/loopstructural/processing/algorithms/thickness_calculator.py +++ b/loopstructural/processing/algorithms/thickness_calculator.py @@ -260,7 +260,7 @@ def processAlgorithm( thickness_type = ['InterpolatedStructure', 'StructuralPoint'][thickness_type_index] dtm_data = self.parameterAsRasterLayer(parameters, self.INPUT_DTM, context) bounding_box_type = self.parameterAsEnum(parameters, self.INPUT_BOUNDING_BOX_TYPE, context) - max_line_length = self.parameterAsSource(parameters, self.INPUT_MAX_LINE_LENGTH, context) + max_line_length = self.parameterAsDouble(parameters, self.INPUT_MAX_LINE_LENGTH, context) basal_contacts = self.parameterAsSource(parameters, self.INPUT_BASAL_CONTACTS, context) geology_source = self.parameterAsSource(parameters, self.INPUT_GEOLOGY, context) if geology_source is None: diff --git a/tests/qgis/test_basal_contacts.py b/tests/qgis/test_basal_contacts.py index bbf2950..a955824 100644 --- a/tests/qgis/test_basal_contacts.py +++ b/tests/qgis/test_basal_contacts.py @@ -1,3 +1,4 @@ +import inspect import unittest from pathlib import Path from unittest.mock import patch @@ -222,6 +223,20 @@ def test_extraction_error_is_raised(self): with self.assertRaisesRegex(RuntimeError, 'bad geology'): self._extract(['A', 'B', 'C']) + def test_the_list_of_ignored_units_does_not_change(self): + ignore_units = ['D'] + for _ in range(2): + extract_basal_contacts( + geology=self._geology(), + stratigraphic_order=['A', 'B', 'C'], + ignore_units=ignore_units, + unit_name_field='unitname', + ) + self.assertEqual(ignore_units, ['D']) + # the default is not a shared list that grows at each call + default = inspect.signature(extract_basal_contacts).parameters['ignore_units'].default + self.assertIsNone(default) + if __name__ == '__main__': unittest.main()