diff --git a/loopstructural/main/m2l_api.py b/loopstructural/main/m2l_api.py index 6bb22a4..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, @@ -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: @@ -315,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 @@ -400,10 +406,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/loopstructural/processing/algorithms/thickness_calculator.py b/loopstructural/processing/algorithms/thickness_calculator.py index 299a6c5..f280ec8 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, @@ -52,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' @@ -132,9 +134,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( @@ -144,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( @@ -248,9 +260,12 @@ 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_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 @@ -261,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) + # 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: - geology_layer = self.parameterAsVectorLayer(parameters, self.INPUT_GEOLOGY, context) - 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: @@ -313,9 +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) - 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..." + ) + # 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_basal_contacts.py b/tests/qgis/test_basal_contacts.py index 385a8de..a955824 100644 --- a/tests/qgis/test_basal_contacts.py +++ b/tests/qgis/test_basal_contacts.py @@ -1,5 +1,7 @@ +import inspect import unittest from pathlib import Path +from unittest.mock import patch from qgis.core import ( QgsFeature, QgsGeometry, @@ -212,6 +214,29 @@ 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']) + + 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() 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()