qgis-tem-loader

tem_loader.py (41054B)

---

 import math
 from pathlib import Path
 
 from qgis.PyQt.QtCore import QMetaType
 from qgis.PyQt.QtWidgets import (
     QAction,
     QCheckBox,
     QDialog,
     QDialogButtonBox,
     QDoubleSpinBox,
     QFileDialog,
     QFormLayout,
     QLabel,
     QListWidget,
     QListWidgetItem,
     QMessageBox,
     QProgressDialog,
     QPushButton,
     QSpinBox,
     QVBoxLayout,
 )
 from qgis.core import (
     Qgis,
     QgsCoordinateReferenceSystem,
     QgsCoordinateTransform,
     QgsCsException,
     QgsFeature,
     QgsField,
     QgsFields,
     QgsGeometry,
     QgsPointXY,
     QgsProject,
     QgsRasterLayer,
     QgsVectorFileWriter,
     QgsVectorLayer,
 )
 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'
 GEOMETRY_FIELD = 'Geometry'
 STRING_FIELDS = {'Line', 'StationNo', 'Color'}
 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 _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}{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
     if field_name in INTEGER_FIELDS:
         return QMetaType.Type.Int
     return QMetaType.Type.Double
 
 
 def _build_fields(rows):
     fields = QgsFields()
     if not rows:
         return fields
     for field_name in rows[0]:
         if field_name == GEOMETRY_FIELD:
             continue
         fields.append(QgsField(field_name, _field_type(field_name)))
     return fields
 
 
 def _rows_to_features(rows, fields):
     field_names = [field.name() for field in fields]
     features = []
     for row in rows:
         feature = QgsFeature(fields)
         feature.setGeometry(QgsGeometry.fromWkt(row[GEOMETRY_FIELD]))
         feature.setAttributes([row.get(name) for name in field_names])
         features.append(feature)
     return features
 
 
 def _writer_no_error():
     writer_error = getattr(QgsVectorFileWriter, 'WriterError', QgsVectorFileWriter)
     return writer_error.NoError
 
 
 def _dialog_button(name):
     standard_button = getattr(
         QDialogButtonBox,
         'StandardButton',
         QDialogButtonBox,
     )
     return getattr(standard_button, name)
 
 
 def _dialog_buttons():
     return _dialog_button('Ok') | _dialog_button('Cancel')
 
 
 def _dialog_code(name):
     dialog_code = getattr(QDialog, 'DialogCode', QDialog)
     return getattr(dialog_code, name)
 
 
 def _layer_filter(name):
     layer_filter = getattr(Qgis, 'LayerFilter', Qgis)
     return getattr(layer_filter, name)
 
 
 def _warn_dem_sample_failed(raster_layer, x, y):
     print(
         'TEM Loader warning: '
         f'point ({x}, {y}) is outside DEM raster {raster_layer.name()}'
     )
 
 
 def _sample_dem_elevation(raster_layer, source_crs, project, x, y):
     transform = QgsCoordinateTransform(source_crs, raster_layer.crs(), project)
     try:
         raster_point = transform.transform(QgsPointXY(x, y))
     except QgsCsException as exc:
         print(
             'TEM Loader warning: '
             f'could not transform point ({x}, {y}) to DEM raster '
             f'{raster_layer.name()}: {exc}'
         )
         return None
 
     value, ok = raster_layer.dataProvider().sample(raster_point, 1)
     try:
         elevation = float(value)
     except (TypeError, ValueError):
         elevation = math.nan
     if not ok or math.isnan(elevation):
         _warn_dem_sample_failed(raster_layer, x, y)
         return None
     return elevation
 
 
 def _point_wkt(x, y, z):
     return f'POINT Z ({x} {y} {z})'
 
 
 def _line_wkt(x, y, z_top, z_bottom):
     return f'LINESTRING Z ({x} {y} {z_top}, {x} {y} {z_bottom})'
 
 
 def _adjust_rows_to_dem(points, doi_points, layers, raster_layer, source_crs,
                         project):
     if raster_layer is None:
         return
 
     dem_z_by_xy = {}
     for point in points:
         x = point['X']
         y = point['Y']
         dem_z = _sample_dem_elevation(raster_layer, source_crs, project, x, y)
         if dem_z is None:
             continue
         dem_z_by_xy[(x, y)] = dem_z
         point['Z'] = dem_z
         point['Geometry'] = _point_wkt(x, y, dem_z)
 
     for doi_point in doi_points:
         x = doi_point['X']
         y = doi_point['Y']
         dem_z = dem_z_by_xy.get((x, y))
         if dem_z is None:
             continue
         z_doi = dem_z - doi_point['DOI']
         doi_point['Z'] = z_doi
         doi_point['ZDOI'] = z_doi
         doi_point['Geometry'] = _point_wkt(x, y, z_doi)
 
     for layer in layers:
         x = layer['X']
         y = layer['Y']
         dem_z = dem_z_by_xy.get((x, y))
         if dem_z is None:
             continue
         z_top = dem_z - layer['DepthTop']
         z_bottom = dem_z - layer['DepthBottom']
         z_mid = (z_top + z_bottom) / 2
         layer['Z'] = dem_z
         layer['ZTop'] = z_top
         layer['ZMid'] = z_mid
         layer['ZBottom'] = z_bottom
         layer['Geometry'] = _line_wkt(x, y, z_top, z_bottom)
 
 
 def _write_geopackage_layer(rows, gpkg_path, layer_name, wkb_type, crs,
                             transform_context, action):
     fields = _build_fields(rows)
     options = QgsVectorFileWriter.SaveVectorOptions()
     options.driverName = 'GPKG'
     options.fileEncoding = 'UTF-8'
     options.layerName = layer_name
     options.actionOnExistingFile = action
 
     writer = QgsVectorFileWriter.create(
         str(gpkg_path),
         fields,
         wkb_type,
         crs,
         transform_context,
         options,
     )
     if writer is None:
         raise ValueError(f'failed to create GeoPackage layer {layer_name}')
 
     try:
         if writer.hasError() != _writer_no_error():
             message = writer.errorMessage()
             details = f': {message}' if message else ''
             raise ValueError(
                 f'failed to create GeoPackage layer {layer_name}{details}'
             )
         features = _rows_to_features(rows, fields)
         if features and not writer.addFeatures(features):
             raise ValueError(f'failed to write GeoPackage layer {layer_name}')
     finally:
         del writer
 
 
 class _ImportOptionsDialog(QDialog):
     def __init__(self, parent=None):
         super().__init__(parent)
         self.setWindowTitle('TEM Loader Options')
 
         self._clip_checkbox = QCheckBox(
             'Hide layers below depth of interest (DOI, 2D and 3D views)'
         )
         self._clip_checkbox.setChecked(True)
 
         self._mask_checkbox = QCheckBox(
             'Make layers below DOI partially transparent (2D view only)'
         )
         self._mask_checkbox.setChecked(True)
 
         self._opacity_spinbox = QSpinBox()
         self._opacity_spinbox.setRange(0, 100)
         self._opacity_spinbox.setSuffix('%')
         self._opacity_spinbox.setValue(core.BELOW_DOI_OPACITY)
 
         self._clip_checkbox.toggled.connect(self._update_below_doi_controls)
         self._mask_checkbox.toggled.connect(self._update_below_doi_controls)
         self._update_below_doi_controls()
 
         self._dem_checkbox = QCheckBox(
             'Adjust vertical position to digital elevation model'
         )
         self._dem_checkbox.setChecked(False)
 
         self._dem_raster_combo = QgsMapLayerComboBox()
         self._dem_raster_combo.setProject(QgsProject.instance())
         self._dem_raster_combo.setFilters(_layer_filter('RasterLayer'))
         self._dem_raster_combo.setAllowEmptyLayer(True, 'No elevation raster')
         self._dem_raster_combo.setEnabled(False)
         self._dem_checkbox.toggled.connect(self._dem_raster_combo.setEnabled)
 
         opacity_form = QFormLayout()
         opacity_form.addRow('Opacity', self._opacity_spinbox)
 
         dem_form = QFormLayout()
         dem_form.addRow('Elevation raster', self._dem_raster_combo)
 
         buttons = QDialogButtonBox(_dialog_buttons())
         buttons.accepted.connect(self.accept)
         buttons.rejected.connect(self.reject)
 
         layout = QVBoxLayout()
         layout.addWidget(self._clip_checkbox)
         layout.addWidget(self._mask_checkbox)
         layout.addLayout(opacity_form)
         layout.addWidget(self._dem_checkbox)
         layout.addLayout(dem_form)
         layout.addWidget(buttons)
         self.setLayout(layout)
 
     def _update_below_doi_controls(self, *_args):
         clipping = self._clip_checkbox.isChecked()
         self._mask_checkbox.setEnabled(not clipping)
         self._opacity_spinbox.setEnabled(
             not clipping and self._mask_checkbox.isChecked()
         )
 
     def options(self):
         elevation_raster_layer = None
         if self._dem_checkbox.isChecked():
             elevation_raster_layer = self._dem_raster_combo.currentLayer()
         return {
             'clip_below_doi': self._clip_checkbox.isChecked(),
             'mask_below_doi': self._mask_checkbox.isChecked(),
             'below_doi_opacity': self._opacity_spinbox.value(),
             'elevation_raster_layer': elevation_raster_layer,
         }
 
 
 def _exec_dialog(dialog):
     execute = getattr(dialog, 'exec', None)
     if execute is None:
         execute = dialog.exec_
     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()
         for sb in (self._dx, self._dy, self._dz):
             sb.setRange(0.001, 1e6)
             sb.setDecimals(3)
             sb.setValue(100.0)
         self._dz.setValue(10.0)
         # Cell counts are derived from extent and cell size, not editable.
         self._nx = QLabel('—')
         self._ny = QLabel('—')
         self._nz = QLabel('—')
         self._count_x = self._count_y = self._count_z = None
 
         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._radius = QDoubleSpinBox()
         self._radius.setRange(0.0, 1e9)
         self._radius.setDecimals(2)
         self._radius.setValue(0.0)
         self._radius.setSpecialValueText('Auto')
 
         self._radius_z = QDoubleSpinBox()
         self._radius_z.setRange(0.0, 1e9)
         self._radius_z.setDecimals(2)
         self._radius_z.setValue(0.0)
         self._radius_z.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._updating = False
 
         for sb in (self._dx, self._dy, self._dz,
                    self._xmin, self._xmax, self._ymin, self._ymax,
                    self._zmin, self._zmax):
             sb.valueChanged.connect(self._recompute_counts)
 
         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('Horizontal search radius (m)', self._radius)
         param_form.addRow('Vertical search radius (m)', self._radius_z)
         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._recompute_counts()
 
     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._recompute_counts()
 
     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
 
     @staticmethod
     def _axis_count(extent, cell_size):
         try:
             return interp.axis_count(extent, cell_size)
         except (ValueError, ZeroDivisionError):
             return None
 
     def _recompute_counts(self):
         if self._updating:
             return
         self._count_x = self._axis_count(
             self._xmax.value() - self._xmin.value(), self._dx.value())
         self._count_y = self._axis_count(
             self._ymax.value() - self._ymin.value(), self._dy.value())
         self._count_z = self._axis_count(
             self._zmax.value() - self._zmin.value(), self._dz.value())
         self._nx.setText('—' if self._count_x is None else str(self._count_x))
         self._ny.setText('—' if self._count_y is None else str(self._count_y))
         self._nz.setText('—' if self._count_z is None else str(self._count_z))
         self._update_summary()
 
     def _update_summary(self):
         if None in (self._count_x, self._count_y, self._count_z):
             self._summary_label.setText(
                 'Set a non-empty extent to compute the grid.'
             )
             return
         nx, ny, nz = self._count_x, self._count_y, self._count_z
         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_size(
             bounds,
             self._dx.value(),
             self._dy.value(),
             self._dz.value(),
         )
         radius = self._radius.value()
         radius_z = self._radius_z.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(),
             'radius': radius if radius > 0 else None,
             'radius_z': radius_z if radius_z > 0 else None,
             'min_neighbors': self._min_neighbors.value(),
             'max_neighbors': max_neighbors if max_neighbors > 0 else None,
         }
 
 
 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())
         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)
 
         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(
             self.iface.mainWindow(),
             'Select TEM XYZ files',
             '',
             'XYZ files (*.xyz);;All files (*)',
         )
         if not paths:
             return
 
         dialog = _ImportOptionsDialog(self.iface.mainWindow())
         if _exec_dialog(dialog) != _dialog_code('Accepted'):
             return
 
         options = dialog.options()
         failed = []
         for path in paths:
             filepath = Path(path)
             try:
                 self._load_xyz(filepath, **options)
             except Exception as exc:
                 failed.append(f'{filepath.name}: {exc}')
         if failed:
             QMessageBox.warning(
                 self.iface.mainWindow(),
                 'TEM Loader',
                 '\n'.join(failed),
             )
 
     def _resolve_crs(self, filepath, project):
         crs = None
         source_epsg = core.detect_source_epsg(filepath)
         if source_epsg:
             candidate = QgsCoordinateReferenceSystem()
             if candidate.createFromString(source_epsg) and candidate.isValid():
                 crs = candidate
 
         if crs is None:
             crs = project.crs()
         if not crs.isValid():
             crs = QgsCoordinateReferenceSystem()
             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 = _output_path(input_layers, '_voxels.gpkg')
         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'],
             radius_z=opts['radius_z'],
             min_neighbors=opts['min_neighbors'],
             max_neighbors=opts['max_neighbors'],
             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 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,
         clip_below_doi=True,
         mask_below_doi=True,
         below_doi_opacity=core.BELOW_DOI_OPACITY,
         elevation_raster_layer=None,
     ):
         points, doi_points, layers = core.process_xyz(
             filepath,
             mask_below_doi=mask_below_doi,
             below_doi_opacity=below_doi_opacity,
             clip_below_doi=clip_below_doi,
         )
 
         project = QgsProject.instance()
         crs = self._resolve_crs(filepath, project)
         _adjust_rows_to_dem(
             points,
             doi_points,
             layers,
             elevation_raster_layer,
             crs,
             project,
         )
         gpkg_path = filepath.with_suffix('.gpkg')
         transform_context = project.transformContext()
 
         output_layers = [
             ('layers', layers, Qgis.WkbType.LineStringZ),
         ]
         if doi_points:
             output_layers.append(('doi', doi_points, Qgis.WkbType.PointZ))
         output_layers.append(('points', points, Qgis.WkbType.PointZ))
 
         written_layer_names = []
         first_layer = True
         for name, rows, wkb_type in output_layers:
             if not rows:
                 continue
             action = (
                 QgsVectorFileWriter.CreateOrOverwriteFile
                 if first_layer
                 else QgsVectorFileWriter.CreateOrOverwriteLayer
             )
             _write_geopackage_layer(
                 rows,
                 gpkg_path,
                 name,
                 wkb_type,
                 crs,
                 transform_context,
                 action,
             )
             written_layer_names.append(name)
             first_layer = False
 
         loaded_layers = {}
         for name in written_layer_names:
             uri = _build_geopackage_uri(gpkg_path, name)
             layer = QgsVectorLayer(uri, name, 'ogr')
             if not layer.isValid():
                 continue
 
             qml = STYLES_DIR / f'{name}.qml'
             if qml.exists():
                 layer.loadNamedStyle(str(qml))
 
             project.addMapLayer(layer, False)
             loaded_layers[name] = layer
 
         if not loaded_layers:
             raise ValueError('failed to load any layers')
 
         group_name = filepath.stem
         root = project.layerTreeRoot()
         group = root.insertGroup(0, group_name)
         insert_index = 0
         for name in ('points', 'doi', 'layers'):
             layer = loaded_layers.get(name)
             if layer is None:
                 continue
             group.insertLayer(insert_index, layer)
             insert_index += 1
 
         failed = [name for name in written_layer_names if name not in loaded_layers]
         if failed:
             QMessageBox.warning(
                 self.iface.mainWindow(),
                 'TEM Loader',
                 f'{filepath.name}: failed to load layers: {", ".join(failed)}',
             )