commit bc1454ab9f9ecc1060fe1377e0d921d92a40cadf
parent bf942038f18bed7953af99b6c2054dbebdb5c7d2
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Sat, 27 Jun 2026 07:55:40 +0200
fix(qgis): handle empty interpolation extent without crashing
Diffstat:
1 file changed, 452 insertions(+), 1 deletion(-)
diff --git a/tem_loader/tem_loader.py b/tem_loader/tem_loader.py
@@ -7,9 +7,15 @@ from qgis.PyQt.QtWidgets import (
QCheckBox,
QDialog,
QDialogButtonBox,
+ QDoubleSpinBox,
QFileDialog,
QFormLayout,
+ QLabel,
+ QListWidget,
+ QListWidgetItem,
QMessageBox,
+ QProgressDialog,
+ QPushButton,
QSpinBox,
QVBoxLayout,
)
@@ -30,18 +36,31 @@ from qgis.core import (
from qgis.gui import QgsMapLayerComboBox
from . import core
+from . import interpolate as interp
STYLES_DIR = Path(__file__).parent / 'styles'
GEOMETRY_FIELD = 'Geometry'
STRING_FIELDS = {'Line', 'StationNo', 'Color'}
-INTEGER_FIELDS = {'NumLayers', 'Layer', 'Opacity'}
+INTEGER_FIELDS = {'NumLayers', 'Layer', 'Opacity', 'NumSamples', 'CellX', 'CellY', 'CellZ'}
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."""
+ 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 None
+
+
def _field_type(field_name):
if field_name in STRING_FIELDS:
return QMetaType.Type.QString
@@ -296,19 +315,302 @@ def _exec_dialog(dialog):
return execute()
+class _InterpolationDialog(QDialog):
+ def __init__(self, parent=None):
+ super().__init__(parent)
+ self.setWindowTitle('3D IDW Interpolation')
+
+ project = QgsProject.instance()
+ all_layers = list(project.mapLayers().values()) if hasattr(project, 'mapLayers') else []
+
+ self._input_list = QListWidget()
+ self._input_list.setSelectionMode(QListWidget.SelectionMode.MultiSelection
+ if hasattr(QListWidget, 'SelectionMode')
+ else QListWidget.MultiSelection)
+ linestring_z_wkb = Qgis.WkbType.LineStringZ
+ for lyr in all_layers:
+ if hasattr(lyr, 'wkbType') and lyr.wkbType() == linestring_z_wkb:
+ item = QListWidgetItem(lyr.name())
+ item.setData(32, lyr)
+ self._input_list.addItem(item)
+
+ self._doi_list = QListWidget()
+ self._doi_list.setSelectionMode(QListWidget.SelectionMode.MultiSelection
+ if hasattr(QListWidget, 'SelectionMode')
+ else QListWidget.MultiSelection)
+ point_z_wkb = Qgis.WkbType.PointZ
+ for lyr in all_layers:
+ if hasattr(lyr, 'wkbType') and lyr.wkbType() == point_z_wkb:
+ item = QListWidgetItem(lyr.name())
+ item.setData(32, lyr)
+ self._doi_list.addItem(item)
+ if 'doi' in lyr.name().lower():
+ item.setSelected(True)
+
+ self._xmin = QDoubleSpinBox()
+ self._xmax = QDoubleSpinBox()
+ self._ymin = QDoubleSpinBox()
+ self._ymax = QDoubleSpinBox()
+ self._zmin = QDoubleSpinBox()
+ self._zmax = QDoubleSpinBox()
+ for sb in (self._xmin, self._xmax, self._ymin, self._ymax,
+ self._zmin, self._zmax):
+ sb.setRange(-1e9, 1e9)
+ sb.setDecimals(3)
+
+ self._dx = QDoubleSpinBox()
+ self._dy = QDoubleSpinBox()
+ self._dz = QDoubleSpinBox()
+ self._nx = QSpinBox()
+ self._ny = QSpinBox()
+ self._nz = QSpinBox()
+ for sb in (self._dx, self._dy, self._dz):
+ sb.setRange(0.001, 1e6)
+ sb.setDecimals(3)
+ sb.setValue(100.0)
+ for sb in (self._nx, self._ny, self._nz):
+ sb.setRange(1, 10000)
+ sb.setValue(10)
+
+ self._summary_label = QLabel()
+
+ self._power = QDoubleSpinBox()
+ self._power.setRange(0.1, 10.0)
+ self._power.setValue(2.0)
+ self._power.setDecimals(1)
+
+ self._log10_checkbox = QCheckBox('Interpolate in log10 space')
+ self._log10_checkbox.setChecked(True)
+
+ self._clip_doi_checkbox = QCheckBox('Clip voxels below per-sounding DOI')
+ self._clip_doi_checkbox.setChecked(True)
+
+ self._anisotropy = QDoubleSpinBox()
+ self._anisotropy.setRange(0.01, 100.0)
+ self._anisotropy.setValue(1.0)
+ self._anisotropy.setDecimals(2)
+
+ self._radius = QDoubleSpinBox()
+ self._radius.setRange(0.0, 1e9)
+ self._radius.setDecimals(2)
+ self._radius.setValue(0.0)
+ self._radius.setSpecialValueText('Auto')
+
+ self._min_neighbors = QSpinBox()
+ self._min_neighbors.setRange(1, 1000)
+ self._min_neighbors.setValue(1)
+
+ self._max_neighbors = QSpinBox()
+ self._max_neighbors.setRange(0, 10000)
+ self._max_neighbors.setValue(0)
+ self._max_neighbors.setSpecialValueText('Unlimited')
+
+ use_extent_btn = QPushButton('Use extent of selected input layers')
+ use_extent_btn.clicked.connect(self._fill_extent_from_inputs)
+
+ self._dx.valueChanged.connect(self._on_cell_size_changed)
+ self._dy.valueChanged.connect(self._on_cell_size_changed)
+ self._dz.valueChanged.connect(self._on_cell_size_changed)
+ self._nx.valueChanged.connect(self._on_cell_count_changed)
+ self._ny.valueChanged.connect(self._on_cell_count_changed)
+ self._nz.valueChanged.connect(self._on_cell_count_changed)
+ self._xmin.valueChanged.connect(self._update_summary)
+ self._xmax.valueChanged.connect(self._update_summary)
+ self._ymin.valueChanged.connect(self._update_summary)
+ self._ymax.valueChanged.connect(self._update_summary)
+ self._zmin.valueChanged.connect(self._update_summary)
+ self._zmax.valueChanged.connect(self._update_summary)
+
+ self._updating = False
+
+ buttons = QDialogButtonBox(_dialog_buttons())
+ buttons.accepted.connect(self.accept)
+ buttons.rejected.connect(self.reject)
+
+ extent_form = QFormLayout()
+ extent_form.addRow('X min', self._xmin)
+ extent_form.addRow('X max', self._xmax)
+ extent_form.addRow('Y min', self._ymin)
+ extent_form.addRow('Y max', self._ymax)
+ extent_form.addRow('Z min', self._zmin)
+ extent_form.addRow('Z max', self._zmax)
+
+ size_form = QFormLayout()
+ size_form.addRow('Cell size X (m)', self._dx)
+ size_form.addRow('Cell size Y (m)', self._dy)
+ size_form.addRow('Cell size Z (m)', self._dz)
+ size_form.addRow('Cell count X', self._nx)
+ size_form.addRow('Cell count Y', self._ny)
+ size_form.addRow('Cell count Z', self._nz)
+ size_form.addRow('', self._summary_label)
+
+ param_form = QFormLayout()
+ param_form.addRow('IDW power', self._power)
+ param_form.addRow('Vertical anisotropy', self._anisotropy)
+ param_form.addRow('Search radius (m)', self._radius)
+ param_form.addRow('Min neighbors', self._min_neighbors)
+ param_form.addRow('Max neighbors', self._max_neighbors)
+
+ 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_extent_btn)
+ layout.addLayout(extent_form)
+ layout.addLayout(size_form)
+ layout.addWidget(self._log10_checkbox)
+ layout.addWidget(self._clip_doi_checkbox)
+ layout.addLayout(param_form)
+ layout.addWidget(buttons)
+ self.setLayout(layout)
+
+ self._update_summary()
+
+ def _fill_extent_from_inputs(self):
+ layers = self._selected_input_layers()
+ if not layers:
+ return
+ xs, ys, zs = [], [], []
+ for lyr in layers:
+ ext = lyr.extent()
+ xs += [ext.xMinimum(), ext.xMaximum()]
+ ys += [ext.yMinimum(), ext.yMaximum()]
+ for feat in lyr.getFeatures():
+ geom = feat.geometry()
+ if geom is None:
+ continue
+ for vertex in geom.vertices():
+ z = vertex.z()
+ if z == z: # exclude NaN
+ zs.append(z)
+ self._updating = True
+ self._xmin.setValue(min(xs))
+ self._xmax.setValue(max(xs))
+ self._ymin.setValue(min(ys))
+ self._ymax.setValue(max(ys))
+ if zs:
+ self._zmin.setValue(min(zs))
+ self._zmax.setValue(max(zs))
+ self._updating = False
+ self._update_summary()
+
+ def _selected_input_layers(self):
+ layers = []
+ for i in range(self._input_list.count()):
+ item = self._input_list.item(i)
+ if item.isSelected():
+ lyr = item.data(32)
+ if lyr is not None:
+ layers.append(lyr)
+ return layers
+
+ def _selected_doi_layers(self):
+ layers = []
+ for i in range(self._doi_list.count()):
+ item = self._doi_list.item(i)
+ if item.isSelected():
+ lyr = item.data(32)
+ if lyr is not None:
+ layers.append(lyr)
+ return layers
+
+ def _on_cell_size_changed(self):
+ if self._updating:
+ return
+ self._updating = True
+ try:
+ xext = self._xmax.value() - self._xmin.value()
+ yext = self._ymax.value() - self._ymin.value()
+ zext = self._zmax.value() - self._zmin.value()
+ if xext > 0:
+ self._nx.setValue(interp.axis_count(xext, self._dx.value()))
+ if yext > 0:
+ self._ny.setValue(interp.axis_count(yext, self._dy.value()))
+ if zext > 0:
+ self._nz.setValue(interp.axis_count(zext, self._dz.value()))
+ except (ValueError, ZeroDivisionError):
+ pass
+ finally:
+ self._updating = False
+ self._update_summary()
+
+ def _on_cell_count_changed(self):
+ if self._updating:
+ return
+ self._updating = True
+ try:
+ xext = self._xmax.value() - self._xmin.value()
+ yext = self._ymax.value() - self._ymin.value()
+ zext = self._zmax.value() - self._zmin.value()
+ if xext > 0:
+ self._dx.setValue(interp.axis_cell_size(xext, self._nx.value()))
+ if yext > 0:
+ self._dy.setValue(interp.axis_cell_size(yext, self._ny.value()))
+ if zext > 0:
+ self._dz.setValue(interp.axis_cell_size(zext, self._nz.value()))
+ except (ValueError, ZeroDivisionError):
+ pass
+ finally:
+ self._updating = False
+ self._update_summary()
+
+ def _update_summary(self):
+ nx, ny, nz = self._nx.value(), self._ny.value(), self._nz.value()
+ total = nx * ny * nz
+ self._summary_label.setText(
+ f'{nx} × {ny} × {nz} = {total:,} cells'
+ )
+
+ def options(self):
+ bounds = (
+ self._xmin.value(), self._xmax.value(),
+ self._ymin.value(), self._ymax.value(),
+ self._zmin.value(), self._zmax.value(),
+ )
+ grid = interp.GridSpec.from_cell_count(
+ bounds,
+ self._nx.value(),
+ self._ny.value(),
+ self._nz.value(),
+ )
+ radius = self._radius.value()
+ max_neighbors = self._max_neighbors.value()
+ return {
+ 'input_layers': self._selected_input_layers(),
+ 'doi_layers': self._selected_doi_layers(),
+ 'grid': grid,
+ 'power': self._power.value(),
+ 'log10': self._log10_checkbox.isChecked(),
+ 'clip_below_doi': self._clip_doi_checkbox.isChecked(),
+ 'anisotropy_z': self._anisotropy.value(),
+ 'radius': radius if radius > 0 else None,
+ 'min_neighbors': self._min_neighbors.value(),
+ 'max_neighbors': max_neighbors if max_neighbors > 0 else None,
+ }
+
+
class TEMLoaderPlugin:
def __init__(self, iface):
self.iface = iface
self._action = None
+ self._interp_action = None
def initGui(self):
self._action = QAction('Load TEM XYZ files\u2026', self.iface.mainWindow())
self._action.triggered.connect(self.run)
self.iface.addPluginToMenu('TEM Loader', self._action)
+ self._interp_action = QAction('3D interpolation\u2026', self.iface.mainWindow())
+ self._interp_action.triggered.connect(self.run_interpolation)
+ self.iface.addPluginToMenu('TEM Loader', self._interp_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
def run(self):
paths, _ = QFileDialog.getOpenFileNames(
@@ -354,6 +656,155 @@ class TEMLoaderPlugin:
crs.createFromString('EPSG:4326')
return crs
+ def run_interpolation(self):
+ dialog = _InterpolationDialog(self.iface.mainWindow())
+ if _exec_dialog(dialog) != _dialog_code('Accepted'):
+ return
+
+ try:
+ opts = dialog.options()
+ except ValueError as exc:
+ QMessageBox.warning(
+ self.iface.mainWindow(),
+ 'TEM Loader',
+ f'Invalid grid: {exc}. Set a non-empty extent (e.g. use '
+ '"Use extent of selected input layers").',
+ )
+ return
+ input_layers = opts['input_layers']
+ if not input_layers:
+ QMessageBox.warning(
+ self.iface.mainWindow(),
+ 'TEM Loader',
+ 'No input layers selected.',
+ )
+ return
+
+ gpkg_path = _voxel_output_path(input_layers)
+ if gpkg_path is None:
+ QMessageBox.warning(
+ self.iface.mainWindow(),
+ 'TEM Loader',
+ 'Could not determine an output location: the selected input '
+ 'layers are not backed by a GeoPackage on disk.',
+ )
+ return
+
+ try:
+ self._run_interpolation(gpkg_path, opts, input_layers)
+ except ValueError as exc:
+ QMessageBox.warning(
+ self.iface.mainWindow(),
+ 'TEM Loader',
+ str(exc),
+ )
+
+ def _run_interpolation(self, gpkg_path, opts, input_layers):
+ required = ('X', 'Y', 'ZMid', '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.'
+ )
+
+ layer_rows = []
+ for lyr in input_layers:
+ for feat in lyr.getFeatures():
+ attrs = {field.name(): feat[field.name()] for field in lyr.fields()}
+ geom = feat.geometry()
+ layer_rows.append({
+ 'X': attrs.get('X'),
+ 'Y': attrs.get('Y'),
+ 'ZMid': attrs.get('ZMid'),
+ 'Resistivity': attrs.get('Resistivity'),
+ 'Geometry': geom.asWkt() if geom else '',
+ })
+
+ doi_rows = []
+ for lyr in opts['doi_layers']:
+ for feat in lyr.getFeatures():
+ attrs = {field.name(): feat[field.name()] for field in lyr.fields()}
+ doi_rows.append({
+ 'X': attrs.get('X'),
+ 'Y': attrs.get('Y'),
+ 'ZDOI': attrs.get('ZDOI'),
+ })
+
+ progress_dialog = QProgressDialog(
+ 'Interpolating voxels…', 'Cancel', 0, 100,
+ self.iface.mainWindow(),
+ )
+ progress_dialog.setWindowTitle('TEM Loader')
+ progress_dialog.setMinimumDuration(500)
+
+ last_pct = [-1]
+
+ def on_progress(done, total):
+ pct = int(100 * done / total) if total else 100
+ # Only touch the widget on whole-percent changes; the callback fires
+ # once per cell (up to several million times).
+ if pct != last_pct[0]:
+ last_pct[0] = pct
+ progress_dialog.setValue(pct)
+ return progress_dialog.wasCanceled()
+
+ grid = opts['grid']
+ voxels = interp.interpolate_voxels(
+ layer_rows,
+ grid,
+ doi_points=doi_rows if doi_rows else None,
+ clip_below_doi=opts['clip_below_doi'],
+ power=opts['power'],
+ radius=opts['radius'],
+ min_neighbors=opts['min_neighbors'],
+ max_neighbors=opts['max_neighbors'],
+ anisotropy_z=opts['anisotropy_z'],
+ log10=opts['log10'],
+ progress=on_progress,
+ )
+ progress_dialog.setValue(100)
+
+ if progress_dialog.wasCanceled():
+ raise ValueError('Interpolation cancelled.')
+
+ if not voxels:
+ raise ValueError('Interpolation produced no voxels. Check extent and input layers.')
+
+ project = QgsProject.instance()
+ crs = input_layers[0].crs()
+ transform_context = project.transformContext()
+
+ _write_geopackage_layer(
+ voxels,
+ gpkg_path,
+ 'voxels',
+ Qgis.WkbType.PointZ,
+ crs,
+ transform_context,
+ QgsVectorFileWriter.CreateOrOverwriteFile,
+ )
+
+ uri = _build_geopackage_uri(gpkg_path, 'voxels')
+ layer = QgsVectorLayer(uri, 'voxels', 'ogr')
+ if not layer.isValid():
+ raise ValueError(f'Failed to load voxel layer from {gpkg_path}')
+
+ qml = STYLES_DIR / 'voxels.qml'
+ if qml.exists():
+ layer.loadNamedStyle(str(qml))
+
+ project.addMapLayer(layer, False)
+ root = project.layerTreeRoot()
+ group = root.insertGroup(0, 'Interpolation')
+ group.insertLayer(0, layer)
+
+ self.iface.messageBar().pushInfo(
+ 'TEM Loader', f'Voxel layer saved to {gpkg_path}'
+ )
+
def _load_xyz(
self,
filepath,
