mean_resistivity.py (7995B)
1 import math 2 from collections import namedtuple 3 4 import numpy as np 5 6 from .interpolate import _default_radius, doi_by_xy_from_points 7 8 MAP_CELL_CAP = 10_000_000 9 TILE_SIZE = 32 10 11 MapGrid = namedtuple('MapGrid', ['xmin', 'ymax', 'cell_size', 'nx', 'ny']) 12 Sounding = namedtuple('Sounding', ['x', 'y', 'ztop', 'zbottom', 'res']) 13 14 15 def map_grid(xs, ys, cell_size): 16 """North-up grid covering all soundings, with its extent snapped outward 17 to whole multiples of ``cell_size`` so that maps of the same survey align.""" 18 if cell_size <= 0: 19 raise ValueError('cell size must be positive') 20 if len(xs) == 0: 21 raise ValueError('no soundings to grid') 22 xmin = math.floor(min(xs) / cell_size) * cell_size 23 xmax = math.ceil(max(xs) / cell_size) * cell_size 24 ymin = math.floor(min(ys) / cell_size) * cell_size 25 ymax = math.ceil(max(ys) / cell_size) * cell_size 26 nx = max(1, round((xmax - xmin) / cell_size)) 27 ny = max(1, round((ymax - ymin) / cell_size)) 28 if nx * ny > MAP_CELL_CAP: 29 raise ValueError( 30 f'map grid has {nx * ny} cells, exceeding the safety cap of ' 31 f'{MAP_CELL_CAP}; increase the cell size' 32 ) 33 return MapGrid(xmin, ymax, cell_size, nx, ny) 34 35 36 def elevation_intervals(top, bottom, thickness): 37 """Intervals ``(top, bottom)`` from ``top`` downward; the last one is 38 truncated at ``bottom`` when the range is not a multiple of ``thickness``.""" 39 if thickness <= 0: 40 raise ValueError('interval thickness must be positive') 41 if top <= bottom: 42 raise ValueError('top elevation must be above bottom elevation') 43 count = math.ceil(round((top - bottom) / thickness, 9)) 44 intervals = [] 45 for i in range(count): 46 upper = top - i * thickness 47 intervals.append((upper, max(bottom, upper - thickness))) 48 return intervals 49 50 51 def soundings_from_layers(layer_rows, doi_by_xy=None): 52 """Group layer rows into soundings by exact (X, Y). Layers are clipped to 53 above the sounding's DOI elevation when ``doi_by_xy`` provides one; layers 54 without a positive resistivity are dropped.""" 55 by_xy = {} 56 for row in layer_rows: 57 try: 58 x = float(row['X']) 59 y = float(row['Y']) 60 ztop = float(row['ZTop']) 61 zbottom = float(row['ZBottom']) 62 res = float(row['Resistivity']) 63 except (KeyError, ValueError, TypeError): 64 continue 65 if any(math.isnan(v) for v in (x, y, ztop, zbottom, res)) or res <= 0: 66 continue 67 if doi_by_xy is not None: 68 zdoi = doi_by_xy.get((x, y)) 69 if zdoi is not None and not math.isnan(zdoi): 70 zbottom = max(zbottom, zdoi) 71 if ztop <= zbottom: 72 continue 73 by_xy.setdefault((x, y), []).append((ztop, zbottom, res)) 74 75 soundings = [] 76 for (x, y), layers in by_xy.items(): 77 arr = np.asarray(layers, dtype=float) 78 soundings.append(Sounding(x, y, arr[:, 0], arr[:, 1], arr[:, 2])) 79 return soundings 80 81 82 def interval_mean(ztop, zbottom, res, top, bottom, min_coverage): 83 """Thickness-weighted mean of log10(res) over the part of each layer inside 84 ``[bottom, top]``. Returns ``None`` when the layers cover less than 85 ``min_coverage`` of the interval.""" 86 overlap = np.clip(np.minimum(ztop, top) - np.maximum(zbottom, bottom), 87 0.0, None) 88 covered = float(overlap.sum()) 89 if covered <= 0 or covered < min_coverage * (top - bottom) - 1e-9: 90 return None 91 return float((overlap * np.log10(res)).sum() / covered) 92 93 94 def idw_grid(xs, ys, values, grid, power=2.0, radius=None, min_neighbors=1): 95 """Inverse-distance weighting of point values onto the grid nodes. Nodes with 96 fewer than ``min_neighbors`` points within ``radius`` are NaN; ``radius=None`` 97 uses every point. Returns an array of shape ``(ny, nx)``, row 0 north.""" 98 if min_neighbors < 1: 99 raise ValueError('min neighbors must be >= 1') 100 xs = np.asarray(xs, dtype=float) 101 ys = np.asarray(ys, dtype=float) 102 values = np.asarray(values, dtype=float) 103 out = np.full((grid.ny, grid.nx), np.nan) 104 if xs.size == 0: 105 return out 106 107 cx = grid.xmin + (np.arange(grid.nx) + 0.5) * grid.cell_size 108 cy = grid.ymax - (np.arange(grid.ny) + 0.5) * grid.cell_size 109 110 # Square tiles bound the distance matrix to tile nodes x nearby points. 111 for r0 in range(0, grid.ny, TILE_SIZE): 112 qy = cy[r0:r0 + TILE_SIZE] 113 for c0 in range(0, grid.nx, TILE_SIZE): 114 qx = cx[c0:c0 + TILE_SIZE] 115 if radius is None: 116 sx, sy, sv = xs, ys, values 117 else: 118 near = ((xs >= qx[0] - radius) & (xs <= qx[-1] + radius) 119 & (ys >= qy[-1] - radius) & (ys <= qy[0] + radius)) 120 if near.sum() < min_neighbors: 121 continue 122 sx, sy, sv = xs[near], ys[near], values[near] 123 124 dx = qx[None, :, None] - sx 125 dy = qy[:, None, None] - sy 126 d2 = dx * dx + dy * dy 127 if radius is None: 128 inside = np.ones(d2.shape, dtype=bool) 129 else: 130 inside = d2 <= radius * radius 131 exact = inside & (d2 == 0.0) 132 n_exact = exact.sum(axis=-1) 133 134 with np.errstate(divide='ignore', invalid='ignore'): 135 weights = np.where(inside & ~exact, d2 ** (-power / 2.0), 0.0) 136 estimate = (weights * sv).sum(axis=-1) / weights.sum(axis=-1) 137 exact_mean = (exact * sv).sum(axis=-1) / n_exact 138 estimate = np.where(n_exact > 0, exact_mean, estimate) 139 estimate = np.where(inside.sum(axis=-1) >= min_neighbors, 140 estimate, np.nan) 141 out[r0:r0 + qy.size, c0:c0 + qx.size] = estimate 142 return out 143 144 145 def mean_resistivity_maps(layer_rows, intervals, cell_size, *, doi_points=None, 146 power=2.0, radius=None, min_neighbors=1, 147 min_coverage=0.5, progress=None): 148 """Mean resistivity map per elevation interval: the thickness-weighted log10 149 mean of each sounding within the interval, gridded with IDW in log10 space. 150 151 Returns ``(grid, maps)`` where each map is ``(top, bottom, array, count)``: 152 ``array`` holds resistivity in ohm-m (NaN where blank) and ``count`` is the 153 number of soundings that contributed. Intervals without data are omitted.""" 154 if not 0.0 <= min_coverage <= 1.0: 155 raise ValueError('minimum coverage must be between 0 and 1') 156 if radius is not None and radius <= 0: 157 raise ValueError('search radius must be positive') 158 159 doi_by_xy = None 160 if doi_points is not None: 161 doi_by_xy = doi_by_xy_from_points(doi_points) 162 soundings = soundings_from_layers(layer_rows, doi_by_xy) 163 if not soundings: 164 raise ValueError('no layers with positive resistivity to map') 165 166 if radius is None: 167 radius = _default_radius(soundings) 168 grid = map_grid([s.x for s in soundings], [s.y for s in soundings], 169 cell_size) 170 171 maps = [] 172 total = len(intervals) 173 for done, (top, bottom) in enumerate(intervals, 1): 174 xs, ys, logs = [], [], [] 175 for s in soundings: 176 value = interval_mean(s.ztop, s.zbottom, s.res, top, bottom, 177 min_coverage) 178 if value is not None: 179 xs.append(s.x) 180 ys.append(s.y) 181 logs.append(value) 182 if logs: 183 grid_logs = idw_grid(xs, ys, logs, grid, power=power, 184 radius=radius, min_neighbors=min_neighbors) 185 if np.isfinite(grid_logs).any(): 186 array = (10.0 ** grid_logs).astype(np.float32) 187 maps.append((top, bottom, array, len(logs))) 188 # A truthy return from progress requests cancellation. 189 if progress is not None and progress(done, total): 190 break 191 return grid, maps