commit f689cd6bd1d37813ce21b1d777a145cdc498dcea
parent 99c3ac3a0f150a3271b3b16e357d3efcd8048b80
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 22 Oct 2026 09:08:41 +0000
feat(qgis): add mean resistivity maps dialog
Diffstat:
2 files changed, 352 insertions(+), 6 deletions(-)
diff --git a/tem_loader/tem_loader.py b/tem_loader/tem_loader.py
@@ -30,6 +30,7 @@ from qgis.core import (
QgsGeometry,
QgsPointXY,
QgsProject,
+ QgsRasterLayer,
QgsVectorFileWriter,
QgsVectorLayer,
)
@@ -37,6 +38,8 @@ from qgis.gui import QgsMapLayerComboBox
from . import core
from . import interpolate as interp
+from . import mean_resistivity as mres
+from .geotiff import write_geotiff
STYLES_DIR = Path(__file__).parent / 'styles'
@@ -49,18 +52,26 @@ def _build_geopackage_uri(gpkg_path, layer_name):
return f'{gpkg_path.resolve()}|layername={layer_name}'
-def _voxel_output_path(input_layers):
- """Auto-pick the voxel GeoPackage path next to the import's GeoPackage,
- named after it with a ``_voxels`` suffix. Returns ``None`` if no input layer
- is backed by a GeoPackage on disk."""
+def _output_path(input_layers, suffix):
+ """Auto-pick an output path next to the import's GeoPackage, named after it
+ with ``suffix`` appended to its stem. Returns ``None`` if no input layer is
+ backed by a GeoPackage on disk."""
for lyr in input_layers:
source = lyr.source().split('|', 1)[0]
src = Path(source)
if src.suffix.lower() == '.gpkg' and src.parent.is_dir():
- return src.with_name(f'{src.stem}_voxels.gpkg')
+ return src.with_name(f'{src.stem}{suffix}')
return None
+def _mean_resistivity_file_name(top, bottom):
+ return f'mean_resistivity_{top:g}_{bottom:g}.tif'
+
+
+def _mean_resistivity_layer_name(top, bottom):
+ return f'Mean resistivity {top:g} to {bottom:g} m'
+
+
def _field_type(field_name):
if field_name in STRING_FIELDS:
return QMetaType.Type.QString
@@ -570,11 +581,186 @@ class _InterpolationDialog(QDialog):
}
+def _layer_list(layers, wkb_type, preselect=None):
+ widget = QListWidget()
+ widget.setSelectionMode(QListWidget.SelectionMode.MultiSelection
+ if hasattr(QListWidget, 'SelectionMode')
+ else QListWidget.MultiSelection)
+ for lyr in layers:
+ if hasattr(lyr, 'wkbType') and lyr.wkbType() == wkb_type:
+ item = QListWidgetItem(lyr.name())
+ item.setData(32, lyr)
+ widget.addItem(item)
+ if preselect and preselect in lyr.name().lower():
+ item.setSelected(True)
+ return widget
+
+
+def _selected_list_layers(widget):
+ layers = []
+ for i in range(widget.count()):
+ item = widget.item(i)
+ if item.isSelected() and item.data(32) is not None:
+ layers.append(item.data(32))
+ return layers
+
+
+class _MeanResistivityDialog(QDialog):
+ def __init__(self, parent=None):
+ super().__init__(parent)
+ self.setWindowTitle('Mean Resistivity Maps')
+
+ all_layers = list(QgsProject.instance().mapLayers().values())
+ self._input_list = _layer_list(all_layers, Qgis.WkbType.LineStringZ)
+ self._input_list.itemSelectionChanged.connect(self._update_summary)
+ self._doi_list = _layer_list(all_layers, Qgis.WkbType.PointZ, 'doi')
+
+ self._top = QDoubleSpinBox()
+ self._bottom = QDoubleSpinBox()
+ for sb in (self._top, self._bottom):
+ sb.setRange(-1e5, 1e5)
+ sb.setDecimals(1)
+ sb.setSuffix(' m')
+ self._top.setValue(50.0)
+ self._bottom.setValue(-50.0)
+
+ self._thickness = QDoubleSpinBox()
+ self._thickness.setRange(0.1, 1e4)
+ self._thickness.setDecimals(1)
+ self._thickness.setSuffix(' m')
+ self._thickness.setValue(10.0)
+
+ self._cell_size = QDoubleSpinBox()
+ self._cell_size.setRange(0.1, 1e5)
+ self._cell_size.setDecimals(1)
+ self._cell_size.setSuffix(' m')
+ self._cell_size.setValue(50.0)
+
+ self._radius = QDoubleSpinBox()
+ self._radius.setRange(0.0, 1e6)
+ self._radius.setDecimals(1)
+ self._radius.setValue(0.0)
+ self._radius.setSpecialValueText('Auto')
+
+ self._power = QDoubleSpinBox()
+ self._power.setRange(0.1, 10.0)
+ self._power.setDecimals(1)
+ self._power.setValue(2.0)
+
+ self._min_neighbors = QSpinBox()
+ self._min_neighbors.setRange(1, 1000)
+ self._min_neighbors.setValue(1)
+
+ self._min_coverage = QSpinBox()
+ self._min_coverage.setRange(0, 100)
+ self._min_coverage.setSuffix('%')
+ self._min_coverage.setValue(50)
+
+ self._summary_label = QLabel()
+
+ for sb in (self._top, self._bottom, self._thickness, self._cell_size):
+ sb.valueChanged.connect(self._update_summary)
+
+ use_range_btn = QPushButton(
+ 'Use elevation range of selected input layers'
+ )
+ use_range_btn.clicked.connect(self._fill_range_from_inputs)
+
+ buttons = QDialogButtonBox(_dialog_buttons())
+ buttons.accepted.connect(self.accept)
+ buttons.rejected.connect(self.reject)
+
+ interval_form = QFormLayout()
+ interval_form.addRow('Top elevation', self._top)
+ interval_form.addRow('Bottom elevation', self._bottom)
+ interval_form.addRow('Interval thickness', self._thickness)
+
+ param_form = QFormLayout()
+ param_form.addRow('Cell size', self._cell_size)
+ param_form.addRow('Search radius (m)', self._radius)
+ param_form.addRow('IDW power', self._power)
+ param_form.addRow('Min neighbors', self._min_neighbors)
+ param_form.addRow('Minimum interval coverage', self._min_coverage)
+ param_form.addRow('', self._summary_label)
+
+ layout = QVBoxLayout()
+ layout.addWidget(QLabel('Input layers (LineStringZ):'))
+ layout.addWidget(self._input_list)
+ layout.addWidget(QLabel('DOI layers (PointZ, optional):'))
+ layout.addWidget(self._doi_list)
+ layout.addWidget(use_range_btn)
+ layout.addLayout(interval_form)
+ layout.addLayout(param_form)
+ layout.addWidget(buttons)
+ self.setLayout(layout)
+
+ self._update_summary()
+
+ def _fill_range_from_inputs(self):
+ tops, bottoms = [], []
+ for lyr in _selected_list_layers(self._input_list):
+ fields = lyr.fields()
+ for name, values, pick in (('ZTop', tops, lyr.maximumValue),
+ ('ZBottom', bottoms, lyr.minimumValue)):
+ index = fields.indexOf(name)
+ if index < 0:
+ continue
+ value = pick(index)
+ if isinstance(value, (int, float)) and value == value:
+ values.append(float(value))
+ if not tops or not bottoms:
+ return
+ step = self._thickness.value()
+ self._top.setValue(math.ceil(max(tops) / step) * step)
+ self._bottom.setValue(math.floor(min(bottoms) / step) * step)
+
+ def _update_summary(self, *_args):
+ try:
+ intervals = mres.elevation_intervals(
+ self._top.value(), self._bottom.value(),
+ self._thickness.value(),
+ )
+ except ValueError as exc:
+ self._summary_label.setText(str(exc).capitalize() + '.')
+ return
+ text = f'{len(intervals)} intervals'
+ layers = _selected_list_layers(self._input_list)
+ if layers:
+ xs, ys = [], []
+ for lyr in layers:
+ ext = lyr.extent()
+ xs += [ext.xMinimum(), ext.xMaximum()]
+ ys += [ext.yMinimum(), ext.yMaximum()]
+ try:
+ grid = mres.map_grid(xs, ys, self._cell_size.value())
+ text += f', {grid.nx} × {grid.ny} cells'
+ except ValueError as exc:
+ text += f'; {exc}'
+ self._summary_label.setText(text)
+
+ def options(self):
+ radius = self._radius.value()
+ return {
+ 'input_layers': _selected_list_layers(self._input_list),
+ 'doi_layers': _selected_list_layers(self._doi_list),
+ 'intervals': mres.elevation_intervals(
+ self._top.value(), self._bottom.value(),
+ self._thickness.value(),
+ ),
+ 'cell_size': self._cell_size.value(),
+ 'radius': radius if radius > 0 else None,
+ 'power': self._power.value(),
+ 'min_neighbors': self._min_neighbors.value(),
+ 'min_coverage': self._min_coverage.value() / 100.0,
+ }
+
+
class TEMLoaderPlugin:
def __init__(self, iface):
self.iface = iface
self._action = None
self._interp_action = None
+ self._mean_res_action = None
def initGui(self):
self._action = QAction('Load TEM XYZ files\u2026', self.iface.mainWindow())
@@ -585,12 +771,21 @@ class TEMLoaderPlugin:
self._interp_action.triggered.connect(self.run_interpolation)
self.iface.addPluginToMenu('TEM Loader', self._interp_action)
+ self._mean_res_action = QAction(
+ 'Mean resistivity maps\u2026', self.iface.mainWindow()
+ )
+ self._mean_res_action.triggered.connect(self.run_mean_resistivity_maps)
+ self.iface.addPluginToMenu('TEM Loader', self._mean_res_action)
+
def unload(self):
self.iface.removePluginMenu('TEM Loader', self._action)
self._action = None
if self._interp_action is not None:
self.iface.removePluginMenu('TEM Loader', self._interp_action)
self._interp_action = None
+ if self._mean_res_action is not None:
+ self.iface.removePluginMenu('TEM Loader', self._mean_res_action)
+ self._mean_res_action = None
def run(self):
paths, _ = QFileDialog.getOpenFileNames(
@@ -660,7 +855,7 @@ class TEMLoaderPlugin:
)
return
- gpkg_path = _voxel_output_path(input_layers)
+ gpkg_path = _output_path(input_layers, '_voxels.gpkg')
if gpkg_path is None:
QMessageBox.warning(
self.iface.mainWindow(),
@@ -785,6 +980,101 @@ class TEMLoaderPlugin:
'TEM Loader', f'Voxel layer saved to {gpkg_path}'
)
+ def run_mean_resistivity_maps(self):
+ dialog = _MeanResistivityDialog(self.iface.mainWindow())
+ if _exec_dialog(dialog) != _dialog_code('Accepted'):
+ return
+ try:
+ self._run_mean_resistivity_maps(dialog.options())
+ except ValueError as exc:
+ QMessageBox.warning(self.iface.mainWindow(), 'TEM Loader', str(exc))
+
+ def _run_mean_resistivity_maps(self, opts):
+ input_layers = opts['input_layers']
+ if not input_layers:
+ raise ValueError('No input layers selected.')
+ required = ('X', 'Y', 'ZTop', 'ZBottom', 'Resistivity')
+ for lyr in input_layers:
+ names = {field.name() for field in lyr.fields()}
+ missing = [name for name in required if name not in names]
+ if missing:
+ raise ValueError(
+ f'Layer "{lyr.name()}" is missing required field(s): '
+ f'{", ".join(missing)}. Select TEM "layers" outputs.'
+ )
+ out_dir = _output_path(input_layers, '_mean_resistivity')
+ if out_dir is None:
+ raise ValueError(
+ 'Could not determine an output location: the selected input '
+ 'layers are not backed by a GeoPackage on disk.'
+ )
+
+ layer_rows = [
+ {name: feat[name] for name in required}
+ for lyr in input_layers for feat in lyr.getFeatures()
+ ]
+ doi_rows = [
+ {name: feat[name] for name in ('X', 'Y', 'ZDOI')}
+ for lyr in opts['doi_layers'] for feat in lyr.getFeatures()
+ ]
+
+ progress_dialog = QProgressDialog(
+ 'Computing mean resistivity maps…', 'Cancel', 0, 100,
+ self.iface.mainWindow(),
+ )
+ progress_dialog.setWindowTitle('TEM Loader')
+ progress_dialog.setMinimumDuration(500)
+
+ def on_progress(done, total):
+ progress_dialog.setValue(int(100 * done / total))
+ return progress_dialog.wasCanceled()
+
+ grid, maps = mres.mean_resistivity_maps(
+ layer_rows,
+ opts['intervals'],
+ opts['cell_size'],
+ doi_points=doi_rows or None,
+ power=opts['power'],
+ radius=opts['radius'],
+ min_neighbors=opts['min_neighbors'],
+ min_coverage=opts['min_coverage'],
+ progress=on_progress,
+ )
+ progress_dialog.setValue(100)
+ if progress_dialog.wasCanceled():
+ raise ValueError('Mean resistivity maps cancelled.')
+ if not maps:
+ raise ValueError(
+ 'No mean resistivity maps produced. Check the elevation range '
+ 'and input layers.'
+ )
+
+ out_dir.mkdir(exist_ok=True)
+ crs_wkt = input_layers[0].crs().toWkt()
+ paths = []
+ for top, bottom, array, _ in maps:
+ path = out_dir / _mean_resistivity_file_name(top, bottom)
+ write_geotiff(path, grid, array, crs_wkt)
+ paths.append((path, _mean_resistivity_layer_name(top, bottom)))
+
+ project = QgsProject.instance()
+ group = project.layerTreeRoot().insertGroup(0, 'Mean resistivity')
+ qml = STYLES_DIR / 'mean_resistivity.qml'
+ for index, (path, name) in enumerate(paths):
+ layer = QgsRasterLayer(str(path), name, 'gdal')
+ if not layer.isValid():
+ raise ValueError(f'Failed to load mean resistivity map {path}')
+ layer.loadNamedStyle(str(qml))
+ project.addMapLayer(layer, False)
+ node = group.insertLayer(index, layer)
+ if index > 0:
+ node.setItemVisibilityChecked(False)
+
+ self.iface.messageBar().pushInfo(
+ 'TEM Loader',
+ f'{len(paths)} mean resistivity maps saved to {out_dir}',
+ )
+
def _load_xyz(
self,
filepath,
diff --git a/test/test_core.py b/test/test_core.py
@@ -1088,6 +1088,9 @@ class PluginTests(unittest.TestCase):
qgis_gui = types.ModuleType("qgis.gui")
qgis_gui.QgsMapLayerComboBox = FakeMapLayerComboBox
+ osgeo = types.ModuleType("osgeo")
+ osgeo.gdal = types.ModuleType("osgeo.gdal")
+
module_map = {
"qgis": types.ModuleType("qgis"),
"qgis.PyQt": types.ModuleType("qgis.PyQt"),
@@ -1095,6 +1098,8 @@ class PluginTests(unittest.TestCase):
"qgis.PyQt.QtWidgets": qtwidgets,
"qgis.core": qgis_core,
"qgis.gui": qgis_gui,
+ "osgeo": osgeo,
+ "osgeo.gdal": osgeo.gdal,
}
with patch.dict(sys.modules, module_map):
@@ -2189,6 +2194,57 @@ class PluginTests(unittest.TestCase):
self.assertEqual(active.attrib["value"], "true")
self.assertEqual(expr.attrib["value"], '"Color"')
+ def test_output_path_appends_suffix_to_geopackage_stem(self):
+ module, _, _ = self._import_plugin_module()
+ with TemporaryDirectory() as tmp:
+ gpkg = Path(tmp) / "survey.gpkg"
+ layers = [Mock(**{"source.return_value": "/missing/a.csv"}),
+ Mock(**{"source.return_value":
+ f"{gpkg}|layername=layers"})]
+
+ self.assertEqual(
+ module._output_path(layers, "_mean_resistivity"),
+ Path(tmp) / "survey_mean_resistivity",
+ )
+ self.assertIsNone(module._output_path(layers[:1], "_voxels.gpkg"))
+
+ def test_mean_resistivity_names(self):
+ module, _, _ = self._import_plugin_module()
+ self.assertEqual(module._mean_resistivity_file_name(0.0, -10.0),
+ "mean_resistivity_0_-10.tif")
+ self.assertEqual(module._mean_resistivity_layer_name(12.5, 2.5),
+ "Mean resistivity 12.5 to 2.5 m")
+
+ def test_mean_resistivity_requires_input_layers(self):
+ module, _, _ = self._import_plugin_module()
+ plugin = module.TEMLoaderPlugin(Mock())
+ with self.assertRaisesRegex(ValueError, "No input layers"):
+ plugin._run_mean_resistivity_maps({"input_layers": []})
+
+ def test_mean_resistivity_requires_layer_fields(self):
+ module, _, _ = self._import_plugin_module()
+ plugin = module.TEMLoaderPlugin(Mock())
+ layer = Mock()
+ layer.name.return_value = "points"
+ layer.fields.return_value = [module.QgsField("X", None),
+ module.QgsField("Y", None)]
+ with self.assertRaisesRegex(ValueError, "ZTop, ZBottom, Resistivity"):
+ plugin._run_mean_resistivity_maps({"input_layers": [layer]})
+
+ def test_run_mean_resistivity_maps_shows_errors(self):
+ module, _, message_box = self._import_plugin_module()
+ iface = Mock()
+ dialog = Mock()
+ dialog.options.return_value = {"input_layers": []}
+ module._MeanResistivityDialog = Mock(return_value=dialog)
+ module._exec_dialog = Mock(return_value=module.QDialog.Accepted)
+ plugin = module.TEMLoaderPlugin(iface)
+
+ plugin.run_mean_resistivity_maps()
+
+ self.assertEqual(len(message_box.warnings), 1)
+ self.assertIn("No input layers", message_box.warnings[0][2])
+
def test_integer_fields_includes_voxel_attrs(self):
module, _, _ = self._import_plugin_module()
for name in ('NumSamples', 'CellX', 'CellY', 'CellZ'):
