commit 200121e8731ef36cb89dccfd6a7f04176bcc96b0
parent c5d68d0c09afe3d3e2ecc1e4f94f822285456253
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 22 Oct 2026 09:04:29 +0000
feat(maps): compute interval mean resistivity grids
Diffstat:
2 files changed, 459 insertions(+), 0 deletions(-)
diff --git a/tem_loader/mean_resistivity.py b/tem_loader/mean_resistivity.py
@@ -0,0 +1,191 @@
+import math
+from collections import namedtuple
+
+import numpy as np
+
+from .interpolate import _default_radius, doi_by_xy_from_points
+
+MAP_CELL_CAP = 10_000_000
+TILE_SIZE = 32
+
+MapGrid = namedtuple('MapGrid', ['xmin', 'ymax', 'cell_size', 'nx', 'ny'])
+Sounding = namedtuple('Sounding', ['x', 'y', 'ztop', 'zbottom', 'res'])
+
+
+def map_grid(xs, ys, cell_size):
+ """North-up grid covering all soundings, with its extent snapped outward
+ to whole multiples of ``cell_size`` so that maps of the same survey align."""
+ if cell_size <= 0:
+ raise ValueError('cell size must be positive')
+ if len(xs) == 0:
+ raise ValueError('no soundings to grid')
+ xmin = math.floor(min(xs) / cell_size) * cell_size
+ xmax = math.ceil(max(xs) / cell_size) * cell_size
+ ymin = math.floor(min(ys) / cell_size) * cell_size
+ ymax = math.ceil(max(ys) / cell_size) * cell_size
+ nx = max(1, round((xmax - xmin) / cell_size))
+ ny = max(1, round((ymax - ymin) / cell_size))
+ if nx * ny > MAP_CELL_CAP:
+ raise ValueError(
+ f'map grid has {nx * ny} cells, exceeding the safety cap of '
+ f'{MAP_CELL_CAP}; increase the cell size'
+ )
+ return MapGrid(xmin, ymax, cell_size, nx, ny)
+
+
+def elevation_intervals(top, bottom, thickness):
+ """Intervals ``(top, bottom)`` from ``top`` downward; the last one is
+ truncated at ``bottom`` when the range is not a multiple of ``thickness``."""
+ if thickness <= 0:
+ raise ValueError('interval thickness must be positive')
+ if top <= bottom:
+ raise ValueError('top elevation must be above bottom elevation')
+ count = math.ceil(round((top - bottom) / thickness, 9))
+ intervals = []
+ for i in range(count):
+ upper = top - i * thickness
+ intervals.append((upper, max(bottom, upper - thickness)))
+ return intervals
+
+
+def soundings_from_layers(layer_rows, doi_by_xy=None):
+ """Group layer rows into soundings by exact (X, Y). Layers are clipped to
+ above the sounding's DOI elevation when ``doi_by_xy`` provides one; layers
+ without a positive resistivity are dropped."""
+ by_xy = {}
+ for row in layer_rows:
+ try:
+ x = float(row['X'])
+ y = float(row['Y'])
+ ztop = float(row['ZTop'])
+ zbottom = float(row['ZBottom'])
+ res = float(row['Resistivity'])
+ except (KeyError, ValueError, TypeError):
+ continue
+ if any(math.isnan(v) for v in (x, y, ztop, zbottom, res)) or res <= 0:
+ continue
+ if doi_by_xy is not None:
+ zdoi = doi_by_xy.get((x, y))
+ if zdoi is not None and not math.isnan(zdoi):
+ zbottom = max(zbottom, zdoi)
+ if ztop <= zbottom:
+ continue
+ by_xy.setdefault((x, y), []).append((ztop, zbottom, res))
+
+ soundings = []
+ for (x, y), layers in by_xy.items():
+ arr = np.asarray(layers, dtype=float)
+ soundings.append(Sounding(x, y, arr[:, 0], arr[:, 1], arr[:, 2]))
+ return soundings
+
+
+def interval_mean(ztop, zbottom, res, top, bottom, min_coverage):
+ """Thickness-weighted mean of log10(res) over the part of each layer inside
+ ``[bottom, top]``. Returns ``None`` when the layers cover less than
+ ``min_coverage`` of the interval."""
+ overlap = np.clip(np.minimum(ztop, top) - np.maximum(zbottom, bottom),
+ 0.0, None)
+ covered = float(overlap.sum())
+ if covered <= 0 or covered < min_coverage * (top - bottom) - 1e-9:
+ return None
+ return float((overlap * np.log10(res)).sum() / covered)
+
+
+def idw_grid(xs, ys, values, grid, power=2.0, radius=None, min_neighbors=1):
+ """Inverse-distance weighting of point values onto the grid nodes. Nodes with
+ fewer than ``min_neighbors`` points within ``radius`` are NaN; ``radius=None``
+ uses every point. Returns an array of shape ``(ny, nx)``, row 0 north."""
+ if min_neighbors < 1:
+ raise ValueError('min neighbors must be >= 1')
+ xs = np.asarray(xs, dtype=float)
+ ys = np.asarray(ys, dtype=float)
+ values = np.asarray(values, dtype=float)
+ out = np.full((grid.ny, grid.nx), np.nan)
+ if xs.size == 0:
+ return out
+
+ cx = grid.xmin + (np.arange(grid.nx) + 0.5) * grid.cell_size
+ cy = grid.ymax - (np.arange(grid.ny) + 0.5) * grid.cell_size
+
+ # Square tiles bound the distance matrix to tile nodes x nearby points.
+ for r0 in range(0, grid.ny, TILE_SIZE):
+ qy = cy[r0:r0 + TILE_SIZE]
+ for c0 in range(0, grid.nx, TILE_SIZE):
+ qx = cx[c0:c0 + TILE_SIZE]
+ if radius is None:
+ sx, sy, sv = xs, ys, values
+ else:
+ near = ((xs >= qx[0] - radius) & (xs <= qx[-1] + radius)
+ & (ys >= qy[-1] - radius) & (ys <= qy[0] + radius))
+ if near.sum() < min_neighbors:
+ continue
+ sx, sy, sv = xs[near], ys[near], values[near]
+
+ dx = qx[None, :, None] - sx
+ dy = qy[:, None, None] - sy
+ d2 = dx * dx + dy * dy
+ if radius is None:
+ inside = np.ones(d2.shape, dtype=bool)
+ else:
+ inside = d2 <= radius * radius
+ exact = inside & (d2 == 0.0)
+ n_exact = exact.sum(axis=-1)
+
+ with np.errstate(divide='ignore', invalid='ignore'):
+ weights = np.where(inside & ~exact, d2 ** (-power / 2.0), 0.0)
+ estimate = (weights * sv).sum(axis=-1) / weights.sum(axis=-1)
+ exact_mean = (exact * sv).sum(axis=-1) / n_exact
+ estimate = np.where(n_exact > 0, exact_mean, estimate)
+ estimate = np.where(inside.sum(axis=-1) >= min_neighbors,
+ estimate, np.nan)
+ out[r0:r0 + qy.size, c0:c0 + qx.size] = estimate
+ return out
+
+
+def mean_resistivity_maps(layer_rows, intervals, cell_size, *, doi_points=None,
+ power=2.0, radius=None, min_neighbors=1,
+ min_coverage=0.5, progress=None):
+ """Mean resistivity map per elevation interval: the thickness-weighted log10
+ mean of each sounding within the interval, gridded with IDW in log10 space.
+
+ Returns ``(grid, maps)`` where each map is ``(top, bottom, array, count)``:
+ ``array`` holds resistivity in ohm-m (NaN where blank) and ``count`` is the
+ number of soundings that contributed. Intervals without data are omitted."""
+ if not 0.0 <= min_coverage <= 1.0:
+ raise ValueError('minimum coverage must be between 0 and 1')
+ if radius is not None and radius <= 0:
+ raise ValueError('search radius must be positive')
+
+ doi_by_xy = None
+ if doi_points is not None:
+ doi_by_xy = doi_by_xy_from_points(doi_points)
+ soundings = soundings_from_layers(layer_rows, doi_by_xy)
+ if not soundings:
+ raise ValueError('no layers with positive resistivity to map')
+
+ if radius is None:
+ radius = _default_radius(soundings)
+ grid = map_grid([s.x for s in soundings], [s.y for s in soundings],
+ cell_size)
+
+ maps = []
+ total = len(intervals)
+ for done, (top, bottom) in enumerate(intervals, 1):
+ xs, ys, logs = [], [], []
+ for s in soundings:
+ value = interval_mean(s.ztop, s.zbottom, s.res, top, bottom,
+ min_coverage)
+ if value is not None:
+ xs.append(s.x)
+ ys.append(s.y)
+ logs.append(value)
+ if logs:
+ grid_logs = idw_grid(xs, ys, logs, grid, power=power,
+ radius=radius, min_neighbors=min_neighbors)
+ if np.isfinite(grid_logs).any():
+ array = (10.0 ** grid_logs).astype(np.float32)
+ maps.append((top, bottom, array, len(logs)))
+ # A truthy return from progress requests cancellation.
+ if progress is not None and progress(done, total):
+ break
+ return grid, maps
diff --git a/test/test_mean_resistivity.py b/test/test_mean_resistivity.py
@@ -0,0 +1,268 @@
+import math
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from tem_loader.core import process_xyz
+from tem_loader.mean_resistivity import (
+ MapGrid,
+ elevation_intervals,
+ idw_grid,
+ interval_mean,
+ map_grid,
+ mean_resistivity_maps,
+ soundings_from_layers,
+)
+
+FIXTURE_DIR = Path(__file__).parent / 'data'
+
+
+def _layer(x, y, ztop, zbottom, res):
+ return {'X': x, 'Y': y, 'ZTop': ztop, 'ZBottom': zbottom,
+ 'Resistivity': res}
+
+
+class ElevationIntervalTests(unittest.TestCase):
+ def test_whole_multiples(self):
+ self.assertEqual(
+ elevation_intervals(20.0, -10.0, 10.0),
+ [(20.0, 10.0), (10.0, 0.0), (0.0, -10.0)],
+ )
+
+ def test_last_interval_truncated(self):
+ self.assertEqual(
+ elevation_intervals(20.0, 5.0, 10.0),
+ [(20.0, 10.0), (10.0, 5.0)],
+ )
+
+ def test_rejects_invalid_input(self):
+ with self.assertRaises(ValueError):
+ elevation_intervals(20.0, 0.0, 0.0)
+ with self.assertRaises(ValueError):
+ elevation_intervals(0.0, 0.0, 10.0)
+ with self.assertRaises(ValueError):
+ elevation_intervals(-10.0, 0.0, 10.0)
+
+
+class MapGridTests(unittest.TestCase):
+ def test_extent_snapped_to_cell_size(self):
+ grid = map_grid([1012.0, 1188.0], [2049.0, 2101.0], 50.0)
+ self.assertEqual(grid, MapGrid(1000.0, 2150.0, 50.0, 4, 3))
+
+ def test_covers_all_points(self):
+ xs = [3.3, 97.1, 51.0]
+ ys = [-12.0, 40.5, 7.7]
+ grid = map_grid(xs, ys, 10.0)
+ self.assertLessEqual(grid.xmin, min(xs))
+ self.assertGreaterEqual(grid.xmin + grid.nx * grid.cell_size, max(xs))
+ self.assertGreaterEqual(grid.ymax, max(ys))
+ self.assertLessEqual(grid.ymax - grid.ny * grid.cell_size, min(ys))
+
+ def test_single_point(self):
+ grid = map_grid([100.0], [200.0], 25.0)
+ self.assertEqual((grid.nx, grid.ny), (1, 1))
+
+ def test_rejects_invalid_input(self):
+ with self.assertRaises(ValueError):
+ map_grid([0.0], [0.0], 0.0)
+ with self.assertRaises(ValueError):
+ map_grid([], [], 10.0)
+ with self.assertRaises(ValueError):
+ map_grid([0.0, 1e6], [0.0, 1e6], 1.0)
+
+
+class SoundingsFromLayersTests(unittest.TestCase):
+ def test_groups_by_xy(self):
+ rows = [
+ _layer(0.0, 0.0, 10.0, 5.0, 100.0),
+ _layer(0.0, 0.0, 5.0, 0.0, 10.0),
+ _layer(50.0, 0.0, 10.0, 0.0, 30.0),
+ ]
+ soundings = sorted(soundings_from_layers(rows))
+ self.assertEqual(len(soundings), 2)
+ self.assertEqual(soundings[0].res.tolist(), [100.0, 10.0])
+ self.assertEqual(soundings[1].ztop.tolist(), [10.0])
+
+ def test_skips_invalid_resistivity(self):
+ rows = [
+ _layer(0.0, 0.0, 10.0, 5.0, None),
+ _layer(0.0, 0.0, 10.0, 5.0, math.nan),
+ _layer(0.0, 0.0, 10.0, 5.0, 0.0),
+ _layer(0.0, 0.0, 10.0, 5.0, -5.0),
+ ]
+ self.assertEqual(soundings_from_layers(rows), [])
+
+ def test_clips_to_doi(self):
+ rows = [
+ _layer(0.0, 0.0, 10.0, 0.0, 100.0),
+ _layer(0.0, 0.0, 0.0, -20.0, 10.0),
+ ]
+ (sounding,) = soundings_from_layers(rows, {(0.0, 0.0): 4.0})
+ self.assertEqual(sounding.ztop.tolist(), [10.0])
+ self.assertEqual(sounding.zbottom.tolist(), [4.0])
+
+
+class IntervalMeanTests(unittest.TestCase):
+ def _mean(self, layers, top, bottom, min_coverage=0.5):
+ arr = np.asarray(layers, dtype=float)
+ return interval_mean(arr[:, 0], arr[:, 1], arr[:, 2], top, bottom,
+ min_coverage)
+
+ def test_single_covering_layer(self):
+ value = self._mean([(50.0, -50.0, 40.0)], 10.0, 0.0)
+ self.assertAlmostEqual(10.0 ** value, 40.0)
+
+ def test_equal_thickness_layers_give_geometric_mean(self):
+ value = self._mean([(10.0, 5.0, 10.0), (5.0, 0.0, 1000.0)], 10.0, 0.0)
+ self.assertAlmostEqual(10.0 ** value, 100.0)
+
+ def test_partial_overlap_weighted_by_thickness(self):
+ layers = [(20.0, 8.0, 10.0), (8.0, -20.0, 1000.0)]
+ value = self._mean(layers, 10.0, 0.0)
+ self.assertAlmostEqual(value, (2.0 * 1.0 + 8.0 * 3.0) / 10.0)
+
+ def test_low_coverage_returns_none(self):
+ self.assertIsNone(self._mean([(10.0, 6.0, 10.0)], 10.0, 0.0))
+ self.assertIsNotNone(self._mean([(10.0, 6.0, 10.0)], 10.0, 0.0,
+ min_coverage=0.4))
+
+ def test_no_overlap_returns_none(self):
+ self.assertIsNone(self._mean([(30.0, 20.0, 10.0)], 10.0, 0.0,
+ min_coverage=0.0))
+
+
+class IdwGridTests(unittest.TestCase):
+ def setUp(self):
+ self.grid = MapGrid(0.0, 10.0, 10.0, 3, 1)
+
+ def test_node_at_point_returns_its_value(self):
+ out = idw_grid([5.0, 25.0], [5.0, 5.0], [1.0, 3.0], self.grid,
+ radius=100.0)
+ self.assertAlmostEqual(out[0, 0], 1.0)
+ self.assertAlmostEqual(out[0, 2], 3.0)
+
+ def test_midpoint_is_mean(self):
+ out = idw_grid([5.0, 25.0], [5.0, 5.0], [1.0, 3.0], self.grid,
+ radius=100.0)
+ self.assertAlmostEqual(out[0, 1], 2.0)
+
+ def test_nodes_beyond_radius_are_blank(self):
+ out = idw_grid([5.0], [5.0], [1.0], self.grid, radius=12.0)
+ self.assertAlmostEqual(out[0, 0], 1.0)
+ self.assertAlmostEqual(out[0, 1], 1.0)
+ self.assertTrue(math.isnan(out[0, 2]))
+
+ def test_min_neighbors_enforced(self):
+ out = idw_grid([5.0, 25.0], [5.0, 5.0], [1.0, 3.0], self.grid,
+ radius=12.0, min_neighbors=2)
+ self.assertTrue(math.isnan(out[0, 0]))
+ self.assertAlmostEqual(out[0, 1], 2.0)
+ self.assertTrue(math.isnan(out[0, 2]))
+
+ def test_unlimited_radius_uses_all_points(self):
+ out = idw_grid([5.0], [5.0], [1.0], self.grid)
+ self.assertTrue(np.all(out == 1.0))
+
+ def test_rows_are_north_up(self):
+ grid = MapGrid(0.0, 20.0, 10.0, 1, 2)
+ out = idw_grid([5.0, 5.0], [15.0, 5.0], [1.0, 2.0], grid,
+ radius=1.0)
+ self.assertEqual(out[:, 0].tolist(), [1.0, 2.0])
+
+ def test_large_grid_spans_tiles(self):
+ grid = MapGrid(0.0, 1000.0, 10.0, 100, 100)
+ out = idw_grid([5.0, 995.0], [995.0, 5.0], [1.0, 2.0], grid,
+ radius=15.0)
+ self.assertAlmostEqual(out[0, 0], 1.0)
+ self.assertAlmostEqual(out[99, 99], 2.0)
+ self.assertEqual(int(np.isfinite(out).sum()), 2 * 4)
+
+
+class MeanResistivityMapsTests(unittest.TestCase):
+ def setUp(self):
+ self.rows = []
+ for x, res_top, res_bottom in ((0.0, 10.0, 100.0),
+ (100.0, 1000.0, 100.0)):
+ self.rows.append(_layer(x, 0.0, 10.0, 0.0, res_top))
+ self.rows.append(_layer(x, 0.0, 0.0, -10.0, res_bottom))
+
+ def test_maps_per_interval(self):
+ grid, maps = mean_resistivity_maps(
+ self.rows, [(10.0, 0.0), (0.0, -10.0)], 50.0, radius=500.0,
+ )
+ self.assertEqual((grid.nx, grid.ny), (2, 1))
+ self.assertEqual([(m[0], m[1]) for m in maps],
+ [(10.0, 0.0), (0.0, -10.0)])
+ top_map = maps[0][2]
+ self.assertEqual(top_map.shape, (1, 2))
+ self.assertTrue(np.all(top_map > 10.0))
+ self.assertTrue(np.all(top_map < 1000.0))
+ np.testing.assert_allclose(maps[1][2], 100.0, rtol=1e-6)
+ self.assertEqual(maps[0][3], 2)
+
+ def test_intervals_without_data_are_omitted(self):
+ _, maps = mean_resistivity_maps(
+ self.rows, [(50.0, 40.0), (10.0, 0.0), (-50.0, -60.0)], 50.0,
+ )
+ self.assertEqual([(m[0], m[1]) for m in maps], [(10.0, 0.0)])
+
+ def test_doi_points_clip_layers(self):
+ doi = [{'X': 0.0, 'Y': 0.0, 'ZDOI': 0.0},
+ {'X': 100.0, 'Y': 0.0, 'ZDOI': 0.0}]
+ _, maps = mean_resistivity_maps(
+ self.rows, [(10.0, 0.0), (0.0, -10.0)], 50.0, doi_points=doi,
+ )
+ self.assertEqual([(m[0], m[1]) for m in maps], [(10.0, 0.0)])
+
+ def test_progress_can_cancel(self):
+ calls = []
+
+ def progress(done, total):
+ calls.append((done, total))
+ return True
+
+ _, maps = mean_resistivity_maps(
+ self.rows, [(10.0, 0.0), (0.0, -10.0)], 50.0, progress=progress,
+ )
+ self.assertEqual(calls, [(1, 2)])
+ self.assertEqual(len(maps), 1)
+
+ def test_rejects_invalid_parameters(self):
+ with self.assertRaises(ValueError):
+ mean_resistivity_maps(self.rows, [(10.0, 0.0)], 50.0,
+ min_coverage=1.5)
+ with self.assertRaises(ValueError):
+ mean_resistivity_maps(self.rows, [(10.0, 0.0)], 50.0, radius=0.0)
+ with self.assertRaises(ValueError):
+ mean_resistivity_maps([], [(10.0, 0.0)], 50.0)
+
+ def test_fixtures(self):
+ for name in ('stem_80x80_workbench_2024_2_0_0.xyz',
+ 'stem_temimager_4_0_4_6.xyz'):
+ with self.subTest(name=name):
+ _, doi_points, layers = process_xyz(FIXTURE_DIR / name)
+ tops = [row['ZTop'] for row in layers]
+ bottoms = [row['ZBottom'] for row in layers]
+ intervals = elevation_intervals(
+ math.ceil(max(tops) / 10.0) * 10.0,
+ math.floor(min(bottoms) / 10.0) * 10.0,
+ 10.0,
+ )
+ _, maps = mean_resistivity_maps(
+ layers, intervals, 100.0,
+ doi_points=doi_points or None,
+ )
+ self.assertTrue(maps)
+ self.assertLess(len(maps), len(intervals) + 1)
+ res = [row['Resistivity'] for row in layers
+ if row['Resistivity'] > 0]
+ for top, bottom, array, count in maps:
+ finite = array[np.isfinite(array)]
+ self.assertGreater(count, 0)
+ self.assertGreaterEqual(finite.min(), min(res) * 0.999)
+ self.assertLessEqual(finite.max(), max(res) * 1.001)
+
+
+if __name__ == '__main__':
+ unittest.main()
