diff --git a/docs/development/usability-plan.md b/docs/development/usability-plan.md index d0f3278..c199fb1 100644 --- a/docs/development/usability-plan.md +++ b/docs/development/usability-plan.md @@ -334,8 +334,218 @@ model, the model uses basal contacts and thicknesses that match the new order. ### Phase 5: View and export -- [ ] Show "Open 3D view" as the next action after a successful build. -- [ ] Add export of surfaces, block model and cross-sections to step 5. +- [x] Show "Open 3D view" as the next action after a successful build. +- [x] Add export of surfaces, block model and cross-sections to step 5. + +Files: `gui/modelling/steps/export_panel.py`, `main/model_export.py`, +`gui/modelling/modelling_widget.py`. + +Acceptance: after a successful build, the dock tells the user to open the 3D +view. From step 5, a user can write the surfaces, the block model and a +cross-section to files without the 3D view. + +### Phase 6: Follow-up changes + +Four changes that come from use of phases 3 to 5. Do them in this order. Each +one is a separate pull request. + +#### 6.1 Two separate workflows + +Problem: with "Interpolate surfaces from constraints", the build still uses the +data of the map workflow. The build reads the basal contacts, the structural +orientations, the fault traces and the stratigraphic column that the user set +in the other steps. The "constraints" choice only hides steps 2 and 3. It does +not stop the model from using their data. + +Rules: + +- With "Interpolate surfaces from constraints", the model has only the features + that the user adds in step 4 (foliations, unconformities, parametric faults). + It has no feature from the column, and no fault from a trace layer. +- With "Build from a geological map", the model has the generated features, and + the features that the user adds. +- The data of one workflow is kept when the user changes the choice. It is not + used, and it is not deleted. + +Tasks: + +- [ ] Give the model manager the workflow mode. `update_model` skips + `update_fault_features` (trace faults only), `update_foliation_features` + and the generated-feature data in the "constraints" mode. It still builds + the parametric faults and the manual foliations. +- [ ] `refresh_feature_data` and `model_state` ignore the column, the contacts + and the fault traces in the "constraints" mode. A change to the column + does not make the model stale. +- [ ] The data manager does not watch, reload or check the map input layers + (contacts, structure, fault traces) in the "constraints" mode. This stops + "Update data and solve" and the "Input layers changed" problem for layers + that the model does not use. `get_layers_outside_bounding_box` uses the + same rule. +- [ ] `sync_extra_constraints` and the read-only "processed" rows are empty in + the "constraints" mode. +- [ ] Step checks: `check_data` does not ask for the geology and structure + layers in the "constraints" mode. Move these two items to + `check_stratigraphy` (see 6.2). +- [ ] The primary button text and tooltip name the workflow ("Build model from + constraints"). +- [ ] Save the choice with the state (done in phase 4). Load old state files + with the "map" mode. + +Acceptance: a user has a map project with a column, contacts and faults. The +user changes the choice to "constraints" and adds one foliation. The build gives +one feature. The solved model does not change when the user edits the column or +the contacts layer. When the user changes back to "map", the generated features +are made again from the same data. + +Tests: unit tests of the model manager with a fake column and fake contacts: the +"constraints" mode builds only the manual features; the "map" mode builds both. +Unit tests for the checks and for `choose_primary_action` in each mode. + +#### 6.2 Stratigraphic layers move to step 2 + +Problem: the "Stratigraphic Layers" group (basal contacts and structural +orientations) is in step 1. It belongs to the column workflow. Its layer picker +shows only line and point layers. With "Calculate from geology polygons", the +input is the polygon geology layer, so the user cannot select it there. + +Tasks: + +- [ ] Remove `StratigraphicLayersWidget` from `ModelDefinitionTab`. Step 1 has + only the bounding box, the CRS and the DEM. +- [ ] Add the widget to `StratigraphyStep`, as a collapsible group (see 6.3). +- [ ] The first layer picker depends on the contacts source: + - "Calculate from geology polygons": the picker shows polygon layers. It + reads and writes the `geology` and `geology_unit_field` roles. It does not + call `set_basal_contacts`. The Z-coordinate check box is hidden. + - "Use a contacts layer": the picker shows line and point layers, as now. It + reads and writes the `basal_contacts` role. + - Change the group title and the label to match ("Geology layer" or + "Contacts layer"). +- [ ] When the user changes the source, do not write the old selection to the + other role. Keep one selection for each source. +- [ ] The `set_basal_contacts` callback (a tool or a build made a new contacts + layer) does not change the picker in the "geology" source. The layer is for + display. +- [ ] Keep the geology picker of the column group in sync with the new picker + (both use the same roles). Decide in the review if one of them must be + removed (see the open questions). +- [ ] Move the "Select the geology layer" and "Select the structure layer" + items from `check_data` to `check_stratigraphy`. +- [ ] Update the text that says "in step 1" for these layers + (`derived_refresh.py`, `checks.py`, `pages.py`, the docs). +- [ ] Keep the saved widget settings key `stratigraphic_layers_widget`, so old + state files load. + +Acceptance: with "Calculate from geology polygons", the user selects the +polygon layer in the group in step 2, and the unit name field list shows the +fields of that layer. After the user changes to "Use a contacts layer", the +picker shows only line and point layers. The selection of each source is kept. + +Tests: unit tests for the role writes of each source. A Qt test of the picker +filter for each source (QGIS test job). + +#### 6.3 Limit the vertical stack of widgets + +Problem: the pages put many widgets one above the other. Step 2 has the +geology group, the button row, the column list, the style group, the derived +data panel and the "Derive from map" row. Adding the stratigraphic layers (6.2) +makes it worse. The group boxes, the feature details panel (its own scroll +area) and the dock header, step bar and footer use the height. On a laptop +screen the part with the content is very small. + +Rules: + +- A page has at most one scroll area, at the page level. A widget inside a page + does not make its own scroll area. +- A page shows at most two expanded sections at one time. When the user + expands a third section, the section that was expanded first collapses. +- The main widget of a page (the column list, the feature list) is not in a + collapsible section. Its minimum height is 120 px. It gets the extra height. +- Do not nest a group box in a group box more than one level. +- A section that the user needs only sometimes starts collapsed ("Style map + layer", "Derived data", "Stratigraphic layers" after the first setup). + +Tasks: + +- [ ] Add a small `SectionStack` widget in `gui/modelling/steps/`. It holds + collapsible sections, the maximum number of open sections, and the + collapsed state of each section. It saves the state in the widget + settings. +- [ ] Use it in steps 1, 2, 3 and in the feature details panel of step 4 + (`Data Layers`, `Interpolator Settings`, `Preview`, `Export Feature`). +- [ ] Remove the scroll areas that are inside other scroll areas + (`BaseTab(scrollable=True)`, the scroll area of + `feature_details_panel/_base.py`). +- [ ] Give the page a header summary for each collapsed section, for example + "Geology layer: Geology, UNITNAME", so the user does not need to expand + it to read the value. +- [ ] Reduce the vertical use of the dock: the header and the footer use one + row each. The footer text is one line with a tooltip. +- [ ] In step 4, the problems list is collapsed to one line ("3 problems") with + the list in a tooltip or a popup. + +Acceptance: at a window height of 700 px, the main widget of each step has at +least 200 px. No page has a scroll area inside a scroll area. + +Tests: unit test of the open-section limit (it does not need QGIS if the logic +is in a plain class). Manual check at 700 px and 1080 px. + +#### 6.4 Choose the number of elements from the data + +Problem: each interpolator has the fixed number of elements of the settings +(default 50 000). A feature with 10 points and a feature with 10 000 points +get the same number. A small number gives a coarse surface. A large number +makes the solve slow. In addition, the model manager reads the default value +of `PlgSettingsStructure` (the class), not the value that the user saved. The +setting of the user has no effect on a build. User-added foliations do not set +the number at all. + +Design: + +- Add a pure function `suggest_nelements(summary)` in + `main/interpolation_size.py` (no QGIS import). The `summary` has the number + of value constraints, the number of orientation constraints, the number of + different values (surfaces) and a measure of the spread of the orientations + (0 for parallel planes, 1 for all directions). +- Rule of thumb: elements = equations x 25 x surface factor x spread factor, + where an orientation counts as two equations, the surface factor is + 1 + 0.1 for each surface after the first (at most 10 surfaces), and the + spread factor is 1 + spread. Round to 1 000. Limit to 5 000 .. 250 000. A + feature with no data gets the minimum. These numbers are a first guess: check + them with real models (see the open questions). +- Add the setting `interpolator_nelements_auto` (default: on). With it on, + the settings page shows the number as "Automatic" and the spin box is + disabled. With it off, the fixed number is used. +- Add one method in the model manager that gives the interpolator arguments + (`nelements`, `npw`, `cpw`, `regularisation`) for a data frame. It reads the + saved settings, not the class defaults. Use it for the generated features, + the faults, the domain faults, the parametric faults and the foliations that + the user added. +- Show the number that was used in the feature details panel ("Elements: 24 000 + (automatic)"). If the user changes it there, the feature keeps the number of + the user until the user selects "Automatic" again. + +Tasks: + +- [ ] `main/interpolation_size.py` with `summarise_data` and `suggest_nelements`. +- [ ] Setting, settings page and preference test. +- [ ] Model manager: one method for the interpolator arguments, used in all + build paths. Fix the use of the class defaults. +- [ ] Feature panel: show the number, and an "Automatic" check box. +- [ ] Docs: say how the number is chosen, and what the user can change. + +Acceptance: a model with 20 contact points and a model with 5 000 points get +different numbers of elements, both inside the limits. A saved fixed number is +used when "Automatic" is off. A user-added foliation gets a number. + +Tests: unit tests for the function: more data gives more or equal elements; the +limits; no data; parallel and spread orientations; orientation signs do not +matter. + +Order and links between the parts: 6.1 first, because it defines what the model +reads in each workflow. 6.2 depends on the same checks, so do it next. 6.3 +comes after 6.2, because it must lay out the final content of the pages. 6.4 +does not depend on the others. ## Risks @@ -355,6 +565,15 @@ model, the model uses basal contacts and thicknesses that match the new order. overwrites these edits. Show a warning before the overwrite, and offer to change the contacts source to "Use a contacts layer". +- **Saved fixed number of elements (phase 6.4).** Now the saved setting has no + effect on a build. After the fix, a user who saved a small or large number + sees a change in the result. Make "Automatic" the default, and tell the user + in the change log. +- **Hidden data (phase 6.1).** The "constraints" mode keeps the map data but + does not use it. Show a short message in step 4 ("The column and the map + layers are not used in this mode"), so the user knows why a unit is not in + the model. + ## Open questions 1. Must the steps be strict (the user cannot go to step 4 before step 2 is @@ -369,3 +588,12 @@ model, the model uses basal contacts and thicknesses that match the new order. 5. Must the Processing algorithms also use the derived-data record, or only the dock? Recommendation: only the dock. Processing runs are single runs with explicit inputs. +6. (6.2) After the move, step 2 has two geology pickers (the column group and + the stratigraphic layers group). Must one of them be removed? + Recommendation: keep the picker in the stratigraphic layers group, and show + the geology layer in the column group as a read-only summary. +7. (6.4) Are the numbers of the element rule right? Test it with 3 or 4 real + models (small and large, folded and simple) before the default is "Automatic". + Does a fault need a different rule from a foliation? +8. (6.1) Must a manual unconformity or fold in the "constraints" mode use the + column? Recommendation: no. It uses only the features of that mode. diff --git a/docs/usage/interface.md b/docs/usage/interface.md index a333f5b..83b247e 100644 --- a/docs/usage/interface.md +++ b/docs/usage/interface.md @@ -7,7 +7,18 @@ The LoopStructural dock has five steps. The step buttons are at the top of the d 2. **Stratigraphy**: the stratigraphic column. The **Build column** menu has the sorters and Paint Order. The **Derive from map** menu has Basal Contacts, Thickness and Sampler. 3. **Faults**: the fault layer, **Calculate topology...** and the adjacency tables. This step is optional. 4. **Model**: the features and the build of the model. -5. **View**: the 3D view. With the setting "separate dock widgets", this step has a button that opens the 3D view dock. +5. **View**: the 3D view and the **Export** panel. With the setting "separate dock widgets", this step has a button that opens the 3D view dock. + +When a build finishes, a message tells you to go to step 5, and the **Next** button in step 4 changes to **Open 3D view**. + +### Export +The **Export** panel in step 5 is open when the model is solved. It has three tabs: + +- **Surfaces**: one file for each stratigraphic surface and fault surface, in a folder. Formats: VTK (.vtk, .vtp), PLY and STL. A surface with no geometry is skipped. +- **Block model**: a grid that fills the bounding box. Set the number of cells in X, Y and Z. Each cell has the fields `stratigraphy_id` and `unit`. Formats: VTK (.vtk, .vti) and CSV. +- **Cross-section**: a vertical section under a line layer (the selected line, or the first line), or a plane from an origin and a normal. Formats: VTK and CSV. + +The exports run in the background and show a message when they finish. The header of the dock has **Save**, **Open**, **Reset** and **Settings**. They apply to all of the plugin. 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 3550252..6686e39 100644 --- a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py +++ b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py @@ -86,6 +86,9 @@ def _build_status_icon(color: str, *, filled: bool, mark: str = None) -> QIcon: class GeologicalModelTab(QWidget): + # Emitted when a solve finished without an error. The dock offers the 3D view. + model_solved = pyqtSignal() + def __init__(self, parent=None, *, model_manager=None, data_manager=None): super().__init__(parent) self.model_manager = model_manager @@ -586,7 +589,9 @@ def _on_task_progress(self, message, current, total): @pyqtSlot() def _on_task_finished(self): - on_success = None if self._task_failed else self._task_on_success + failed = self._task_failed + solving = self._task_title in ("Solving Model", "Building Model") + on_success = None if failed else self._task_on_success try: # notify observers now on the GUI thread try: @@ -601,6 +606,8 @@ def _on_task_finished(self): self._finish_task() if on_success is not None: on_success() + if not failed and solving and self.model_manager.model_state == 'solved': + self.model_solved.emit() @pyqtSlot(str) def _on_task_cancelled(self, message): diff --git a/loopstructural/gui/modelling/modelling_widget.py b/loopstructural/gui/modelling/modelling_widget.py index 0e12f65..cf1b4be 100644 --- a/loopstructural/gui/modelling/modelling_widget.py +++ b/loopstructural/gui/modelling/modelling_widget.py @@ -8,6 +8,7 @@ QWidget, ) +from loopstructural.gui.messages import push_success from loopstructural.gui.modelling.steps import build_plan, checks from loopstructural.gui.modelling.steps.header import DockHeader from loopstructural.gui.modelling.steps.pages import ( @@ -106,6 +107,7 @@ def __init__( self.step_bar.currentChanged.connect(self.go_to) self._last_checks = {} self.model_step.tab.set_problems_provider(self.problems) + self.model_step.tab.model_solved.connect(self._on_model_solved) if self.data_manager is not None: self.data_manager.add_workflow_mode_callback(self._apply_workflow_mode) self._timer = QTimer(self) @@ -156,6 +158,26 @@ def _apply_workflow_mode(self, mode): else: self.refresh_status() + def _on_model_solved(self): + """After a successful build, offer the 3D view as the next action.""" + self.refresh_status() + push_success("Model solved", "Go to step 5 to see the model in 3D and to export it.") + + def _next_text(self, index): + """The text of the Next button. After a solve, it names the 3D view.""" + next_index = self._neighbour(1, index) + if ( + next_index != index + and self.pages[index].key == checks.STEP_MODEL + and self.pages[next_index].key == checks.STEP_VIEW + and self._is_solved() + ): + return "Open 3D view >" + return "Next >" + + def _is_solved(self): + return self.model_manager is not None and self.model_manager.model_state == 'solved' + def problems(self): """Return the problems of the steps that show, as ``(step_key, message)`` pairs.""" visible = self._visible_indices() @@ -179,6 +201,8 @@ def refresh_status(self): self._last_checks = results self.step_bar.set_checks(results) self.model_step.tab.refresh_primary_action() + self.view_step.refresh() + self.next_button.setText(self._next_text(self.current_index)) page = self.pages[self.current_index] check = results[page.key] if check.messages: diff --git a/loopstructural/gui/modelling/steps/export_panel.py b/loopstructural/gui/modelling/steps/export_panel.py new file mode 100644 index 0000000..e9f8959 --- /dev/null +++ b/loopstructural/gui/modelling/steps/export_panel.py @@ -0,0 +1,288 @@ +"""The export panel of step 5: surfaces, block model and cross-sections.""" + +from pathlib import Path + +from qgis.core import QgsMapLayerProxyModel +from qgis.gui import QgsCollapsibleGroupBox, QgsFileWidget, QgsMapLayerComboBox +from qgis.PyQt.QtWidgets import ( + QCheckBox, + QComboBox, + QDoubleSpinBox, + QFormLayout, + QLabel, + QMessageBox, + QPushButton, + QSpinBox, + QTabWidget, + QVBoxLayout, + QWidget, +) + +from loopstructural.main import model_export + +from ...background_task import finish_background_task, start_background_task +from ...compatibility import configure_layer_combo +from ...messages import push_success +from ...visualisation.line_layer import line_xy_from_layer + + +def _format_combo(parent, formats): + combo = QComboBox(parent) + for extension, label in formats.items(): + combo.addItem(label, extension) + return combo + + +def _spin(parent, value, minimum=1, maximum=2000): + spin = QSpinBox(parent) + spin.setRange(minimum, maximum) + spin.setValue(value) + return spin + + +def _double_spin(parent, value=0.0): + spin = QDoubleSpinBox(parent) + spin.setRange(-1e9, 1e9) + spin.setDecimals(3) + spin.setValue(value) + return spin + + +class ExportPanel(QgsCollapsibleGroupBox): + """Write the results of a solved model to files. + + The exports run on a background task. The panel is enabled when the model + is solved (see `set_model_solved`). + """ + + def __init__(self, parent=None, *, model_manager=None, data_manager=None): + super().__init__("Export", parent) + self.model_manager = model_manager + self.data_manager = data_manager + self._task = None # (thread, worker, progress dialog) of the run + self._model_solved = False + + layout = QVBoxLayout(self) + self.hint = QLabel(self) + self.hint.setWordWrap(True) + layout.addWidget(self.hint) + self.tabs = QTabWidget(self) + self.tabs.addTab(self._surfaces_tab(), "Surfaces") + self.tabs.addTab(self._block_model_tab(), "Block model") + self.tabs.addTab(self._cross_section_tab(), "Cross-section") + layout.addWidget(self.tabs) + self.set_model_solved(False) + self.setCollapsed(True) + + # -- tabs --------------------------------------------------------------- + + def _surfaces_tab(self): + tab = QWidget(self) + form = QFormLayout(tab) + self.surface_format = _format_combo(tab, model_export.SURFACE_FORMATS) + self.surface_strat = QCheckBox("Stratigraphic surfaces", tab) + self.surface_strat.setChecked(True) + self.surface_faults = QCheckBox("Fault surfaces", tab) + self.surface_faults.setChecked(True) + self.surface_folder = QgsFileWidget(tab) + self.surface_folder.setStorageMode(QgsFileWidget.StorageMode.GetDirectory) + self.surface_button = QPushButton("Export surfaces", tab) + self.surface_button.clicked.connect(self.export_surfaces) + form.addRow("Format", self.surface_format) + form.addRow(self.surface_strat) + form.addRow(self.surface_faults) + form.addRow("Folder", self.surface_folder) + form.addRow(self.surface_button) + return tab + + def _block_model_tab(self): + tab = QWidget(self) + form = QFormLayout(tab) + self.block_format = _format_combo(tab, model_export.BLOCK_MODEL_FORMATS) + self.block_nx = _spin(tab, 50) + self.block_ny = _spin(tab, 50) + self.block_nz = _spin(tab, 25) + self.block_file = QgsFileWidget(tab) + self.block_file.setStorageMode(QgsFileWidget.StorageMode.SaveFile) + self.block_button = QPushButton("Export block model", tab) + self.block_button.clicked.connect(self.export_block_model) + form.addRow("Format", self.block_format) + form.addRow("Cells in X", self.block_nx) + form.addRow("Cells in Y", self.block_ny) + form.addRow("Cells in Z", self.block_nz) + form.addRow("File", self.block_file) + form.addRow(self.block_button) + return tab + + def _cross_section_tab(self): + tab = QWidget(self) + form = QFormLayout(tab) + self.section_format = _format_combo(tab, model_export.CROSS_SECTION_FORMATS) + self.section_kind = QComboBox(tab) + self.section_kind.addItem("Vertical section under a line layer", 'line') + self.section_kind.addItem("Plane (origin and normal)", 'plane') + self.section_kind.currentIndexChanged.connect(self._update_section_inputs) + self.section_line = QgsMapLayerComboBox(tab) + configure_layer_combo(self.section_line, QgsMapLayerProxyModel.Filter.LineLayer) + self.section_origin = [_double_spin(tab) for _ in range(3)] + self.section_normal = [_double_spin(tab, v) for v in (1.0, 0.0, 0.0)] + self.section_size = _double_spin(tab, 1000.0) + self.section_resolution = _spin(tab, 100) + self.section_file = QgsFileWidget(tab) + self.section_file.setStorageMode(QgsFileWidget.StorageMode.SaveFile) + self.section_button = QPushButton("Export cross-section", tab) + self.section_button.clicked.connect(self.export_cross_section) + form.addRow("Format", self.section_format) + form.addRow("Type", self.section_kind) + form.addRow("Line layer", self.section_line) + for label, widgets in (("Origin", self.section_origin), ("Normal", self.section_normal)): + for axis, widget in zip("XYZ", widgets): + form.addRow(f"{label} {axis}", widget) + form.addRow("Plane size", self.section_size) + form.addRow("Resolution", self.section_resolution) + form.addRow("File", self.section_file) + form.addRow(self.section_button) + self._section_plane_rows = ( + self.section_origin + self.section_normal + [self.section_size] + ) + self._section_form = form + self._update_section_inputs() + return tab + + def _update_section_inputs(self, *_args): + is_line = self.section_kind.currentData() == 'line' + self.section_line.setVisible(is_line) + self._section_form.labelForField(self.section_line).setVisible(is_line) + for widget in self._section_plane_rows: + widget.setVisible(not is_line) + self._section_form.labelForField(widget).setVisible(not is_line) + + # -- state -------------------------------------------------------------- + + def set_model_solved(self, solved): + """Enable the exports when the model is solved.""" + self._model_solved = bool(solved) + self.hint.setText( + "Write the results of the model to files." + if solved + else "Build the model in step 4 to export its results." + ) + for tab_index in range(self.tabs.count()): + self.tabs.widget(tab_index).setEnabled(self._model_solved and self._task is None) + + # -- actions ------------------------------------------------------------ + + def _with_extension(self, file_widget, format_combo): + """Return the output path with the extension of the format, or None if empty.""" + text = file_widget.filePath().strip() + if not text: + return None + return Path(text).with_suffix(format_combo.currentData()) + + def _warn(self, text): + QMessageBox.warning(self, "Export", text) + + def export_surfaces(self): + folder = self.surface_folder.filePath().strip() + if not folder: + return self._warn("Select a folder for the surfaces.") + if not (self.surface_strat.isChecked() or self.surface_faults.isChecked()): + return self._warn("Select the stratigraphic surfaces, the fault surfaces, or both.") + extension = self.surface_format.currentData() + strat, faults = self.surface_strat.isChecked(), self.surface_faults.isChecked() + + def target(progress): + files = model_export.export_surfaces( + self.model_manager, + folder, + extension, + stratigraphic=strat, + faults=faults, + progress=progress, + ) + return f"{len(files)} surfaces written to {folder}" + + self._run(target, "Export surfaces") + + def export_block_model(self): + path = self._with_extension(self.block_file, self.block_format) + if path is None: + return self._warn("Select a file for the block model.") + ncells = (self.block_nx.value(), self.block_ny.value(), self.block_nz.value()) + + def target(progress): + model_export.export_block_model(self.model_manager, path, ncells, progress=progress) + return f"Block model written to {path}" + + self._run(target, "Export block model") + + def export_cross_section(self): + path = self._with_extension(self.section_file, self.section_format) + if path is None: + return self._warn("Select a file for the cross-section.") + resolution = self.section_resolution.value() + if self.section_kind.currentData() == 'line': + layer = self.section_line.currentLayer() + if layer is None: + return self._warn("Select a line layer.") + target_crs = self.data_manager.get_model_crs() if self.data_manager else None + line_xy = line_xy_from_layer(layer, target_crs) + if line_xy is None: + return self._warn("The line layer has no usable line.") + options = {'line_xy': line_xy, 'z_resolution': resolution} + else: + options = { + 'origin': [w.value() for w in self.section_origin], + 'normal': [w.value() for w in self.section_normal], + 'size': self.section_size.value(), + } + + def target(progress): + model_export.export_cross_section( + self.model_manager, path, resolution=resolution, progress=progress, **options + ) + return f"Cross-section written to {path}" + + self._run(target, "Export cross-section") + + # -- background task ---------------------------------------------------- + + def _run(self, target, title): + if not self._model_solved or self._task is not None: + return + self._task = start_background_task( + self, + target, + title=title, + initial_label=f"{title}...", + on_progress=self._on_progress, + on_finished=self._on_finished, + on_error=self._on_error, + ) + self.set_model_solved(self._model_solved) + + def _end_task(self): + if self._task is not None: + finish_background_task(*self._task) + self._task = None + self.set_model_solved(self._model_solved) + + def _on_progress(self, message): + try: + self._task[2].setLabelText(message) + except Exception: + pass + + def _on_finished(self, message): + self._end_task() + push_success("Export", message) + + def _on_error(self, traceback_text): + self._end_task() + reason = traceback_text.strip().splitlines()[-1] if traceback_text.strip() else "unknown" + box = QMessageBox(self) + box.setIcon(QMessageBox.Icon.Critical) + box.setWindowTitle("Export failed") + box.setText(reason) + box.setDetailedText(traceback_text) + box.exec() diff --git a/loopstructural/gui/modelling/steps/pages.py b/loopstructural/gui/modelling/steps/pages.py index 847fc84..5adaa8b 100644 --- a/loopstructural/gui/modelling/steps/pages.py +++ b/loopstructural/gui/modelling/steps/pages.py @@ -28,6 +28,7 @@ from loopstructural.main.workflow_mode import WORKFLOW_MODE_LABELS, WORKFLOW_MODES from . import checks +from .export_panel import ExportPanel class StepPage(QWidget): @@ -193,7 +194,7 @@ def __init__(self, parent=None, **kwargs): class ViewStep(StepPage): - """Step 5: the 3D view. + """Step 5: the 3D view and the export of the results. In one dock, the visualisation widget is in this page. With separate docks, the page has a button that shows the visualisation dock. @@ -218,3 +219,12 @@ def __init__(self, parent=None, *, view_widget=None, **kwargs): layout.addWidget(label) layout.addWidget(button) layout.addStretch(1) + self.export_panel = ExportPanel( + self, model_manager=self.model_manager, data_manager=self.data_manager + ) + layout.addWidget(self.export_panel) + + def refresh(self): + """Enable the exports when the model is solved.""" + state = self.model_manager.model_state if self.model_manager is not None else 'empty' + self.export_panel.set_model_solved(state == 'solved') diff --git a/loopstructural/gui/visualisation/feature_list_widget.py b/loopstructural/gui/visualisation/feature_list_widget.py index e1ba2c8..a174c7f 100644 --- a/loopstructural/gui/visualisation/feature_list_widget.py +++ b/loopstructural/gui/visualisation/feature_list_widget.py @@ -4,14 +4,7 @@ import numpy as np import pyvista as pv from LoopStructural.datatypes import VectorPoints -from qgis.core import ( - QgsApplication, - QgsCoordinateTransform, - QgsGeometry, - QgsMapLayerProxyModel, - QgsProject, - QgsWkbTypes, -) +from qgis.core import QgsApplication, QgsMapLayerProxyModel from qgis.gui import QgsMapLayerComboBox from qgis.PyQt.QtCore import QSize from qgis.PyQt.QtGui import QIcon @@ -43,6 +36,7 @@ build_line_extrusion_mesh, build_plane_mesh, ) +from .line_layer import line_xy_from_layer from .mesh_scalar_utils import stratigraphic_ids_to_rgb logger = logging.getLogger(__name__) @@ -963,41 +957,13 @@ def _extract_line_xy(self, layer) -> Optional[np.ndarray]: selected feature (or first feature, if nothing is selected), reprojected into the model's CRS if one is available. """ - features = ( - list(layer.getSelectedFeatures()) - if layer.selectedFeatureCount() > 0 - else list(layer.getFeatures()) - ) - if not features: - return None - geom = features[0].geometry() - if geom is None or geom.isEmpty(): - return None - target_crs = None if self.data_manager is not None: try: target_crs = self.data_manager.get_model_crs() 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 - ): - geom = QgsGeometry(geom) - geom.transform(QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance())) - - if QgsWkbTypes.isMultiType(geom.wkbType()): - parts = geom.asMultiPolyline() - polyline = parts[0] if parts else [] - else: - polyline = geom.asPolyline() - if len(polyline) < 2: - return None - return np.array([[pt.x(), pt.y()] for pt in polyline]) + return line_xy_from_layer(layer, target_crs) def add_line_cross_section(self): """Extrude the selected QGIS line layer vertically across the model's diff --git a/loopstructural/gui/visualisation/line_layer.py b/loopstructural/gui/visualisation/line_layer.py new file mode 100644 index 0000000..064ebf1 --- /dev/null +++ b/loopstructural/gui/visualisation/line_layer.py @@ -0,0 +1,44 @@ +"""Read the line of a cross-section from a QGIS line layer.""" + +from typing import Optional + +import numpy as np +from qgis.core import QgsCoordinateTransform, QgsGeometry, QgsProject, QgsWkbTypes + + +def line_xy_from_layer(layer, target_crs=None) -> Optional[np.ndarray]: + """Return an (M, 2) array of the vertices of the first line of a layer. + + The line is the first selected feature, or the first feature if nothing is + selected. It is reprojected to ``target_crs`` if this is valid. The result + is None if the layer has no usable line. + """ + features = ( + list(layer.getSelectedFeatures()) + if layer.selectedFeatureCount() > 0 + else list(layer.getFeatures()) + ) + if not features: + return None + geom = features[0].geometry() + if geom is None or geom.isEmpty(): + return None + + source_crs = layer.sourceCrs() + 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())) + + if QgsWkbTypes.isMultiType(geom.wkbType()): + parts = geom.asMultiPolyline() + polyline = parts[0] if parts else [] + else: + polyline = geom.asPolyline() + if len(polyline) < 2: + return None + return np.array([[pt.x(), pt.y()] for pt in polyline]) diff --git a/loopstructural/main/model_export.py b/loopstructural/main/model_export.py new file mode 100644 index 0000000..5602fad --- /dev/null +++ b/loopstructural/main/model_export.py @@ -0,0 +1,212 @@ +"""Export of the results of a solved model to files. + +The functions do not import QGIS, so the unit tests can run them. They read +the model through the model manager (`model`, `evaluate_stratigraphy_on_points` +and `get_stratigraphic_unit_names`) and write files with pyvista. +""" + +import re +from pathlib import Path + +import numpy as np + +from loopstructural.gui.visualisation.cross_section_utils import ( + build_block_model_mesh, + build_line_extrusion_mesh, + build_plane_mesh, +) + +# The file formats of each kind of result: extension -> label +SURFACE_FORMATS = {'.vtk': 'VTK legacy (.vtk)', '.vtp': 'VTK PolyData (.vtp)', '.ply': 'PLY (.ply)', '.stl': 'STL (.stl)'} +BLOCK_MODEL_FORMATS = {'.vtk': 'VTK legacy (.vtk)', '.vti': 'VTK ImageData (.vti)', '.csv': 'CSV (.csv)'} +CROSS_SECTION_FORMATS = {'.vtk': 'VTK legacy (.vtk)', '.csv': 'CSV (.csv)'} + +STRATIGRAPHY_ID_FIELD = 'stratigraphy_id' +UNIT_FIELD = 'unit' + + +class ExportError(Exception): + """The export cannot run, for example because the model is not solved.""" + + +def safe_file_name(name: str) -> str: + """Return ``name`` with only letters, digits, ``-``, ``_`` and ``.``.""" + cleaned = re.sub(r'[^A-Za-z0-9_.-]+', '_', str(name)).strip('._') + return cleaned or 'surface' + + +def unique_names(names): + """Return a file-safe name for each name, with a suffix for the repeats.""" + seen = {} + result = [] + for name in names: + base = safe_file_name(name) + count = seen.get(base, 0) + seen[base] = count + 1 + result.append(base if count == 0 else f"{base}_{count + 1}") + return result + + +def _require_model(model_manager): + model = getattr(model_manager, 'model', None) + if model is None: + raise ExportError("There is no model to export. Build the model in step 4 first.") + return model + + +def _check_extension(extension, formats): + extension = extension.lower() + if extension not in formats: + raise ExportError(f"Unknown file format: {extension}") + return extension + + +def unit_names_for_ids(model_manager, ids) -> np.ndarray: + """Return the unit name for each stratigraphic id. An id of -1 has no unit.""" + names = list(model_manager.get_stratigraphic_unit_names()) + ids = np.asarray(ids, dtype=int) + return np.array( + [names[i] if 0 <= i < len(names) else '' for i in ids], dtype=object + ) + + +def add_stratigraphy(model_manager, mesh, *, cells=False): + """Evaluate the stratigraphic column on a mesh and add the result to it. + + Adds the fields ``stratigraphy_id`` and ``unit``, to the cells or to the points. + """ + points = mesh.cell_centers().points if cells else mesh.points + ids = np.asarray(model_manager.evaluate_stratigraphy_on_points(points)) + data = mesh.cell_data if cells else mesh.point_data + data[STRATIGRAPHY_ID_FIELD] = ids + return ids + + +def _write_csv(path, points, ids, names): + points = np.asarray(points) + with open(path, 'w', encoding='utf-8', newline='') as handle: + handle.write(f"x,y,z,{STRATIGRAPHY_ID_FIELD},{UNIT_FIELD}\n") + for (x, y, z), unit_id, name in zip(points, ids, names): + handle.write(f"{x:.6f},{y:.6f},{z:.6f},{int(unit_id)},\"{name}\"\n") + + +def export_surfaces( + model_manager, + folder, + extension='.vtk', + *, + stratigraphic=True, + faults=True, + progress=None, +): + """Write one file for each surface of the model into a folder. + + Parameters + ---------- + folder : str or Path + The output folder. It is made if it does not exist. + extension : str + One of `SURFACE_FORMATS`. + stratigraphic, faults : bool + Which surfaces to write. + progress : callable, optional + Called with a message before each surface. + + Returns + ------- + list of Path + The files that were written. A surface with no geometry is skipped: a + unit with no data of its own can have an isovalue that the field does + not reach. + """ + model = _require_model(model_manager) + extension = _check_extension(extension, SURFACE_FORMATS) + surfaces = [] + if stratigraphic: + surfaces += [(f"strat_{s.name}", s) for s in model.get_stratigraphic_surfaces()] + if faults: + surfaces += [(f"fault_{s.name}", s) for s in model.get_fault_surfaces()] + if not surfaces: + raise ExportError("The model has no surfaces to export.") + folder = Path(folder) + folder.mkdir(parents=True, exist_ok=True) + written = [] + for file_name, (_label, surface) in zip(unique_names(l for l, _ in surfaces), surfaces): + if progress is not None: + progress(f"Writing {surface.name}...") + mesh = surface.vtk() + if mesh.n_points == 0: + continue + path = folder / f"{file_name}{extension}" + mesh.save(str(path)) + written.append(path) + if not written: + raise ExportError("All surfaces are empty. Check the bounding box and the data.") + return written + + +def export_block_model(model_manager, path, ncells, *, progress=None): + """Write the block model: a grid that fills the bounding box. + + Each cell has the stratigraphic id and the unit name at its centre. The + format comes from the extension of ``path`` (see `BLOCK_MODEL_FORMATS`). + """ + model = _require_model(model_manager) + path = Path(path) + extension = _check_extension(path.suffix, BLOCK_MODEL_FORMATS) + if progress is not None: + progress("Making the block model grid...") + mesh = build_block_model_mesh(model.bounding_box.origin, model.bounding_box.maximum, ncells) + if progress is not None: + progress("Evaluating the stratigraphy in the blocks...") + ids = add_stratigraphy(model_manager, mesh, cells=True) + path.parent.mkdir(parents=True, exist_ok=True) + if extension == '.csv': + _write_csv(path, mesh.cell_centers().points, ids, unit_names_for_ids(model_manager, ids)) + else: + mesh.save(str(path)) + return path + + +def export_cross_section( + model_manager, + path, + *, + origin=None, + normal=None, + size=None, + line_xy=None, + resolution=100, + z_resolution=100, + progress=None, +): + """Write a cross-section: a plane (``origin``, ``normal``, ``size``) or a + vertical section under a line (``line_xy``, an (M, 2) array). + """ + model = _require_model(model_manager) + path = Path(path) + extension = _check_extension(path.suffix, CROSS_SECTION_FORMATS) + if line_xy is not None: + bounding_box = model.bounding_box + mesh = build_line_extrusion_mesh( + line_xy, + float(bounding_box.origin[2]), + float(bounding_box.maximum[2]), + resolution=resolution, + z_resolution=z_resolution, + ) + elif origin is not None and normal is not None and size: + if np.allclose(normal, 0.0): + raise ExportError("The normal of the cross-section cannot be the zero vector.") + mesh = build_plane_mesh(origin, normal, size, resolution) + else: + raise ExportError("Give a line, or an origin, a normal and a size.") + if progress is not None: + progress("Evaluating the stratigraphy on the cross-section...") + ids = add_stratigraphy(model_manager, mesh) + path.parent.mkdir(parents=True, exist_ok=True) + if extension == '.csv': + _write_csv(path, mesh.points, ids, unit_names_for_ids(model_manager, ids)) + else: + mesh.save(str(path)) + return path diff --git a/tests/unit/test_model_export.py b/tests/unit/test_model_export.py new file mode 100644 index 0000000..0f546eb --- /dev/null +++ b/tests/unit/test_model_export.py @@ -0,0 +1,123 @@ +"""Pytest tests for `main.model_export`. + +A fake model manager gives surfaces and a stratigraphy, so only pyvista and +numpy are needed. +""" + +from types import SimpleNamespace + +import numpy as np +import pyvista as pv +import pytest + +from loopstructural.main import model_export + + +class FakeSurface: + def __init__(self, name, empty=False): + self.name = name + self._empty = empty + + def vtk(self): + return pv.PolyData() if self._empty else pv.Sphere(radius=1.0) + + +class FakeModel: + def __init__(self, strat=(), faults=()): + self.bounding_box = SimpleNamespace(origin=(0, 0, 0), maximum=(10, 10, 10)) + self._strat, self._faults = list(strat), list(faults) + + def get_stratigraphic_surfaces(self): + return self._strat + + def get_fault_surfaces(self): + return self._faults + + +class FakeManager: + def __init__(self, model): + self.model = model + + def evaluate_stratigraphy_on_points(self, points): + # unit 0 below z = 5, unit 1 above + return (np.asarray(points)[:, 2] > 5).astype(int) + + def get_stratigraphic_unit_names(self): + return ['lower', 'upper'] + + +def test_safe_file_name_and_unique_names(): + assert model_export.safe_file_name('a b/c') == 'a_b_c' + assert model_export.safe_file_name('///') == 'surface' + assert model_export.unique_names(['a b', 'a_b', 'a_b']) == ['a_b', 'a_b_2', 'a_b_3'] + + +def test_export_surfaces_writes_one_file_each_and_skips_empty(tmp_path): + model = FakeModel( + strat=[FakeSurface('unit a'), FakeSurface('top', empty=True)], faults=[FakeSurface('F1')] + ) + files = model_export.export_surfaces(FakeManager(model), tmp_path / 'out', '.vtk') + assert sorted(f.name for f in files) == ['fault_F1.vtk', 'strat_unit_a.vtk'] + assert all(f.exists() for f in files) + + +def test_export_surfaces_can_select_kinds(tmp_path): + model = FakeModel(strat=[FakeSurface('a')], faults=[FakeSurface('F1')]) + files = model_export.export_surfaces( + FakeManager(model), tmp_path, '.vtp', stratigraphic=False + ) + assert [f.name for f in files] == ['fault_F1.vtp'] + + +def test_export_surfaces_errors(tmp_path): + with pytest.raises(model_export.ExportError): + model_export.export_surfaces(FakeManager(FakeModel()), tmp_path) + with pytest.raises(model_export.ExportError): + model_export.export_surfaces(FakeManager(None), tmp_path) + with pytest.raises(model_export.ExportError): + model_export.export_surfaces( + FakeManager(FakeModel(strat=[FakeSurface('a', empty=True)])), tmp_path + ) + with pytest.raises(model_export.ExportError): + model_export.export_surfaces(FakeManager(FakeModel(strat=[FakeSurface('a')])), tmp_path, '.xyz') + + +def test_export_block_model_vtk_has_stratigraphy(tmp_path): + path = model_export.export_block_model( + FakeManager(FakeModel()), tmp_path / 'block.vtk', (2, 2, 2) + ) + mesh = pv.read(str(path)) + assert mesh.n_cells == 8 + assert sorted(set(mesh.cell_data[model_export.STRATIGRAPHY_ID_FIELD])) == [0, 1] + + +def test_export_block_model_csv(tmp_path): + path = model_export.export_block_model( + FakeManager(FakeModel()), tmp_path / 'block.csv', (1, 1, 2) + ) + lines = path.read_text(encoding='utf-8').splitlines() + assert lines[0] == 'x,y,z,stratigraphy_id,unit' + assert lines[1].endswith('0,"lower"') + assert lines[2].endswith('1,"upper"') + + +def test_export_cross_section_plane_and_line(tmp_path): + manager = FakeManager(FakeModel()) + plane = model_export.export_cross_section( + manager, tmp_path / 'plane.vtk', origin=(5, 5, 5), normal=(0, 1, 0), size=10, resolution=5 + ) + assert pv.read(str(plane)).n_cells == 25 + line = model_export.export_cross_section( + manager, tmp_path / 'line.csv', line_xy=[(0, 0), (10, 10)], resolution=4, z_resolution=3 + ) + assert len(line.read_text(encoding='utf-8').splitlines()) == 1 + 4 * 3 + + +def test_export_cross_section_needs_input(tmp_path): + manager = FakeManager(FakeModel()) + with pytest.raises(model_export.ExportError): + model_export.export_cross_section(manager, tmp_path / 'x.vtk') + with pytest.raises(model_export.ExportError): + model_export.export_cross_section( + manager, tmp_path / 'x.vtk', origin=(0, 0, 0), normal=(0, 0, 0), size=1 + )