From 589d71176e9fcc6c9b168f2c4519eda63e3369f4 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Wed, 7 Oct 2026 10:23:56 -0300 Subject: [PATCH] feat: model step with one build button and direct interpolation Phase 4 of the usability plan. - One primary button replaces Initialize Model, Solve Model and Update Model Data. Its text and action come from the model state. The problems of all steps show above the feature list. - A build calculates out-of-date basal contacts and thicknesses first, on a background task with progress. With "Calculate from geology polygons" the model reads the extracted contacts directly; the project layer is for display only. - Start choice: build from a geological map, or interpolate surfaces from constraints (hides steps 2 and 3). Saved with the state. - Constraint types: value, interface, gradient/normal, tangent, form line, inequality and pairwise inequality, each with a weight and a Z source (layer, DEM or constant). - Detach a generated feature to keep a fixed copy of its data. Rows added to a generated feature are now used in the build. - Preview one feature as isolines on the map. - Add Fault from a centre, strike, dip and size. --- docs/development/usability-plan.md | 27 +- .../geological_model_tab/add_fault_dialog.py | 153 +++++++- .../geological_model_tab/add_fault_dialog.ui | 141 ------- .../feature_details_panel/_base.py | 214 ++++++++++ .../geological_model_tab.py | 371 ++++++++++++------ .../layer_selection_table.py | 141 ++++++- .../gui/modelling/modelling_widget.py | 55 ++- .../gui/modelling/steps/build_plan.py | 111 ++++++ loopstructural/gui/modelling/steps/checks.py | 5 +- loopstructural/gui/modelling/steps/pages.py | 29 ++ .../gui/modelling/steps/step_bar.py | 8 + loopstructural/main/constraints.py | 186 +++++++++ loopstructural/main/data_manager.py | 93 ++++- loopstructural/main/derived_refresh.py | 316 +++++++++++++++ loopstructural/main/model_manager.py | 267 +++++++++++-- loopstructural/main/parametric_fault.py | 168 ++++++++ loopstructural/main/preview.py | 106 +++++ loopstructural/main/workflow_mode.py | 19 + tests/qgis/test_derived_data_state.py | 35 ++ tests/qgis/test_detach_feature.py | 128 ++++++ tests/qgis/test_manual_foliations.py | 57 +++ tests/qgis/test_parametric_fault.py | 97 +++++ tests/unit/test_build_plan.py | 174 ++++++++ tests/unit/test_constraints.py | 173 ++++++++ tests/unit/test_parametric_fault.py | 113 ++++++ tests/unit/test_preview.py | 79 ++++ 26 files changed, 2930 insertions(+), 336 deletions(-) delete mode 100644 loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.ui create mode 100644 loopstructural/gui/modelling/steps/build_plan.py create mode 100644 loopstructural/main/constraints.py create mode 100644 loopstructural/main/derived_refresh.py create mode 100644 loopstructural/main/parametric_fault.py create mode 100644 loopstructural/main/preview.py create mode 100644 loopstructural/main/workflow_mode.py create mode 100644 tests/qgis/test_detach_feature.py create mode 100644 tests/qgis/test_parametric_fault.py create mode 100644 tests/unit/test_build_plan.py create mode 100644 tests/unit/test_constraints.py create mode 100644 tests/unit/test_parametric_fault.py create mode 100644 tests/unit/test_preview.py diff --git a/docs/development/usability-plan.md b/docs/development/usability-plan.md index b3a9270..d0f3278 100644 --- a/docs/development/usability-plan.md +++ b/docs/development/usability-plan.md @@ -302,27 +302,27 @@ Acceptance: a new user can go from an empty project to a solved model with the ### Phase 4: Model step and direct interpolation -- [ ] Replace "Initialize Model", "Solve Model" and "Update Model Data" with one +- [x] Replace "Initialize Model", "Solve Model" and "Update Model Data" with one primary button. Its text and action come from the model state. -- [ ] Show the problems from all steps before the build. -- [ ] Before the build, calculate all out-of-date derived data again: basal +- [x] Show the problems from all steps before the build. +- [x] Before the build, calculate all out-of-date derived data again: basal contacts first, then calculated thicknesses. Run this in a background task with progress. If the calculation fails, stop the build and show the error. -- [ ] With "Calculate from geology polygons", give the extracted contacts +- [x] With "Calculate from geology polygons", give the extracted contacts directly to the model. Update the project contacts layer for display only. -- [ ] Add the start choice: "Build from a geological map" or "Interpolate +- [x] Add the start choice: "Build from a geological map" or "Interpolate surfaces from constraints". Save the choice with the state. -- [ ] Constraint list for each feature: source layer, constraint type, field +- [x] Constraint list for each feature: source layer, constraint type, field mapping, weight, Z source (layer Z, DEM or constant). -- [ ] Add the constraint types: value, interface, gradient/normal, tangent, +- [x] Add the constraint types: value, interface, gradient/normal, tangent, inequality, pairwise inequality. -- [ ] Show generated constraints as read-only rows. Add "Detach" to make a +- [x] Show generated constraints as read-only rows. Add "Detach" to make a generated feature editable. -- [ ] Build and preview one feature: an isoline on the map canvas, or a surface +- [x] Build and preview one feature: an isoline on the map canvas, or a surface in the 3D view. -- [ ] Implement "Add Fault" in the model step. +- [x] Implement "Add Fault" in the model step. Files: `geological_model_tab/*.py`, `layer_selection_table.py`, `feature_details_panel/*.py`, `main/model_manager.py`. @@ -361,9 +361,10 @@ model, the model uses basal contacts and thicknesses that match the new order. done), or only a guide? Recommendation: only a guide. Show problems, but do not lock steps. 2. Is the Data Conversion dialog part of step 1, or a separate tool? -3. In direct mode, which constraint types does each interpolator (FDI, PLI, - surfe) support? The UI must hide types that the selected interpolator does - not support. +3. ~~In direct mode, which constraint types does each interpolator (FDI, PLI, + surfe) support?~~ Resolved: all interpolators accept the same constraint + types. LoopStructural converts them as needed. The UI shows the same types + for every interpolator and does not hide any. 4. Is there a minimum QGIS version for the new widgets? 5. Must the Processing algorithms also use the derived-data record, or only the dock? Recommendation: only the dock. Processing runs are single runs diff --git a/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.py b/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.py index 7d35d46..64aac59 100644 --- a/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.py +++ b/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.py @@ -1,33 +1,144 @@ -import os +from qgis.PyQt.QtWidgets import ( + QDialog, + QDialogButtonBox, + QDoubleSpinBox, + QFormLayout, + QGroupBox, + QHBoxLayout, + QLabel, + QLineEdit, + QVBoxLayout, +) -from qgis.PyQt.QtWidgets import QDialog -from qgis.PyQt.uic import loadUi +from ....main import parametric_fault + + +def _spin(low, high, decimals=2, step=1.0): + spin = QDoubleSpinBox() + spin.setRange(low, high) + spin.setDecimals(decimals) + spin.setSingleStep(step) + return spin class AddFaultDialog(QDialog): - def __init__(self, parent=None): + """Ask for a fault that is given by numbers: a centre, a strike, a dip and a size. + + Faults that have a trace on the map come from the fault layer in step 3. + The start values come from the bounding box of the model. + """ + + def __init__(self, parent=None, *, model_manager=None): super().__init__(parent) - ui_path = os.path.join(os.path.dirname(__file__), 'add_fault_dialog.ui') - loadUi(ui_path, self) - self.setWindowTitle('Add Fault Feature') - # You can access widgets by their objectName from the .ui file - # Example: self.strike_input, self.dip_input, etc. + self.model_manager = model_manager + self.setWindowTitle('Add Fault') + + layout = QVBoxLayout(self) + name_form = QFormLayout() + self.name_input = QLineEdit() + name_form.addRow("Name:", self.name_input) + layout.addLayout(name_form) + + orientation = QGroupBox("Orientation") + form = QFormLayout(orientation) + self.strike_input = _spin(0, 360) + self.strike_input.setToolTip("Degrees clockwise from north. The fault dips to the right.") + self.dip_input = _spin(0.01, 90) + self.pitch_input = _spin(-180, 180) + self.pitch_input.setToolTip( + "The angle of the slip in the fault plane, from the strike direction towards " + "the down-dip direction. 0 is strike-slip. 90 is dip-slip." + ) + self.displacement_input = _spin(-1e9, 1e9) + self.displacement_input.setToolTip("The size of the slip, in model units.") + form.addRow("Strike (°):", self.strike_input) + form.addRow("Dip (°):", self.dip_input) + form.addRow("Pitch (°):", self.pitch_input) + form.addRow("Displacement:", self.displacement_input) + layout.addWidget(orientation) + + position = QGroupBox("Position") + form = QFormLayout(position) + centre_row = QHBoxLayout() + self.centre_inputs = [_spin(-1e9, 1e9) for _ in range(3)] + for label, spin in zip("XYZ", self.centre_inputs): + centre_row.addWidget(QLabel(label)) + centre_row.addWidget(spin) + form.addRow("Centre:", centre_row) + layout.addWidget(position) + + size = QGroupBox("Size") + form = QFormLayout(size) + self.major_input = _spin(0.01, 1e9) + self.major_input.setToolTip("How far the fault extends along its strike.") + self.intermediate_input = _spin(0.01, 1e9) + self.intermediate_input.setToolTip("How far the fault extends down its dip.") + self.minor_input = _spin(0.01, 1e9) + self.minor_input.setToolTip("How far the fault changes the field, across the fault.") + form.addRow("Length along strike:", self.major_input) + form.addRow("Extent down dip:", self.intermediate_input) + form.addRow("Influence distance:", self.minor_input) + layout.addWidget(size) + + self.problem_label = QLabel() + self.problem_label.setWordWrap(True) + layout.addWidget(self.problem_label) + self.button_box = QDialogButtonBox( + QDialogButtonBox.StandardButton.Ok | QDialogButtonBox.StandardButton.Cancel + ) + self.button_box.accepted.connect(self.accept) + self.button_box.rejected.connect(self.reject) + layout.addWidget(self.button_box) + + self._set_start_values() + for spin in ( + self.strike_input, + self.dip_input, + self.pitch_input, + self.displacement_input, + self.major_input, + self.intermediate_input, + self.minor_input, + *self.centre_inputs, + ): + spin.valueChanged.connect(self._validate) + self.name_input.textChanged.connect(self._validate) + self._validate() + + def _set_start_values(self): + model = getattr(self.model_manager, 'model', None) + box = getattr(model, 'bounding_box', None) + if box is None: + spec = parametric_fault.default_spec([0, 0, 0], [1000, 1000, 1000]) + else: + spec = parametric_fault.default_spec(box.origin, box.maximum) + self.strike_input.setValue(spec['strike']) + self.dip_input.setValue(spec['dip']) + self.pitch_input.setValue(spec['pitch']) + self.displacement_input.setValue(spec['displacement']) + for spin, value in zip(self.centre_inputs, spec['centre']): + spin.setValue(value) + self.major_input.setValue(spec['major_axis']) + self.intermediate_input.setValue(spec['intermediate_axis']) + self.minor_input.setValue(spec['minor_axis']) + + def _taken_names(self): + return self.model_manager.used_names() if self.model_manager is not None else set() + + def _validate(self, *args): + found = parametric_fault.problems(self.get_fault_data(), self._taken_names()) + self.problem_label.setText("\n".join(found)) + self.button_box.button(QDialogButtonBox.StandardButton.Ok).setEnabled(not found) def get_fault_data(self): return { + 'name': self.name_input.text().strip(), 'strike': self.strike_input.value(), 'dip': self.dip_input.value(), - 'centre': ( - self.centre_x_input.value(), - self.centre_y_input.value(), - self.centre_z_input.value(), - ), - 'ellipsoid_extents': ( - self.extent_x_input.value(), - self.extent_y_input.value(), - self.extent_z_input.value(), - ), - 'displacement': self.displacement_input.value(), 'pitch': self.pitch_input.value(), - 'name': self.name_input.text(), + 'displacement': self.displacement_input.value(), + 'centre': tuple(spin.value() for spin in self.centre_inputs), + 'major_axis': self.major_input.value(), + 'intermediate_axis': self.intermediate_input.value(), + 'minor_axis': self.minor_input.value(), } diff --git a/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.ui b/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.ui deleted file mode 100644 index 1d93ee4..0000000 --- a/loopstructural/gui/modelling/geological_model_tab/add_fault_dialog.ui +++ /dev/null @@ -1,141 +0,0 @@ - - - AddFaultDialog - - - Add Fault Feature - - - - - - - - Name - - - - - - - - - - - - Orientation - - - - - - Strike (°) - - - - - - - - - - Dip (°) - - - - - - - - - - Pitch (°) - - - - - - - - - - Displacement - - - - - - - - - - - - - Position - - - - - - Centre (X, Y, Z) - - - - - - - - - - - - - - - - - - - - - - - Ellipsoid Size - - - - - - Extents (X, Y, Z) - - - - - - - - - - - - - - - - - - - - - - - QDialogButtonBox::Ok|QDialogButtonBox::Cancel - - - - - - - - diff --git a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py index c528973..3bce57d 100644 --- a/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py +++ b/loopstructural/gui/modelling/geological_model_tab/feature_details_panel/_base.py @@ -1,3 +1,4 @@ +import numpy as np from LoopStructural.modelling.features import StructuralFrame from LoopStructural.utils import normal_vector_to_strike_and_dip from qgis.gui import QgsCollapsibleGroupBox, QgsMapLayerComboBox @@ -6,17 +7,22 @@ QComboBox, QDoubleSpinBox, QFormLayout, + QHBoxLayout, QLabel, QMessageBox, QPushButton, QScrollArea, + QSpinBox, QVBoxLayout, QWidget, ) from qgis.utils import plugins from LoopStructural import getLogger +from loopstructural.main import preview +from ....background_task import finish_background_task, start_background_task +from ....messages import push_info, push_warning from ..bounding_box_widget import BoundingBoxWidget from ..layer_selection_table import LayerSelectionTable @@ -149,6 +155,7 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage name_validator=lambda: (True, ''), # Always valid in this context ) table_layout = QVBoxLayout() + table_layout.addWidget(self._build_detach_widget()) table_layout.addWidget(self.layer_table) self.view_constraint_data_button = QPushButton("View Data Used by Interpolator") self.view_constraint_data_button.setToolTip( @@ -170,11 +177,218 @@ def __init__(self, parent=None, *, feature=None, model_manager=None, data_manage group_box = QgsCollapsibleGroupBox('Interpolator Settings') group_box.setLayout(form_layout) self.layout.addWidget(group_box) + self.layout.addWidget(self._build_preview_widget()) self.layout.addWidget(table_group_box) # this will call the addMidBlock and addExportBlock methods self.addMidBlock() self.addExportBlock() + def _build_preview_widget(self): + """Return the group with the button that shows the isolines of this feature on the map.""" + group = QgsCollapsibleGroupBox('Preview') + row = QHBoxLayout(group) + self.preview_levels_spin = QSpinBox() + self.preview_levels_spin.setRange(1, 100) + self.preview_levels_spin.setValue(preview.DEFAULT_LEVEL_COUNT) + self.preview_levels_spin.setPrefix("Lines: ") + self.preview_levels_spin.setToolTip("The number of isolines.") + self.preview_button = QPushButton("Preview on map") + self.preview_button.setToolTip( + "Solve this feature only, and add its isolines on the ground (the DEM) to the " + "project as a temporary layer. The other features are not solved." + ) + self.preview_button.clicked.connect(self.preview_on_map) + row.addWidget(self.preview_levels_spin) + row.addWidget(self.preview_button, 1) + self._preview_layer_id = None + return group + + def preview_on_map(self): + """Solve this feature and show its isolines on the map canvas. + + The grid is read on the GUI thread (the DEM can be a raster layer). The + solve and the lines are made on a background thread. + """ + manager = self.model_manager + name = self.feature.name + model = getattr(manager, 'model', None) + if model is None or model.bounding_box is None: + QMessageBox.warning(self, "Preview", "Build the model first.") + return + bounding_box = model.bounding_box + x, y, points = preview.grid_points( + bounding_box.origin, + bounding_box.maximum, + preview.DEFAULT_RESOLUTION, + getattr(manager, 'dem_function', None), + ) + count = self.preview_levels_spin.value() + + def target(progress_callback): + progress_callback("Solving the feature...") + values = np.asarray(manager.evaluate_feature_on_points(name, points), dtype=float) + if values.ndim != 1: + raise ValueError("The feature has no scalar field to draw.") + values = values.reshape(len(y), len(x)) + progress_callback("Making the lines...") + levels = preview.default_levels(values, count) + return preview.isolines(x, y, values, levels), levels + + self.preview_button.setEnabled(False) + self._preview_task = start_background_task( + self, + target, + title="Preview", + initial_label="Solving the feature...", + on_progress=self._on_preview_progress, + on_finished=self._on_preview_finished, + on_error=self._on_preview_error, + ) + + def _on_preview_progress(self, message): + try: + self._preview_task[2].setLabelText(message) + except Exception: + pass + + def _end_preview_task(self): + finish_background_task(*self._preview_task) + self.preview_button.setEnabled(True) + + def _on_preview_error(self, traceback_text): + self._end_preview_task() + lines = traceback_text.strip().splitlines() + QMessageBox.critical( + self, "Preview failed", lines[-1] if lines else "The preview stopped with an error." + ) + + def _on_preview_finished(self, result): + self._end_preview_task() + lines, levels = result + if not lines: + push_warning( + "Preview", + f"'{self.feature.name}' has no isolines on the map. The field has no range " + "or no values in the model area.", + ) + return + self._show_preview_lines(lines) + # the feature is solved now, so the ticks of the feature list change + self.model_manager.notify('model_updated') + push_info( + "Preview", + f"{len(lines)} lines for '{self.feature.name}', from {levels.min():g} to {levels.max():g}.", + ) + + def _show_preview_lines(self, lines): + """Add the isolines as a temporary layer. An earlier preview of this feature is replaced.""" + from qgis.core import ( + QgsFeature, + QgsField, + QgsGeometry, + QgsPointXY, + QgsProject, + QgsVectorLayer, + ) + + from loopstructural.gui.compatibility import QVariantCompat + + project = QgsProject.instance() + previous = getattr(self, '_preview_layer_id', None) + if previous and project.mapLayer(previous) is not None: + project.removeMapLayer(previous) + crs = self.data_manager.get_model_crs() if self.data_manager is not None else None + crs_text = crs.authid() if crs is not None and crs.isValid() else project.crs().authid() + layer = QgsVectorLayer( + f"LineString?crs={crs_text}", f"Preview: {self.feature.name}", "memory" + ) + provider = layer.dataProvider() + provider.addAttributes([QgsField("value", QVariantCompat.Double)]) + layer.updateFields() + features = [] + for level, coordinates in lines: + feature = QgsFeature(layer.fields()) + feature.setGeometry( + QgsGeometry.fromPolylineXY([QgsPointXY(px, py) for px, py in coordinates]) + ) + feature.setAttribute("value", level) + features.append(feature) + provider.addFeatures(features) + layer.updateExtents() + project.addMapLayer(layer) + self._preview_layer_id = layer.id() + + def _build_detach_widget(self): + """Return the row that tells where the data of a generated feature comes from. + + A generated feature is one that the stratigraphic column makes. Its + contact and orientation rows are read-only. The user can add rows to + them, or detach the feature to keep a fixed copy of the data. + """ + self.detach_widget = QWidget() + row = QHBoxLayout(self.detach_widget) + row.setContentsMargins(0, 0, 0, 0) + self.detach_label = QLabel() + self.detach_label.setWordWrap(True) + self.detach_button = QPushButton() + self.detach_button.clicked.connect(self._toggle_detached) + row.addWidget(self.detach_label, 1) + row.addWidget(self.detach_button) + self._refresh_detach_widget() + return self.detach_widget + + def _refresh_detach_widget(self): + manager = self.model_manager + name = self.feature.name + if manager is None or not manager.is_generated(name): + self.detach_widget.hide() + return + self.detach_widget.show() + if manager.is_detached(name): + self.detach_label.setText( + "Detached. The rows from the column are a fixed copy. They do not change " + "when the column or the contacts change." + ) + self.detach_button.setText("Follow the column again") + self.detach_button.setToolTip( + "Use the data of the stratigraphic column again. The model needs a rebuild." + ) + else: + self.detach_label.setText( + "The read-only rows come from the stratigraphic column. Rows that you " + "add are used together with them." + ) + self.detach_button.setText("Detach") + self.detach_button.setToolTip( + "Keep a fixed copy of the data of this feature, so that it does not " + "follow the column." + ) + + def _toggle_detached(self): + manager = self.model_manager + name = self.feature.name + if manager.is_detached(name): + reply = QMessageBox.question( + self, + "Follow the column again", + f"'{name}' will use the data of the stratigraphic column again. " + "The fixed copy is removed. Continue?", + ) + if reply != QMessageBox.StandardButton.Yes: + return + manager.attach_feature(name) + push_info("Feature", f"'{name}' follows the stratigraphic column again. Rebuild the model.") + elif manager.detach_feature(name): + push_info("Feature", f"'{name}' is detached. Its data is a fixed copy.") + else: + QMessageBox.warning( + self, "Cannot detach", f"'{name}' has no data from the column to keep." + ) + return + self.data_manager._sync_processed_feature_data() + self.layer_table.restore_table_state() + self._refresh_detach_widget() + def addMidBlock(self): """Base mid block is intentionally empty now — bounding-box controls were moved into the export/evaluation section so they appear alongside diff --git a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py index 2c2a78b..3550252 100644 --- a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py +++ b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py @@ -1,3 +1,5 @@ +import html + from LoopStructural.modelling.features import FeatureType from qgis.PyQt.QtCore import QObject, Qt, QThread, pyqtSignal, pyqtSlot from qgis.PyQt.QtGui import QColor, QIcon, QPainter, QPen, QPixmap @@ -16,7 +18,12 @@ QWidget, ) +from ....main.derived_refresh import DerivedRefresh, DerivedRefreshError, names_to_refresh from ....main.model_manager import ModelSolveCancelled +from ....main.workflow_mode import DEFAULT_WORKFLOW_MODE, WORKFLOW_MODE_CONSTRAINTS +from ...messages import push_info +from ..steps import build_plan +from .add_fault_dialog import AddFaultDialog from .add_foliation_dialog import AddFoliationDialog from .add_unconformity_dialog import AddUnconformityDialog from .feature_details_panel import ( @@ -67,14 +74,14 @@ def _build_status_icon(color: str, *, filled: bool, mark: str = None) -> QIcon: 'solved': "Model status: solved", 'stale': ( "Model status: faults, stratigraphic column or input data changed — " - "re-run Initialize Model" + "rebuild the model" ), } -# Solve Model only rebuilds interpolators for features that already exist; it -# can't apply a fault topology edit (that requires re-running Initialize -# Model, see GeologicalModelManager._on_fault_topology_changed), so it stays -# disabled outside these two states. +# Solving only rebuilds interpolators for features that already exist; it +# can't apply a fault topology edit (that requires a build, see +# GeologicalModelManager._on_fault_topology_changed), so the primary button +# builds in the other states. _SOLVABLE_STATES = {'initialized', 'solved'} @@ -153,17 +160,18 @@ def __init__(self, parent=None, *, model_manager=None, data_manager=None): splitter.setStretchFactor(1, 0) # Feature details panel splitter.setOrientation(Qt.Orientation.Horizontal) # Add horizontal slider - # Initialize / Solve Model buttons + a status summary of where the - # model currently is: empty -> initialized (unsolved) -> solved. - self.initializeModelButton = QPushButton("Initialize Model") - self.solveModelButton = QPushButton("Solve Model") - self.solveModelButton.setEnabled(False) # nothing to solve until initialized + # One primary button. Its text and action come from the model state + # (see build_plan.choose_primary_action). The status label shows where + # the model is: empty -> initialized (unsolved) -> solved. + self.primaryButton = QPushButton("Build model") + self.primaryButton.clicked.connect(self.on_primary_clicked) + self._primary_action = build_plan.PrimaryAction(build_plan.ACTION_BUILD, "Build model") + self._task_running = False self.modelStatusLabel = QLabel(_MODEL_STATE_LABELS['empty']) buttonRow = QHBoxLayout() buttonRow.setContentsMargins(0, 0, 0, 0) - buttonRow.addWidget(self.initializeModelButton) - buttonRow.addWidget(self.solveModelButton) + buttonRow.addWidget(self.primaryButton) buttonRow.addStretch(1) buttonRow.addWidget(self.modelStatusLabel) buttonRowWidget = QWidget() @@ -173,36 +181,17 @@ def __init__(self, parent=None, *, model_manager=None, data_manager=None): buttonRowWidget.setSizePolicy(QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Fixed) mainLayout.insertWidget(0, buttonRowWidget, 0) - # Shown when an input layer changed after the model data was read - # from it. "Update Model Data" puts the new data into the existing - # features, so the changes made to them are kept (Initialize Model - # builds every feature again). - self.layerChangedLabel = QLabel() - self.layerChangedLabel.setWordWrap(True) - self.updateModelDataButton = QPushButton("Update Model Data") - self.updateModelDataButton.setToolTip( - "Read the changed layers again and put the new data into the current " - "features. Solve Model then uses the new data." - ) - self.updateModelDataButton.clicked.connect(self.update_model_data) - layerChangedRow = QHBoxLayout() - layerChangedRow.setContentsMargins(0, 0, 0, 0) - layerChangedRow.addWidget(self.layerChangedLabel, 1) - layerChangedRow.addWidget(self.updateModelDataButton) - self.layerChangedWidget = QWidget() - self.layerChangedWidget.setLayout(layerChangedRow) - self.layerChangedWidget.setSizePolicy( - QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Fixed - ) - mainLayout.insertWidget(1, self.layerChangedWidget, 0) - self.layerChangedWidget.hide() + # The problems of all steps. They show before the build, so the user + # sees why a build can fail or can give a poor model. + self.problemsLabel = QLabel() + self.problemsLabel.setWordWrap(True) + self.problemsLabel.setTextFormat(Qt.TextFormat.RichText) + self.problemsLabel.setSizePolicy(QSizePolicy.Policy.Preferred, QSizePolicy.Policy.Fixed) + mainLayout.insertWidget(1, self.problemsLabel, 0) + self.problemsLabel.hide() + self._problems_provider = None if self.data_manager is not None: - self.data_manager.add_layer_data_changed_callback(self._refresh_changed_layers) - - # Action buttons - - self.initializeModelButton.clicked.connect(self.initialize_model) - self.solveModelButton.clicked.connect(self.solve_model) + self.data_manager.add_layer_data_changed_callback(self.refresh_primary_action) # Connect feature selection to update details panel self.featureList.itemClicked.connect(self.on_feature_selected) @@ -218,24 +207,37 @@ def __init__(self, parent=None, *, model_manager=None, data_manager=None): def show_add_feature_menu(self, *args): menu = QMenu(self) - add_fault = menu.addAction("Add Fault (not yet implemented)") - # Unlike Add Foliation/Add Unconformity, there's no model_manager entry - # point yet for a parametric (strike/dip/centre) fault -- faults are - # currently only created from trace data via update_fault_points. Keep - # the menu entry visible (so it's discoverable) but disabled, rather - # than silently accepting input and doing nothing with it. - add_fault.setEnabled(False) - add_fault.setToolTip("Adding a fault from parameters isn't implemented yet.") + add_fault = menu.addAction("Add Fault") + add_fault.setToolTip( + "Add a fault from a centre, a strike, a dip and a size. " + "Faults from a trace layer are in step 3." + ) add_foliaton = menu.addAction("Add Foliation") add_unconformity = menu.addAction("Add Unconformity") buttonPosition = self.sender().mapToGlobal(self.sender().rect().bottomLeft()) action = menu.exec(buttonPosition) - if action == add_foliaton: + if action == add_fault: + self.open_add_fault_dialog() + elif action == add_foliaton: self.open_add_foliation_dialog() elif action == add_unconformity: self.open_add_unconformity_dialog() + def open_add_fault_dialog(self): + dialog = AddFaultDialog(self, model_manager=self.model_manager) + if dialog.exec() != dialog.Accepted: + return + try: + self.model_manager.add_parametric_fault(dialog.get_fault_data()) + except ValueError as err: + QMessageBox.critical(self, "Cannot add the fault", str(err)) + return + push_info( + "Fault", f"'{dialog.get_fault_data()['name']}' is added. Rebuild the model to use it." + ) + self.refresh_primary_action() + def open_add_foliation_dialog(self): dialog = AddFoliationDialog( self, data_manager=self.data_manager, model_manager=self.model_manager @@ -250,39 +252,85 @@ def open_add_unconformity_dialog(self): if dialog.exec() == dialog.Accepted: pass - def initialize_model(self): - # Run update_model in a background thread to avoid blocking the UI. - if not self.model_manager: - return - if self.data_manager is not None: - if not self.data_manager.is_bounding_box_set(): - QMessageBox.critical( - self, - "Bounding box required", - "Please set the bounding box before initializing the model.", - ) - return + def set_problems_provider(self, provider): + """Set the function that gives the problems of all steps. - # Validate model CRS - if not self.data_manager.is_model_crs_valid(): - crs = self.data_manager.get_model_crs() - if crs is None or not crs.isValid(): - msg = "Model CRS is not set or invalid. Please select a valid projected CRS in the Model Definition tab." - else: - # Safely get CRS description - try: - crs_desc = crs.description() or crs.authid() or "Unknown" - except Exception: - crs_desc = crs.authid() if hasattr(crs, 'authid') else "Unknown" - msg = f"Model CRS must be projected (in meters), not geographic.\nSelected CRS: {crs_desc}\n\nPlease select a valid projected CRS in the Model Definition tab." + ``provider()`` returns ``(step_key, message)`` pairs. The dock sets it, + because the tab does not know the other steps. + """ + self._problems_provider = provider + self.refresh_primary_action() + + def _blocked_reason(self): + """Return why no build can start now, or None.""" + if self.model_manager is None or self.data_manager is None: + return None + if not self.data_manager.is_bounding_box_set(): + return "Set the bounding box in step 1." + if not self.data_manager.is_model_crs_valid(): + return "Select a projected model CRS in step 1." + return None + + def _workflow_mode(self): + return getattr(self.data_manager, 'workflow_mode', DEFAULT_WORKFLOW_MODE) + + def _derived_to_refresh(self): + """Return the derived results that the build calculates again. + + With "Interpolate surfaces from constraints", there is no column, so + there is nothing to calculate. + """ + if self.data_manager is None or self._workflow_mode() == WORKFLOW_MODE_CONSTRAINTS: + return [] + return names_to_refresh(self.data_manager) - QMessageBox.critical( - self, - "Invalid Model CRS", - msg, - ) - return + def refresh_primary_action(self, *args, **kwargs): + """Show the action of the primary button and the problems of the steps.""" + if self.model_manager is None: + return + changed = self.data_manager.get_changed_layers() if self.data_manager else [] + self._primary_action = build_plan.choose_primary_action( + self.model_manager.model_state, + derived_out_of_date=self._derived_to_refresh(), + layers_changed=bool(changed), + blocked_reason=self._blocked_reason(), + ) + action = self._primary_action + self.primaryButton.setText(action.text) + self.primaryButton.setToolTip(action.tooltip) + if not self._task_running: + self.primaryButton.setEnabled(action.enabled) + + problems = self._problems_provider() if self._problems_provider is not None else [] + if problems: + lines = "".join(f"
  • {html.escape(message)}
    • " for _key, message in problems) + self.problemsLabel.setText( + f"Check before the build:" + ) + self.problemsLabel.show() + else: + self.problemsLabel.hide() + def on_primary_clicked(self): + """Do the action of the primary button.""" + if self.model_manager is None or self._task_running: + return + self.refresh_primary_action() + action = self._primary_action + if action.action == build_plan.ACTION_BUILD: + self.build_model() + elif action.action == build_plan.ACTION_SOLVE: + self.solve_model() + + def build_model(self): + """Calculate the out-of-date derived data, make the features, and solve them. + + The derived data is calculated first, because a build must not use + out-of-date contacts or thicknesses. If a calculation fails, the build + stops. + """ + if not self.model_manager or not self._validate_before_build(): + return if not self._confirm_bounding_box_contains_data(): return @@ -290,55 +338,123 @@ def initialize_model(self): # build from the current layer data, not the data read before # the layers changed self.data_manager.reload_changed_layers() + # the rows that the user added to generated features + self.data_manager.sync_extra_constraints() - self._run_model_task( - lambda progress_callback: self.model_manager.update_model( + refresh = None + names = self._derived_to_refresh() + if names: + try: + refresh = DerivedRefresh(self.data_manager, self.model_manager, names) + except DerivedRefreshError as err: + QMessageBox.critical(self, "Cannot build the model", str(err)) + return + + if refresh is not None: + self._run_model_task( + lambda progress_callback: refresh.run(progress_callback), + title="Updating derived data", + initial_label="Calculating the data that comes from the column...", + cancellable=False, + on_success=lambda: self._after_derived_refresh(refresh), + ) + else: + self._run_build_task() + + def _after_derived_refresh(self, refresh): + """Put the new derived data in place, then make the features.""" + try: + skipped = refresh.finish() + except Exception as err: + QMessageBox.critical(self, "Cannot build the model", f"{type(err).__name__}: {err}") + return + if skipped: + push_info( + "Thickness", + "These units keep the thickness that you typed: " + ", ".join(skipped) + ".", + ) + self._run_build_task() + + def _run_build_task(self): + def target(progress_callback): + self.model_manager.update_model( notify_observers=False, progress_callback=progress_callback - ), - title="Updating Model", - initial_label="Updating geological model...", + ) + self.model_manager.update_all_features( + progress_callback=progress_callback, notify_observers=False + ) + + self._run_model_task( + target, + title="Building Model", + initial_label="Building geological model...", ) - def update_model_data(self): + def _validate_before_build(self): + if self.data_manager is None: + return True + if not self.data_manager.is_bounding_box_set(): + QMessageBox.critical( + self, + "Bounding box required", + "Please set the bounding box before building the model.", + ) + return False + if not self.data_manager.is_model_crs_valid(): + crs = self.data_manager.get_model_crs() + if crs is None or not crs.isValid(): + msg = ( + "Model CRS is not set or invalid. " + "Please select a valid projected CRS in step 1." + ) + else: + # Safely get CRS description + try: + crs_desc = crs.description() or crs.authid() or "Unknown" + except Exception: + crs_desc = crs.authid() if hasattr(crs, 'authid') else "Unknown" + msg = ( + "Model CRS must be projected (in meters), not geographic.\n" + f"Selected CRS: {crs_desc}\n\nPlease select a valid projected CRS in step 1." + ) + QMessageBox.critical(self, "Invalid Model CRS", msg) + return False + return True + + def _update_model_data(self): """Put the data of the changed input layers into the current - features, without Initialize Model.""" - if self.data_manager is None or self.model_manager is None: - return + features, without a build. Returns False if the caller must not solve now.""" try: result = self.data_manager.refresh_model_data() except Exception as e: QMessageBox.critical(self, "Update model data failed", str(e)) - return - self._refresh_model_status() - if result['needs_initialize']: - names = "\n".join(f" - {name}" for name in result['needs_initialize']) - QMessageBox.information( - self, - "Initialize Model needed", - "The new data for these features cannot be put into the current " - f"model:\n{names}\n\n" - "Run Initialize Model to use it. Initialize Model builds all " - "features again.", - ) - - def _refresh_changed_layers(self): - names = self.data_manager.get_changed_layers() if self.data_manager else [] - if not names: - self.layerChangedWidget.hide() - return - self.layerChangedLabel.setText( - "Input layers changed after the model data was read: " + ", ".join(names) + return False + self.refresh_primary_action() + if not result['needs_initialize']: + return True + names = "\n".join(f" - {name}" for name in result['needs_initialize']) + reply = QMessageBox.question( + self, + "Rebuild needed", + "The new data for these features cannot be put into the current " + f"model:\n{names}\n\n" + "Rebuild the model to use it? A rebuild makes all features again.", ) - self.layerChangedWidget.show() + if reply == QMessageBox.StandardButton.Yes: + self.build_model() + return False def solve_model(self): - # Build/interpolate every feature already added to the model. Only - # meaningful once Initialize Model has created some features, and not - # while a fault topology edit is pending re-Initialize. + # Solve every feature already added to the model. Only meaningful + # once a build has created some features, and not while a fault + # topology edit is pending a build. if not self.model_manager or self.model_manager.model_state not in _SOLVABLE_STATES: return if not self._confirm_bounding_box_contains_data(): return + if self.data_manager is not None and self.data_manager.get_changed_layers(): + if not self._update_model_data(): + return self._run_model_task( lambda progress_callback: self.model_manager.update_all_features( progress_callback=progress_callback, notify_observers=False @@ -378,7 +494,9 @@ def _confirm_bounding_box_contains_data(self): ) return reply == QMessageBox.StandardButton.Yes - def _run_model_task(self, target, *, title, initial_label): + def _run_model_task( + self, target, *, title, initial_label, cancellable=True, on_success=None + ): """Run `target(progress_callback)` on a background QThread with a non-modal progress dialog, so the rest of QGIS stays usable. Both Initialize Model and Solve Model share this: they disable each other @@ -398,11 +516,18 @@ def _run_model_task(self, target, *, title, initial_label): progress.setWindowModality(Qt.WindowModality.NonModal) progress.setWindowTitle(title) progress.setMinimumDuration(0) - progress.canceled.connect(self._on_task_cancel_requested) + if cancellable: + progress.canceled.connect(self._on_task_cancel_requested) + else: + # the calculation of derived data cannot stop part way + progress.setCancelButton(None) progress.show() - self.initializeModelButton.setEnabled(False) - self.solveModelButton.setEnabled(False) + self._task_running = True + self._task_failed = False + self._task_on_success = on_success + self._task_cancellable = cancellable + self.primaryButton.setEnabled(False) # Only one task runs at a time (buttons are disabled above for the # duration), so it's safe to stash the per-run state needed by the @@ -461,6 +586,7 @@ def _on_task_progress(self, message, current, total): @pyqtSlot() def _on_task_finished(self): + on_success = None if self._task_failed else self._task_on_success try: # notify observers now on the GUI thread try: @@ -473,6 +599,8 @@ def _on_task_finished(self): self._debug.log_error("Error notifying observer", e) finally: self._finish_task() + if on_success is not None: + on_success() @pyqtSlot(str) def _on_task_cancelled(self, message): @@ -480,10 +608,12 @@ def _on_task_cancelled(self, message): # right after this, from the worker's `finally`) run its normal # cleanup/refresh -- the feature list will reflect whatever was # actually built before the cancellation took effect. + self._task_failed = True print(f"{self._task_title} cancelled: {message}") @pyqtSlot(str, str) def _on_task_error(self, reason, tb): + self._task_failed = True try: box = QMessageBox(self) box.setIcon(QMessageBox.Icon.Critical) @@ -497,8 +627,8 @@ def _on_task_error(self, reason, tb): self._finish_task() def _finish_task(self): - self.initializeModelButton.setEnabled(True) - self.solveModelButton.setEnabled(self.model_manager.model_state in _SOLVABLE_STATES) + self._task_running = False + self.refresh_primary_action() try: # QProgressDialog.close() emits canceled() itself (same as # clicking the Cancel button), so disconnect first -- otherwise @@ -566,7 +696,7 @@ def _status_icon(self, built): def _refresh_model_status(self): state = self.model_manager.model_state if self.model_manager is not None else 'empty' self.modelStatusLabel.setText(_MODEL_STATE_LABELS.get(state, "Model status: unknown")) - self.solveModelButton.setEnabled(state in _SOLVABLE_STATES) + self.refresh_primary_action() def on_feature_selected(self, item): feature_name = item.text(0) @@ -663,6 +793,7 @@ def delete_feature(self, item): feature_name = item.text(0) # so that Initialize Model does not build it again self.model_manager.remove_manual_foliation(feature_name) + self.model_manager.remove_parametric_fault(feature_name) # Attempt to remove from the underlying model in a few ways try: # Try model's __delitem__ if supported diff --git a/loopstructural/gui/modelling/geological_model_tab/layer_selection_table.py b/loopstructural/gui/modelling/geological_model_tab/layer_selection_table.py index 826aa8c..0a8c3a5 100644 --- a/loopstructural/gui/modelling/geological_model_tab/layer_selection_table.py +++ b/loopstructural/gui/modelling/geological_model_tab/layer_selection_table.py @@ -1,5 +1,6 @@ from qgis.core import QgsMapLayerProxyModel from qgis.gui import QgsFieldComboBox, QgsMapLayerComboBox +from qgis.PyQt.QtCore import Qt from qgis.PyQt.QtWidgets import ( QCheckBox, QComboBox, @@ -16,6 +17,7 @@ ) from loopstructural.gui.compatibility import configure_layer_combo +from loopstructural.main import constraints class LayerSelectionTable(QWidget): @@ -173,7 +175,13 @@ def add_item_row(self): def _create_type_combo(self): """Create type selection combo box.""" combo = QComboBox() - combo.addItems(["Value", "Form Line", "Orientation", "Inequality"]) + for layer_type in constraints.CONSTRAINT_TYPES: + combo.addItem(layer_type) + combo.setItemData( + combo.count() - 1, + constraints.DESCRIPTIONS[layer_type], + Qt.ItemDataRole.ToolTipRole, + ) return combo def _create_select_layer_button(self, row, type_combo): @@ -444,6 +452,114 @@ def _setup_type_specific_fields(self, layout): self._setup_inequality_fields(layout) elif self.layer_type == "Form Line": self._setup_form_line_fields(layout) + elif self.layer_type == constraints.INTERFACE: + self._setup_interface_fields(layout) + elif self.layer_type in (constraints.GRADIENT_NORMAL, constraints.TANGENT): + self._setup_vector_fields(layout) + elif self.layer_type == constraints.PAIRWISE_INEQUALITY: + self._setup_pairwise_fields(layout) + self._setup_common_fields(layout) + + def _field_combo(self, existing_key, allow_empty=False): + """Return a field combo of the layer, set to the field of an existing row.""" + combo = QgsFieldComboBox() + if allow_empty: + combo.setAllowEmptyFieldName(True) + combo.setLayer(self.layer_combo.currentLayer()) + self.layer_combo.layerChanged.connect(combo.setLayer) + if self.existing_data.get(existing_key): + combo.setField(self.existing_data[existing_key]) + return combo + + def _setup_interface_fields(self, layout): + """Setup fields for the interface type. + + The group field is optional. Without it, each feature of the layer + (for example each line) is one surface. + """ + form = QFormLayout() + self.group_field_combo = self._field_combo('group_field', allow_empty=True) + self.group_field_combo.setToolTip( + "Points with the same value in this field are on one surface. " + "If it is empty, each feature of the layer is one surface." + ) + form.addRow("Group field (optional):", self.group_field_combo) + layout.addLayout(form) + self.field_combos = {'group_field': self.group_field_combo} + + def _setup_vector_fields(self, layout): + """Setup fields for the gradient/normal and tangent types.""" + form = QFormLayout() + self.field_combos = {} + if self.layer_type == constraints.GRADIENT_NORMAL: + self.vector_kind_combo = QComboBox() + self.vector_kind_combo.addItem("Gradient", constraints.KIND_GRADIENT) + self.vector_kind_combo.addItem("Normal", constraints.KIND_NORMAL) + index = self.vector_kind_combo.findData( + self.existing_data.get('vector_kind', constraints.KIND_GRADIENT) + ) + self.vector_kind_combo.setCurrentIndex(max(index, 0)) + self.vector_kind_combo.setToolTip( + "A gradient has a direction and a size. A normal has a direction only." + ) + form.addRow("Vector:", self.vector_kind_combo) + self.field_combos['vector_kind'] = self.vector_kind_combo + for key, label in ( + ('vector_x_field', "X component:"), + ('vector_y_field', "Y component:"), + ('vector_z_field', "Z component:"), + ): + combo = self._field_combo(key) + form.addRow(label, combo) + self.field_combos[key] = combo + layout.addLayout(form) + + def _setup_pairwise_fields(self, layout): + """Setup fields for the pairwise inequality type.""" + form = QFormLayout() + self.pair_field_combo = self._field_combo('pair_field') + self.pair_field_combo.setToolTip( + "A number for each group of points. The groups are ordered by this number." + ) + form.addRow("Group number field:", self.pair_field_combo) + layout.addLayout(form) + self.field_combos = {'pair_field': self.pair_field_combo} + + def _setup_common_fields(self, layout): + """Setup the weight and the source of Z. Every type has them.""" + form = QFormLayout() + self.weight_spin = QDoubleSpinBox() + self.weight_spin.setDecimals(3) + self.weight_spin.setRange(0.001, 1000.0) + self.weight_spin.setSingleStep(0.1) + self.weight_spin.setValue(self.existing_data.get('weight', constraints.DEFAULT_WEIGHT)) + self.weight_spin.setToolTip("How much the constraint counts in the interpolation.") + form.addRow("Weight:", self.weight_spin) + + self.z_source_combo = QComboBox() + for source in constraints.Z_SOURCES: + self.z_source_combo.addItem(constraints.Z_LABELS[source], source) + self.z_source_combo.setCurrentIndex( + max(self.z_source_combo.findData(self.existing_data.get('z_source', constraints.Z_LAYER)), 0) + ) + self.z_source_combo.setToolTip( + "Where the Z of the points comes from: the geometry of the layer, the DEM " + "of the project, or one number for all points." + ) + form.addRow("Z source:", self.z_source_combo) + + self.z_value_spin = QDoubleSpinBox() + self.z_value_spin.setDecimals(3) + self.z_value_spin.setRange(-1e9, 1e9) + self.z_value_spin.setValue(float(self.existing_data.get('z_value', 0.0))) + form.addRow("Z value:", self.z_value_spin) + layout.addLayout(form) + + def update_z_value_state(): + self.z_value_spin.setEnabled(self.z_source_combo.currentData() == constraints.Z_CONSTANT) + + self.z_source_combo.currentIndexChanged.connect(update_z_value_state) + update_z_value_state() def _setup_orientation_fields(self, layout): """Setup fields for orientation data type.""" @@ -673,7 +789,11 @@ def _on_accepted(self): 'layer': self.layer_combo.currentLayer(), 'layer_name': self._unique_layer_name(self.layer_combo.currentLayer()), 'type': self.layer_type, + 'weight': self.weight_spin.value(), + 'z_source': self.z_source_combo.currentData(), } + if self.layer_data['z_source'] == constraints.Z_CONSTANT: + self.layer_data['z_value'] = self.z_value_spin.value() # Add type-specific data if self.layer_type == "Orientation": @@ -713,6 +833,25 @@ def _on_accepted(self): 'dip_weight_spin' ].value() + elif self.layer_type == constraints.INTERFACE: + # the group field is optional + group_field = self.field_combos['group_field'].currentField() + if group_field: + self.layer_data['group_field'] = group_field + + elif self.layer_type in ( + constraints.GRADIENT_NORMAL, + constraints.TANGENT, + constraints.PAIRWISE_INEQUALITY, + ): + for key, combo in self.field_combos.items(): + if key == 'vector_kind': + self.layer_data[key] = combo.currentData() + else: + self.layer_data[key] = combo.currentField() + if constraints.missing_fields(self.layer_data): + return + self.accept() def get_layer_data(self): diff --git a/loopstructural/gui/modelling/modelling_widget.py b/loopstructural/gui/modelling/modelling_widget.py index 9b57d3d..0e12f65 100644 --- a/loopstructural/gui/modelling/modelling_widget.py +++ b/loopstructural/gui/modelling/modelling_widget.py @@ -8,7 +8,7 @@ QWidget, ) -from loopstructural.gui.modelling.steps import checks +from loopstructural.gui.modelling.steps import build_plan, checks from loopstructural.gui.modelling.steps.header import DockHeader from loopstructural.gui.modelling.steps.pages import ( DataStep, @@ -90,8 +90,8 @@ def __init__( self.footer_label.setWordWrap(True) self.back_button = QPushButton("< Back", self) self.next_button = QPushButton("Next >", self) - self.back_button.clicked.connect(lambda _checked=False: self.go_to(self.current_index - 1)) - self.next_button.clicked.connect(lambda _checked=False: self.go_to(self.current_index + 1)) + self.back_button.clicked.connect(lambda _checked=False: self.go_to(self._neighbour(-1))) + self.next_button.clicked.connect(lambda _checked=False: self.go_to(self._neighbour(1))) footer = QHBoxLayout() footer.addWidget(self.footer_label, 1) footer.addWidget(self.back_button) @@ -104,9 +104,14 @@ def __init__( mainLayout.addLayout(footer) self.step_bar.currentChanged.connect(self.go_to) + self._last_checks = {} + self.model_step.tab.set_problems_provider(self.problems) + if self.data_manager is not None: + self.data_manager.add_workflow_mode_callback(self._apply_workflow_mode) self._timer = QTimer(self) self._timer.setInterval(CHECK_INTERVAL_MS) self._timer.timeout.connect(self.refresh_status) + self._apply_workflow_mode(getattr(self.data_manager, 'workflow_mode', None)) self.go_to(0) @property @@ -116,12 +121,50 @@ def current_index(self): def go_to(self, index): """Show a step. The index is limited to the steps that exist.""" index = max(0, min(index, len(self.pages) - 1)) + if index not in self._visible_indices(): + # a hidden step: show the next visible step, or the one before it + after = self._neighbour(1, index) + index = after if after != index else self._neighbour(-1, index) self.stack.setCurrentIndex(index) self.step_bar.set_current_index(index) - self.back_button.setEnabled(index > 0) - self.next_button.setEnabled(index < len(self.pages) - 1) + self.back_button.setEnabled(self._neighbour(-1, index) != index) + self.next_button.setEnabled(self._neighbour(1, index) != index) self.refresh_status() + def _visible_indices(self): + return [i for i, page in enumerate(self.pages) if self.step_bar.is_step_visible(page.key)] + + def _neighbour(self, direction, index=None): + """Return the index of the next visible step, or ``index`` if there is none.""" + index = self.current_index if index is None else index + candidates = [ + i for i in self._visible_indices() if (i - index) * direction > 0 + ] + if not candidates: + return index + return min(candidates, key=lambda i: abs(i - index)) + + def _apply_workflow_mode(self, mode): + """Show the steps that the start choice needs.""" + mode = mode or getattr(self.data_manager, 'workflow_mode', None) + visible = build_plan.steps_for_mode(mode, [page.key for page in self.pages]) + for page in self.pages: + self.step_bar.set_step_visible(page.key, page.key in visible) + self.data_step.set_workflow_mode(mode) + if self.stack.currentIndex() >= 0 and self.pages[self.stack.currentIndex()].key not in visible: + self.go_to(self.stack.currentIndex()) + else: + self.refresh_status() + + def problems(self): + """Return the problems of the steps that show, as ``(step_key, message)`` pairs.""" + visible = self._visible_indices() + return build_plan.collect_problems( + (page.key, self._last_checks[page.key]) + for i, page in enumerate(self.pages) + if i in visible and page.key in self._last_checks + ) + def show_step(self, key): """Show the step with this key, for example `checks.STEP_VIEW`.""" self.go_to([page.key for page in self.pages].index(key)) @@ -133,7 +176,9 @@ def check_all(self): def refresh_status(self): """Run the checks again and show the results in the bar and the footer.""" results = self.check_all() + self._last_checks = results self.step_bar.set_checks(results) + self.model_step.tab.refresh_primary_action() page = self.pages[self.current_index] check = results[page.key] if check.messages: diff --git a/loopstructural/gui/modelling/steps/build_plan.py b/loopstructural/gui/modelling/steps/build_plan.py new file mode 100644 index 0000000..34899a0 --- /dev/null +++ b/loopstructural/gui/modelling/steps/build_plan.py @@ -0,0 +1,111 @@ +"""The primary action of the model step. + +One button replaces "Initialize Model", "Solve Model" and "Update Model +Data". The text and the action of the button come from the state of the model +and from the derived data. This module does not import QGIS, so the unit tests +can run it. +""" + +from dataclasses import dataclass + +from loopstructural.main.workflow_mode import ( + CONSTRAINT_MODE_HIDDEN_STEPS, + WORKFLOW_MODE_CONSTRAINTS, +) + +# The actions of the primary button +ACTION_BUILD = 'build' # make the features again, then solve them +ACTION_SOLVE = 'solve' # solve the features that exist +ACTION_NONE = 'none' # nothing can be done now + +@dataclass(frozen=True) +class PrimaryAction: + """The text, the action and the enabled state of the primary button.""" + + action: str + text: str + enabled: bool = True + tooltip: str = "" + + +def choose_primary_action( + model_state, + *, + derived_out_of_date=(), + layers_changed=False, + blocked_reason=None, +) -> PrimaryAction: + """Return the action of the primary button. + + Parameters + ---------- + model_state : str + `GeologicalModelManager.model_state`: empty, initialized, solved or stale. + derived_out_of_date : iterable of str + Names of the derived results that the build calculates again. + layers_changed : bool + True if an input layer changed after the model data was read. + blocked_reason : str, optional + A problem that stops every build, for example no bounding box. + """ + if blocked_reason: + return PrimaryAction(ACTION_NONE, "Build model", False, blocked_reason) + derived_out_of_date = list(derived_out_of_date) + if model_state == 'empty': + return PrimaryAction( + ACTION_BUILD, "Build model", tooltip="Make the features and solve them." + ) + if model_state == 'stale' or derived_out_of_date: + if derived_out_of_date: + reason = "The data that comes from the column is out of date." + else: + reason = "The inputs changed after the model was built." + return PrimaryAction( + ACTION_BUILD, + "Rebuild model", + tooltip=f"{reason} Make the features again and solve them.", + ) + if layers_changed: + return PrimaryAction( + ACTION_SOLVE, + "Update data and solve", + tooltip=( + "An input layer changed. Read the layers again, put the new data into " + "the features, and solve them." + ), + ) + if model_state == 'initialized': + return PrimaryAction(ACTION_SOLVE, "Solve model", tooltip="Solve the features.") + return PrimaryAction( + ACTION_SOLVE, "Solve again", tooltip="Solve the features again with their current settings." + ) + + +def collect_problems(step_checks): + """Return the problems of all steps as ``(step_key, message)`` pairs. + + Parameters + ---------- + step_checks : iterable of (str, StepCheck) + The step keys and the results of their checks, in the order of the steps. + """ + problems = [] + seen = set() + for key, check in step_checks: + for message in check.problems: + # The same problem can show in more than one step + if message in seen: + continue + seen.add(message) + problems.append((key, message)) + return problems + + +def steps_for_mode(mode, all_steps): + """Return the step keys that the start choice shows. + + With "Interpolate surfaces from constraints", the stratigraphy and fault + steps are hidden. The user can show them later. + """ + hidden = CONSTRAINT_MODE_HIDDEN_STEPS if mode == WORKFLOW_MODE_CONSTRAINTS else () + return [key for key in all_steps if key not in hidden] diff --git a/loopstructural/gui/modelling/steps/checks.py b/loopstructural/gui/modelling/steps/checks.py index 4606291..8159e0c 100644 --- a/loopstructural/gui/modelling/steps/checks.py +++ b/loopstructural/gui/modelling/steps/checks.py @@ -6,6 +6,7 @@ """ from loopstructural.main import derived_data, layer_roles +from loopstructural.main.workflow_mode import WORKFLOW_MODE_CONSTRAINTS from .status import StepCheck @@ -106,7 +107,9 @@ def check_model(data_manager, model_manager=None) -> StepCheck: problems.append("Set the bounding box in step 1.") if not data_manager.is_model_crs_valid(): problems.append("The model CRS must be a projected CRS (in metres).") - problems.extend(_out_of_date_messages(data_manager)) + # With constraints only, there is no column, so no result comes from it + if getattr(data_manager, 'workflow_mode', None) != WORKFLOW_MODE_CONSTRAINTS: + problems.extend(_out_of_date_messages(data_manager)) changed = data_manager.get_changed_layers() if changed: problems.append("Input layers changed after the model was built: " + ", ".join(changed)) diff --git a/loopstructural/gui/modelling/steps/pages.py b/loopstructural/gui/modelling/steps/pages.py index 7c7ec72..847fc84 100644 --- a/loopstructural/gui/modelling/steps/pages.py +++ b/loopstructural/gui/modelling/steps/pages.py @@ -8,6 +8,7 @@ from qgis.gui import QgsCollapsibleGroupBox from qgis.PyQt.QtCore import Qt, pyqtSignal from qgis.PyQt.QtWidgets import ( + QComboBox, QHBoxLayout, QLabel, QMenu, @@ -24,6 +25,8 @@ from loopstructural.gui.modelling.model_definition import ModelDefinitionTab from loopstructural.gui.modelling.model_definition.fault_layers import FaultLayersWidget +from loopstructural.main.workflow_mode import WORKFLOW_MODE_LABELS, WORKFLOW_MODES + from . import checks @@ -76,10 +79,36 @@ def __init__(self, parent=None, **kwargs): row.addWidget(QLabel("Select the area and the source layers.", self), 1) row.addWidget(convert) layout.addLayout(row) + + # The start choice. The "constraints" choice hides steps 2 and 3. + self.mode_combo = QComboBox(self) + for mode in WORKFLOW_MODES: + self.mode_combo.addItem(WORKFLOW_MODE_LABELS[mode], mode) + self.mode_combo.setToolTip( + "Steps 2 and 3 make features and constraints from a geological map. " + "Interpolate from constraints hides them. Change the choice to show them again." + ) + self.mode_combo.currentIndexChanged.connect(self._on_mode_changed) + mode_row = QHBoxLayout() + mode_row.addWidget(QLabel("Method:", self)) + mode_row.addWidget(self.mode_combo, 1) + layout.addLayout(mode_row) self.tab = ModelDefinitionTab(self, data_manager=self.data_manager) layout.addWidget(self.tab, 1) + def set_workflow_mode(self, mode): + index = self.mode_combo.findData(mode) + if index >= 0 and index != self.mode_combo.currentIndex(): + self.mode_combo.blockSignals(True) + self.mode_combo.setCurrentIndex(index) + self.mode_combo.blockSignals(False) + + def _on_mode_changed(self, _index): + if self.data_manager is not None: + self.data_manager.set_workflow_mode(self.mode_combo.currentData()) + + class StratigraphyStep(StepPage): """Step 2: the stratigraphic column, and the results that come from the map.""" diff --git a/loopstructural/gui/modelling/steps/step_bar.py b/loopstructural/gui/modelling/steps/step_bar.py index e8d96c0..811ddf1 100644 --- a/loopstructural/gui/modelling/steps/step_bar.py +++ b/loopstructural/gui/modelling/steps/step_bar.py @@ -79,6 +79,14 @@ def __init__(self, steps, parent=None): def count(self): return len(self._buttons) + def set_step_visible(self, key, visible): + """Show or hide the button of a step.""" + self._buttons[self._keys.index(key)].setVisible(visible) + + def is_step_visible(self, key): + # `isVisibleTo` does not depend on the visibility of the parent + return self._buttons[self._keys.index(key)].isVisibleTo(self) + def current_index(self): return self._group.checkedId() diff --git a/loopstructural/main/constraints.py b/loopstructural/main/constraints.py new file mode 100644 index 0000000..951d317 --- /dev/null +++ b/loopstructural/main/constraints.py @@ -0,0 +1,186 @@ +"""The constraint types of a feature that the user builds from layers. + +Each row of the constraint list of a feature is a dictionary (see +`LayerSelectionTable`). This module knows the types, the keys that each type +needs, and how to turn the points of a layer into the data frame that +LoopStructural reads. All interpolators (FDI, PLI and surfe) take the same +types. LoopStructural converts them when an interpolator needs another form. + +This module does not import QGIS, so the unit tests can run it. +""" + +import numpy as np +import pandas as pd + +VALUE = 'Value' +INTERFACE = 'Interface' +ORIENTATION = 'Orientation' +GRADIENT_NORMAL = 'Gradient/Normal' +TANGENT = 'Tangent' +FORM_LINE = 'Form Line' +INEQUALITY = 'Inequality' +PAIRWISE_INEQUALITY = 'Pairwise Inequality' + +# In the order of the type list of the user interface +CONSTRAINT_TYPES = ( + VALUE, + INTERFACE, + ORIENTATION, + GRADIENT_NORMAL, + TANGENT, + FORM_LINE, + INEQUALITY, + PAIRWISE_INEQUALITY, +) + +DESCRIPTIONS = { + VALUE: "The scalar field has a known value at the points.", + INTERFACE: "The points are on one surface. The value of the surface is not known.", + ORIENTATION: "Strike and dip of the surface at the points.", + GRADIENT_NORMAL: "A vector at the points, from three fields (x, y, z). A gradient has " + "a direction and a size. A normal has a direction only.", + TANGENT: "A vector in the surface at the points, from three fields (x, y, z).", + FORM_LINE: "Lines on the surface. The line gives a constant value or the strike.", + INEQUALITY: "The value at the points is above a lower limit and below an upper limit.", + PAIRWISE_INEQUALITY: "Groups of points. The values of one group are ordered against the " + "values of the other groups.", +} + +Z_LAYER = 'layer' +Z_DEM = 'dem' +Z_CONSTANT = 'constant' +Z_SOURCES = (Z_LAYER, Z_DEM, Z_CONSTANT) +Z_LABELS = { + Z_LAYER: "Z of the layer", + Z_DEM: "DEM", + Z_CONSTANT: "Constant", +} + +KIND_GRADIENT = 'gradient' +KIND_NORMAL = 'normal' + +DEFAULT_WEIGHT = 1.0 + +# The keys of a row that name a field of the layer, for each type. A row +# without these keys is not complete. +REQUIRED_FIELDS = { + VALUE: ('value_field',), + INTERFACE: (), + ORIENTATION: ('strike_field', 'dip_field'), + GRADIENT_NORMAL: ('vector_x_field', 'vector_y_field', 'vector_z_field'), + TANGENT: ('vector_x_field', 'vector_y_field', 'vector_z_field'), + FORM_LINE: (), + INEQUALITY: ('lower_field', 'upper_field'), + PAIRWISE_INEQUALITY: ('pair_field',), +} + + +def missing_fields(layer_data): + """Return the keys that a row needs and does not have.""" + return [key for key in REQUIRED_FIELDS.get(layer_data.get('type'), ()) if not layer_data.get(key)] + + +def solver_for(layer_type): + """Return the solver that a type needs, or None for the default. + + Inequality constraints need the ADMM solver. + """ + return 'admm' if layer_type in (INEQUALITY, PAIRWISE_INEQUALITY) else None + + +def sample_layer(sampler, layer_data, dem_function, default_use_z=False): + """Return the points of a row, with the Z that the row asks for. + + Parameters + ---------- + sampler : callable + Gives a data frame of points: ``sampler(df, dem_function, use_z)``. + layer_data : dict + The row. ``z_source`` is ``layer``, ``dem`` or ``constant``. For a + constant, ``z_value`` is the Z. A row with no ``z_source`` (made by an + older version) uses ``default_use_z``. + """ + source = layer_data.get('z_source') + use_z = default_use_z if source is None else source == Z_LAYER + points = sampler(layer_data['df'], dem_function, use_z) + if source == Z_CONSTANT: + points = points.copy() + points['Z'] = float(layer_data.get('z_value', 0.0)) + return points + + +def add_weight(rows, points, layer_data): + """Add the weight of the row to ``rows``, a slice of ``points``. + + Rows that already have a weight keep it. Without a weight in the row of the + constraint list, LoopStructural uses a weight of 1. + """ + weight = layer_data.get('weight') + if weight is None or 'w' in rows: + return rows + rows = rows.copy() + rows['w'] = float(weight) + return rows + + +def _vector_rows(points, layer_data, names): + fields = [layer_data[key] for key in ('vector_x_field', 'vector_y_field', 'vector_z_field')] + rows = points[['X', 'Y', 'Z']].copy() + for name, field in zip(names, fields): + rows[name] = pd.to_numeric(points[field], errors='coerce') + rows = rows.dropna(subset=list(names)) + # A zero vector has no direction + zero = (rows[list(names)] == 0).all(axis=1) + return rows[~zero] + + +def interface_rows(points, layer_data, feature_name, offset=0): + """Return the interface rows, and the offset for the next interface row. + + With a group field, the points that have the same value are on one + surface. Without it, each feature of the layer (for example each line) is + one surface. + """ + field = layer_data.get('group_field') + if field: + codes, _ = pd.factorize(points[field]) + points = points[codes >= 0].copy() + interface = codes[codes >= 0].astype(float) + else: + interface = points['feature_id'].to_numpy(float) + count = len(np.unique(interface)) + rows = points[['X', 'Y', 'Z']].copy() + # Different groups need different numbers, also between the rows of a feature + _, interface = np.unique(interface, return_inverse=True) + rows['interface'] = interface.astype(float) + offset + rows['feature_name'] = feature_name + return rows, offset + count + + +def constraint_rows(points, layer_data, feature_name): + """Return the data frame rows of the types that this module builds. + + Types: gradient/normal, tangent, pairwise inequality. Value, orientation, + form line and inequality are in `ModelManager._foliation_data`. Interface + rows come from `interface_rows`. + + Returns + ------- + pandas.DataFrame + The rows, with the weight of the row of the constraint list. + """ + layer_type = layer_data['type'] + if layer_type == GRADIENT_NORMAL: + kind = layer_data.get('vector_kind', KIND_GRADIENT) + names = ('nx', 'ny', 'nz') if kind == KIND_NORMAL else ('gx', 'gy', 'gz') + rows = _vector_rows(points, layer_data, names) + elif layer_type == TANGENT: + rows = _vector_rows(points, layer_data, ('tx', 'ty', 'tz')) + elif layer_type == PAIRWISE_INEQUALITY: + rows = points[['X', 'Y', 'Z']].copy() + rows['pair_id'] = pd.to_numeric(points[layer_data['pair_field']], errors='coerce') + rows = rows.dropna(subset=['pair_id']) + else: + raise ValueError(f"Unknown layer type: {layer_type}") + rows['feature_name'] = feature_name + return add_weight(rows, points, layer_data) diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py index cea7c0d..fceda2d 100644 --- a/loopstructural/main/data_manager.py +++ b/loopstructural/main/data_manager.py @@ -33,6 +33,7 @@ from .layer_roles import LayerRoles from .m2l_api import paint_stratigraphic_order from .vectorLayerWrapper import qgsLayerToGeoDataFrame +from .workflow_mode import WORKFLOW_MODE_MAP, WORKFLOW_MODES def _lookup_colour_ramp(ramp_name): @@ -169,6 +170,10 @@ def __init__(self, *, project=None, mapCanvas=None, logger=None): # For each derived result, the inputs of the last run. See # `derived_data` for how an out-of-date result is found. self.derived = DerivedData() + # How the user builds the model: from a geological map, or by + # interpolating surfaces from constraints. See `gui.modelling.steps.build_plan`. + self.workflow_mode = WORKFLOW_MODE_MAP + self._workflow_mode_callbacks = [] self.derived.register(derived_data.BASAL_CONTACTS, self.basal_contacts_inputs) self.derived.register(derived_data.THICKNESS, self.thickness_inputs) self.derived.register(derived_data.STYLED_FIELDS, self.styled_fields_inputs) @@ -301,6 +306,20 @@ def set_debug_manager(self, debug_manager): """Set the debug manager for the tools that the data manager can start.""" self.debug_manager = debug_manager + def set_workflow_mode(self, mode): + """Set the start choice of the user and tell the listeners if it changed.""" + if mode not in WORKFLOW_MODES: + raise ValueError(f"Unknown workflow mode '{mode}'.") + if mode != self.workflow_mode: + self.workflow_mode = mode + for callback in list(self._workflow_mode_callbacks): + callback(mode) + + def add_workflow_mode_callback(self, callback): + """Call ``callback(mode)`` each time the start choice changes.""" + if callback not in self._workflow_mode_callbacks: + self._workflow_mode_callbacks.append(callback) + def _on_layer_role_changed(self, role, value): """A layer role changed: the derived results can be out of date.""" self.derived.refresh() @@ -1326,6 +1345,35 @@ def get_changed_layers(self): model data was last read from them.""" return sorted(self._watched_layers[i][0].name() for i in self._changed_layer_ids) + def sync_extra_constraints(self): + """Read the layers of the constraints that the user added to generated features. + + A generated feature is one that the stratigraphic column makes. Its + table of data layers can have rows of the user. The model manager adds + them to the data of the column when it builds the feature. A feature + that the user added (a manual foliation) has its own build, so it is + not in this list. + """ + if self._model_manager is None: + return + manual = self._model_manager.manual_foliations + generated = set(self._model_manager.generated_feature_names()) + model_crs = self.get_model_crs() + extras = {} + for name, entries in self.feature_data.items(): + if name in manual or name not in generated: + continue + rows = {} + for key, entry in entries.items(): + if entry.get('processed') or entry.get('layer') is None: + continue + row = dict(entry) + row['df'] = qgsLayerToGeoDataFrame(entry['layer'], target_crs=model_crs) + rows[key] = row + if rows: + extras[name] = {'data': rows, 'use_z_coordinate': True} + self._model_manager.extra_constraints = extras + def reload_changed_layers(self): """Read the model data again from the input layers that changed. @@ -1371,6 +1419,7 @@ def refresh_model_data(self): if self._model_manager is None: raise RuntimeError("Model manager is not set.") self.reload_changed_layers() + self.sync_extra_constraints() return self._model_manager.refresh_feature_data() def get_layers_outside_bounding_box(self): @@ -1505,12 +1554,28 @@ def get_stratigraphic_column(self): def clear_stratigraphic_column(self): self._stratigraphic_column.clear() - def update_stratigraphy(self): - """Update the foliation features in the model manager.""" + def update_stratigraphy(self, basal_contacts=None, unit_name_field=None): + """Update the foliation features in the model manager. + + Parameters + ---------- + basal_contacts : geopandas.GeoDataFrame, optional + Contacts that the model uses in place of the basal contacts layer. + A build with "Calculate from geology polygons" gives the extracted + contacts here, so that the model does not read the project layer. + unit_name_field : str, optional + The unit name field of ``basal_contacts``. + """ self.logger(message="Updating stratigraphy...", log_level=4) if self._model_manager is not None: model_crs = self.get_model_crs() - if self._basal_contacts is not None: + if basal_contacts is not None: + self._model_manager.update_contact_traces( + basal_contacts, + unit_name_field=unit_name_field, + use_z_coordinate=False, + ) + elif self._basal_contacts is not None: self._model_manager.update_contact_traces( qgsLayerToGeoDataFrame(self._basal_contacts['layer'], target_crs=model_crs), unit_name_field=self._basal_contacts['unitname_field'], @@ -1599,6 +1664,7 @@ def _sync_processed_feature_data(self): unit_data = self._model_manager.get_stratigraphy_entry(unit.name) if not unit_data: continue + suffix = 'detached' if self._model_manager.is_detached(group.name) else 'auto' contact = unit_data.get('contact') if ( contact is not None @@ -1608,8 +1674,9 @@ def _sync_processed_feature_data(self): self._add_processed_feature_row( group.name, self._basal_contacts.get('layer'), - 'Contact (auto)', + f'Contact ({suffix})', unit.name, + suffix, ) orientations = unit_data.get('orientations') if ( @@ -1620,8 +1687,9 @@ def _sync_processed_feature_data(self): self._add_processed_feature_row( group.name, self._structural_orientations.get('layer'), - 'Orientation (auto)', + f'Orientation ({suffix})', unit.name, + suffix, ) if self._fault_traces is not None: @@ -1635,7 +1703,9 @@ def _sync_processed_feature_data(self): fault_name, ) - def _add_processed_feature_row(self, feature_name, layer, type_label, source_name): + def _add_processed_feature_row( + self, feature_name, layer, type_label, source_name, suffix='auto' + ): """Add a single read-only, workflow-derived row to `feature_data`. Keyed on a string distinct from a plain layer name so a processed row @@ -1644,7 +1714,7 @@ def _add_processed_feature_row(self, feature_name, layer, type_label, source_nam """ if layer is None: return - display_name = f"{source_name} ({layer.name()}, auto)" + display_name = f"{source_name} ({layer.name()}, {suffix})" self.feature_data[feature_name][display_name] = { 'layer': layer, 'layer_name': display_name, @@ -1701,6 +1771,7 @@ def reset(self): self.layer_roles.clear() self.thickness_sources.clear() self.derived.clear() + self.set_workflow_mode(WORKFLOW_MODE_MAP) self.set_dem_layer(None) self.use_dem = True @@ -1741,6 +1812,8 @@ def save_state(self, filepath): self._model_manager.save_model(str(model_path)) state['model_file'] = model_path.name state['manual_foliations'] = self._model_manager.manual_foliations_to_dict() + state['detached_features'] = self._model_manager.detached_to_dict() + state['parametric_faults'] = self._model_manager.parametric_faults_to_dict() with open(path, 'w') as f: json.dump(state, f, indent=2) @@ -1773,6 +1846,8 @@ def load_state(self, filepath): # the pickled model already has these features; this lets # Initialize Model build them again self._model_manager.manual_foliations_from_dict(state.get('manual_foliations', {})) + self._model_manager.detached_from_dict(state.get('detached_features', {})) + self._model_manager.parametric_faults_from_dict(state.get('parametric_faults', {})) # the data was just read from the layers self._changed_layer_ids.clear() self.refresh_layer_watchers() @@ -1856,6 +1931,7 @@ def to_dict(self): 'use_project_crs': self._use_project_crs, 'layer_roles': self.layer_roles.to_dict(), 'derived_data': self.derived.to_dict(), + 'workflow_mode': self.workflow_mode, 'thickness_sources': self.thickness_sources.to_dict(), } @@ -1876,6 +1952,9 @@ def _restore_derived_state(self, data): self.layer_roles.contacts_source = layer_roles.CONTACTS_FROM_GEOLOGY self.thickness_sources.from_dict(data.get('thickness_sources')) self.derived.from_dict(data.get('derived_data')) + # A state file of an older version has no start choice + mode = data.get('workflow_mode') + self.set_workflow_mode(mode if mode in WORKFLOW_MODES else WORKFLOW_MODE_MAP) def from_dict(self, data): """Load data from a dictionary.""" diff --git a/loopstructural/main/derived_refresh.py b/loopstructural/main/derived_refresh.py new file mode 100644 index 0000000..66dde38 --- /dev/null +++ b/loopstructural/main/derived_refresh.py @@ -0,0 +1,316 @@ +"""Calculate the out-of-date derived data again before a model build. + +A model build must not use out-of-date basal contacts or thicknesses. This +module finds the results to calculate, reads their inputs on the GUI thread, +calculates them (on a background thread) and records the new inputs (on the +GUI thread). + +The order is fixed: the basal contacts first, then the thicknesses, because +the thickness calculation reads the contacts. + +With "Calculate from geology polygons", the contacts go directly to the model. +The project layer of the contacts is updated afterwards, for display only. +""" + +from typing import Callable, List, Optional + +from . import derived_data, layer_roles + +BASAL_CONTACT_UNIT_FIELD = 'basal_unit' + + +class DerivedRefreshError(RuntimeError): + """A derived result could not be calculated. The build must stop.""" + + +def names_to_refresh(data_manager) -> List[str]: + """Return the derived results that a build calculates, in the order of calculation. + + - Basal contacts: with "Calculate from geology polygons", when they are + out of date, or never calculated. With "Use a contacts layer", never, + because the layer is an input of the user. + - Thickness: only when it is out of date. A thickness that was never + calculated can be typed, so the build does not calculate it. + """ + derived = data_manager.derived + names = [] + if ( + data_manager.layer_roles.contacts_source == layer_roles.CONTACTS_FROM_GEOLOGY + and data_manager.get_stratigraphic_unit_names() + and data_manager.get_layer_role(layer_roles.GEOLOGY) is not None + and data_manager.get_layer_role(layer_roles.GEOLOGY_UNIT_FIELD) + and derived.status(derived_data.BASAL_CONTACTS) != derived_data.CURRENT + ): + names.append(derived_data.BASAL_CONTACTS) + if derived.is_out_of_date(derived_data.THICKNESS): + names.append(derived_data.THICKNESS) + return names + + +def _check_contacts_inputs(inputs): + if not inputs.get('geology') or not inputs.get('unit_field'): + raise DerivedRefreshError( + "The basal contacts cannot be calculated: select the geology layer and the " + "unit name field in step 1." + ) + + +class DerivedRefresh: + """One refresh of the out-of-date derived data. + + Make the object on the GUI thread (it reads the layers and the settings). + Call `run` on the background thread. Call `finish` on the GUI thread. + """ + + def __init__(self, data_manager, model_manager, names=None): + self.data_manager = data_manager + self.model_manager = model_manager + self.names = list(names if names is not None else names_to_refresh(data_manager)) + self._contacts = None + self._thickness = None + self._contacts_result = None + self._thickness_result = None + if derived_data.BASAL_CONTACTS in self.names: + self._contacts = self._read_contacts_job() + if derived_data.THICKNESS in self.names: + self._thickness = self._read_thickness_job() + + def __bool__(self): + return bool(self.names) + + # -- the reads on the GUI thread ------------------------------------- + + def _read_contacts_job(self): + dm = self.data_manager + settings = dm.get_widget_settings('basal_contacts_widget', {}) or {} + geology = dm.get_layer_role(layer_roles.GEOLOGY) + faults = ( + dm.find_layer_by_name(settings['faults_layer']) + if settings.get('faults_layer') + else dm.get_layer_role(layer_roles.FAULT_TRACES) + ) + unit_field = dm.get_layer_role(layer_roles.GEOLOGY_UNIT_FIELD) + inputs = dm.basal_contacts_inputs( + geology=geology, + unit_field=unit_field, + faults=faults, + ignore_units=settings.get('ignore_units', []), + override_units=settings.get('basal_override_units', []), + ) + _check_contacts_inputs(inputs) + target_crs = dm.get_model_crs() + if target_crs is None or not target_crs.isValid(): + target_crs = geology.crs() + return { + 'inputs': inputs, + 'params': dict( + geology=geology, + stratigraphic_order=dm.get_stratigraphic_unit_names(), + faults=faults, + ignore_units=list(settings.get('ignore_units', [])), + unit_name_field=unit_field, + target_crs=target_crs, + unit_colours=dm.get_stratigraphic_unit_colours(), + basal_override_units=list(settings.get('basal_override_units', [])), + debug_manager=getattr(dm, 'debug_manager', None), + ), + } + + def _read_thickness_job(self): + dm = self.data_manager + settings = dm.get_widget_settings('thickness_calculator_widget', {}) or {} + calculator_type = settings.get('calculator_type') + if not calculator_type: + raise DerivedRefreshError( + "The thicknesses are out of date, but the thickness calculator has no settings. " + "Run Calculate thickness in step 2 one time." + ) + + def layer(key, role=None): + name = settings.get(key) + if name: + return dm.find_layer_by_name(name) + return dm.get_layer_role(role) if role else None + + geology = layer('geology_layer', layer_roles.GEOLOGY) + structure = layer('structure_layer', layer_roles.STRUCTURE) + cross_sections = layer('cross_sections_layer') + unit_field = settings.get('unit_name_field') or dm.get_layer_role( + layer_roles.GEOLOGY_UNIT_FIELD + ) + if geology is None or not unit_field: + raise DerivedRefreshError( + "The thicknesses cannot be calculated: select the geology layer and the " + "unit name field in step 1." + ) + if calculator_type == 'StructuralPoint' and structure is None: + raise DerivedRefreshError( + "The thicknesses cannot be calculated: select the structure layer in step 1." + ) + if calculator_type == 'AlongSection' and cross_sections is None: + raise DerivedRefreshError( + "The thicknesses cannot be calculated: select the cross-sections layer in " + "the thickness calculator." + ) + contacts_layer = None + if dm.layer_roles.contacts_source == layer_roles.CONTACTS_FROM_LAYER: + contacts_layer = dm.get_layer_role(layer_roles.BASAL_CONTACTS) + elif settings.get('basal_contacts_layer'): + contacts_layer = dm.find_layer_by_name(settings['basal_contacts_layer']) + inputs = dm.thickness_inputs( + geology=geology, + unit_field=unit_field, + contacts_layer=contacts_layer, + calculator_type=calculator_type, + structure=structure, + cross_sections=cross_sections, + ) + orientation_types = ("Dip Direction", "Strike") + index = settings.get('orientation_type_index', 0) + orientation_type = orientation_types[index] if 0 <= index < len(orientation_types) else None + params = dict( + calculator_type=calculator_type, + dtm=layer('dtm_layer', layer_roles.DEM), + geology=geology, + basal_contacts=contacts_layer, + sampling_frequency=settings.get('sampling_frequency', 200), + structure=structure, + cross_sections=cross_sections, + unit_name_field=unit_field, + dip_field=settings.get('dip_field') or 'DIP', + dipdir_field=settings.get('dipdir_field') or 'DIPDIR', + basal_contacts_unit_name=settings.get('basal_unit_field'), + max_line_length=( + settings.get('max_line_length') + if calculator_type in ("StructuralPoint", "InterpolatedStructure") + else None + ), + stratigraphic_order=dm.get_stratigraphic_unit_names(), + debug_manager=getattr(dm, 'debug_manager', None), + ) + if orientation_type is not None: + params['orientation_type'] = orientation_type + return {'inputs': inputs, 'params': params} + + # -- the calculation on the background thread ------------------------ + + def run(self, progress: Optional[Callable[[str, int, int], None]] = None): + """Calculate the results. Raise `DerivedRefreshError` if one fails.""" + total = len(self.names) + + def report(message, current): + if progress is not None: + progress(message, current, total) + + for index, name in enumerate(self.names): + label = derived_data.DESCRIPTIONS[name].lower() + report(f"Calculating the {label} ({index + 1} of {total})...", index) + + def updater(message, _label=label, _index=index): + report(f"Calculating the {_label}: {message}", _index) + + try: + if name == derived_data.BASAL_CONTACTS: + self._run_contacts(updater) + elif name == derived_data.THICKNESS: + self._run_thickness(updater) + except DerivedRefreshError: + raise + except Exception as err: + raise DerivedRefreshError( + f"The {label} could not be calculated: {type(err).__name__}: {err}" + ) from err + report("Derived data is up to date.", total) + + def _run_contacts(self, updater): + from .m2l_api import extract_basal_contacts + + job = self._contacts + result = extract_basal_contacts(updater=updater, **job['params']) + contacts = result['basal_contacts'] + if contacts is None or contacts.empty: + raise DerivedRefreshError( + "No basal contacts were found with the geology layer and the stratigraphic " + "column. Check the unit name field and the order of the column." + ) + self._contacts_result = contacts + # The model reads the contacts directly. The project layer is only for display. + self.data_manager.update_stratigraphy( + basal_contacts=contacts, unit_name_field=BASAL_CONTACT_UNIT_FIELD + ) + + def _run_thickness(self, updater): + from .m2l_api import calculate_thickness + + result = calculate_thickness(updater=updater, **self._thickness['params']) + if not isinstance(result, dict) or result.get('thicknesses') is None: + raise DerivedRefreshError("The thickness calculation gave no result.") + self._thickness_result = result + + # -- the changes on the GUI thread ---------------------------------- + + def finish(self): + """Put the results into the data manager and record their inputs. + + Returns + ------- + list of str + The names of the units that keep a thickness that the user typed. + """ + dm = self.data_manager + skipped = [] + if self._contacts_result is not None: + self._update_contacts_layer(self._contacts_result) + dm.derived.record(derived_data.BASAL_CONTACTS, inputs=self._contacts['inputs']) + if self._thickness_result is not None: + values = thickness_values(self._thickness_result['thicknesses']) + _, skipped = dm.apply_calculated_thicknesses(values) + dm.derived.record(derived_data.THICKNESS, inputs=self._thickness['inputs']) + return skipped + + def _update_contacts_layer(self, contacts): + """Show the contacts in the project. The model does not read this layer.""" + from qgis.core import QgsProject + + from .vectorLayerWrapper import addGeoDataFrameToproject + + dm = self.data_manager + old = dm.get_layer_role(layer_roles.BASAL_CONTACTS) + layer = addGeoDataFrameToproject(contacts, "Basal contacts") + dm.apply_stratigraphic_colours_to_layer(layer, BASAL_CONTACT_UNIT_FIELD) + dm.set_basal_contacts(layer, unitname_field=BASAL_CONTACT_UNIT_FIELD) + # An earlier layer of the plugin is replaced. A layer of the user is not removed. + try: + if ( + old is not None + and old.id() != layer.id() + and old.dataProvider().name() == 'memory' + and old.name() == "Basal contacts" + ): + QgsProject.instance().removeMapLayer(old.id()) + except RuntimeError: + pass + + +def thickness_values(thicknesses): + """Return unit name -> thickness from the table of the thickness calculator. + + The median is used if it is there, then the mean. A unit with no result + (map2loop uses -1) is left out. + """ + columns = getattr(thicknesses, 'columns', []) + if 'ThicknessMedian' in columns: + column = 'ThicknessMedian' + elif 'ThicknessMean' in columns: + column = 'ThicknessMean' + else: + return {} + values = {} + for _, row in thicknesses.iterrows(): + name = row.get('name') or row.get('UNITNAME') + if not name: + continue + value = row.get(column) + if derived_data.thickness_is_set(value): + values[name] = float(value) + return values diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 903b81c..3a79f3b 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -29,6 +29,7 @@ from LoopStructural import GeologicalModel from loopstructural.toolbelt.preferences import PlgSettingsStructure +from ..main import constraints, parametric_fault from ..main.data_types import FaultEntry, StratigraphyEntry from ..main.helpers import qgisAttributeIsNone @@ -193,6 +194,16 @@ def __init__(self, debug_manager=None): # features, so it builds these again from here; see # `_build_manual_foliations`. Insertion order is the build order. self.manual_foliations: Dict[str, dict] = {} + # Features that the stratigraphic column makes (one for each group of + # units). `detached` keeps a copy of the data of a feature that the + # user detached: the data no longer follows the column and the + # contacts. `extra_constraints` are the rows that the user added to a + # generated feature, in the form of `manual_foliations[name]`. + self.detached: Dict[str, pd.DataFrame] = {} + # Faults that the user gave by numbers (see `parametric_fault`), by name. + # They have no fault trace, so they are not in `faults`. + self.parametric_faults: Dict[str, dict] = {} + self.extra_constraints: Dict[str, dict] = {} # Observers managed by Observable base class self.dem_function = lambda x, y: 0 # internal flag to temporarily suppress notifications (used when @@ -271,6 +282,9 @@ def reset(self): self.faults = defaultdict(dict) self.stratigraphy = defaultdict(dict) self.manual_foliations = {} + self.detached = {} + self.extra_constraints = {} + self.parametric_faults = {} self.dem_function = lambda x, y: 0 self._topology_dirty = False self._data_dirty = False @@ -1020,6 +1034,7 @@ def update_foliation_features(self): groupname, data=data, force_constrained=True, + **self._extra_kwargs(groupname), nelements=PlgSettingsStructure.interpolator_nelements, npw=PlgSettingsStructure.interpolator_npw, cpw=PlgSettingsStructure.interpolator_cpw, @@ -1035,12 +1050,32 @@ def update_foliation_features(self): # foliation features were rebuilt; let observers know self._emit('foliation_features_updated') - def _group_data(self, group, isovalues) -> Optional[pd.DataFrame]: + def _group_data( + self, group, isovalues, include_extra=True, use_detached=True + ) -> Optional[pd.DataFrame]: """Return the contact and orientation data of the units in `group` as one data frame for its foliation, or None if there is no data. - `isovalues` is `stratigraphic_column.get_isovalues()`. + `isovalues` is `stratigraphic_column.get_isovalues()`. A detached + group gives its kept data, not the data of the column. With + `include_extra`, the rows that the user added to the feature are in the + result. With `use_detached` False, the data of the column is used also + for a detached group. """ + if use_detached and group.name in self.detached: + data = [self.detached[group.name].copy()] + else: + data = self._column_group_data(group, isovalues) + extra = self._extra_data(group.name) if include_extra else None + if extra is not None: + data.append(extra) + if len(data) == 0: + return None + return pd.concat(data, ignore_index=True) + + def _column_group_data(self, group, isovalues) -> list: + """Return the data frames of the contacts and the orientations of the + units in `group`.""" data = [] for u in group.units: val = isovalues[u.name]['value'] @@ -1059,9 +1094,108 @@ def _group_data(self, group, isovalues) -> Optional[pd.DataFrame]: orientations['val'] = np.nan orientations['feature_name'] = group.name data.append(orientations) - if len(data) == 0: + return data + + def _extra_data(self, name) -> Optional[pd.DataFrame]: + """Return the rows that the user added to the generated feature `name`.""" + spec = self.extra_constraints.get(name) + if not spec or not spec.get('data'): return None - return pd.concat(data, ignore_index=True) + try: + data, _kwargs = self._foliation_data( + name, spec['data'], AllSampler(), spec.get('use_z_coordinate', True) + ) + except Exception as e: + if self._debug_manager is not None: + self._debug_manager.log( + f"Could not read the added constraints of '{name}': {e}", log_level=2 + ) + return None + return data + + def _extra_kwargs(self, name) -> dict: + """Return the `create_and_add_foliation` arguments that the rows added to `name` need.""" + spec = self.extra_constraints.get(name) + if not spec or not spec.get('data'): + return {} + kwargs = {} + for layer_data in spec['data'].values(): + solver = constraints.solver_for(layer_data.get('type')) + if solver: + kwargs['solver'] = solver + return kwargs + + # -- detach a generated feature --------------------------------------- + + def generated_feature_names(self): + """Return the names of the features that the stratigraphic column makes.""" + if self.stratigraphic_column is None: + return [] + return [group.name for group in self.stratigraphic_column.get_groups()] + + def is_generated(self, name) -> bool: + return name in self.generated_feature_names() + + def is_detached(self, name) -> bool: + return name in self.detached + + def detach_feature(self, name) -> bool: + """Keep the data of a generated feature, so that the user can edit it. + + The data stays the same as it is now. It no longer changes when the + column, the contacts or the orientations change. + + Returns + ------- + bool + True if the feature is detached. False if it is not a generated + feature, or it has no data. + """ + group = self._group_by_name(name) + if group is None: + return False + data = self._group_data( + group, self.stratigraphic_column.get_isovalues(), include_extra=False + ) + if data is None: + return False + self.detached[name] = data.copy() + self._emit('model_updated') + return True + + def attach_feature(self, name) -> bool: + """Make a detached feature follow the column again. + + The model needs a build to use the data of the column. + """ + if self.detached.pop(name, None) is None: + return False + self._data_dirty = True + self._emit('model_updated') + return True + + def _group_by_name(self, name): + if self.stratigraphic_column is None: + return None + for group in self.stratigraphic_column.get_groups(): + if group.name == name: + return group + return None + + def detached_to_dict(self) -> dict: + """Return `detached` in a form that `json.dump` can write.""" + return { + name: { + 'columns': {str(c): [_json_value(v) for v in data[c]] for c in data.columns} + } + for name, data in self.detached.items() + } + + def detached_from_dict(self, detached: dict): + """Replace `detached` with the data from `detached_to_dict`.""" + self.detached = { + name: pd.DataFrame(entry['columns']) for name, entry in (detached or {}).items() + } def _strip_spurious_regions_from_domain_faults(self): """Work around a LoopStructural core gap that corrupts a domain @@ -1252,8 +1386,68 @@ def update_fault_features(self): cpw=PlgSettingsStructure.interpolator_cpw, regularisation=PlgSettingsStructure.interpolator_regularisation, ) + self._build_parametric_faults() self.apply_fault_abutting_relationships() + # -- faults that the user gives by numbers --------------------------- + + def used_names(self) -> set: + """Return the names that a new feature or fault cannot use.""" + names = {f.name for f in self.features()} + names |= set(self.faults) | set(self.parametric_faults) | set(self.manual_foliations) + return names + + def add_parametric_fault(self, spec: dict): + """Add a fault that is given by a centre, a strike, a dip and a size. + + The fault is built the next time the model is built, together with the + other faults, so that the features of the column use it. Until then + the model is stale. + + Raises + ------ + ValueError + If the spec has a problem, see `parametric_fault.problems`. + """ + spec = parametric_fault.clean_spec(spec) + found = parametric_fault.problems(spec, self.used_names()) + if found: + raise ValueError(" ".join(found)) + self.parametric_faults[spec['name']] = spec + self._data_dirty = True + self._emit('model_updated') + + def remove_parametric_fault(self, name: str): + """Stop building the parametric fault `name`. The model is stale until the next build.""" + if self.parametric_faults.pop(name, None) is not None: + self._data_dirty = True + + def _build_parametric_faults(self): + for name, spec in self.parametric_faults.items(): + self._report_progress(f"Building fault '{name}'") + try: + self.model.create_and_add_fault( + name, + data=parametric_fault.frame_data(spec), + nelements=PlgSettingsStructure.interpolator_nelements, + npw=PlgSettingsStructure.interpolator_npw, + cpw=PlgSettingsStructure.interpolator_cpw, + regularisation=PlgSettingsStructure.interpolator_regularisation, + **parametric_fault.fault_arguments(spec), + ) + except Exception as e: + raise ValueError(f"Could not build fault '{name}': {e}") from e + + def parametric_faults_to_dict(self) -> dict: + """Return `parametric_faults` in a form that `json.dump` can write.""" + return {name: parametric_fault.clean_spec(spec) for name, spec in self.parametric_faults.items()} + + def parametric_faults_from_dict(self, faults: dict): + """Replace `parametric_faults` with the specs from `parametric_faults_to_dict`.""" + self.parametric_faults = { + name: parametric_fault.spec_from_dict(spec) for name, spec in (faults or {}).items() + } + def _get_feature_by_name_or_none(self, name): """Non-raising counterpart to `GeologicalModel.get_feature_by_name`. @@ -1460,7 +1654,12 @@ def update_model( ) self._progress_callback = progress_callback displacement_fault_count = len(set(self.faults) - set(self.fault_boundaries.values())) - self._progress_total = displacement_fault_count + group_count + len(self.manual_foliations) + self._progress_total = ( + displacement_fault_count + + len(self.parametric_faults) + + group_count + + len(self.manual_foliations) + ) self._progress_current = 0 dbg = getattr(self, '_debug_manager', None) if dbg is not None: @@ -1741,10 +1940,12 @@ def add_foliation( Name for the new foliation feature. data : dict Mapping of layer identifiers to dicts describing each layer. Each - layer dict must include a 'type' key (one of 'Orientation', - 'Form Line', 'Value', 'Inequality') and the fields required by - that type (e.g. 'strike_field', 'dip_field', 'value_field', - 'form_line_constraint', ...). + layer dict must include a 'type' key (one of + `constraints.CONSTRAINT_TYPES`) and the fields required by that + type (e.g. 'strike_field', 'dip_field', 'value_field', + 'form_line_constraint', ...). Each layer dict can also have a + 'weight' and a 'z_source' ('layer', 'dem' or 'constant', with + 'z_value'); see `constraints`. folded_feature_name : str or None Optional name of a feature to which the foliation should be associated/converted (currently unused in this helper). @@ -1880,21 +2081,25 @@ def _foliation_data(self, name, data, sampler=AllSampler(), use_z_coordinate=Fal # auto-synced rows from the map2loop workflow (type values # like 'Contact (auto)') aren't meant for this manual path continue - if layer_data['type'] == 'Orientation': - df = sampler(layer_data['df'], self.dem_function, use_z_coordinate) + layer_type = layer_data['type'] + if layer_type not in constraints.CONSTRAINT_TYPES: + raise ValueError(f"Unknown layer type: {layer_type}") + df = constraints.sample_layer(sampler, layer_data, self.dem_function, use_z_coordinate) + if layer_type == 'Orientation': df['strike'] = df[layer_data['strike_field']] if layer_data.get('orientation_format') == 'Dip Direction': df['strike'] = df['strike'] - 90 df['dip'] = df[layer_data['dip_field']] df['feature_name'] = name - dfs.append(df[['X', 'Y', 'Z', 'strike', 'dip', 'feature_name']]) - elif layer_data['type'] == 'Form Line': - df = sampler(layer_data['df'], self.dem_function, use_z_coordinate) + rows = df[['X', 'Y', 'Z', 'strike', 'dip', 'feature_name']] + dfs.append(constraints.add_weight(rows, df, layer_data)) + elif layer_type == 'Form Line': df['feature_name'] = name if layer_data.get('form_line_constraint') == 'strike': df[['tx', 'ty', 'tz']] = _form_line_tangent_vectors(df) df = df.dropna(subset=['tx', 'ty', 'tz']) - dfs.append(df[['X', 'Y', 'Z', 'tx', 'ty', 'tz', 'feature_name']]) + rows = df[['X', 'Y', 'Z', 'tx', 'ty', 'tz', 'feature_name']] + dfs.append(constraints.add_weight(rows, df, layer_data)) dip = layer_data.get('form_line_dip') if dip is not None and not df.empty: # The line's own direction is a precisely known hard @@ -1919,24 +2124,32 @@ def _foliation_data(self, name, data, sampler=AllSampler(), use_z_coordinate=Fal dip_df['w'] = layer_data.get('form_line_dip_weight', 0.1) dfs.append(dip_df[['X', 'Y', 'Z', 'strike', 'dip', 'w', 'feature_name']]) else: - df['interface'] = df['feature_id'].astype(float) + interface_offset - interface_offset += df['feature_id'].nunique() - dfs.append(df[['X', 'Y', 'Z', 'interface', 'feature_name']]) - elif layer_data['type'] == 'Value': - df = sampler(layer_data['df'], self.dem_function, use_z_coordinate) + rows, interface_offset = constraints.interface_rows( + df, {}, name, interface_offset + ) + dfs.append(constraints.add_weight(rows, df, layer_data)) + elif layer_type == 'Interface': + rows, interface_offset = constraints.interface_rows( + df, layer_data, name, interface_offset + ) + dfs.append(constraints.add_weight(rows, df, layer_data)) + elif layer_type == 'Value': df['val'] = df[layer_data['value_field']] df['feature_name'] = name - dfs.append(df[['X', 'Y', 'Z', 'val', 'feature_name']]) - - elif layer_data['type'] == 'Inequality': - df = sampler(layer_data['df'], self.dem_function, use_z_coordinate) + dfs.append( + constraints.add_weight(df[['X', 'Y', 'Z', 'val', 'feature_name']], df, layer_data) + ) + elif layer_type == 'Inequality': df['l'] = df[layer_data['lower_field']] df['u'] = df[layer_data['upper_field']] df['feature_name'] = name - dfs.append(df[['X', 'Y', 'Z', 'l', 'u', 'feature_name']]) - kwargs['solver'] = 'admm' + rows = df[['X', 'Y', 'Z', 'l', 'u', 'feature_name']] + dfs.append(constraints.add_weight(rows, df, layer_data)) else: - raise ValueError(f"Unknown layer type: {layer_data['type']}") + dfs.append(constraints.constraint_rows(df, layer_data, name)) + solver = constraints.solver_for(layer_type) + if solver: + kwargs['solver'] = solver return pd.concat(dfs, ignore_index=True), kwargs def add_unconformity( diff --git a/loopstructural/main/parametric_fault.py b/loopstructural/main/parametric_fault.py new file mode 100644 index 0000000..c275b65 --- /dev/null +++ b/loopstructural/main/parametric_fault.py @@ -0,0 +1,168 @@ +"""A fault that the user gives by numbers: a centre, a strike, a dip and a size. + +Most faults come from the fault trace layer. A parametric fault has no +trace. It is one ellipsoid in space. This module makes the vectors and the +data that LoopStructural needs from the numbers of the user. + +This module does not import QGIS, so the unit tests can run it. +""" + +import math +from typing import Iterable, List, Optional, Tuple + +import numpy as np +import pandas as pd + +FIELDS = ( + 'name', + 'strike', + 'dip', + 'pitch', + 'displacement', + 'centre', + 'major_axis', + 'intermediate_axis', + 'minor_axis', +) + + +def plane_vectors(strike: float, dip: float, pitch: float) -> Tuple[np.ndarray, np.ndarray]: + """Return the normal vector and the slip vector of a fault plane. + + Parameters + ---------- + strike : float + Strike in degrees clockwise from north. The plane dips to the right of + the strike direction (right-hand rule). + dip : float + Dip in degrees, from 0 (flat) to 90 (vertical). + pitch : float + The angle of the slip in the plane, in degrees from the strike + direction towards the down-dip direction. 0 is a strike-slip fault and + 90 is a dip-slip fault. + + Returns + ------- + normal, slip : numpy.ndarray + Unit vectors (x east, y north, z up). The normal points up. The slip + vector is in the plane, so it is at right angles to the normal. + """ + strike_r, dip_r, pitch_r = (math.radians(a) for a in (strike, dip, pitch)) + dip_direction = strike_r + math.pi / 2 + normal = np.array( + [ + math.sin(dip_r) * math.sin(dip_direction), + math.sin(dip_r) * math.cos(dip_direction), + math.cos(dip_r), + ] + ) + along_strike = np.array([math.sin(strike_r), math.cos(strike_r), 0.0]) + down_dip = np.array( + [ + math.cos(dip_r) * math.sin(dip_direction), + math.cos(dip_r) * math.cos(dip_direction), + -math.sin(dip_r), + ] + ) + slip = math.cos(pitch_r) * along_strike + math.sin(pitch_r) * down_dip + return normal, slip / np.linalg.norm(slip) + + +def default_spec(origin: Iterable[float], maximum: Iterable[float]) -> dict: + """Return the start values for a fault in the model area (the centre and the sizes).""" + origin = np.asarray(list(origin), dtype=float) + maximum = np.asarray(list(maximum), dtype=float) + size = maximum - origin + major = float(max(size[0], size[1])) + return { + 'name': '', + 'strike': 0.0, + 'dip': 90.0, + 'pitch': 90.0, + 'displacement': 100.0, + 'centre': tuple(float(v) for v in (origin + maximum) / 2), + 'major_axis': major, + 'intermediate_axis': major, + 'minor_axis': major / 2, + } + + +def problems(spec: dict, taken_names: Iterable[str] = ()) -> List[str]: + """Return what is wrong with a fault, as text for the user. An empty list is a good fault.""" + found = [] + name = str(spec.get('name', '')).strip() + if not name: + found.append("Give the fault a name.") + elif name in set(taken_names): + found.append(f"The name '{name}' is already used by a feature or a fault.") + if not 0 < float(spec.get('dip', 0)) <= 90: + found.append("The dip must be above 0 and not more than 90 degrees.") + for key, label in ( + ('major_axis', "The length along strike"), + ('intermediate_axis', "The extent down dip"), + ('minor_axis', "The influence distance"), + ): + if not float(spec.get(key) or 0) > 0: + found.append(f"{label} must be more than 0.") + centre = spec.get('centre') + if centre is None or len(centre) != 3 or not all(math.isfinite(float(v)) for v in centre): + found.append("The centre needs X, Y and Z.") + return found + + +def frame_data(spec: dict) -> pd.DataFrame: + """Return the data frame with the one point that LoopStructural needs. + + The point is the centre of the fault, with the value 0 of the fault + surface. The vectors and the sizes go to LoopStructural as separate + arguments, see `fault_arguments`. + """ + x, y, z = (float(v) for v in spec['centre']) + return pd.DataFrame( + { + 'X': [x], + 'Y': [y], + 'Z': [z], + 'feature_name': [spec['name']], + 'val': [0.0], + 'coord': [0], + } + ) + + +def fault_arguments(spec: dict) -> dict: + """Return the arguments of ``GeologicalModel.create_and_add_fault`` for a fault. + + ``name`` and ``data`` are not in the result: pass the name, and + ``frame_data(spec)`` as the data. + """ + normal, slip = plane_vectors(spec['strike'], spec['dip'], spec.get('pitch', 90.0)) + return { + 'displacement': float(spec['displacement']), + 'fault_normal_vector': normal, + 'fault_slip_vector': slip, + 'fault_center': np.array([float(v) for v in spec['centre']]), + 'major_axis': float(spec['major_axis']), + 'intermediate_axis': float(spec['intermediate_axis']), + 'minor_axis': float(spec['minor_axis']), + 'fault_dip': float(spec['dip']), + } + + +def clean_spec(spec: dict) -> dict: + """Return the spec as plain numbers and text, safe for JSON.""" + result = {key: spec[key] for key in FIELDS if key in spec} + result['name'] = str(result['name']).strip() + for key in ('strike', 'dip', 'pitch', 'displacement', 'major_axis', 'intermediate_axis', + 'minor_axis'): + if key in result: + result[key] = float(result[key]) + result['centre'] = [float(v) for v in result['centre']] + return result + + +def spec_from_dict(data: Optional[dict]) -> dict: + """Read a spec that was saved. The centre is a tuple again.""" + spec = clean_spec(data) + spec['centre'] = tuple(spec['centre']) + return spec diff --git a/loopstructural/main/preview.py b/loopstructural/main/preview.py new file mode 100644 index 0000000..8451457 --- /dev/null +++ b/loopstructural/main/preview.py @@ -0,0 +1,106 @@ +"""The map preview of one feature: isolines of its scalar field. + +The scalar field is evaluated on a regular grid over the model area. Each +point has the Z of the DEM, so the lines are the trace of the surfaces on the +ground, as on a geological map. + +This module does not import QGIS, so the unit tests can run it. +""" + +from typing import Callable, List, Optional, Tuple + +import numpy as np + +DEFAULT_RESOLUTION = 150 +DEFAULT_LEVEL_COUNT = 10 + + +def grid_points( + origin, maximum, resolution: int, dem_function: Optional[Callable] = None +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return the grid of the preview. + + Parameters + ---------- + origin, maximum : sequence of float + The lower and the upper corner of the model area. Only X and Y are used. + resolution : int + The number of points along each side. + dem_function : callable, optional + ``dem_function(x, y)`` gives the Z of the ground. Without it, Z is 0. + + Returns + ------- + x, y : numpy.ndarray + The coordinates along each side, with shape ``(resolution,)``. + points : numpy.ndarray + The points, with shape ``(resolution * resolution, 3)``. The order is + row by row: the first ``resolution`` points have the first Y. + """ + if resolution < 2: + raise ValueError("The resolution must be 2 or more.") + x = np.linspace(origin[0], maximum[0], resolution) + y = np.linspace(origin[1], maximum[1], resolution) + xx, yy = np.meshgrid(x, y) + if dem_function is None: + zz = np.zeros_like(xx) + else: + zz = np.vectorize(dem_function, otypes=[float])(xx, yy) + points = np.column_stack([xx.ravel(), yy.ravel(), zz.ravel()]) + return x, y, points + + +def default_levels(values, count: int = DEFAULT_LEVEL_COUNT) -> np.ndarray: + """Return ``count`` levels, spaced evenly between the lowest and the highest value. + + The lowest and the highest value are not levels, because a line there is a + point or an edge. An empty array means that the field has no range. + """ + values = np.asarray(values, dtype=float) + values = values[np.isfinite(values)] + if values.size == 0 or count < 1: + return np.array([]) + low, high = values.min(), values.max() + if not high > low: + return np.array([]) + return np.linspace(low, high, count + 2)[1:-1] + + +def isolines(x, y, values, levels) -> List[Tuple[float, np.ndarray]]: + """Return the lines of the field ``values`` at each level. + + Parameters + ---------- + x, y : array_like + The coordinates along each side of the grid. + values : array_like + The field, with shape ``(len(y), len(x))``. A value that is not finite + is a gap: no line crosses it. + levels : iterable of float + The levels of the lines. + + Returns + ------- + list of tuple + ``(level, coordinates)`` for each line. ``coordinates`` has shape + ``(n, 2)`` with X and Y. + """ + try: + import contourpy + except ImportError as err: # pragma: no cover - contourpy comes with matplotlib + raise RuntimeError("The preview needs the contourpy package (part of matplotlib).") from err + values = np.asarray(values, dtype=float) + if values.shape != (len(y), len(x)): + raise ValueError("values must have the shape (len(y), len(x)).") + generator = contourpy.contour_generator( + np.asarray(x, dtype=float), + np.asarray(y, dtype=float), + np.ma.masked_invalid(values), + ) + lines = [] + for level in levels: + for line in generator.lines(float(level)): + line = np.asarray(line, dtype=float) + if len(line) >= 2: + lines.append((float(level), line)) + return lines diff --git a/loopstructural/main/workflow_mode.py b/loopstructural/main/workflow_mode.py new file mode 100644 index 0000000..db1aea1 --- /dev/null +++ b/loopstructural/main/workflow_mode.py @@ -0,0 +1,19 @@ +"""The start choice of the user: how the model is built. + +This module does not import QGIS, so the unit tests can run it. +""" + +# Steps 2 (stratigraphy) and 3 (faults) make features and constraints for step 4. +WORKFLOW_MODE_MAP = 'map' +# Steps 2 and 3 are hidden. The user gives the constraints in step 4. +WORKFLOW_MODE_CONSTRAINTS = 'constraints' +WORKFLOW_MODES = (WORKFLOW_MODE_MAP, WORKFLOW_MODE_CONSTRAINTS) +DEFAULT_WORKFLOW_MODE = WORKFLOW_MODE_MAP + +WORKFLOW_MODE_LABELS = { + WORKFLOW_MODE_MAP: "Build from a geological map", + WORKFLOW_MODE_CONSTRAINTS: "Interpolate surfaces from constraints", +} + +# The steps that the "constraints" choice hides. The user can show them later. +CONSTRAINT_MODE_HIDDEN_STEPS = ('stratigraphy', 'faults') diff --git a/tests/qgis/test_derived_data_state.py b/tests/qgis/test_derived_data_state.py index 4f4ca42..d4732bb 100644 --- a/tests/qgis/test_derived_data_state.py +++ b/tests/qgis/test_derived_data_state.py @@ -267,3 +267,38 @@ def test_the_sources_are_saved_and_loaded(self, data_manager): data_manager.thickness_sources.clear() data_manager.update_from_dict(state) assert data_manager.get_thickness_source(unit.uuid) == derived_data.CALCULATED + + +class TestWorkflowMode: + def test_the_default_is_the_map(self, data_manager): + assert data_manager.workflow_mode == 'map' + + def test_the_choice_is_saved_and_loaded(self, data_manager): + data_manager.set_workflow_mode('constraints') + state = json.loads(json.dumps(data_manager.to_dict())) + data_manager.set_workflow_mode('map') + data_manager.update_from_dict(state) + assert data_manager.workflow_mode == 'constraints' + + def test_a_state_file_of_an_older_version_gives_the_map(self, data_manager): + data_manager.set_workflow_mode('constraints') + state = json.loads(json.dumps(data_manager.to_dict())) + del state['workflow_mode'] + data_manager.update_from_dict(state) + assert data_manager.workflow_mode == 'map' + + def test_a_change_calls_the_listeners_one_time(self, data_manager): + seen = [] + data_manager.add_workflow_mode_callback(seen.append) + data_manager.set_workflow_mode('constraints') + data_manager.set_workflow_mode('constraints') + assert seen == ['constraints'] + + def test_an_unknown_choice_is_an_error(self, data_manager): + with pytest.raises(ValueError): + data_manager.set_workflow_mode('other') + + def test_reset_gives_the_map(self, data_manager): + data_manager.set_workflow_mode('constraints') + data_manager.reset() + assert data_manager.workflow_mode == 'map' diff --git a/tests/qgis/test_detach_feature.py b/tests/qgis/test_detach_feature.py new file mode 100644 index 0000000..1232f32 --- /dev/null +++ b/tests/qgis/test_detach_feature.py @@ -0,0 +1,128 @@ +"""Pytest tests for detaching a feature that the stratigraphic column makes. + +A detached feature keeps a copy of its contact and orientation data. Rows that +the user adds to a generated feature are used together with the data of the +column. +""" + +import geopandas as gpd +import pandas as pd +import pytest +from LoopStructural import StratigraphicColumn + +from loopstructural.main.model_manager import GeologicalModelManager + + +def _contact(x): + return pd.DataFrame({'X': [x], 'Y': [0.0], 'Z': [0.0]}) + + +@pytest.fixture +def manager(monkeypatch): + manager = GeologicalModelManager() + captured = [] + + def fake_create_and_add_foliation(name, data=None, **kwargs): + captured.append((name, data, kwargs)) + return object() + + monkeypatch.setattr(manager.model, 'create_and_add_foliation', fake_create_and_add_foliation) + monkeypatch.setattr(manager.model, 'add_unconformity', lambda *a, **k: None) + manager._captured = captured + + column = StratigraphicColumn() + column.clear(basement=False) + column.add_unit(name='oldest', thickness=100.0, where='top') + column.add_unit(name='youngest', thickness=200.0, where='top') + manager.stratigraphic_column = column + manager.stratigraphy['oldest']['contact'] = _contact(1.0) + manager.stratigraphy['youngest']['contact'] = _contact(2.0) + return manager + + +def _group_name(manager): + return manager.stratigraphic_column.get_groups()[0].name + + +def _built_x(manager): + manager.update_foliation_features() + return sorted(manager._captured[-1][1]['X']) + + +class TestDetach: + def test_a_group_of_the_column_is_a_generated_feature(self, manager): + assert manager.is_generated(_group_name(manager)) + assert not manager.is_generated('other') + + def test_a_detached_feature_keeps_its_data(self, manager): + name = _group_name(manager) + assert manager.detach_feature(name) + manager.stratigraphy['oldest']['contact'] = _contact(99.0) + assert _built_x(manager) == [1.0, 2.0] + + def test_a_feature_that_is_not_detached_follows_the_contacts(self, manager): + manager.stratigraphy['oldest']['contact'] = _contact(99.0) + assert _built_x(manager) == [2.0, 99.0] + + def test_a_feature_without_data_cannot_be_detached(self, manager): + manager.stratigraphy.clear() + assert not manager.detach_feature(_group_name(manager)) + + def test_a_feature_that_is_not_generated_cannot_be_detached(self, manager): + assert not manager.detach_feature('other') + + def test_attach_follows_the_column_again_and_the_model_is_stale(self, manager): + name = _group_name(manager) + manager.detach_feature(name) + manager.stratigraphy['oldest']['contact'] = _contact(99.0) + assert manager.attach_feature(name) + assert _built_x(manager) == [2.0, 99.0] + assert manager._data_dirty + + def test_the_detached_data_is_saved_and_loaded(self, manager): + name = _group_name(manager) + manager.detach_feature(name) + written = manager.detached_to_dict() + manager.detached = {} + manager.detached_from_dict(written) + assert manager.is_detached(name) + assert sorted(manager.detached[name]['X']) == [1.0, 2.0] + + def test_reset_forgets_detached_features(self, manager): + manager.detach_feature(_group_name(manager)) + manager.reset() + assert manager.detached == {} + + +class TestAddedRows: + def _extra(self, x): + layer = { + 'layer_name': 'values', + 'type': 'Value', + 'value_field': 'v', + 'df': gpd.GeoDataFrame({'v': [5.0]}, geometry=gpd.points_from_xy([x], [0.0])), + } + return {'data': {'values': layer}, 'use_z_coordinate': False} + + def test_added_rows_are_used_with_the_data_of_the_column(self, manager): + name = _group_name(manager) + manager.extra_constraints[name] = self._extra(7.0) + assert _built_x(manager) == [1.0, 2.0, 7.0] + + def test_added_rows_stay_when_a_feature_is_detached(self, manager): + name = _group_name(manager) + manager.extra_constraints[name] = self._extra(7.0) + manager.detach_feature(name) + # the copy has the data of the column only + assert sorted(manager.detached[name]['X']) == [1.0, 2.0] + assert _built_x(manager) == [1.0, 2.0, 7.0] + + def test_an_inequality_row_gives_the_solver(self, manager): + name = _group_name(manager) + spec = self._extra(7.0) + spec['data']['values'].update( + {'type': 'Inequality', 'lower_field': 'v', 'upper_field': 'v'} + ) + manager.extra_constraints[name] = spec + manager.update_foliation_features() + assert manager._captured[-1][2]['solver'] == 'admm' diff --git a/tests/qgis/test_manual_foliations.py b/tests/qgis/test_manual_foliations.py index f699818..8e27583 100644 --- a/tests/qgis/test_manual_foliations.py +++ b/tests/qgis/test_manual_foliations.py @@ -136,3 +136,60 @@ def test_loaded_foliation_is_built_again_by_update_model(self, manager): manager.update_model(notify_observers=False) assert 's1' in _names(manager) + + +def _points(**fields): + points = [Point(x, y, 0.0) for x in (20.0, 50.0, 80.0) for y in (20.0, 50.0, 80.0)] + columns = {name: [fn(p) for p in points] for name, fn in fields.items()} + return gpd.GeoDataFrame(columns, geometry=points) + + +class TestConstraintTypes: + """A feature built from constraint layers only, with no stratigraphic column data.""" + + def _rows(self, manager, layer): + data, kwargs = manager._foliation_data('f', {'layer': layer}, use_z_coordinate=True) + return data, kwargs + + def test_interface_rows_use_the_group_field(self, manager): + df = _points(group=lambda p: 'a' if p.x < 50 else 'b') + layer = {'type': 'Interface', 'group_field': 'group', 'df': df} + data, _ = self._rows(manager, layer) + assert sorted(data['interface'].unique()) == [0.0, 1.0] + + def test_weight_goes_to_the_rows(self, manager): + layer = _value_layer() + layer['weight'] = 0.25 + data, _ = self._rows(manager, layer) + assert (data['w'] == 0.25).all() + + def test_a_constant_z_source_gives_one_z(self, manager): + layer = _value_layer() + layer.update({'z_source': 'constant', 'z_value': 12.5}) + data, _ = self._rows(manager, layer) + assert (data['Z'] == 12.5).all() + + def test_a_pairwise_inequality_needs_the_admm_solver(self, manager): + df = _points(order=lambda p: 1 if p.x < 50 else 2) + layer = {'type': 'Pairwise Inequality', 'pair_field': 'order', 'df': df} + data, kwargs = self._rows(manager, layer) + assert kwargs == {'solver': 'admm'} + assert 'pair_id' in data + + def test_a_feature_from_value_and_normal_layers_solves_without_a_column(self, manager): + values = _value_layer() + normals = { + 'layer_name': 'normals', + 'type': 'Gradient/Normal', + 'vector_kind': 'normal', + 'vector_x_field': 'nx', + 'vector_y_field': 'ny', + 'vector_z_field': 'nz', + 'df': _points(nx=lambda p: 1.0, ny=lambda p: 0.0, nz=lambda p: 0.0), + } + manager.add_foliation( + 's1', {'values': values, 'normals': normals}, use_z_coordinate=True + ) + manager.update_all_features(notify_observers=False) + result = manager.model['s1'].evaluate_value(np.array([[50.0, 50.0, 0.0]])) + assert not np.any(np.isnan(result)) diff --git a/tests/qgis/test_parametric_fault.py b/tests/qgis/test_parametric_fault.py new file mode 100644 index 0000000..d8fd937 --- /dev/null +++ b/tests/qgis/test_parametric_fault.py @@ -0,0 +1,97 @@ +"""Pytest tests for the faults that the user gives by numbers (Add Fault). + +The fault has no trace layer. `update_model` builds it together with the other +faults, so a fault that was added must be built again at each build. +""" + +import json + +import numpy as np +import pytest +from LoopStructural import StratigraphicColumn +from LoopStructural.datatypes import BoundingBox + +from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.toolbelt.preferences import PlgSettingsStructure + + +class _DebugManager: + def log(self, *args, **kwargs): + pass + + +@pytest.fixture +def manager(monkeypatch): + monkeypatch.setattr(PlgSettingsStructure, 'interpolator_nelements', 200) + manager = GeologicalModelManager(debug_manager=_DebugManager()) + manager.update_bounding_box(BoundingBox(origin=[0, 0, -50], maximum=[100, 100, 50])) + manager.stratigraphic_column = StratigraphicColumn() + return manager + + +def _spec(**changes): + spec = { + 'name': 'F1', + 'strike': 0.0, + 'dip': 90.0, + 'pitch': 90.0, + 'displacement': 10.0, + 'centre': (50.0, 50.0, 0.0), + 'major_axis': 100.0, + 'intermediate_axis': 100.0, + 'minor_axis': 30.0, + } + spec.update(changes) + return spec + + +def _names(manager): + return [f.name for f in manager.model.features] + + +class TestParametricFault: + def test_the_fault_is_built_by_update_model(self, manager): + manager.add_parametric_fault(_spec()) + manager.update_model(notify_observers=False) + assert 'F1' in _names(manager) + + def test_the_fault_surface_is_where_the_fault_is(self, manager): + manager.add_parametric_fault(_spec()) + manager.update_model(notify_observers=False) + manager.update_all_features(notify_observers=False) + # the fault strikes north through x=50, so the field changes sign across it + west = manager.model['F1'][0].evaluate_value(np.array([[30.0, 50.0, 0.0]])) + east = manager.model['F1'][0].evaluate_value(np.array([[70.0, 50.0, 0.0]])) + assert west[0] * east[0] < 0 + + def test_adding_a_fault_makes_the_model_stale(self, manager): + manager.add_parametric_fault(_spec()) + assert manager._data_dirty + + def test_a_name_that_is_used_is_an_error(self, manager): + manager.add_parametric_fault(_spec()) + with pytest.raises(ValueError): + manager.add_parametric_fault(_spec()) + + def test_a_bad_fault_is_an_error(self, manager): + with pytest.raises(ValueError): + manager.add_parametric_fault(_spec(dip=0.0)) + + def test_a_removed_fault_is_not_built_again(self, manager): + manager.add_parametric_fault(_spec()) + manager.remove_parametric_fault('F1') + manager.update_model(notify_observers=False) + assert 'F1' not in _names(manager) + + def test_the_faults_are_saved_and_loaded(self, manager): + manager.add_parametric_fault(_spec()) + written = json.loads(json.dumps(manager.parametric_faults_to_dict())) + manager.parametric_faults = {} + manager.parametric_faults_from_dict(written) + manager.update_model(notify_observers=False) + assert 'F1' in _names(manager) + + def test_reset_forgets_the_faults(self, manager): + manager.add_parametric_fault(_spec()) + manager.reset() + assert manager.parametric_faults == {} diff --git a/tests/unit/test_build_plan.py b/tests/unit/test_build_plan.py new file mode 100644 index 0000000..72d6a48 --- /dev/null +++ b/tests/unit/test_build_plan.py @@ -0,0 +1,174 @@ +"""Pytest tests for the primary action of the model step and for the choice of +the derived data that a build calculates again. + +The modules do not import QGIS, so the tests use fake managers and run in the +fast tests/unit/ job. +""" + +from types import SimpleNamespace + +import pytest + +from loopstructural.gui.modelling.steps import build_plan +from loopstructural.gui.modelling.steps.status import StepCheck +from loopstructural.main import derived_data, layer_roles +from loopstructural.main.derived_data import DerivedData +from loopstructural.main.derived_refresh import names_to_refresh, thickness_values +from loopstructural.main.workflow_mode import WORKFLOW_MODE_CONSTRAINTS, WORKFLOW_MODE_MAP + + +class TestPrimaryAction: + def test_an_empty_model_is_built(self): + action = build_plan.choose_primary_action('empty') + assert (action.action, action.text) == (build_plan.ACTION_BUILD, "Build model") + + def test_an_initialized_model_is_solved(self): + action = build_plan.choose_primary_action('initialized') + assert (action.action, action.text) == (build_plan.ACTION_SOLVE, "Solve model") + + def test_a_solved_model_can_be_solved_again(self): + action = build_plan.choose_primary_action('solved') + assert action.action == build_plan.ACTION_SOLVE + assert action.text == "Solve again" + + def test_a_stale_model_is_rebuilt(self): + action = build_plan.choose_primary_action('stale') + assert (action.action, action.text) == (build_plan.ACTION_BUILD, "Rebuild model") + + @pytest.mark.parametrize('state', ['initialized', 'solved']) + def test_out_of_date_derived_data_makes_a_rebuild(self, state): + action = build_plan.choose_primary_action( + state, derived_out_of_date=[derived_data.BASAL_CONTACTS] + ) + assert action.action == build_plan.ACTION_BUILD + assert 'out of date' in action.tooltip + + def test_a_changed_layer_updates_the_data_and_solves(self): + action = build_plan.choose_primary_action('solved', layers_changed=True) + assert action.action == build_plan.ACTION_SOLVE + assert action.text == "Update data and solve" + + def test_a_blocked_build_is_disabled_with_the_reason(self): + action = build_plan.choose_primary_action('solved', blocked_reason="Set the bounding box.") + assert action.action == build_plan.ACTION_NONE + assert not action.enabled + assert action.tooltip == "Set the bounding box." + + +class TestProblems: + def test_problems_of_all_steps_in_the_order_of_the_steps(self): + problems = build_plan.collect_problems( + [ + ('data', StepCheck(problems=('Bad CRS.',))), + ('stratigraphy', StepCheck(todo=('Add units.',))), + ('faults', StepCheck(problems=('No faults.',))), + ] + ) + assert problems == [('data', 'Bad CRS.'), ('faults', 'No faults.')] + + def test_the_same_problem_shows_one_time(self): + problems = build_plan.collect_problems( + [ + ('stratigraphy', StepCheck(problems=('Contacts are out of date.',))), + ('model', StepCheck(problems=('Contacts are out of date.',))), + ] + ) + assert problems == [('stratigraphy', 'Contacts are out of date.')] + + +class TestModes: + ALL = ['data', 'stratigraphy', 'faults', 'model', 'view'] + + def test_the_map_choice_shows_all_steps(self): + assert build_plan.steps_for_mode(WORKFLOW_MODE_MAP, self.ALL) == self.ALL + + def test_the_constraints_choice_hides_steps_2_and_3(self): + shown = build_plan.steps_for_mode(WORKFLOW_MODE_CONSTRAINTS, self.ALL) + assert shown == ['data', 'model', 'view'] + + +class FakeRoles: + def __init__(self): + self.values = {layer_roles.GEOLOGY: object(), layer_roles.GEOLOGY_UNIT_FIELD: 'UNITNAME'} + self.contacts_source = layer_roles.CONTACTS_FROM_GEOLOGY + + +class FakeDataManager: + """The parts of the data manager that `names_to_refresh` reads.""" + + def __init__(self): + self.units = ['a', 'b'] + self.layer_roles = FakeRoles() + self.inputs = {'unit_order': ['a', 'b']} + self.derived = DerivedData() + for name in (derived_data.BASAL_CONTACTS, derived_data.THICKNESS): + self.derived.register(name, lambda: self.inputs) + + def get_stratigraphic_unit_names(self): + return list(self.units) + + def get_layer_role(self, role): + return self.layer_roles.values.get(role) + + +@pytest.fixture +def dm(): + return FakeDataManager() + + +class TestNamesToRefresh: + def test_contacts_that_were_never_calculated_are_calculated(self, dm): + assert names_to_refresh(dm) == [derived_data.BASAL_CONTACTS] + + def test_current_contacts_are_not_calculated_again(self, dm): + dm.derived.record(derived_data.BASAL_CONTACTS) + assert names_to_refresh(dm) == [] + + def test_a_reorder_makes_contacts_and_thickness_out_of_date_in_order(self, dm): + dm.derived.record(derived_data.BASAL_CONTACTS) + dm.derived.record(derived_data.THICKNESS) + dm.inputs = {'unit_order': ['b', 'a']} + assert names_to_refresh(dm) == [derived_data.BASAL_CONTACTS, derived_data.THICKNESS] + + def test_a_reorder_and_its_undo_calculate_nothing(self, dm): + dm.derived.record(derived_data.BASAL_CONTACTS) + dm.inputs = {'unit_order': ['b', 'a']} + dm.inputs = {'unit_order': ['a', 'b']} + assert names_to_refresh(dm) == [] + + def test_a_thickness_that_was_never_calculated_can_be_typed(self, dm): + dm.derived.record(derived_data.BASAL_CONTACTS) + assert derived_data.THICKNESS not in names_to_refresh(dm) + + def test_a_contacts_layer_of_the_user_is_never_calculated(self, dm): + dm.layer_roles.contacts_source = layer_roles.CONTACTS_FROM_LAYER + assert derived_data.BASAL_CONTACTS not in names_to_refresh(dm) + + def test_no_column_calculates_no_contacts(self, dm): + dm.units = [] + assert names_to_refresh(dm) == [] + + def test_no_geology_layer_calculates_no_contacts(self, dm): + dm.layer_roles.values[layer_roles.GEOLOGY] = None + assert names_to_refresh(dm) == [] + + +class TestThicknessValues: + def test_the_median_is_used_and_a_unit_with_no_result_is_left_out(self): + pd = pytest.importorskip('pandas') + table = pd.DataFrame( + { + 'name': ['a', 'b', 'c'], + 'ThicknessMedian': [10.0, -1.0, float('nan')], + 'ThicknessMean': [11.0, 12.0, 13.0], + } + ) + assert thickness_values(table) == {'a': 10.0} + + def test_the_mean_is_used_without_a_median(self): + pd = pytest.importorskip('pandas') + table = pd.DataFrame({'name': ['a'], 'ThicknessMean': [5.0]}) + assert thickness_values(table) == {'a': 5.0} + + def test_a_table_without_thickness_columns_gives_nothing(self): + assert thickness_values(SimpleNamespace(columns=['name'])) == {} diff --git a/tests/unit/test_constraints.py b/tests/unit/test_constraints.py new file mode 100644 index 0000000..4a28caa --- /dev/null +++ b/tests/unit/test_constraints.py @@ -0,0 +1,173 @@ +"""Pytest tests for the constraint types of a feature. + +`loopstructural.main.constraints` does not import QGIS, so the tests run in +the fast tests/unit/ job. +""" + +import numpy as np +import pandas as pd +import pytest + +from loopstructural.main import constraints + + +def sampler_of(points): + """A sampler that ignores the layer and gives fixed points.""" + calls = [] + + def sampler(df, dem, use_z): + calls.append(use_z) + return points.copy() + + sampler.calls = calls + return sampler + + +@pytest.fixture +def points(): + return pd.DataFrame( + { + 'X': [0.0, 1.0, 2.0, 3.0], + 'Y': [0.0, 1.0, 2.0, 3.0], + 'Z': [10.0, 11.0, 12.0, 13.0], + 'feature_id': [0, 0, 1, 1], + 'gx': [0.0, 0.0, 1.0, np.nan], + 'gy': [0.0, 0.0, 0.0, 1.0], + 'gz': [1.0, 0.0, 0.0, 1.0], + 'group': ['a', 'b', 'a', None], + 'order': [1, 1, 2, 2], + } + ) + + +class TestTypes: + def test_every_type_has_a_description_and_a_field_rule(self): + for layer_type in constraints.CONSTRAINT_TYPES: + assert constraints.DESCRIPTIONS[layer_type] + assert layer_type in constraints.REQUIRED_FIELDS + + def test_the_six_types_of_the_plan_are_there(self): + for name in ( + 'Value', + 'Interface', + 'Gradient/Normal', + 'Tangent', + 'Inequality', + 'Pairwise Inequality', + ): + assert name in constraints.CONSTRAINT_TYPES + + def test_missing_fields(self): + row = {'type': constraints.GRADIENT_NORMAL, 'vector_x_field': 'gx'} + assert constraints.missing_fields(row) == ['vector_y_field', 'vector_z_field'] + assert constraints.missing_fields({'type': constraints.INTERFACE}) == [] + + def test_only_inequality_types_need_the_admm_solver(self): + assert constraints.solver_for(constraints.INEQUALITY) == 'admm' + assert constraints.solver_for(constraints.PAIRWISE_INEQUALITY) == 'admm' + assert constraints.solver_for(constraints.VALUE) is None + + +class TestZSource: + def test_a_row_of_an_older_version_uses_the_default(self, points): + sampler = sampler_of(points) + constraints.sample_layer(sampler, {'df': None}, None, default_use_z=True) + assert sampler.calls == [True] + + @pytest.mark.parametrize( + 'source, expected', [('layer', True), ('dem', False), ('constant', False)] + ) + def test_the_source_decides_if_the_layer_z_is_used(self, points, source, expected): + sampler = sampler_of(points) + constraints.sample_layer(sampler, {'df': None, 'z_source': source}, None, True) + assert sampler.calls == [expected] + + def test_a_constant_replaces_all_z(self, points): + rows = constraints.sample_layer( + sampler_of(points), {'df': None, 'z_source': 'constant', 'z_value': -50}, None + ) + assert (rows['Z'] == -50.0).all() + + +class TestWeight: + def test_no_weight_adds_no_column(self, points): + assert 'w' not in constraints.add_weight(points[['X', 'Y', 'Z']], points, {}) + + def test_the_weight_is_added(self, points): + rows = constraints.add_weight(points[['X', 'Y', 'Z']], points, {'weight': 0.5}) + assert (rows['w'] == 0.5).all() + + def test_a_weight_in_the_rows_is_kept(self, points): + rows = points[['X', 'Y', 'Z']].assign(w=0.1) + assert (constraints.add_weight(rows, points, {'weight': 5})['w'] == 0.1).all() + + +class TestInterface: + def test_each_feature_is_a_surface_without_a_group_field(self, points): + rows, offset = constraints.interface_rows(points, {}, 'f') + assert list(rows['interface']) == [0, 0, 1, 1] + assert offset == 2 + + def test_the_offset_keeps_the_surfaces_of_two_rows_apart(self, points): + rows, offset = constraints.interface_rows(points, {}, 'f', offset=5) + assert list(rows['interface']) == [5, 5, 6, 6] + assert offset == 7 + + def test_a_group_field_gives_the_surfaces_and_drops_empty_values(self, points): + rows, offset = constraints.interface_rows(points, {'group_field': 'group'}, 'f') + assert list(rows['interface']) == [0, 1, 0] + assert offset == 2 + assert (rows['feature_name'] == 'f').all() + + +class TestVectors: + def test_a_gradient_gives_the_g_columns_and_drops_bad_rows(self, points): + row = { + 'type': constraints.GRADIENT_NORMAL, + 'vector_x_field': 'gx', + 'vector_y_field': 'gy', + 'vector_z_field': 'gz', + } + rows = constraints.constraint_rows(points, row, 'f') + # row 1 is a zero vector and row 3 has no x + assert list(rows.index) == [0, 2] + assert {'gx', 'gy', 'gz'} <= set(rows.columns) + + def test_a_normal_gives_the_n_columns(self, points): + row = { + 'type': constraints.GRADIENT_NORMAL, + 'vector_kind': constraints.KIND_NORMAL, + 'vector_x_field': 'gx', + 'vector_y_field': 'gy', + 'vector_z_field': 'gz', + 'weight': 2.0, + } + rows = constraints.constraint_rows(points, row, 'f') + assert {'nx', 'ny', 'nz'} <= set(rows.columns) + assert (rows['w'] == 2.0).all() + + def test_a_tangent_gives_the_t_columns(self, points): + row = { + 'type': constraints.TANGENT, + 'vector_x_field': 'gx', + 'vector_y_field': 'gy', + 'vector_z_field': 'gz', + } + assert {'tx', 'ty', 'tz'} <= set(constraints.constraint_rows(points, row, 'f').columns) + + +class TestPairwise: + def test_the_group_number_gives_the_pair_id(self, points): + row = {'type': constraints.PAIRWISE_INEQUALITY, 'pair_field': 'order'} + rows = constraints.constraint_rows(points, row, 'f') + assert list(rows['pair_id']) == [1, 1, 2, 2] + + def test_a_row_without_a_number_is_dropped(self, points): + points.loc[0, 'order'] = np.nan + row = {'type': constraints.PAIRWISE_INEQUALITY, 'pair_field': 'order'} + assert len(constraints.constraint_rows(points, row, 'f')) == 3 + + +def test_an_unknown_type_is_an_error(points): + with pytest.raises(ValueError): + constraints.constraint_rows(points, {'type': 'Other'}, 'f') diff --git a/tests/unit/test_parametric_fault.py b/tests/unit/test_parametric_fault.py new file mode 100644 index 0000000..b6699c9 --- /dev/null +++ b/tests/unit/test_parametric_fault.py @@ -0,0 +1,113 @@ +"""Pytest tests for the faults that the user gives by numbers. + +`loopstructural.main.parametric_fault` does not import QGIS, so the tests run +in the fast tests/unit/ job. +""" + +import json + +import numpy as np +import pytest + +from loopstructural.main import parametric_fault + + +def good_spec(**changes): + spec = { + 'name': 'F1', + 'strike': 0.0, + 'dip': 60.0, + 'pitch': 90.0, + 'displacement': 50.0, + 'centre': (5.0, 6.0, 7.0), + 'major_axis': 100.0, + 'intermediate_axis': 80.0, + 'minor_axis': 20.0, + } + spec.update(changes) + return spec + + +class TestPlaneVectors: + def test_a_vertical_plane_that_strikes_north_has_an_east_normal(self): + normal, _ = parametric_fault.plane_vectors(0, 90, 0) + assert np.allclose(normal, [1, 0, 0], atol=1e-9) + + def test_a_flat_plane_has_an_up_normal(self): + normal, _ = parametric_fault.plane_vectors(123, 0, 0) + assert np.allclose(normal, [0, 0, 1], atol=1e-9) + + def test_the_plane_dips_to_the_right_of_the_strike(self): + # strike east, so the dip direction is south: the normal tilts to the south + normal, _ = parametric_fault.plane_vectors(90, 45, 0) + assert normal[1] < 0 < normal[2] + + @pytest.mark.parametrize('strike', [0, 37, 90, 210]) + @pytest.mark.parametrize('dip', [10, 45, 90]) + @pytest.mark.parametrize('pitch', [0, 30, 90, -45]) + def test_the_slip_is_a_unit_vector_in_the_plane(self, strike, dip, pitch): + normal, slip = parametric_fault.plane_vectors(strike, dip, pitch) + assert np.isclose(np.linalg.norm(normal), 1) + assert np.isclose(np.linalg.norm(slip), 1) + assert np.isclose(normal @ slip, 0, atol=1e-9) + + def test_pitch_0_is_along_the_strike_and_pitch_90_is_down_the_dip(self): + _, along = parametric_fault.plane_vectors(0, 90, 0) + _, down = parametric_fault.plane_vectors(0, 90, 90) + assert np.allclose(along, [0, 1, 0], atol=1e-9) + assert np.allclose(down, [0, 0, -1], atol=1e-9) + + +class TestProblems: + def test_a_good_fault_has_no_problem(self): + assert parametric_fault.problems(good_spec()) == [] + + def test_the_name_is_needed_and_must_be_new(self): + assert parametric_fault.problems(good_spec(name=' ')) + assert 'already used' in parametric_fault.problems(good_spec(), ['F1'])[0] + + @pytest.mark.parametrize('dip', [0, -5, 91]) + def test_the_dip_must_be_in_range(self, dip): + assert parametric_fault.problems(good_spec(dip=dip)) + + @pytest.mark.parametrize('axis', ['major_axis', 'intermediate_axis', 'minor_axis']) + def test_a_size_must_be_above_zero(self, axis): + assert parametric_fault.problems(good_spec(**{axis: 0})) + + def test_the_centre_needs_three_numbers(self): + assert parametric_fault.problems(good_spec(centre=(1.0, 2.0))) + assert parametric_fault.problems(good_spec(centre=(1.0, 2.0, float('nan')))) + + +class TestLoopStructuralArguments: + def test_the_frame_data_is_the_centre_on_the_fault_surface(self): + data = parametric_fault.frame_data(good_spec()) + row = data.iloc[0] + assert (row['X'], row['Y'], row['Z']) == (5.0, 6.0, 7.0) + assert row['val'] == 0 and row['coord'] == 0 and row['feature_name'] == 'F1' + + def test_the_arguments(self): + arguments = parametric_fault.fault_arguments(good_spec()) + assert arguments['displacement'] == 50.0 + assert arguments['major_axis'] == 100.0 + assert arguments['intermediate_axis'] == 80.0 + assert arguments['minor_axis'] == 20.0 + assert list(arguments['fault_center']) == [5.0, 6.0, 7.0] + assert np.isclose(arguments['fault_normal_vector'] @ arguments['fault_slip_vector'], 0) + + +class TestSaving: + def test_a_spec_survives_json(self): + spec = good_spec(name=' F1 ', centre=(np.float64(5.0), 6, 7)) + written = json.loads(json.dumps(parametric_fault.clean_spec(spec))) + loaded = parametric_fault.spec_from_dict(written) + assert loaded['name'] == 'F1' + assert loaded['centre'] == (5.0, 6.0, 7.0) + assert loaded['dip'] == 60.0 + + def test_the_default_spec_is_in_the_middle_of_the_area(self): + spec = parametric_fault.default_spec([0, 0, -100], [1000, 400, 100]) + assert spec['centre'] == (500.0, 200.0, 0.0) + assert spec['major_axis'] == 1000.0 + assert spec['minor_axis'] == 500.0 + assert parametric_fault.problems({**spec, 'name': 'F'}) == [] diff --git a/tests/unit/test_preview.py b/tests/unit/test_preview.py new file mode 100644 index 0000000..4665200 --- /dev/null +++ b/tests/unit/test_preview.py @@ -0,0 +1,79 @@ +"""Pytest tests for the map preview of a feature (grid, levels, isolines). + +`loopstructural.main.preview` does not import QGIS, so the tests run in the +fast tests/unit/ job. +""" + +import numpy as np +import pytest + +pytest.importorskip('contourpy') + +from loopstructural.main import preview + + +class TestGrid: + def test_the_points_go_row_by_row(self): + x, y, points = preview.grid_points([0, 10, -5], [4, 20, 5], 3) + assert list(x) == [0, 2, 4] + assert list(y) == [10, 15, 20] + assert points.shape == (9, 3) + assert list(points[:3, 1]) == [10, 10, 10] + assert list(points[:3, 0]) == [0, 2, 4] + + def test_z_is_the_dem(self): + _, _, points = preview.grid_points([0, 0], [1, 1], 2, lambda x, y: 100 + x + 2 * y) + assert sorted(points[:, 2]) == [100, 101, 102, 103] + + def test_without_a_dem_z_is_zero(self): + _, _, points = preview.grid_points([0, 0], [1, 1], 2) + assert (points[:, 2] == 0).all() + + def test_a_grid_needs_two_points_on_each_side(self): + with pytest.raises(ValueError): + preview.grid_points([0, 0], [1, 1], 1) + + +class TestLevels: + def test_the_levels_are_inside_the_range(self): + levels = preview.default_levels([0.0, 10.0, np.nan], 4) + assert len(levels) == 4 + assert levels.min() > 0 and levels.max() < 10 + assert np.allclose(np.diff(levels), np.diff(levels)[0]) + + @pytest.mark.parametrize('values', [[], [np.nan], [5.0, 5.0]]) + def test_a_field_without_a_range_has_no_levels(self, values): + assert len(preview.default_levels(values)) == 0 + + +class TestIsolines: + def setup_method(self): + self.x = np.linspace(0, 10, 11) + self.y = np.linspace(0, 10, 11) + # a field that grows with X: the lines are straight and run along Y + self.values = np.tile(self.x, (len(self.y), 1)) + + def test_a_line_is_at_the_level(self): + lines = preview.isolines(self.x, self.y, self.values, [2.5]) + assert len(lines) == 1 + level, coordinates = lines[0] + assert level == 2.5 + assert np.allclose(coordinates[:, 0], 2.5) + assert coordinates[:, 1].min() == 0 and coordinates[:, 1].max() == 10 + + def test_a_level_outside_the_range_has_no_line(self): + assert preview.isolines(self.x, self.y, self.values, [50.0]) == [] + + def test_a_gap_stops_the_line(self): + values = self.values.copy() + values[:, 2:4] = np.nan + assert preview.isolines(self.x, self.y, values, [2.5]) == [] + + def test_the_shape_must_fit_the_grid(self): + with pytest.raises(ValueError): + preview.isolines(self.x, self.y, self.values.T[:, :5], [1.0]) + + def test_the_points_of_the_grid_make_the_same_lines(self): + x, y, points = preview.grid_points([0, 0], [10, 10], 11) + values = points[:, 0].reshape(len(y), len(x)) + assert len(preview.isolines(x, y, values, preview.default_levels(values, 3))) == 3