commit 25990b06d886078b7ad2350ff8098316bf49c222
parent c554d2c95bfd8d2b3902f2baac65955e674b4fb4
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Sat, 27 Jun 2026 08:50:33 +0200
feat(interpolate): replace vertical anisotropy with anisotropic search ellipsoid
Horizontal and vertical search radii now define a search ellipsoid that governs both sample culling and IDW weighting, preserving the vertical layering of the soundings.
Diffstat:
5 files changed, 121 insertions(+), 45 deletions(-)
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,7 +3,7 @@
## [0.3.0] - 2026-06-27
### Added
-- "3D interpolation…" menu item: IDW interpolation of loaded TEM layer features into a voxel PointZ layer written to a GeoPackage. The interpolation dialog lets the user select input LineStringZ layers and optional DOI PointZ layers, define a 3D grid extent with cell sizes and a live derived cell-count preview, and control IDW parameters (power, log10 mode, vertical anisotropy, per-sounding DOI clipping). Output voxels are styled with the resistivity color ramp and rendered as cubes in the 3D Map View.
+- "3D interpolation…" menu item: IDW interpolation of loaded TEM layer features into a voxel PointZ layer written to a GeoPackage. The interpolation dialog lets the user select input LineStringZ layers and optional DOI PointZ layers, define a 3D grid extent with cell sizes and a live derived cell-count preview, and control IDW parameters (power, log10 mode, horizontal and vertical search radii, per-sounding DOI clipping). The horizontal and vertical search radii define an anisotropic search ellipsoid (after GeoScene3D) that governs both which samples a voxel uses and how they are weighted, preserving the vertical layering of the soundings. Output voxels are styled with the resistivity color ramp and rendered as cubes in the 3D Map View.
## [0.2.2] - 2026-05-28
diff --git a/README.md b/README.md
@@ -43,7 +43,7 @@ After loading one or more TEM datasets, the plugin can interpolate their resisti
2. Go to **Plugins > TEM Loader > 3D interpolation…**.
3. Select the input `layers` from the list (multi-select). Select any `doi` layers to enable per-sounding DOI clipping.
4. Set the 3D grid extent (X, Y, Z) and cell size. The cell counts are derived from the extent and cell size and shown as a live `nx × ny × nz = N cells` preview that updates as you edit. Use **Use extent of selected input layers** to fill X/Y from the layer bounding box.
-5. Adjust IDW parameters (power, log10 mode, vertical anisotropy, search radius, and min/max neighbors) and optionally disable DOI clipping. Leave the search radius at **Auto** to derive it from the median sounding spacing, or **Max neighbors** at **Unlimited** to weight every sample in range.
+5. Adjust IDW parameters (power, log10 mode, horizontal and vertical search radii, and min/max neighbors) and optionally disable DOI clipping. The horizontal and vertical search radii define an anisotropic search ellipsoid: a wide horizontal reach bridges sparse soundings while a tight vertical radius preserves the layering of the resistivity models. Leave **Horizontal search radius** at **Auto** to derive it from the median sounding spacing, **Vertical search radius** at **Auto** to derive it from the median layer thickness, or **Max neighbors** at **Unlimited** to weight every sample in range.
6. Click **OK** and choose an output `.gpkg` path.
7. The resulting `voxels` PointZ layer is added under an **Interpolation** group, styled with the resistivity color ramp, and visible as colored cubes in the 3D Map View.
diff --git a/tem_loader/interpolate.py b/tem_loader/interpolate.py
@@ -157,6 +157,29 @@ def _default_radius(samples):
return float(np.median(nearest) * 3.0)
+def _default_radius_z(samples):
+ """Vertical search window derived from the typical layer thickness: 1.5x the
+ median gap between consecutive sample depths within a single sounding
+ (roughly a couple of layers). Operating per sounding keeps the estimate
+ independent of how the soundings are laid out horizontally. Returns ``None``
+ when no sounding has two or more distinct depths."""
+ by_xy = {}
+ for s in samples:
+ by_xy.setdefault((s.x, s.y), []).append(s.z)
+ gaps = []
+ for zs in by_xy.values():
+ if len(zs) < 2:
+ continue
+ zs = sorted(zs)
+ for a, b in zip(zs, zs[1:]):
+ gap = b - a
+ if gap > 0:
+ gaps.append(gap)
+ if not gaps:
+ return None
+ return float(np.median(gaps) * 1.5)
+
+
def _bucket_index(arrays, radius):
"""Map integer (bx, by) cells of side ``radius`` to the sample indices that
fall inside them, so a query only weighs samples in its 3x3 neighborhood."""
@@ -187,14 +210,22 @@ def _candidate_indices(buckets, bucket_size, x, y):
return np.concatenate(parts)
-def _interpolate_at(arrays, idx, qx, qy, qz, *, power, radius, max_neighbors,
- min_neighbors, anisotropy_z, log10, gate_doi):
+def _interpolate_at(arrays, idx, qx, qy, qz, *, power, radius, radius_z,
+ max_neighbors, min_neighbors, log10, gate_doi):
"""Inverse-distance interpolation at a single query point.
+ ``radius`` (horizontal) and ``radius_z`` (vertical) define a search
+ ellipsoid that governs both which samples are used and how they are
+ weighted. The ellipsoid is a sphere in vertically-scaled coordinates with
+ exaggeration ``k = radius / radius_z``, so the scaled squared distance
+ ``d2 = horizontal2 + (dz * k)**2`` drives the cull, neighbor selection, and
+ weighting consistently. When either radius is ``None`` weighting is isotropic
+ (``k = 1``) and only the axes whose radius is set are culled.
+
Returns ``(value, num_valid)``. ``value`` is ``None`` when fewer than
``min_neighbors`` usable samples remain; ``num_valid`` is always the count of
- samples that passed the DOI gate and radius cull. Shared by :func:`idw` and
- :func:`interpolate_voxels` so both honor identical weighting rules."""
+ samples that passed the DOI gate and ellipsoid cull. Shared by :func:`idw`
+ and :func:`interpolate_voxels` so both honor identical rules."""
if idx is None:
sx, sy, sz, sres, sdoi = (arrays.x, arrays.y, arrays.z, arrays.res,
arrays.doi)
@@ -204,26 +235,33 @@ def _interpolate_at(arrays, idx, qx, qy, qz, *, power, radius, max_neighbors,
dx = sx - qx
dy = sy - qy
+ dz = sz - qz
horizontal2 = dx * dx + dy * dy
+ if radius is not None and radius_z is not None:
+ k = radius / radius_z
+ else:
+ k = 1.0
+ vdz = dz * k
+ d2 = horizontal2 + vdz * vdz
+
keep = np.ones(horizontal2.shape, dtype=bool)
if gate_doi:
# NaN DOI (no DOI for that sounding) never gates; only drop samples
# whose DOI elevation sits above the query depth.
keep &= ~((~np.isnan(sdoi)) & (qz < sdoi))
- if radius is not None:
- keep &= horizontal2 <= radius * radius
+ if radius is not None and radius_z is not None:
+ # Ellipsoid: d2 <= radius**2 is horizontal2/radius**2 + dz**2/radius_z**2 <= 1.
+ keep &= d2 <= radius * radius
+ else:
+ if radius is not None:
+ keep &= horizontal2 <= radius * radius
+ if radius_z is not None:
+ keep &= dz * dz <= radius_z * radius_z
- horizontal2 = horizontal2[keep]
- sz = sz[keep]
+ d2 = d2[keep]
res = sres[keep]
- if anisotropy_z != 0:
- vdz = (sz - qz) / anisotropy_z
- d2 = horizontal2 + vdz * vdz
- else:
- d2 = horizontal2
-
if max_neighbors is not None and d2.shape[0] > max_neighbors:
nearest = np.argpartition(d2, max_neighbors)[:max_neighbors]
d2 = d2[nearest]
@@ -262,8 +300,8 @@ def _interpolate_at(arrays, idx, qx, qy, qz, *, power, radius, max_neighbors,
return (10.0 ** avg if log10 else avg), num_valid
-def idw(samples, queries, power=2.0, radius=None, max_neighbors=None,
- min_neighbors=1, anisotropy_z=1.0, log10=True, gate_doi=True):
+def idw(samples, queries, power=2.0, radius=None, radius_z=None,
+ max_neighbors=None, min_neighbors=1, log10=True, gate_doi=True):
if radius is None:
radius = _default_radius(samples)
@@ -275,8 +313,8 @@ def idw(samples, queries, power=2.0, radius=None, max_neighbors=None,
idx = _candidate_indices(buckets, bucket_size, qx, qy)
value, _ = _interpolate_at(
arrays, idx, qx, qy, qz,
- power=power, radius=radius, max_neighbors=max_neighbors,
- min_neighbors=min_neighbors, anisotropy_z=anisotropy_z,
+ power=power, radius=radius, radius_z=radius_z,
+ max_neighbors=max_neighbors, min_neighbors=min_neighbors,
log10=log10, gate_doi=gate_doi,
)
results.append(value)
@@ -285,7 +323,7 @@ def idw(samples, queries, power=2.0, radius=None, max_neighbors=None,
def interpolate_voxels(layer_rows, grid, *, doi_points=None,
clip_below_doi=True, power=2.0, radius=None,
- max_neighbors=None, min_neighbors=1, anisotropy_z=1.0,
+ radius_z=None, max_neighbors=None, min_neighbors=1,
log10=True, progress=None):
if grid.cell_count() > VOXEL_SAFETY_CAP:
raise ValueError(
@@ -301,6 +339,8 @@ def interpolate_voxels(layer_rows, grid, *, doi_points=None,
if radius is None:
radius = _default_radius(samples)
+ if radius_z is None:
+ radius_z = _default_radius_z(samples)
arrays = _SampleArrays(samples)
buckets, bucket_size = _bucket_index(arrays, radius)
@@ -313,8 +353,8 @@ def interpolate_voxels(layer_rows, grid, *, doi_points=None,
idx = _candidate_indices(buckets, bucket_size, cx, cy)
res_val, num_valid = _interpolate_at(
arrays, idx, cx, cy, cz,
- power=power, radius=radius, max_neighbors=max_neighbors,
- min_neighbors=min_neighbors, anisotropy_z=anisotropy_z,
+ power=power, radius=radius, radius_z=radius_z,
+ max_neighbors=max_neighbors, min_neighbors=min_neighbors,
log10=log10, gate_doi=clip_below_doi,
)
done += 1
diff --git a/tem_loader/tem_loader.py b/tem_loader/tem_loader.py
@@ -385,17 +385,18 @@ class _InterpolationDialog(QDialog):
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._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)
@@ -438,8 +439,8 @@ class _InterpolationDialog(QDialog):
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('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)
@@ -553,6 +554,7 @@ class _InterpolationDialog(QDialog):
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(),
@@ -561,8 +563,8 @@ class _InterpolationDialog(QDialog):
'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,
+ '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,
}
@@ -737,9 +739,9 @@ class TEMLoaderPlugin:
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'],
- anisotropy_z=opts['anisotropy_z'],
log10=opts['log10'],
progress=on_progress,
)
diff --git a/test/test_interpolate.py b/test/test_interpolate.py
@@ -5,6 +5,7 @@ from tem_loader.interpolate import (
GridSpec,
Sample,
_default_radius,
+ _default_radius_z,
axis_cell_size,
axis_count,
doi_by_xy_from_points,
@@ -132,21 +133,21 @@ class IDWTests(unittest.TestCase):
result = idw(samples, [(100.0, 0.0, 0.0)], radius=5.0, min_neighbors=2)
self.assertIsNone(result[0])
- def test_anisotropy_z_changes_vertical_weighting(self):
- # Sample A is 1 unit away horizontally; sample B is 10 units away
- # vertically. With anisotropy_z=1 the vertical distance dominates so A
- # wins; with anisotropy_z=10 the vertical gap shrinks so B gains weight.
+ def test_radius_z_excludes_distant_depths(self):
+ # Sample A sits 5 units away horizontally; sample B sits 50 units away
+ # vertically. A horizontal radius of 100 keeps both, but a tight vertical
+ # radius of 10 culls B, leaving A's value only.
samples = [
- self._s(1.0, 0.0, 0.0, 10.0),
- self._s(0.0, 0.0, 10.0, 100.0),
+ self._s(5.0, 0.0, 0.0, 10.0),
+ self._s(0.0, 0.0, 50.0, 100.0),
]
- result_tight = idw(samples, [(0.0, 0.0, 0.0)], anisotropy_z=1.0,
- log10=False, radius=None)
- result_loose = idw(samples, [(0.0, 0.0, 0.0)], anisotropy_z=10.0,
- log10=False, radius=None)
- self.assertIsNotNone(result_tight[0])
- self.assertIsNotNone(result_loose[0])
- self.assertNotAlmostEqual(result_tight[0], result_loose[0], places=1)
+ with_both = idw(samples, [(0.0, 0.0, 0.0)], radius=100.0, radius_z=100.0,
+ log10=False)
+ tight_z = idw(samples, [(0.0, 0.0, 0.0)], radius=100.0, radius_z=10.0,
+ log10=False)
+ self.assertIsNotNone(with_both[0])
+ self.assertAlmostEqual(tight_z[0], 10.0, places=6)
+ self.assertNotAlmostEqual(with_both[0], tight_z[0], places=1)
def test_log10_two_samples_geometric_mean(self):
samples = [
@@ -219,6 +220,18 @@ class DefaultRadiusTests(unittest.TestCase):
self.assertAlmostEqual(_default_radius(stacked), 90.0)
+class DefaultRadiusZTests(unittest.TestCase):
+ def test_stacked_column_median_gap(self):
+ # A single sounding stacked at 10 m spacing -> median gap 10 -> 1.5x.
+ column = [Sample(0.0, 0.0, z, 100.0, None)
+ for z in range(0, 100, 10)]
+ self.assertAlmostEqual(_default_radius_z(column), 15.0)
+
+ def test_returns_none_without_two_depths(self):
+ single = [Sample(0.0, 0.0, 0.0, 100.0, None)]
+ self.assertIsNone(_default_radius_z(single))
+
+
class SamplesFromLayersTests(unittest.TestCase):
def _row(self, x, y, zmid, res):
return {
@@ -345,6 +358,27 @@ class InterpolateVoxelsTests(unittest.TestCase):
self.assertIsNotNone(both, 'expected a voxel in both-usable zone at z=-45')
self.assertGreaterEqual(both['NumSamples'], 2)
+ def test_thin_layer_stays_sharp_with_radius_z(self):
+ # A single sounding with a sharp resistivity contrast at two depths.
+ # A tight vertical search radius must keep each voxel pinned to its own
+ # depth instead of bleeding the two values together across CellZ.
+ rows = [
+ self._make_layer_rows(0.0, 0.0, 0.0, 10.0),
+ self._make_layer_rows(0.0, 0.0, -50.0, 1000.0),
+ ]
+ g = GridSpec(-5.0, -5.0, -60.0, 10.0, 10.0, 10.0, 1, 1, 6)
+ voxels_by_z = {v['CellZ']: v for v in interpolate_voxels(
+ rows, g, radius=20.0, radius_z=10.0
+ )}
+ top = voxels_by_z.get(5) # center z = -5, nearest the 10 ohm-m sample
+ bottom = voxels_by_z.get(0) # center z = -55, nearest the 1000 ohm-m
+ self.assertIsNotNone(top)
+ self.assertIsNotNone(bottom)
+ self.assertAlmostEqual(top['Resistivity'], 10.0, places=6)
+ self.assertAlmostEqual(bottom['Resistivity'], 1000.0, places=6)
+ self.assertEqual(top['NumSamples'], 1)
+ self.assertEqual(bottom['NumSamples'], 1)
+
def test_safety_cap_raises(self):
bounds = (0.0, 200.0, 0.0, 200.0, 0.0, 200.0)
g = GridSpec.from_cell_size(bounds, 1.0, 1.0, 1.0)
