From ed17d931a04efd574eb6dec0abe4903dddc683b0 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 20 Aug 2026 14:30:37 +0930 Subject: [PATCH 1/7] test: pin training-value/isovalue direction agreement for stratigraphic columns Guards against the swap fixed in 814be12: model_manager.py's per-unit training value and LoopStructural's get_isovalues() must agree on which direction values increase, or extracted isosurfaces get labelled with the wrong unit while keeping correct geometry. Co-Authored-By: Claude Sonnet 5 --- .../test_stratigraphic_value_consistency.py | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 tests/qgis/test_stratigraphic_value_consistency.py diff --git a/tests/qgis/test_stratigraphic_value_consistency.py b/tests/qgis/test_stratigraphic_value_consistency.py new file mode 100644 index 0000000..5135a20 --- /dev/null +++ b/tests/qgis/test_stratigraphic_value_consistency.py @@ -0,0 +1,102 @@ +"""Regression test for the training-value / isovalue direction bug. + +`GeologicalModelManager.update_foliation_features` assigns a scalar `val` to +each unit's basal contact before handing the data to the interpolator. +`StratigraphicColumn.get_isovalues` (LoopStructural core) later decides which +name to stamp on each extracted isosurface, using its own idea of which +value belongs to which unit. + +These two must agree on direction (does value increase from oldest-to- +youngest, or youngest-to-oldest?), or every extracted surface gets labelled +with the wrong unit while keeping correct geometry -- see the "stratigraphic +column was reversed" fixes in model_manager.py (2025-07-21) and the widget +(2025-08-21, reverted 2025-09-08). This has flipped back and forth as this +plugin and LoopStructural evolved independently; this test pins the +invariant so a future change on either side fails loudly here instead of +silently inverting a user's model. +""" + +import pandas as pd +import pytest +from LoopStructural import StratigraphicColumn + +from loopstructural.main.model_manager import GeologicalModelManager + + +def _contact(unit_name): + """A minimal single-point basal contact, tagged with its unit name so + the test can recover which row came from which unit after the group + DataFrames get concatenated.""" + return pd.DataFrame({'X': [0.0], 'Y': [0.0], 'Z': [0.0], 'source_unit': [unit_name]}) + + +@pytest.fixture +def manager(monkeypatch): + manager = GeologicalModelManager() + + captured_calls = [] + + def fake_create_and_add_foliation(name, data=None, **kwargs): + captured_calls.append(data) + return object() # stand-in foliation, only passed back into add_unconformity + + 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_calls = captured_calls + return manager + + +class TestTrainingValueMatchesIsovalue: + def test_single_group_three_units(self, manager): + column = StratigraphicColumn() + column.clear(basement=False) # single flat group, no unconformities + column.add_unit(name='oldest', thickness=100.0, where='top') + column.add_unit(name='middle', thickness=200.0, where='top') + column.add_unit(name='youngest', thickness=300.0, where='top') + + manager.stratigraphic_column = column + for name in ('oldest', 'middle', 'youngest'): + manager.stratigraphy[name]['contact'] = _contact(name) + + manager.update_foliation_features() + + training_values = self._training_values_by_unit(manager._captured_calls) + expected_values = { + name: entry['value'] for name, entry in column.get_isovalues().items() + } + + for unit_name in ('oldest', 'middle', 'youngest'): + assert training_values[unit_name] == pytest.approx(expected_values[unit_name]), ( + f"'{unit_name}' was trained with val={training_values[unit_name]} but " + f"get_isovalues() will label the value={expected_values[unit_name]} surface " + f"with this unit's name -- the trained field and the isosurface labels " + f"disagree on direction, so extracted surfaces will get the wrong unit name." + ) + + def test_two_groups_split_by_unconformity(self, manager): + column = StratigraphicColumn() + column.clear(basement=False) + column.add_unit(name='basin_floor', thickness=50.0, where='top') + column.add_unit(name='basin_fill', thickness=150.0, where='top') + column.add_unconformity(name='regional_unconformity', where='top') + column.add_unit(name='cover_lower', thickness=80.0, where='top') + column.add_unit(name='cover_upper', thickness=120.0, where='top') + + manager.stratigraphic_column = column + for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): + manager.stratigraphy[name]['contact'] = _contact(name) + + manager.update_foliation_features() + + training_values = self._training_values_by_unit(manager._captured_calls) + expected_values = { + name: entry['value'] for name, entry in column.get_isovalues().items() + } + + for unit_name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): + assert training_values[unit_name] == pytest.approx(expected_values[unit_name]) + + @staticmethod + def _training_values_by_unit(captured_calls): + combined = pd.concat(captured_calls, ignore_index=True) + return dict(zip(combined['source_unit'], combined['val'])) From fbe4ce92ed08c3c128b6fc813d70e8ad47980c9a Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Thu, 20 Aug 2026 16:35:24 +0930 Subject: [PATCH 2/7] feat: allow a fault to act as a stratigraphic domain boundary Lets a stratigraphic-column unconformity be linked to an existing fault instead of a flat isovalue surface, so the fault's own (non-displacing) geometry splits the model into two domains -- built via LoopStructural's create_and_add_domain_fault, reusing the same trace data already ingested for the fault. - Stratigraphic column UI gains a "fault" boundary type with a fault picker; faults used this way are excluded from the fault topology's FAULTED/ABUTTING and fault-stratigraphy tables, since those assume a displacement-modelled fault. - The fault's trace is automatically extended to the model's bounding box edges along its own trend, and given synthetic strike/dip orientation constraints, so the interpolated surface spans and properly varies across the whole domain rather than only being reliable near the digitised trace. - A domain-boundary fault is skipped by the ordinary displacement-fault build loop, and any region a later unconformity incorrectly attaches to it is stripped after each build (defensive; the root cause is fixed upstream in LoopStructural core separately). - Fixes a project-load ordering bug where the model CRS was restored after the layers that get reprojected against it, silently skipping reprojection for any layer already in the project's own CRS. Co-Authored-By: Claude Sonnet 5 --- .../gui/modelling/fault_adjacency_tab.py | 40 ++- .../stratigraphic_column.py | 60 +++- .../stratigraphic_column/unconformity.py | 81 ++++- .../stratigraphic_column/unconformity.ui | 15 + loopstructural/main/data_manager.py | 139 ++++++++- loopstructural/main/model_manager.py | 292 +++++++++++++++++- tests/qgis/test_fault_domain_boundary.py | 260 ++++++++++++++++ 7 files changed, 837 insertions(+), 50 deletions(-) create mode 100644 tests/qgis/test_fault_domain_boundary.py diff --git a/loopstructural/gui/modelling/fault_adjacency_tab.py b/loopstructural/gui/modelling/fault_adjacency_tab.py index 6478f42..ffc910a 100644 --- a/loopstructural/gui/modelling/fault_adjacency_tab.py +++ b/loopstructural/gui/modelling/fault_adjacency_tab.py @@ -80,7 +80,7 @@ def _update(self, observable, event, *args, **kwargs): self.update_fault_adjacency_table() self.update_stratigraphic_units_table() - def change_button_color(self, button, row, col): + def change_button_color(self, button, fault1, fault2): """Cycle the button color and update the fault relationship.""" current_color = button.styleSheet() if "red" in current_color: @@ -94,15 +94,23 @@ def change_button_color(self, button, row, col): relationship = FaultRelationshipType.ABUTTING button.setStyleSheet(f"background-color: {new_color};") - f1 = self.data_manager._fault_topology.faults[row] - f2 = self.data_manager._fault_topology.faults[col] - self.data_manager._fault_topology.update_fault_relationship(f1, f2, relationship) + self.data_manager._fault_topology.update_fault_relationship(fault1, fault2, relationship) + + def _displacement_fault_names(self): + """Fault names to show in these tables, excluding faults used as + stratigraphic-column domain boundaries (see `set_fault_boundary`). + Those are non-displacing splits, not faults that cut/abut other + faults or offset stratigraphic units, so FAULTED/ABUTTING and + fault-stratigraphy relationships don't apply to them. + """ + domain_boundary_faults = self.data_manager.get_fault_boundary_fault_names() + return [ + f for f in self.data_manager._fault_topology.faults if f not in domain_boundary_faults + ] def update_fault_adjacency_table(self): """Update the fault adjacency table with QPushButtons.""" - faults = ( - self.data_manager._fault_topology.faults - ) # Assuming faults is a list of fault names + faults = self._displacement_fault_names() if not faults: self.fault_table_group.hide() return @@ -145,15 +153,15 @@ def update_fault_adjacency_table(self): else: button.setStyleSheet("background-color: white;") button.clicked.connect( - lambda _, b=button, r=row, c=col: self.change_button_color(b, r, c) + lambda _, b=button, f1=faults[row], f2=faults[ + col + ]: self.change_button_color(b, f1, f2) ) self.table.setCellWidget(row, col, button) def update_stratigraphic_units_table(self): """Update the stratigraphic units table with QPushButtons.""" - faults = ( - self.data_manager._fault_topology.faults - ) # Assuming faults is a list of fault names + faults = self._displacement_fault_names() group_units_pairs = self.data_manager._stratigraphic_column.get_group_unit_pairs() if not faults or not group_units_pairs: @@ -185,11 +193,13 @@ def update_stratigraphic_units_table(self): # Default to white if no relationship or not faulted button.setStyleSheet("background-color: white;") button.clicked.connect( - lambda _, b=button, r=row, c=col: self.change_button_colour_binary(b, r, c) + lambda _, b=button, u=units[row], f=faults[ + col + ]: self.change_button_colour_binary(b, u, f) ) self.stratigraphic_table.setCellWidget(row, col, button) - def change_button_colour_binary(self, button, row, col): + def change_button_colour_binary(self, button, unit_name, fault_name): """Cycle the button color between red, green, and black.""" current_color = button.styleSheet() @@ -199,8 +209,6 @@ def change_button_colour_binary(self, button, row, col): else: button.setStyleSheet("background-color: red;") flag = True - fault = self.data_manager._fault_topology.faults[col] - unit = self.data_manager._stratigraphic_column.get_group_unit_pairs()[row] self.data_manager._fault_topology.update_fault_stratigraphy_relationship( - unit[1], fault, flag + unit_name, fault_name, flag ) diff --git a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py index 0f1131b..27c66d6 100644 --- a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py +++ b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py @@ -232,7 +232,11 @@ def update_display(self): widget, _ = self._widget_cache[unit.uuid] # Update widget data without rebuilding if hasattr(widget, 'setData'): - widget.setData(unit.to_dict()) + unit_data = unit.to_dict() + if isinstance(widget, UnconformityWidget): + unit_data = self._enrich_unconformity_data(unit_data) + widget.set_available_faults(self._get_available_fault_names()) + widget.setData(unit_data) return # If order/content differs, do a full rebuild @@ -263,7 +267,31 @@ def _full_rebuild_display(self, current_order): if unit.element_type == StratigraphicColumnElementType.UNIT: self.add_unit(unit_data=unit.to_dict(), create_new=False) elif unit.element_type == StratigraphicColumnElementType.UNCONFORMITY: - self.add_unconformity(unconformity_data=unit.to_dict(), create_new=False) + self.add_unconformity( + unconformity_data=self._enrich_unconformity_data(unit.to_dict()), + create_new=False, + ) + + def _enrich_unconformity_data(self, unconformity_data): + """Merge in the plugin-side fault-boundary link for an unconformity row. + + `StratigraphicUnconformity.to_dict()` (core) only knows `erode`/ + `onlap`; the fault link is tracked separately in the data manager + (see `ModellingDataManager.set_fault_boundary`), so it has to be + folded in here for display. + """ + fault_name = self.data_manager.get_fault_boundary(unconformity_data.get('uuid')) + if fault_name: + unconformity_data = dict(unconformity_data) + unconformity_data['unconformity_type'] = 'fault' + unconformity_data['fault_name'] = fault_name + return unconformity_data + + def _get_available_fault_names(self): + """Fault names offered when marking an unconformity as a domain boundary.""" + if not self.data_manager: + return [] + return list(self.data_manager._fault_topology.faults) def init_stratigraphic_column_from_basal_contacts(self): if self.data_manager: @@ -482,11 +510,13 @@ def add_unconformity(self, *, unconformity_data=None, create_new=True): widget, _ = self._widget_cache[unconformity.uuid] # Just update the data, don't recreate the widget if hasattr(widget, 'setData'): + widget.set_available_faults(self._get_available_fault_names()) widget.setData(unconformity_data) return unconformity_widget = UnconformityWidget(uuid=unconformity.uuid) unconformity_widget.deleteRequested.connect(self.delete_unit) + unconformity_widget.dataChanged.connect(lambda: self.update_element(unconformity_widget)) unconformity_widget.dragHandlePressed.connect( lambda: self._on_drag_start(unconformity_widget) ) @@ -500,6 +530,8 @@ def add_unconformity(self, *, unconformity_data=None, create_new=True): item.setSizeHint(unconformity_widget.sizeHint()) self.unitList.addItem(item) self.unitList.setItemWidget(item, unconformity_widget) + unconformity_widget.set_available_faults(self._get_available_fault_names()) + unconformity_widget.setData(unconformity_data) # Cache the widget for efficient updates self._widget_cache[unconformity.uuid] = (unconformity_widget, item) @@ -605,6 +637,30 @@ def update_element(self, unit_widget): """ if self.data_manager: unit_data = unit_widget.getData() + if isinstance(unit_widget, UnconformityWidget): + fault_name = unit_data.pop('fault_name', None) + is_fault_boundary = unit_data.get('unconformity_type') == 'fault' + if is_fault_boundary: + # The core stratigraphic column only knows erode/onlap -- + # the fault link lives in the data manager's side table + # (see set_fault_boundary), so store it as a plain + # erosional boundary here. + unit_data['unconformity_type'] = 'erode' + if is_fault_boundary and fault_name: + self.data_manager.set_fault_boundary(unit_widget.uuid, fault_name) + if not self.data_manager.fault_spans_model_domain(fault_name): + QMessageBox.information( + self, + "Fault Domain Boundary", + f"Fault '{fault_name}' does not reach every edge of the model " + "bounding box.\n\nA fault used as a domain boundary crops the " + "whole model, so its digitised trace will automatically be " + "extended out to the domain edges along its overall trend when " + "the model is built. For best results the trace should still " + "roughly follow the fault's real direction across the gap.", + ) + else: + self.data_manager.clear_fault_boundary(unit_widget.uuid) self.data_manager._stratigraphic_column.update_element(unit_data) # Trigger callback to notify all listeners of the change if self.data_manager.stratigraphic_column_callback: diff --git a/loopstructural/gui/modelling/stratigraphic_column/unconformity.py b/loopstructural/gui/modelling/stratigraphic_column/unconformity.py index 115b9df..a828edd 100644 --- a/loopstructural/gui/modelling/stratigraphic_column/unconformity.py +++ b/loopstructural/gui/modelling/stratigraphic_column/unconformity.py @@ -10,6 +10,7 @@ class UnconformityWidget(QWidget): deleteRequested = pyqtSignal(QWidget) # Signal to request deletion + dataChanged = pyqtSignal() # Type or fault-name changed dragHandlePressed = pyqtSignal() # Drag handle mouse-down dragHandleMoved = pyqtSignal(QPoint) # Drag handle mouse-move (global pos) dragHandleReleased = pyqtSignal() # Drag handle mouse-up @@ -25,9 +26,9 @@ def __init__( self.buttonDelete.clicked.connect(self.request_delete) self.uuid = uuid self.unconformity_type = 'erode' - # self.comboBoxUnconformityType.currentIndexChanged.connect( - # lambda: setattr(self, 'unconformity_type', self.comboBoxUnconformityType.currentText()) - # ) + self.fault_name = None + self.comboBoxUnconformityType.currentIndexChanged.connect(self._on_type_changed) + self.comboBoxFaultName.currentIndexChanged.connect(self._on_fault_name_changed) # The row's combo box/buttons cover the whole widget, so a QListWidget's # built-in drag-and-drop can never see a mouse press to start a # reorder. Route presses on the dedicated grip label through here instead. @@ -55,21 +56,73 @@ def request_delete(self): self.deleteRequested.emit(self) + def _on_type_changed(self, _index): + self.unconformity_type = self.comboBoxUnconformityType.currentText() + self.comboBoxFaultName.setVisible(self.unconformity_type == 'fault') + if self.unconformity_type == 'fault': + self.fault_name = self.comboBoxFaultName.currentText() or None + else: + self.fault_name = None + self.dataChanged.emit() + + def _on_fault_name_changed(self, _index): + if self.unconformity_type != 'fault': + return + self.fault_name = self.comboBoxFaultName.currentText() or None + self.dataChanged.emit() + + def set_available_faults(self, fault_names): + """Populate the fault-name picker, keeping the current selection if + it is still available (e.g. after the fault trace layer changes). + """ + fault_names = list(fault_names or []) + if [ + self.comboBoxFaultName.itemText(i) for i in range(self.comboBoxFaultName.count()) + ] == fault_names: + return + self.comboBoxFaultName.blockSignals(True) + try: + self.comboBoxFaultName.clear() + self.comboBoxFaultName.addItems(fault_names) + if self.fault_name and self.fault_name in fault_names: + self.comboBoxFaultName.setCurrentText(self.fault_name) + finally: + self.comboBoxFaultName.blockSignals(False) + def setData(self, data: Optional[dict] = None): """Set the data for the unconformity widget. Parameters ---------- data : dict or None - Dictionary containing 'unconformity_type' key. If None, defaults are used. + Dictionary with an 'unconformity_type' key ('erode', 'onlap' or + 'fault'), and a 'fault_name' key when the type is 'fault'. If + None, defaults are used. """ - if data: - self.unconformity_type = data.get("unconformity_type", "") - # self.unconformityTypeComboBox.setCurrentIndex( - # self.unconformityTypeComboBox.findText(self.unconformity_type) - # ) - else: - self.unconformity_type = 'erode' - # self.unconformityTypeComboBox.setCurrentIndex( - # self.unconformityTypeComboBox.findText(self.unconformity_type) - # ) + self.unconformity_type = (data or {}).get("unconformity_type", "erode") + self.fault_name = ( + (data or {}).get("fault_name") if self.unconformity_type == 'fault' else None + ) + + self.comboBoxUnconformityType.blockSignals(True) + self.comboBoxFaultName.blockSignals(True) + try: + index = self.comboBoxUnconformityType.findText(self.unconformity_type) + if index >= 0: + self.comboBoxUnconformityType.setCurrentIndex(index) + self.comboBoxFaultName.setVisible(self.unconformity_type == 'fault') + if self.fault_name: + self.comboBoxFaultName.setCurrentText(self.fault_name) + finally: + self.comboBoxUnconformityType.blockSignals(False) + self.comboBoxFaultName.blockSignals(False) + + def getData(self): + """Return this row's data for the data manager: uuid, unconformity_type + and (when the boundary is fault-linked) fault_name. + """ + return { + 'uuid': self.uuid, + 'unconformity_type': self.unconformity_type, + 'fault_name': self.fault_name, + } diff --git a/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui b/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui index cda9c0f..48239f4 100644 --- a/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui +++ b/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui @@ -61,6 +61,11 @@ onlap + + + fault + + @@ -70,6 +75,16 @@ + + + + Fault whose surface realises this boundary as a domain split + + + false + + + diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py index c91a4c8..d47703d 100644 --- a/loopstructural/main/data_manager.py +++ b/loopstructural/main/data_manager.py @@ -89,6 +89,13 @@ def __init__(self, *, project=None, mapCanvas=None, logger=None): self.logger = logger self._stratigraphic_column = StratigraphicColumn() self._fault_topology = FaultTopology(self._stratigraphic_column) + # Maps a stratigraphic-column unconformity's uuid to the name of an + # existing fault that should realise that boundary, instead of a + # flat isovalue surface -- see `set_fault_boundary`. Kept as a plain + # dict (not part of `StratigraphicColumn` itself) so this stays a + # plugin-side concept for now; passed by reference to the model + # manager the same way `_stratigraphic_column`/`_fault_topology` are. + self._fault_boundaries: dict[str, str] = {} self._model_manager = None self.bounding_box_callback = None self.basal_contacts_callback = None @@ -137,6 +144,7 @@ def set_model_manager(self, model_manager): self._model_manager = model_manager self._model_manager.set_stratigraphic_column(self._stratigraphic_column) self._model_manager.set_fault_topology(self._fault_topology) + self._model_manager.set_fault_boundaries(self._fault_boundaries) self._model_manager.update_bounding_box(self._bounding_box) def set_bounding_box( @@ -557,10 +565,103 @@ def add_to_stratigraphic_column(self, unit_data): def remove_from_stratigraphic_column(self, unit_uuid): """Remove a unit or unconformity from the stratigraphic column.""" self._stratigraphic_column.remove_unit(uuid=unit_uuid) + self._fault_boundaries.pop(unit_uuid, None) self.update_stratigraphy() if self.stratigraphic_column_callback: self.stratigraphic_column_callback() + def set_fault_boundary(self, unconformity_uuid, fault_name): + """Mark a stratigraphic-column unconformity as realised by an + existing fault instead of a flat isovalue surface. + + Parameters + ---------- + unconformity_uuid : str + uuid of the `StratigraphicUnconformity` element in the column. + fault_name : str + Name of an existing fault (as known to `_fault_topology`) whose + surface should be used as the domain boundary at this point in + the column. + """ + self._fault_boundaries[unconformity_uuid] = fault_name + if self.stratigraphic_column_callback: + self.stratigraphic_column_callback() + + def clear_fault_boundary(self, unconformity_uuid): + """Undo `set_fault_boundary`, reverting the unconformity to a plain isovalue boundary.""" + if self._fault_boundaries.pop(unconformity_uuid, None) is not None: + if self.stratigraphic_column_callback: + self.stratigraphic_column_callback() + + def get_fault_boundary(self, unconformity_uuid): + """Return the fault name linked to this unconformity, or None.""" + return self._fault_boundaries.get(unconformity_uuid) + + def get_fault_boundaries(self): + """Return the uuid -> fault_name mapping of all fault-linked boundaries.""" + return dict(self._fault_boundaries) + + def get_fault_boundary_fault_names(self): + """Return the set of fault names currently used as domain boundaries. + + These faults are built as non-displacing domain splits (see + `GeologicalModelManager.update_foliation_features`), so they should + not also be offered in fault-fault (FAULTED/ABUTTING) or + fault-stratigraphy relationship editors, which assume a + displacement-modelled fault. + """ + return set(self._fault_boundaries.values()) + + def fault_spans_model_domain(self, fault_name, *, tolerance=0.0): + """Check whether a fault's trace data reaches every edge of the + model's XY bounding box. + + A fault used as a domain boundary crops the *entire* model on + either side of its interpolated surface (see + `LoopStructural.modelling.core._model_relationships`), so unlike an + ordinary local fault trace it needs to be constrained across the + whole domain -- otherwise the interpolator extrapolates the crop + surface into areas with no supporting data. Returns True if no + fault trace data is available yet (nothing to check against). + + Parameters + ---------- + fault_name : str + Name of the fault to check, as found in `get_fault_traces()`'s layer. + tolerance : float, optional + Allowed gap, in model units, between the trace's extent and the + bounding box edge before it is considered "not spanning". + """ + if self._fault_traces is None or self._fault_traces['layer'] is None: + return True + layer = self._fault_traces['layer'] + name_field = self._fault_traces['fault_name_field'] + trace_extent = None + for feature in layer.getFeatures(): + if name_field is not None and str(feature[name_field]) != str(fault_name): + continue + geom = feature.geometry() + if geom is None or geom.isEmpty(): + continue + bbox = geom.boundingBox() + trace_extent = bbox if trace_extent is None else trace_extent.combineExtentWith(bbox) + if trace_extent is None: + return True + xmin, ymin = self._bounding_box.origin[0], self._bounding_box.origin[1] + xmax, ymax = self._bounding_box.maximum[0], self._bounding_box.maximum[1] + # A boundary only needs to fully cross the domain along one axis (an + # east-west or north-south cut) to split the whole model -- it does + # not need to cover the full bounding box in both directions. + spans_x = ( + trace_extent.xMinimum() <= xmin + tolerance + and trace_extent.xMaximum() >= xmax - tolerance + ) + spans_y = ( + trace_extent.yMinimum() <= ymin + tolerance + and trace_extent.yMaximum() >= ymax - tolerance + ) + return spans_x or spans_y + def update_stratigraphic_column_order(self, new_order): """Update the order of units in the stratigraphic column.""" if not isinstance(new_order, list): @@ -956,6 +1057,7 @@ def to_dict(self): 'stratigraphic_column': ( self._stratigraphic_column.to_dict() if self._stratigraphic_column else None ), + 'fault_boundaries': dict(self._fault_boundaries), 'dem_layer': dem_layer_name if self.dem_layer else None, 'use_dem': self.use_dem, 'elevation': self.elevation, @@ -998,6 +1100,9 @@ def from_dict(self, data): if 'stratigraphic_column' in data: self._stratigraphic_column = StratigraphicColumn.from_dict(data['stratigraphic_column']) self.stratigraphic_column_callback() + self._fault_boundaries.clear() + if data.get('fault_boundaries'): + self._fault_boundaries.update(data['fault_boundaries']) if 'widget_settings' in data: self.widget_settings = data['widget_settings'] @@ -1011,6 +1116,23 @@ def from_dict(self, data): def update_from_dict(self, data): """Update the data manager from a dictionary.""" + # Model CRS must be restored before anything below that reprojects + # a layer against it (basal_contacts/fault_traces/ + # structural_orientations, all via get_model_crs()) -- restoring it + # last meant every reprojection during project load ran against + # whatever `_use_project_crs`/`_model_crs` still held from __init__ + # (use_project_crs=True, i.e. the *project's* CRS) instead of the + # saved model CRS, silently skipping reprojection whenever a layer + # happened to already be in the project's CRS but not the model's. + if 'use_project_crs' in data: + self._use_project_crs = data['use_project_crs'] + else: + self._use_project_crs = True + if 'model_crs' in data and data['model_crs'] is not None: + crs = QgsCoordinateReferenceSystem(data['model_crs']) + if crs.isValid(): + self.set_model_crs(crs, use_project_crs=self._use_project_crs) + if 'bounding_box' in data: self.set_bounding_box( xmin=data['bounding_box']['origin'][0], @@ -1080,22 +1202,17 @@ def update_from_dict(self, data): else: self._stratigraphic_column.clear() + # Mutate in place rather than reassign: `_fault_boundaries` is + # shared by reference with the model manager (see set_model_manager). + self._fault_boundaries.clear() + if data.get('fault_boundaries'): + self._fault_boundaries.update(data['fault_boundaries']) + if 'widget_settings' in data: self.widget_settings = data['widget_settings'] else: self.widget_settings = {} - # Load model CRS settings - if 'use_project_crs' in data: - self._use_project_crs = data['use_project_crs'] - else: - self._use_project_crs = True - - if 'model_crs' in data and data['model_crs'] is not None: - crs = QgsCoordinateReferenceSystem(data['model_crs']) - if crs.isValid(): - self.set_model_crs(crs, use_project_crs=self._use_project_crs) - if self.stratigraphic_column_callback: self.stratigraphic_column_callback() diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 779689f..afa3798 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -18,8 +18,11 @@ import pandas as pd from LoopStructural.datatypes import BoundingBox from LoopStructural.modelling.core.fault_topology import FaultRelationshipType -from LoopStructural.modelling.core.stratigraphic_column import StratigraphicColumn -from LoopStructural.modelling.features import FeatureType, StructuralFrame +from LoopStructural.modelling.core.stratigraphic_column import ( + StratigraphicColumn, + StratigraphicUnconformity, +) +from LoopStructural.modelling.features import FeatureType, StructuralFrame, UnconformityFeature from LoopStructural.modelling.features.fold import FoldFrame from LoopStructural.utils.observer import Observable @@ -111,6 +114,10 @@ def __init__(self, debug_manager=None): self.stratigraphy: Dict[str, StratigraphyEntry] = defaultdict(dict) self.stratigraphic_column = None self.fault_topology = None + # uuid (of a StratigraphicUnconformity in stratigraphic_column) -> + # fault name; see `set_fault_boundaries`. Shared by reference with + # ModellingDataManager, same as stratigraphic_column/fault_topology. + self.fault_boundaries: Dict[str, str] = {} # Observers managed by Observable base class self.dem_function = lambda x, y: 0 # internal flag to temporarily suppress notifications (used when @@ -241,6 +248,12 @@ def set_fault_topology(self, fault_topology): except Exception: pass + def set_fault_boundaries(self, fault_boundaries: Dict[str, str]): + """Set the uuid -> fault_name mapping of fault-linked stratigraphic + column boundaries (see `ModellingDataManager.set_fault_boundary`). + """ + self.fault_boundaries = fault_boundaries + # Topology events that change what actually feeds the interpolator (as # opposed to just which side of an already-solved fault gets cropped # away) and therefore need Initialize Model to be re-run before Solve @@ -531,6 +544,201 @@ def update_stratigraphic_column(self, stratigraphic_column: StratigraphicColumn) # def update_stratigraphic_unit(self, unit_data): # self.data + def _closing_fault_boundary(self, group): + """Return the fault name that closes `group` from above, if any. + + The boundary "closing" a group off from the next (younger) group is + the first `StratigraphicUnconformity` above the group's youngest + unit in `stratigraphic_column.order`. If that boundary has been + linked to a fault (`ModellingDataManager.set_fault_boundary`), the + group should be capped by that fault's surface -- a non-displacing + domain split, see `create_and_add_domain_fault` -- instead of the + flat isovalue-0 surface `add_unconformity` uses. + """ + if not group.units or not self.fault_boundaries or self.stratigraphic_column is None: + return None + order = self.stratigraphic_column.order + youngest_uuid = group.units[0].uuid + start = next((i for i, e in enumerate(order) if e.uuid == youngest_uuid), None) + if start is None: + return None + for element in order[start + 1 :]: + if isinstance(element, StratigraphicUnconformity): + return self.fault_boundaries.get(element.uuid) + return None + + def _clip_line_to_bounding_box(self, centroid, direction): + """Return the (t_min, t_max) range along `centroid + t*direction` + (XY only) that stays within the model's bounding box, or None if + the line never crosses it. Standard slab-method line/box clip. + """ + origin_xy = np.array(self.model.bounding_box.origin[:2], dtype=float) + maximum_xy = np.array(self.model.bounding_box.maximum[:2], dtype=float) + t_min, t_max = -np.inf, np.inf + for axis in (0, 1): + d = direction[axis] + if abs(d) < 1e-12: + if centroid[axis] < origin_xy[axis] or centroid[axis] > maximum_xy[axis]: + return None + continue + t0 = (origin_xy[axis] - centroid[axis]) / d + t1 = (maximum_xy[axis] - centroid[axis]) / d + t0, t1 = min(t0, t1), max(t0, t1) + t_min = max(t_min, t0) + t_max = min(t_max, t1) + if t_min > t_max: + return None + return t_min, t_max + + def _extend_fault_trace_to_domain(self, fault_data): + """Add two synthetic points that extend a fault's trace out to the + edges of the model's bounding box along its own overall trend. + + `create_and_add_domain_fault` interpolates a scalar field only from + the points it is given, over the model's exact bounding box (no + buffer, unlike a displacement fault's mesh) -- so a locally + digitised trace only reliably constrains the surface near itself, + and the domain crop can wander unpredictably further away. Fitting + a line through the existing XY points and adding two constraint + points where that line meets the bounding box edges keeps the + interpolated surface following the fault's actual trend all the + way across the domain, rather than an arbitrary extrapolation -- + this is what makes the fault behave as an "infinite" domain + boundary rather than a locally-anchored patch. + + Z at each synthetic point is extrapolated linearly against + distance along that line, so a dipping trace keeps its dip. + """ + xy = fault_data[['X', 'Y']].to_numpy() + if len(xy) < 2: + return fault_data + centroid = xy.mean(axis=0) + # Principal direction of the trace via SVD -- robust to a + # near-vertical (large-Y-range, small-X-range) trace, unlike a + # simple polyfit of Y against X. + _, _, vt = np.linalg.svd(xy - centroid) + direction = vt[0] + clipped = self._clip_line_to_bounding_box(centroid, direction) + if clipped is None: + return fault_data + t_min, t_max = clipped + projections = (xy - centroid) @ direction + z = fault_data['Z'].to_numpy() + if len(np.unique(projections)) > 1: + z_slope, z_intercept = np.polyfit(projections, z, 1) + else: + z_slope, z_intercept = 0.0, float(np.mean(z)) + new_rows = [] + for t in (t_min, t_max): + point_xy = centroid + direction * t + new_rows.append({'X': point_xy[0], 'Y': point_xy[1], 'Z': z_slope * t + z_intercept}) + return pd.concat([fault_data, pd.DataFrame(new_rows)], ignore_index=True) + + def _domain_fault_orientation_rows(self, points_xyz, fault_entry): + """Build strike/dip orientation constraint rows for a domain-boundary fault. + + A domain fault built only from same-valued (val=0) points has no + information about which direction the field should vary -- the + minimum-curvature solution that exactly satisfies "value=0 along + this line" with nothing else to go on is a trivial constant/flat + field (zero curvature everywhere, zero misfit). That's + geometrically useless: the crop condition + `domain_fault.evaluate_value(pos) > 0` + (see `LoopStructural.modelling.core._model_relationships`) is then + never true anywhere, so everything on the "positive" side reads as + NaN. Confirmed against a live project: a fault built from 31 trace + points and 2 extension points, all val=0, interpolated to an + exactly flat 0.0 field everywhere; adding one orientation + constraint per point produced a properly varying field crossing + zero along the fault's own trend. + + Strike comes from the trace's own principal direction (matching + `_extend_fault_trace_to_domain`'s line fit); dip comes from the + ingested fault trace data's `dip` column if present (matching how + `update_fault_features` picks up dip for a displacement fault), + otherwise defaults to vertical (90 degrees). + """ + xy = points_xyz[['X', 'Y']].to_numpy() + if len(xy) < 2: + return None + centroid = xy.mean(axis=0) + _, _, vt = np.linalg.svd(xy - centroid) + direction = vt[0] + # strikedip2vector's strike is degrees clockwise from North (+Y); + # atan2(dx, dy) matches that convention directly. + strike = float(np.degrees(np.arctan2(direction[0], direction[1])) % 360) + dip = 90.0 + raw_data = fault_entry.get('data') if fault_entry else None + if raw_data is not None and 'dip' in raw_data: + dip_values = raw_data['dip'].dropna() + if not dip_values.empty: + dip = float(dip_values.mean()) + rows = points_xyz[['X', 'Y', 'Z']].copy() + rows['strike'] = strike + rows['dip'] = dip + return rows + + def _build_domain_fault_boundary(self, fault_name, groupname): + """Build `fault_name` as a domain-fault boundary in place of a flat unconformity. + + `GeologicalModel.create_and_add_domain_fault` (unlike + `create_and_add_fault`/`create_and_add_foliation`) has no `data=` + parameter -- it always reads the fault's points from `model.data` + filtered by `feature_name`, so the fault's trace data is registered + there first. Registration replaces any rows already tagged with + this fault name so repeated Initialize Model runs stay idempotent + instead of accumulating duplicate points on every rebuild. + + Returns True if the domain fault was built, False if it fell back + to a flat unconformity for lack of trace data (caller should then + call `self.model.add_unconformity` itself). + """ + fault_entry = self.faults.get(fault_name) + fault_data = fault_entry.get('data') if fault_entry else None + if fault_data is None or fault_data.empty: + self._debug_manager and self._debug_manager.log( + f"Fault '{fault_name}' is linked as a domain boundary for group " + f"'{groupname}' but has no trace data; using a flat unconformity instead.", + log_level=2, + ) + return False + data_for_fault = self._extend_fault_trace_to_domain(fault_data[['X', 'Y', 'Z']].copy()) + orientation_rows = self._domain_fault_orientation_rows(data_for_fault, fault_entry) + data_for_fault['feature_name'] = fault_name + data_for_fault['val'] = 0 + if orientation_rows is not None: + orientation_rows['feature_name'] = fault_name + orientation_rows['val'] = np.nan + data_for_fault = pd.concat([data_for_fault, orientation_rows], ignore_index=True) + # Unlike create_and_add_foliation/create_and_add_fault (which + # normalise their own `data=` argument internally via + # model.prepare_data before building), create_and_add_domain_fault + # reads straight from model.data with no normalisation -- it + # expects every standard column (gx/gy/gz/nx/ny/nz/...) to already + # be present, or the interpolator crashes looking them up. Run it + # through prepare_data ourselves before writing it in. + data_for_fault = self.model.prepare_data(data_for_fault, include_feature_name=True) + existing_data = self.model.data + if existing_data is not None and not existing_data.empty: + existing_data = existing_data.loc[existing_data['feature_name'] != fault_name] + self.model.data = pd.concat([existing_data, data_for_fault], ignore_index=True) + else: + # An empty placeholder frame built with `pd.DataFrame(columns=...)` + # defaults every column to object dtype; concatenating that with + # `data_for_fault`'s float columns can leave the result as + # object dtype too, and `add_data_to_interpolator` then fails + # calling `np.isnan` on an object-dtype column. Assign directly + # instead of concatenating with a dtype-less placeholder. + self.model.data = data_for_fault + self.model.create_and_add_domain_fault( + fault_name, + nelements=PlgSettingsStructure.interpolator_nelements, + npw=PlgSettingsStructure.interpolator_npw, + cpw=PlgSettingsStructure.interpolator_cpw, + regularisation=PlgSettingsStructure.interpolator_regularisation, + ) + return True + def update_foliation_features(self): """Builds the stratigraphic feature from the stratigraphic column data and the basal contacts and structural orientations data. @@ -583,11 +791,51 @@ def update_foliation_features(self): cpw=PlgSettingsStructure.interpolator_cpw, regularisation=PlgSettingsStructure.interpolator_regularisation, ) - self.model.add_unconformity(foliation, 0) + fault_name = self._closing_fault_boundary(group) + if fault_name is None or not self._build_domain_fault_boundary(fault_name, groupname): + self.model.add_unconformity(foliation, 0) + self._strip_spurious_regions_from_domain_faults() self.model.stratigraphic_column = self.stratigraphic_column # foliation features were rebuilt; let observers know self._emit('foliation_features_updated') + def _strip_spurious_regions_from_domain_faults(self): + """Work around a LoopStructural core gap that corrupts a domain + fault's own scalar field. + + `add_unconformity`'s backward crop walk (in LoopStructural's + `_model_relationships.FeatureRelationshipManager.add_unconformity`) + only recognises `FeatureType.FAULT`/`INACTIVEFAULT` as "already + handled, skip" -- it doesn't know about `FeatureType.DOMAINFAULT`. + So whenever a later group in the same Initialize Model run falls + back to a plain `add_unconformity` (no fault linked to its + boundary), that call walks straight through any domain-boundary + fault built earlier and incorrectly adds itself as a region on it. + Since `GeologicalFeature.evaluate_value` returns NaN wherever a + feature's regions don't hold, the domain fault's own field then + reads as NaN on whichever side that unrelated unconformity's + condition fails -- "NaN on one side" of an otherwise valid domain + boundary. + + A domain fault is meant to crop other features, not be cropped + itself -- except by an *earlier* domain fault, which is a + legitimate, intentional cascade (`_add_domain_fault_above` adds + that as a plain lambda region, not an `UnconformityFeature`). So + strip only the `UnconformityFeature`-typed regions injected by the + bug, leaving any real domain-fault-vs-domain-fault crop intact. + """ + for feature in self.model.features: + if getattr(feature, 'type', None) != FeatureType.DOMAINFAULT: + continue + kept = [r for r in feature.regions if not isinstance(r, UnconformityFeature)] + if len(kept) != len(feature.regions): + self._debug_manager and self._debug_manager.log( + f"Removing {len(feature.regions) - len(kept)} unconformity region(s) " + f"incorrectly applied to domain-boundary fault '{feature.name}'.", + log_level=2, + ) + feature.regions = kept + def _report_progress(self, message: str): """Report progress on a long-running model update. @@ -685,7 +933,13 @@ def _cutting_faults_for(self, fault_name): def update_fault_features(self): """Update the fault features in the geological model.""" + domain_boundary_faults = set(self.fault_boundaries.values()) for fault_name in self._fault_build_order(): + if fault_name in domain_boundary_faults: + # Built as a non-displacing domain fault in + # update_foliation_features instead -- see + # `_build_domain_fault_boundary`. + continue fault_data = self.faults[fault_name] self._report_progress(f"Building fault '{fault_name}'") if qgisAttributeIsNone(fault_name): @@ -736,6 +990,19 @@ def update_fault_features(self): ) self.apply_fault_abutting_relationships() + def _get_feature_by_name_or_none(self, name): + """Non-raising counterpart to `GeologicalModel.get_feature_by_name`. + + The wrapped `GeologicalModel` raises `ValueError` for a name that + hasn't been built yet rather than returning None, which several + call sites in this manager treat as "not built yet, skip it" -- + `__contains__` (`name in self.model`) checks `feature_name_index` + directly, so this restores that non-raising lookup. + """ + if name not in self.model: + return None + return self.model.get_feature_by_name(name) + def apply_fault_abutting_relationships(self): """Re-apply fault-fault ABUTTING relationships as region crops on the already-built fault surfaces, for every pair currently in the topology. @@ -756,20 +1023,30 @@ def apply_fault_abutting_relationships(self): """ if self.fault_topology is None: return + # Domain-boundary faults (see `_build_domain_fault_boundary`) are + # non-displacing GeologicalFeatures, not FaultSegments, and are only + # built later in `update_foliation_features` -- calling this before + # that has run would look up a feature that doesn't exist yet. + # ABUTTING relationships don't apply to them either way (that's a + # FaultSegment-specific crop), matching their exclusion from the + # Fault Adjacency tab. + domain_boundary_faults = set(self.fault_boundaries.values()) for f in self.fault_topology.faults: - fault_feature = self.model.get_feature_by_name(f) + if f in domain_boundary_faults: + continue + fault_feature = self._get_feature_by_name_or_none(f) if fault_feature is None or not hasattr(fault_feature, 'abut'): continue coord0 = fault_feature.__getitem__(0) for f2 in self.fault_topology.faults: - if f == f2: + if f == f2 or f2 in domain_boundary_faults: continue relationship = self.fault_topology.get_fault_relationship(f, f2) existing_region = fault_feature.abut.get(f2) if relationship is FaultRelationshipType.ABUTTING: if existing_region is not None: continue # already cropped against f2 - f2_feature = self.model.get_feature_by_name(f2) + f2_feature = self._get_feature_by_name_or_none(f2) if f2_feature is None: continue # Determine which side of f2 to keep ourselves, ignoring any @@ -913,7 +1190,8 @@ def update_model( len(self.stratigraphic_column.get_groups()) if self.stratigraphic_column else 0 ) self._progress_callback = progress_callback - self._progress_total = len(self.faults) + group_count + displacement_fault_count = len(set(self.faults) - set(self.fault_boundaries.values())) + self._progress_total = displacement_fault_count + group_count self._progress_current = 0 dbg = getattr(self, '_debug_manager', None) if dbg is not None: diff --git a/tests/qgis/test_fault_domain_boundary.py b/tests/qgis/test_fault_domain_boundary.py new file mode 100644 index 0000000..df608b4 --- /dev/null +++ b/tests/qgis/test_fault_domain_boundary.py @@ -0,0 +1,260 @@ +"""Regression tests for using a fault as a stratigraphic domain boundary. + +`GeologicalModelManager.update_foliation_features` normally closes each +built stratigraphic group with a flat isovalue unconformity +(`GeologicalModel.add_unconformity`). When the column boundary closing a +group has been linked to a fault (via `fault_boundaries`, populated through +`ModellingDataManager.set_fault_boundary`), the group should instead be +capped by that fault's own surface -- a non-displacing domain split, built +with `GeologicalModel.create_and_add_domain_fault` -- and that fault must be +skipped by the ordinary displacement-fault build loop +(`update_fault_features`), since a domain fault and a displacement fault are +mutually exclusive roles for the same fault name. +""" + +import numpy as np +import pandas as pd +import pytest +from LoopStructural import FaultTopology, StratigraphicColumn +from LoopStructural.datatypes import BoundingBox + +from loopstructural.main.model_manager import GeologicalModelManager +from loopstructural.toolbelt.preferences import PlgSettingsStructure + + +def _contact(unit_name): + return pd.DataFrame({'X': [0.0], 'Y': [0.0], 'Z': [0.0], 'source_unit': [unit_name]}) + + +def _fault_trace(): + return pd.DataFrame({'X': [0.0, 10.0], 'Y': [0.0, 10.0], 'Z': [0.0, 0.0]}) + + +def _two_group_column(): + """basin_floor/basin_fill, then a named unconformity, then cover_lower/cover_upper.""" + column = StratigraphicColumn() + column.clear(basement=False) + column.add_unit(name='basin_floor', thickness=50.0, where='top') + column.add_unit(name='basin_fill', thickness=150.0, where='top') + boundary = column.add_unconformity(name='regional_unconformity', where='top') + column.add_unit(name='cover_lower', thickness=80.0, where='top') + column.add_unit(name='cover_upper', thickness=120.0, where='top') + return column, boundary + + +@pytest.fixture +def manager(monkeypatch): + manager = GeologicalModelManager() + calls = {'foliation': [], 'unconformity': [], 'domain_fault': [], 'fault': []} + + def fake_create_and_add_foliation(name, data=None, **kwargs): + calls['foliation'].append((name, data)) + return object() + + def fake_add_unconformity(feature, value, **kwargs): + calls['unconformity'].append((feature, value)) + + def fake_create_and_add_domain_fault(fault_surface_data, **kwargs): + calls['domain_fault'].append(fault_surface_data) + return object() + + def fake_create_and_add_fault(fault_name, displacement, **kwargs): + calls['fault'].append(fault_name) + return object() + + monkeypatch.setattr(manager.model, 'create_and_add_foliation', fake_create_and_add_foliation) + monkeypatch.setattr(manager.model, 'add_unconformity', fake_add_unconformity) + monkeypatch.setattr( + manager.model, 'create_and_add_domain_fault', fake_create_and_add_domain_fault + ) + monkeypatch.setattr(manager.model, 'create_and_add_fault', fake_create_and_add_fault) + manager._calls = calls + return manager + + +class TestFaultDomainBoundary: + def test_group_boundary_linked_to_fault_builds_domain_fault(self, manager): + column, boundary = _two_group_column() + manager.stratigraphic_column = column + for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): + manager.stratigraphy[name]['contact'] = _contact(name) + manager.faults['boundary_fault']['data'] = _fault_trace() + manager.fault_boundaries[boundary.uuid] = 'boundary_fault' + + manager.update_foliation_features() + + assert manager._calls['domain_fault'] == ['boundary_fault'] + # Only the topmost group (nothing above it in the column) still + # falls back to a flat unconformity. + assert len(manager._calls['unconformity']) == 1 + + registered = manager.model.data + fault_rows = registered.loc[registered['feature_name'] == 'boundary_fault'] + # The original 2 trace points plus 2 synthetic points extending the + # trace to the model's bounding box edges (see + # `_extend_fault_trace_to_domain`). + assert len(fault_rows) == 4 + assert set(fault_rows['val']) == {0} + + def test_group_boundary_without_fault_link_uses_flat_unconformity(self, manager): + column, _boundary = _two_group_column() + manager.stratigraphic_column = column + for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): + manager.stratigraphy[name]['contact'] = _contact(name) + + manager.update_foliation_features() + + assert manager._calls['domain_fault'] == [] + assert len(manager._calls['unconformity']) == 2 + + def test_missing_fault_data_falls_back_to_flat_unconformity(self, manager): + column, boundary = _two_group_column() + manager.stratigraphic_column = column + for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): + manager.stratigraphy[name]['contact'] = _contact(name) + # Linked to a fault, but no trace data was ever ingested for it. + manager.fault_boundaries[boundary.uuid] = 'boundary_fault' + + manager.update_foliation_features() + + assert manager._calls['domain_fault'] == [] + assert len(manager._calls['unconformity']) == 2 + + def test_domain_boundary_fault_skipped_by_displacement_fault_build(self, manager): + _column, boundary = _two_group_column() + manager.faults['boundary_fault']['data'] = _fault_trace() + manager.fault_boundaries[boundary.uuid] = 'boundary_fault' + manager.faults['normal_fault']['data'] = _fault_trace() + + manager.update_fault_features() + + assert manager._calls['fault'] == ['normal_fault'] + + def test_abutting_relationships_skip_unbuilt_domain_boundary_fault(self, manager): + """`update_fault_features` calls `apply_fault_abutting_relationships` + before `update_foliation_features` has run, so a domain-boundary + fault's feature (built there, not here) doesn't exist in the model + yet. `GeologicalModel.get_feature_by_name` raises ValueError rather + than returning None for a name that isn't built -- this must not + propagate out of `apply_fault_abutting_relationships`. + """ + column, boundary = _two_group_column() + manager.fault_topology = FaultTopology(column) + manager.fault_topology.add_fault('boundary_fault') + manager.fault_topology.add_fault('normal_fault') + manager.faults['boundary_fault']['data'] = _fault_trace() + manager.faults['normal_fault']['data'] = _fault_trace() + manager.fault_boundaries[boundary.uuid] = 'boundary_fault' + + # Neither fault has been built into manager.model yet. + manager.apply_fault_abutting_relationships() # must not raise + + +class TestDomainFaultBuildsAndSolves: + """End-to-end regression test against the real LoopStructural build/solve + path (nothing mocked here, unlike the other test classes in this module). + + `create_and_add_domain_fault` reads its points straight from + `model.data` with no column normalisation, unlike + `create_and_add_foliation`/`create_and_add_fault` (which run their + `data=` argument through `GeologicalModel.prepare_data` internally) -- + so registering raw X/Y/Z/feature_name/val rows for a domain-boundary + fault builds fine but crashes later, during Solve Model, with + `KeyError: "None of [Index(['gx', 'gy', 'gz']...` deep inside + `add_data_to_interpolator`. `_build_domain_fault_boundary` must run its + data through `model.prepare_data` itself first. + + It also pins the fix for a second, more subtle bug: LoopStructural's + `add_unconformity` backward crop walk only recognises FAULT/ + INACTIVEFAULT feature types as "already handled, skip" -- not + DOMAINFAULT. `_two_group_column()`'s topmost group (cover) has no fault + link, so it falls back to a plain `add_unconformity` call, which walks + straight through the already-built `boundary_fault` domain fault and + incorrectly adds itself as a region on it. Since + `GeologicalFeature.evaluate_value` returns NaN wherever a feature's + regions don't hold, this reads as the domain fault's own scalar field + being NaN on whichever side that unrelated unconformity's condition + fails. `_strip_spurious_regions_from_domain_faults` removes it. + """ + + def test_domain_fault_builds_and_solves(self, monkeypatch): + monkeypatch.setattr(PlgSettingsStructure, 'interpolator_nelements', 200) + manager = GeologicalModelManager() + manager.update_bounding_box(BoundingBox(origin=[0, 0, -50], maximum=[100, 100, 50])) + column, boundary = _two_group_column() + manager.stratigraphic_column = column + for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): + manager.stratigraphy[name]['contact'] = _contact(name) + manager.faults['boundary_fault']['data'] = pd.DataFrame( + {'X': [10.0, 90.0], 'Y': [10.0, 90.0], 'Z': [0.0, 0.0]} + ) + manager.fault_boundaries[boundary.uuid] = 'boundary_fault' + + manager.update_model(notify_observers=False) + manager.update_all_features(notify_observers=False) # must not raise + + domain_fault = manager.model.get_feature_by_name('boundary_fault') + assert domain_fault is not None + assert domain_fault.regions == [] + + xs, ys = np.meshgrid(np.linspace(0, 100, 11), np.linspace(0, 100, 11)) + pts = np.column_stack([xs.ravel(), ys.ravel(), np.zeros(xs.size)]) + values = domain_fault.evaluate_value(pts) + assert not np.any(np.isnan(values)), ( + "domain fault scalar field is NaN somewhere in the model domain -- " + "an unconformity region was incorrectly left on it" + ) + + +class TestExtendFaultTraceToDomain: + """`_extend_fault_trace_to_domain` is what makes a domain-boundary fault + behave as an "infinite" cut: it adds two synthetic points extending the + trace's own trend out to the model's bounding box edges, so the + interpolated surface isn't only reliable near the digitised trace. + """ + + def test_extends_a_short_trace_to_the_domain_edges(self, manager): + manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100])) + # A short trace entirely inside the domain, running along y=x. + trace = pd.DataFrame({'X': [40.0, 60.0], 'Y': [40.0, 60.0], 'Z': [0.0, 0.0]}) + + extended = manager._extend_fault_trace_to_domain(trace) + + assert len(extended) == 4 + xy = extended[['X', 'Y']].to_numpy() + # Two of the four points must sit on the bounding box boundary. + on_boundary = [ + np.isclose(x, 0) or np.isclose(x, 100) or np.isclose(y, 0) or np.isclose(y, 100) + for x, y in xy + ] + assert sum(on_boundary) == 2 + + def test_extends_a_dipping_trace_keeping_its_trend(self, manager): + manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100])) + # Trace along y=x with Z increasing 1:1 with distance along the line. + trace = pd.DataFrame({'X': [40.0, 60.0], 'Y': [40.0, 60.0], 'Z': [-10.0, 10.0]}) + + extended = manager._extend_fault_trace_to_domain(trace) + + synthetic = extended.iloc[2:] + # The synthetic points should extrapolate the same Z-vs-along-line + # trend rather than reusing the nearest original Z value. + assert synthetic['Z'].max() > 10.0 + assert synthetic['Z'].min() < -10.0 + + def test_single_point_trace_is_left_unchanged(self, manager): + manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100])) + trace = pd.DataFrame({'X': [50.0], 'Y': [50.0], 'Z': [0.0]}) + + extended = manager._extend_fault_trace_to_domain(trace) + + assert len(extended) == 1 + + def test_direction_parallel_to_axis_outside_domain_is_left_unchanged(self, manager): + manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[10, 10, 10])) + # A vertical (constant-X) trace entirely outside the domain in X. + trace = pd.DataFrame({'X': [500.0, 500.0], 'Y': [1.0, 2.0], 'Z': [0.0, 0.0]}) + + extended = manager._extend_fault_trace_to_domain(trace) + + assert len(extended) == 2 From d4cbdc995eb40bea294f44a3f312c0ff99f52ab4 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Fri, 21 Aug 2026 14:23:05 +0930 Subject: [PATCH 3/7] fix: correct stratigraphic value assignment and domain-fault boundary follow-ups Several fixes to make domain-fault-bounded stratigraphic columns build and display correctly: - update_foliation_features now trains each unit's basal-contact data at its own max() (the boundary with the next-older unit, i.e. its true base) instead of min() (the boundary with the next-younger unit, i.e. its top). Basal contacts represent a unit's base, so training at min() anchored every unit's own data to the wrong boundary -- confirmed on a live project where units evaluated into their next-younger neighbour's bracket instead of their own, and a basement unit with no contact data of its own never appeared in the model at all. - Unit thickness now accumulates unconditionally while building that training data, so an undigitised placeholder unit no longer shifts every later unit's value by its own thickness. - Use each fault trace point's own local tangent (rather than one global best-fit line) when extending a domain-boundary fault to the model's bounding box and deriving its orientation constraints, so a curved trace doesn't get flattened into the wrong extrapolation. - Recompute stratigraphic unit value ranges after restoring a column from a saved project (both initial load and reload), matching what a fresh column already gets -- otherwise every restored unit kept the default (0, inf) range and couldn't be told apart from its neighbours. - Show the generic details panel for a domain-fault feature instead of an empty widget. - Skip an isosurface with no geometry when adding stratigraphic surfaces to the 3D viewer instead of crashing, since an undigitised unit can legitimately have no constrained geometry anywhere in the model. Co-Authored-By: Claude Sonnet 5 --- .../geological_model_tab.py | 12 + .../gui/visualisation/feature_list_widget.py | 12 +- loopstructural/main/data_manager.py | 12 + loopstructural/main/model_manager.py | 205 +++++++++++------- tests/qgis/test_fault_domain_boundary.py | 39 +++- .../test_stratigraphic_value_consistency.py | 109 ++++++++-- 6 files changed, 282 insertions(+), 107 deletions(-) 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 13ac391..7dd559b 100644 --- a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py +++ b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py @@ -20,6 +20,7 @@ from .add_foliation_dialog import AddFoliationDialog from .add_unconformity_dialog import AddUnconformityDialog from .feature_details_panel import ( + BaseFeatureDetailsPanel, FaultFeatureDetailsPanel, FoldedFeatureDetailsPanel, FoliationFeatureDetailsPanel, @@ -461,6 +462,17 @@ def on_feature_selected(self, item): self.featureDetailsPanel = FoldedFeatureDetailsPanel( feature=feature, model_manager=self.model_manager, data_manager=self.data_manager ) + elif feature.type == FeatureType.DOMAINFAULT: + # A domain fault is built by the same GeologicalFeatureBuilder + # as a foliation (see create_and_add_domain_fault), just with a + # different .type tag -- the generic base panel (interpolator + # settings, data layers, export/evaluate) already applies to it + # unchanged. Skip FoliationFeatureDetailsPanel's fold-frame + # attachment controls, which don't make sense for a domain + # boundary. + self.featureDetailsPanel = BaseFeatureDetailsPanel( + feature=feature, model_manager=self.model_manager, data_manager=self.data_manager + ) else: self.featureDetailsPanel = QWidget() # Default empty panel diff --git a/loopstructural/gui/visualisation/feature_list_widget.py b/loopstructural/gui/visualisation/feature_list_widget.py index e6e6539..2d699d2 100644 --- a/loopstructural/gui/visualisation/feature_list_widget.py +++ b/loopstructural/gui/visualisation/feature_list_widget.py @@ -458,8 +458,18 @@ def add_stratigraphic_surfaces(self): stratigraphic_surfaces = self.model_manager.model.get_stratigraphic_surfaces() for surface in stratigraphic_surfaces: + mesh = surface.vtk() + if mesh.n_points == 0: + # A unit with no digitised data of its own (e.g. an + # undigitised placeholder like "Top") can have no + # constrained geometry anywhere in the model, so its + # isovalue may not intersect the solved field at all -- + # pyvista refuses to plot an empty mesh, so skip it rather + # than crashing every surface after it in this loop. + logger.info(f"Skipping '{surface.name}': isosurface has no geometry.") + continue self.viewer.add_mesh_object( - surface.vtk(), + mesh, name=surface.name, color=surface.colour, source_feature=surface.name, diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py index d47703d..6ef0f6a 100644 --- a/loopstructural/main/data_manager.py +++ b/loopstructural/main/data_manager.py @@ -1099,6 +1099,8 @@ def from_dict(self, data): self._structural_orientations = data['structural_orientations'] if 'stratigraphic_column' in data: self._stratigraphic_column = StratigraphicColumn.from_dict(data['stratigraphic_column']) + # See the matching call in update_from_dict for why this is needed. + self._stratigraphic_column.update_unit_values() self.stratigraphic_column_callback() self._fault_boundaries.clear() if data.get('fault_boundaries'): @@ -1199,6 +1201,16 @@ def update_from_dict(self, data): ) if 'stratigraphic_column' in data: self._stratigraphic_column.update_from_dict(data['stratigraphic_column']) + # update_from_dict restores elements via add_element, not + # add_unit -- only add_unit computes each unit's min/max + # scalar-field range as a side effect. Without this, every + # restored unit keeps the default (0, inf) range, so + # evaluate_model can't tell any unit in a group apart from any + # other and just labels every point with whichever unit was + # last in the group (see GeologicalModelManager. + # set_stratigraphic_column, which already does this for the + # very first load -- this covers every reload afterwards). + self._stratigraphic_column.update_unit_values() else: self._stratigraphic_column.clear() diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index afa3798..595cb9a 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -592,91 +592,115 @@ def _clip_line_to_bounding_box(self, centroid, direction): def _extend_fault_trace_to_domain(self, fault_data): """Add two synthetic points that extend a fault's trace out to the - edges of the model's bounding box along its own overall trend. + edges of the model's bounding box, and attach a `strike` column + derived from the trace's own *local* tangent at each point. `create_and_add_domain_fault` interpolates a scalar field only from the points it is given, over the model's exact bounding box (no buffer, unlike a displacement fault's mesh) -- so a locally digitised trace only reliably constrains the surface near itself, - and the domain crop can wander unpredictably further away. Fitting - a line through the existing XY points and adding two constraint - points where that line meets the bounding box edges keeps the - interpolated surface following the fault's actual trend all the - way across the domain, rather than an arbitrary extrapolation -- + and the domain crop can wander unpredictably further away. + Extending each end along its own local tangent, out to where it + meets the bounding box edge, keeps the interpolated surface + following the trace's actual trend all the way across the domain -- this is what makes the fault behave as an "infinite" domain boundary rather than a locally-anchored patch. + Using each point's *local* tangent (rather than one global + best-fit line through the whole trace) matters for a genuinely + curved trace: fitting a single global line flattens that curvature + out, and extrapolating along it can land an extension point (or + bias the interpolated field generally) on the wrong side of the + real curve relative to data that's actually near the trace. + Confirmed on a live project: a global-line fit classified a + stratigraphic unit's own contact data as being on the opposite + side of the domain fault from a bounding-box corner that a proper + local (nearest-segment) classification put on the *same* side as + that data -- i.e. the global fit was extrapolating the wrong way. + Z at each synthetic point is extrapolated linearly against - distance along that line, so a dipping trace keeps its dip. + distance along the local end segment, so a dipping trace keeps its + dip at the point it's extended from. """ xy = fault_data[['X', 'Y']].to_numpy() - if len(xy) < 2: - return fault_data - centroid = xy.mean(axis=0) - # Principal direction of the trace via SVD -- robust to a - # near-vertical (large-Y-range, small-X-range) trace, unlike a - # simple polyfit of Y against X. - _, _, vt = np.linalg.svd(xy - centroid) - direction = vt[0] - clipped = self._clip_line_to_bounding_box(centroid, direction) - if clipped is None: - return fault_data - t_min, t_max = clipped - projections = (xy - centroid) @ direction z = fault_data['Z'].to_numpy() - if len(np.unique(projections)) > 1: - z_slope, z_intercept = np.polyfit(projections, z, 1) - else: - z_slope, z_intercept = 0.0, float(np.mean(z)) - new_rows = [] - for t in (t_min, t_max): - point_xy = centroid + direction * t - new_rows.append({'X': point_xy[0], 'Y': point_xy[1], 'Z': z_slope * t + z_intercept}) - return pd.concat([fault_data, pd.DataFrame(new_rows)], ignore_index=True) - - def _domain_fault_orientation_rows(self, points_xyz, fault_entry): - """Build strike/dip orientation constraint rows for a domain-boundary fault. - - A domain fault built only from same-valued (val=0) points has no - information about which direction the field should vary -- the - minimum-curvature solution that exactly satisfies "value=0 along - this line" with nothing else to go on is a trivial constant/flat - field (zero curvature everywhere, zero misfit). That's - geometrically useless: the crop condition - `domain_fault.evaluate_value(pos) > 0` - (see `LoopStructural.modelling.core._model_relationships`) is then - never true anywhere, so everything on the "positive" side reads as - NaN. Confirmed against a live project: a fault built from 31 trace - points and 2 extension points, all val=0, interpolated to an - exactly flat 0.0 field everywhere; adding one orientation - constraint per point produced a properly varying field crossing - zero along the fault's own trend. - - Strike comes from the trace's own principal direction (matching - `_extend_fault_trace_to_domain`'s line fit); dip comes from the - ingested fault trace data's `dip` column if present (matching how - `update_fault_features` picks up dip for a displacement fault), - otherwise defaults to vertical (90 degrees). - """ - xy = points_xyz[['X', 'Y']].to_numpy() - if len(xy) < 2: - return None - centroid = xy.mean(axis=0) - _, _, vt = np.linalg.svd(xy - centroid) - direction = vt[0] + n = len(xy) + result = fault_data.copy() + if n < 2: + result['strike'] = np.nan + return result + + # Local tangent per point: central difference for interior points, + # forward/backward difference at the ends. Assumes points follow + # the digitised line's vertex order (true for AllSampler-derived + # trace data, which walks each LineString's coords in order). + tangents = np.zeros((n, 2)) + tangents[0] = xy[1] - xy[0] + tangents[-1] = xy[-1] - xy[-2] + if n > 2: + tangents[1:-1] = xy[2:] - xy[:-2] # strikedip2vector's strike is degrees clockwise from North (+Y); # atan2(dx, dy) matches that convention directly. - strike = float(np.degrees(np.arctan2(direction[0], direction[1])) % 360) - dip = 90.0 + result['strike'] = np.degrees(np.arctan2(tangents[:, 0], tangents[:, 1])) % 360 + + new_rows = [] + # Extend backward past the first point, continuing on its own + # local tangent (pointing away from the second point). + start_seg = xy[1] - xy[0] + start_len = np.linalg.norm(start_seg) + if start_len > 1e-9: + direction = -start_seg / start_len + clipped = self._clip_line_to_bounding_box(xy[0], direction) + if clipped is not None: + _, t_max = clipped + if t_max > 0: + point_xy = xy[0] + direction * t_max + z_slope = (z[1] - z[0]) / start_len + new_rows.append( + { + 'X': point_xy[0], + 'Y': point_xy[1], + 'Z': z[0] - z_slope * t_max, + 'strike': result['strike'].iloc[0], + } + ) + # Extend forward past the last point, continuing on its own local + # tangent (pointing away from the second-to-last point). + end_seg = xy[-1] - xy[-2] + end_len = np.linalg.norm(end_seg) + if end_len > 1e-9: + direction = end_seg / end_len + clipped = self._clip_line_to_bounding_box(xy[-1], direction) + if clipped is not None: + _, t_max = clipped + if t_max > 0: + point_xy = xy[-1] + direction * t_max + z_slope = (z[-1] - z[-2]) / end_len + new_rows.append( + { + 'X': point_xy[0], + 'Y': point_xy[1], + 'Z': z[-1] + z_slope * t_max, + 'strike': result['strike'].iloc[-1], + } + ) + if not new_rows: + return result + return pd.concat([result, pd.DataFrame(new_rows)], ignore_index=True) + + def _domain_fault_dip(self, fault_entry): + """Dip (degrees from horizontal) for a domain-boundary fault. + + Uses the ingested fault trace data's `dip` column if present + (matching how `update_fault_features` picks up dip for a + displacement fault), otherwise defaults to vertical (90 degrees). + """ raw_data = fault_entry.get('data') if fault_entry else None if raw_data is not None and 'dip' in raw_data: dip_values = raw_data['dip'].dropna() if not dip_values.empty: - dip = float(dip_values.mean()) - rows = points_xyz[['X', 'Y', 'Z']].copy() - rows['strike'] = strike - rows['dip'] = dip - return rows + return float(dip_values.mean()) + return 90.0 def _build_domain_fault_boundary(self, fault_name, groupname): """Build `fault_name` as a domain-fault boundary in place of a flat unconformity. @@ -702,14 +726,19 @@ def _build_domain_fault_boundary(self, fault_name, groupname): log_level=2, ) return False - data_for_fault = self._extend_fault_trace_to_domain(fault_data[['X', 'Y', 'Z']].copy()) - orientation_rows = self._domain_fault_orientation_rows(data_for_fault, fault_entry) - data_for_fault['feature_name'] = fault_name - data_for_fault['val'] = 0 - if orientation_rows is not None: + extended = self._extend_fault_trace_to_domain(fault_data[['X', 'Y', 'Z']].copy()) + + value_rows = extended[['X', 'Y', 'Z']].copy() + value_rows['feature_name'] = fault_name + value_rows['val'] = 0 + + orientation_rows = extended.dropna(subset=['strike'])[['X', 'Y', 'Z', 'strike']].copy() + data_for_fault = value_rows + if not orientation_rows.empty: + orientation_rows['dip'] = self._domain_fault_dip(fault_entry) orientation_rows['feature_name'] = fault_name orientation_rows['val'] = np.nan - data_for_fault = pd.concat([data_for_fault, orientation_rows], ignore_index=True) + data_for_fault = pd.concat([value_rows, orientation_rows], ignore_index=True) # Unlike create_and_add_foliation/create_and_add_fault (which # normalise their own `data=` argument internally via # model.prepare_data before building), create_and_add_domain_fault @@ -758,10 +787,32 @@ def update_foliation_features(self): groupname = group.name stratigraphic_column[groupname] = {} for u in reversed(group.units): + # `reversed(group.units)` walks youngest-to-oldest (matching + # StratigraphicColumn.update_unit_values's own cumulative + # walk), so `val` must accumulate every unit's thickness + # *before* being used as that unit's own training value -- + # regardless of whether the unit has any digitised data -- + # to land on `u.max()`, not `u.min()`. + # + # `u.min()` is the boundary shared with the next *younger* + # neighbour (this unit's top); `u.max()` is the boundary + # shared with the next *older* neighbour (this unit's true + # base). Digitised "basal contact" data represents a unit's + # base, so it belongs at `u.max()`. Using `u.min()` instead + # anchors every unit's own contact points to its top + # boundary rather than its base -- confirmed on a live + # project: every unit's own mapped points evaluated into its + # next-younger neighbour's bracket instead of its own. + # + # Accumulating unconditionally (not skipped for a unit with + # no digitised data, e.g. an undigitised "Top"/basement + # placeholder) also keeps every later unit's value aligned + # with `get_isovalues()`'s own cumulative-thickness bracket + # boundaries, which don't know or care which units were + # actually mapped. + val += u.thickness unit_data = self.stratigraphy.get(u.name, None) - if unit_data is None: - continue - else: + if unit_data is not None: if 'contact' in unit_data: contact = unit_data['contact'] if not contact.empty: @@ -774,8 +825,6 @@ def update_foliation_features(self): orientations['val'] = np.nan orientations['feature_name'] = groupname data.append(orientations) - - val += u.thickness if len(data) == 0: self._debug_manager.log( f"No data found for group {groupname}, skipping.", log_level=2 diff --git a/tests/qgis/test_fault_domain_boundary.py b/tests/qgis/test_fault_domain_boundary.py index df608b4..c4cfb35 100644 --- a/tests/qgis/test_fault_domain_boundary.py +++ b/tests/qgis/test_fault_domain_boundary.py @@ -90,11 +90,16 @@ def test_group_boundary_linked_to_fault_builds_domain_fault(self, manager): registered = manager.model.data fault_rows = registered.loc[registered['feature_name'] == 'boundary_fault'] - # The original 2 trace points plus 2 synthetic points extending the - # trace to the model's bounding box edges (see - # `_extend_fault_trace_to_domain`). + # The default bounding box here (never set explicitly) is smaller + # than the trace itself, so no synthetic edge-extension points get + # added (see TestExtendFaultTraceToDomain for that, with a + # realistic bounding box). What's registered is the 2 trace points + # as value (val=0) constraints, plus the same 2 points again as + # orientation (val=NaN, strike/dip) constraints -- a domain fault + # needs both, see _domain_fault_dip / _extend_fault_trace_to_domain. assert len(fault_rows) == 4 - assert set(fault_rows['val']) == {0} + assert len(fault_rows.loc[fault_rows['val'] == 0]) == 2 + assert fault_rows['val'].isna().sum() == 2 def test_group_boundary_without_fault_link_uses_flat_unconformity(self, manager): column, _boundary = _two_group_column() @@ -258,3 +263,29 @@ def test_direction_parallel_to_axis_outside_domain_is_left_unchanged(self, manag extended = manager._extend_fault_trace_to_domain(trace) assert len(extended) == 2 + + def test_local_tangent_varies_along_a_curved_trace(self, manager): + """Regression test for a real bug: fitting one global best-fit line + through a curved trace (the old approach) flattens its curvature + out, and can extrapolate/orient the interpolated surface on the + wrong side of real nearby data. Confirmed on a live project where + a global-line fit put a stratigraphic unit's own contact data on + the opposite side of its domain-boundary fault from a bounding-box + corner that a correct local (nearest-segment) classification put + on the *same* side. Each point's strike must instead follow its + own local tangent. + """ + manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100])) + # An L-shaped trace: a horizontal leg then a vertical leg. + trace = pd.DataFrame( + {'X': [20.0, 50.0, 50.0], 'Y': [50.0, 50.0, 80.0], 'Z': [0.0, 0.0, 0.0]} + ) + + extended = manager._extend_fault_trace_to_domain(trace) + + strikes = extended['strike'].to_numpy() + # Row 0's local tangent is horizontal (towards row 1); row 2's is + # vertical (away from row 1) -- these must differ substantially. A + # single global best-fit line would instead give every point close + # to the same strike. + assert abs(((strikes[0] - strikes[2] + 180) % 360) - 180) > 45 diff --git a/tests/qgis/test_stratigraphic_value_consistency.py b/tests/qgis/test_stratigraphic_value_consistency.py index 5135a20..8331f80 100644 --- a/tests/qgis/test_stratigraphic_value_consistency.py +++ b/tests/qgis/test_stratigraphic_value_consistency.py @@ -1,19 +1,34 @@ -"""Regression test for the training-value / isovalue direction bug. +"""Regression test for the basal-contact training-value direction. `GeologicalModelManager.update_foliation_features` assigns a scalar `val` to each unit's basal contact before handing the data to the interpolator. -`StratigraphicColumn.get_isovalues` (LoopStructural core) later decides which -name to stamp on each extracted isosurface, using its own idea of which -value belongs to which unit. - -These two must agree on direction (does value increase from oldest-to- -youngest, or youngest-to-oldest?), or every extracted surface gets labelled -with the wrong unit while keeping correct geometry -- see the "stratigraphic -column was reversed" fixes in model_manager.py (2025-07-21) and the widget -(2025-08-21, reverted 2025-09-08). This has flipped back and forth as this -plugin and LoopStructural evolved independently; this test pins the -invariant so a future change on either side fails loudly here instead of -silently inverting a user's model. +`StratigraphicColumn.update_unit_values` (LoopStructural core) computes each +unit's `min()`/`max()` by walking the column youngest-to-oldest, accumulating +thickness from 0 -- forced by the fact that a basement unit's open-ended +range (`thickness=inf`) only works when it's the *last* unit processed in +that walk (an infinite thickness earlier would poison every unit after it). +That makes `u.min()` the boundary shared with the next *younger* neighbour +(a unit's top) and `u.max()` the boundary shared with the next *older* +neighbour (a unit's true base). + +A digitised "basal contact" represents a unit's base, so it must be trained +at `val = u.max()`, not `u.min()`. Training at `u.min()` (the old behaviour) +anchors every unit's own contact points to its top boundary instead of its +base -- confirmed on a live project: every unit's own mapped points +evaluated into its next-younger neighbour's bracket instead of its own, +e.g. "Formacao Betari"'s own contact data landing inside "Formacao +Guaricanga"'s value range. + +Note `get_isovalues()` also reports `u.min()` per unit -- that's a separate +concern (naming which unit an *extracted isosurface* belongs to), not a +statement about which value basal-contact training data should use, so this +test does not compare against it. + +This direction has flipped back and forth as this plugin and LoopStructural +evolved independently -- see the "stratigraphic column was reversed" fixes +in model_manager.py (2025-07-21) and the widget (2025-08-21, reverted +2025-09-08). This test pins the invariant so a future change fails loudly +here instead of silently inverting a user's model. """ import pandas as pd @@ -46,7 +61,7 @@ def fake_create_and_add_foliation(name, data=None, **kwargs): return manager -class TestTrainingValueMatchesIsovalue: +class TestTrainingValueIsUnitsOwnBase: def test_single_group_three_units(self, manager): column = StratigraphicColumn() column.clear(basement=False) # single flat group, no unconformities @@ -61,16 +76,14 @@ def test_single_group_three_units(self, manager): manager.update_foliation_features() training_values = self._training_values_by_unit(manager._captured_calls) - expected_values = { - name: entry['value'] for name, entry in column.get_isovalues().items() - } + expected_values = self._own_base_by_unit(column, ('oldest', 'middle', 'youngest')) for unit_name in ('oldest', 'middle', 'youngest'): assert training_values[unit_name] == pytest.approx(expected_values[unit_name]), ( - f"'{unit_name}' was trained with val={training_values[unit_name]} but " - f"get_isovalues() will label the value={expected_values[unit_name]} surface " - f"with this unit's name -- the trained field and the isosurface labels " - f"disagree on direction, so extracted surfaces will get the wrong unit name." + f"'{unit_name}' was trained with val={training_values[unit_name]} but its " + f"own base (boundary with the next-older neighbour) is " + f"{expected_values[unit_name]} -- basal-contact data must train at a unit's " + f"own max(), not min(), or extracted surfaces get the wrong unit name." ) def test_two_groups_split_by_unconformity(self, manager): @@ -89,13 +102,61 @@ def test_two_groups_split_by_unconformity(self, manager): manager.update_foliation_features() training_values = self._training_values_by_unit(manager._captured_calls) - expected_values = { - name: entry['value'] for name, entry in column.get_isovalues().items() - } + expected_values = self._own_base_by_unit( + column, ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper') + ) for unit_name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'): assert training_values[unit_name] == pytest.approx(expected_values[unit_name]) + def test_undigitised_unit_does_not_shift_later_units_in_group(self, manager): + """Regression test for a real bug: a unit with no digitised contact + or orientation data (e.g. a "Top" unit nobody has mapped points for) + must still contribute its own thickness to `val` for every unit + that follows it in the group -- `update_foliation_features` used to + `continue` past an undigitised unit before accumulating its + thickness, which shifted every later unit's trained value by that + unit's thickness relative to its own true base. + """ + column = StratigraphicColumn() + column.clear(basement=False) + # Added first so it ends up last in `group.units` (each `add_unit` + # prepends) and therefore *first* in the `reversed(group.units)` + # build loop -- matching the live project, where the undigitised + # unit was the one whose skipped increment shifted every unit + # after it. + column.add_unit(name='Top', thickness=999.0, where='top') + column.add_unit(name='basin_fill', thickness=150.0, where='top') + column.add_unit(name='basin_floor', thickness=50.0, where='top') + + manager.stratigraphic_column = column + for name in ('basin_floor', 'basin_fill'): + manager.stratigraphy[name]['contact'] = _contact(name) + # 'Top' deliberately has no entry in manager.stratigraphy at all. + + manager.update_foliation_features() + + training_values = self._training_values_by_unit(manager._captured_calls) + expected_values = self._own_base_by_unit(column, ('basin_floor', 'basin_fill')) + + for unit_name in ('basin_floor', 'basin_fill'): + assert training_values[unit_name] == pytest.approx(expected_values[unit_name]), ( + f"'{unit_name}' was trained with val={training_values[unit_name]} but its " + f"own base is {expected_values[unit_name]} -- an undigitised unit earlier in " + f"the group must still shift later units' trained values by its own thickness." + ) + + @staticmethod + def _own_base_by_unit(column, unit_names): + """Each unit's own base: the boundary with the next-*older* neighbour, + i.e. `u.max()` -- see module docstring for why max() (not min(), which + `get_isovalues()` reports) is the correct target for basal-contact + training data.""" + units_by_name = { + u.name: u for group in column.get_groups() for u in group.units + } + return {name: units_by_name[name].max() for name in unit_names} + @staticmethod def _training_values_by_unit(captured_calls): combined = pd.concat(captured_calls, ignore_index=True) From e09ed0b7ac80db414f58ddaf659411ceb3ab7a77 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Fri, 21 Aug 2026 14:33:41 +0930 Subject: [PATCH 4/7] fix: increase thickness spin box limit + nelements spin box. BUmp default nelements to 50k --- loopstructural/gui/dlg_settings.ui | 6 +++++- .../modelling/stratigraphic_column/stratigraphic_unit.ui | 2 +- loopstructural/toolbelt/preferences.py | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/loopstructural/gui/dlg_settings.ui b/loopstructural/gui/dlg_settings.ui index 998b597..477f844 100644 --- a/loopstructural/gui/dlg_settings.ui +++ b/loopstructural/gui/dlg_settings.ui @@ -274,7 +274,11 @@ - + + + 1000000 + + diff --git a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui index 3f9a3b9..306f738 100644 --- a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui +++ b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui @@ -80,7 +80,7 @@ 0.000000000000000 - 10000.000000000000000 + 100000000.000000000000000 diff --git a/loopstructural/toolbelt/preferences.py b/loopstructural/toolbelt/preferences.py index 99cfb6e..083ce8b 100644 --- a/loopstructural/toolbelt/preferences.py +++ b/loopstructural/toolbelt/preferences.py @@ -31,7 +31,7 @@ class PlgSettingsStructure: debug_directory: str = "" version: str = __version__ interpolator_type: str = 'FDI' - interpolator_nelements: int = 10000 + interpolator_nelements: int = 50000 interpolator_regularisation: float = 1.0 interpolator_cpw: float = 1.0 interpolator_npw: float = 1.0 From c9d3c608a207cfe4c183b8921d5975bff81fda32 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 24 Aug 2026 12:49:25 +0930 Subject: [PATCH 5/7] fix: value should be the basal value of a unit --- loopstructural/main/model_manager.py | 45 ++++++++++++++-------------- 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 595cb9a..4c6bbb3 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -787,30 +787,27 @@ def update_foliation_features(self): groupname = group.name stratigraphic_column[groupname] = {} for u in reversed(group.units): - # `reversed(group.units)` walks youngest-to-oldest (matching - # StratigraphicColumn.update_unit_values's own cumulative - # walk), so `val` must accumulate every unit's thickness - # *before* being used as that unit's own training value -- - # regardless of whether the unit has any digitised data -- - # to land on `u.max()`, not `u.min()`. + # A unit's own `val` is `u.min()` -- the cumulative + # thickness *before* this unit's own thickness is added -- + # matching `StratigraphicColumn.update_unit_values` (a unit + # is only added to the column oldest-first via `where= + # 'top'`, so `min()` is the boundary shared with the + # next-*older* neighbour processed just before it, i.e. + # this unit's own base) and `get_isovalues()` (LoopStructural + # core; both walk `reversed(group.units)` accumulating the + # same way, so a unit's own training value and the isovalue + # `get_isovalues()` later labels with this unit's name + # agree -- see test_stratigraphic_value_consistency.py). # - # `u.min()` is the boundary shared with the next *younger* - # neighbour (this unit's top); `u.max()` is the boundary - # shared with the next *older* neighbour (this unit's true - # base). Digitised "basal contact" data represents a unit's - # base, so it belongs at `u.max()`. Using `u.min()` instead - # anchors every unit's own contact points to its top - # boundary rather than its base -- confirmed on a live - # project: every unit's own mapped points evaluated into its - # next-younger neighbour's bracket instead of its own. - # - # Accumulating unconditionally (not skipped for a unit with - # no digitised data, e.g. an undigitised "Top"/basement - # placeholder) also keeps every later unit's value aligned - # with `get_isovalues()`'s own cumulative-thickness bracket - # boundaries, which don't know or care which units were - # actually mapped. - val += u.thickness + # `val` must accumulate every unit's thickness regardless of + # whether that unit has any digitised data -- get_isovalues() + # assigns each unit's isovalue purely from cumulative + # thickness, with no knowledge of which units were actually + # mapped. Skipping the increment for an unmapped unit (e.g. + # a "Top" placeholder with no contact points) would shift + # every val assigned to units after it in this loop, so + # extracted isosurfaces would get labelled with the wrong + # unit name even though the geometry itself is fine. unit_data = self.stratigraphy.get(u.name, None) if unit_data is not None: if 'contact' in unit_data: @@ -825,6 +822,8 @@ def update_foliation_features(self): orientations['val'] = np.nan orientations['feature_name'] = groupname data.append(orientations) + + val += u.thickness if len(data) == 0: self._debug_manager.log( f"No data found for group {groupname}, skipping.", log_level=2 From f617814b9195b1d71b33971b2410ddb893898178 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 24 Aug 2026 13:19:02 +0930 Subject: [PATCH 6/7] style: fix black formatting and ruff lint (unnecessary generator) Co-Authored-By: Claude Sonnet 5 --- .../gui/map2loop_tools/fault_topology_widget.py | 6 ++---- .../gui/modelling/model_definition/bounding_box.py | 8 ++++---- .../gui/visualisation/feature_list_widget.py | 11 +++++++---- loopstructural/main/m2l_api.py | 4 +++- tests/qgis/test_stratigraphic_value_consistency.py | 4 +--- 5 files changed, 17 insertions(+), 16 deletions(-) diff --git a/loopstructural/gui/map2loop_tools/fault_topology_widget.py b/loopstructural/gui/map2loop_tools/fault_topology_widget.py index 5d786ac..76fa7ff 100644 --- a/loopstructural/gui/map2loop_tools/fault_topology_widget.py +++ b/loopstructural/gui/map2loop_tools/fault_topology_widget.py @@ -174,7 +174,7 @@ def _run_topology(self): # not just the ones map2loop found a relationship for. A fault with # no detected topological relationship is still a real fault and # must not be dropped from the fault topology. - new_faults = set(str(v) for v in gdf['ID'].unique()) + new_faults = {str(v) for v in gdf['ID'].unique()} # Add new faults; never remove existing ones here, so faults # without a detected relationship (or ones the user added @@ -209,9 +209,7 @@ def _run_topology(self): else: f1 = str(row.iloc[0]) f2 = str(row.iloc[1]) - ft.update_fault_relationship( - f1, f2, FaultRelationshipType.ABUTTING - ) + ft.update_fault_relationship(f1, f2, FaultRelationshipType.ABUTTING) except Exception: pass diff --git a/loopstructural/gui/modelling/model_definition/bounding_box.py b/loopstructural/gui/modelling/model_definition/bounding_box.py index d5062f4..4c594f6 100644 --- a/loopstructural/gui/modelling/model_definition/bounding_box.py +++ b/loopstructural/gui/modelling/model_definition/bounding_box.py @@ -32,11 +32,11 @@ def __init__(self, parent=None, data_manager=None): self.selectFromCurrentLayerButton, "mActionZoomToLayer.svg", "Select from Current Layer" ) self._style_tool_button( - self.useCurrentViewExtentButton, "mActionSetToCanvasExtent.svg", "Use Current View Extent" - ) - self._style_tool_button( - self.drawOnMapButton, "mActionAddBasicRectangle.svg", "Draw on Map" + self.useCurrentViewExtentButton, + "mActionSetToCanvasExtent.svg", + "Use Current View Extent", ) + self._style_tool_button(self.drawOnMapButton, "mActionAddBasicRectangle.svg", "Draw on Map") self.drawOnMapButton.setCheckable(True) self.drawOnMapButton.clicked.connect(self.drawOnMap) self._draw_extent_tool = None diff --git a/loopstructural/gui/visualisation/feature_list_widget.py b/loopstructural/gui/visualisation/feature_list_widget.py index 2d699d2..4c22961 100644 --- a/loopstructural/gui/visualisation/feature_list_widget.py +++ b/loopstructural/gui/visualisation/feature_list_widget.py @@ -56,9 +56,7 @@ def __init__(self, parent=None, *, model_manager=None, viewer=None, data_manager self.data_manager = data_manager # Add buttons - self.addBoundingBoxButton = self._make_tool_button( - "extents.svg", "Add Model Bounding Box" - ) + self.addBoundingBoxButton = self._make_tool_button("extents.svg", "Add Model Bounding Box") self.addFaultSurfacesButton = self._make_custom_icon_tool_button( "fault.svg", "Add Fault Surfaces" ) @@ -734,7 +732,12 @@ def _extract_line_xy(self, layer) -> Optional[np.ndarray]: except Exception: target_crs = None source_crs = layer.sourceCrs() - if target_crs is not None and target_crs.isValid() and source_crs.isValid() and source_crs != target_crs: + if ( + target_crs is not None + and target_crs.isValid() + and source_crs.isValid() + and source_crs != target_crs + ): geom = QgsGeometry(geom) geom.transform(QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance())) diff --git a/loopstructural/main/m2l_api.py b/loopstructural/main/m2l_api.py index 2b5f1d8..82aa0a8 100644 --- a/loopstructural/main/m2l_api.py +++ b/loopstructural/main/m2l_api.py @@ -98,7 +98,9 @@ def extract_basal_contacts( unit_name_col = 'UNITNAME' if 'UNITNAME' in geology.columns else unit_name_field if unit_name_col and unit_name_col in geology.columns: geology_unit_names = {str(v).strip() for v in geology[unit_name_col].dropna().unique()} - stratigraphic_names = {str(name).strip() for name in stratigraphic_order if name is not None} + stratigraphic_names = { + str(name).strip() for name in stratigraphic_order if name is not None + } ignored_names = {str(unit).strip() for unit in ignore_units if unit is not None} missing_from_column = sorted(geology_unit_names - stratigraphic_names - ignored_names) if missing_from_column: diff --git a/tests/qgis/test_stratigraphic_value_consistency.py b/tests/qgis/test_stratigraphic_value_consistency.py index 8331f80..ea94dde 100644 --- a/tests/qgis/test_stratigraphic_value_consistency.py +++ b/tests/qgis/test_stratigraphic_value_consistency.py @@ -152,9 +152,7 @@ def _own_base_by_unit(column, unit_names): i.e. `u.max()` -- see module docstring for why max() (not min(), which `get_isovalues()` reports) is the correct target for basal-contact training data.""" - units_by_name = { - u.name: u for group in column.get_groups() for u in group.units - } + units_by_name = {u.name: u for group in column.get_groups() for u in group.units} return {name: units_by_name[name].max() for name in unit_names} @staticmethod From 538a12a7362836144a0d195aa3b40c1af3d28641 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 24 Aug 2026 15:35:48 +0930 Subject: [PATCH 7/7] fix: training values no longer disagree with get_isovalues() direction update_foliation_features re-reversed group.units, which get_groups() already returns in the order get_isovalues() walks, so every basal contact was trained with the wrong scalar value. Co-Authored-By: Claude Sonnet 5 --- loopstructural/main/model_manager.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py index 4c6bbb3..cbd6074 100644 --- a/loopstructural/main/model_manager.py +++ b/loopstructural/main/model_manager.py @@ -786,7 +786,7 @@ def update_foliation_features(self): data = [] groupname = group.name stratigraphic_column[groupname] = {} - for u in reversed(group.units): + for u in group.units: # A unit's own `val` is `u.min()` -- the cumulative # thickness *before* this unit's own thickness is added -- # matching `StratigraphicColumn.update_unit_values` (a unit @@ -794,10 +794,12 @@ def update_foliation_features(self): # 'top'`, so `min()` is the boundary shared with the # next-*older* neighbour processed just before it, i.e. # this unit's own base) and `get_isovalues()` (LoopStructural - # core; both walk `reversed(group.units)` accumulating the - # same way, so a unit's own training value and the isovalue - # `get_isovalues()` later labels with this unit's name - # agree -- see test_stratigraphic_value_consistency.py). + # core). `group.units` (from `get_groups()`) is already in + # that oldest-after-youngest walk order, and `get_isovalues()` + # accumulates over it directly with no extra reversal, so + # this loop must not reverse it either -- doing so trains + # each unit with the wrong scalar value, see + # test_stratigraphic_value_consistency.py. # # `val` must accumulate every unit's thickness regardless of # whether that unit has any digitised data -- get_isovalues()