Source code for pygimli.physics.traveltime.utils

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Tools for traveltime/refraction tomography."""

import numpy as np
import pygimli as pg

[docs] def createCrossholeData(sensors=None, x=None, z=None): """ Create crosshole scheme assuming two boreholes with equal sensor numbers. Parameters ---------- sensors : array (Nx2) Array with position of sensors. Alternatively, specify x and z x : iterable [M] horizontal positions of boreholes z : iterable [M] vertical positions of sensors in each borehole Returns ------- scheme : DataContainer Data container with `sensors` predefined sensor indices 's' and 'g' for shot and receiver numbers. """ from itertools import product if sensors is None: if x is None or z is None: raise Exception("Both x and z must be given") sensors = np.array([(xi, zi) for zi in z for xi in x]) if len(sensors) % 2 > 0: pg.error("createCrossholeData is only defined for an equal number of" " sensors in two boreholes.") n = len(sensors) // 2 numbers = np.arange(n) rays = np.array(list(product(numbers, numbers + n))) # Empty container scheme = pg.DataContainer() # Add sensors for sen in sensors: scheme.createSensor(sen) # Add measurements scheme.resize(len(rays)) scheme["s"] = rays[:, 0] scheme["g"] = rays[:, 1] scheme["valid"] = np.ones(len(rays)) scheme.registerSensorIndex("s") scheme.registerSensorIndex("g") return scheme
[docs] def shotReceiverDistances(data, full=True): """Distance vector between shot and for each 's' and 'g' in data. Parameters ---------- data : pg.DataContainerERT full : bool [True] Get distances between shot and receiver position when full is True or only form x coordinate if full is False Returns ------- dists : array Array of distances """ if full: pos = data.sensors() s, g ="s"),"g") off = [pos[s[i]].distance(pos[g[i]]) for i in range(data.size())] return np.absolute(off) else: px = pg.x(data) gx = np.array([px[g] for g in"g")]) sx = np.array([px[s] for s in"s")]) return np.absolute(gx - sx)
[docs] def createRAData(sensors, shotDistance=1): """Create a refraction data container. Default data container for shot and geophon at every sensor position. Predefined sensor indices's 's' as shot position and 'g' as geophon position. Parameters ---------- sensors: ndarray | R3Vector Geophon and shot positions (same) shotDistances: int [1] Distance between shot indices. Returns ------- data : DataContainer Data container with predefined sensor indices 's' and 'g' for shot and receiver numbers. """ data = pg.DataContainer() data.registerSensorIndex("s") data.registerSensorIndex("g") if isinstance(sensors, np.ndarray): if len(sensors.shape) == 1: for x in sensors: data.createSensor([x, 0.0, 0.0]) else: data.setSensorPositions(sensors) else: data.setSensorPositions(sensors) S, G = [], [] for s in range(0, data.sensorCount(), shotDistance): for g in range(data.sensorCount()): if s is not g: S.append(s) G.append(g) data.resize(len(S)) data.set("s", S) data.set("g", G) data.set("valid", np.abs(np.sign(data("g") - data("s")))) return data
[docs] def createGradientModel2D(data, mesh, vTop, vBot, flat=False): """Create 2D velocity gradient model. Creates a smooth, linear, starting model that takes the slope of the topography into account. This is done by fitting a straight line and using the distance to that as the depth value. Parameters ---------- data: pygimli DataContainer The topography list is in here. mesh: pygimli.Mesh The parametric mesh used for the inversion vTop: float The velocity at the surface of the mesh vBot: float The velocity at the bottom of the mesh Returns ------- model: pygimli Vector, length M A numpy array with slowness values that can be used to start the inversion. """ xVals = pg.x(data) yVals = pg.y(data) if mesh.dim() == 3 or abs(min(yVals)) < 1e-8 and abs(max(yVals)) < 1e-8: yVals = pg.z(data) x = pg.x(mesh.cellCenters()) z = pg.y(mesh.cellCenters()) if mesh.dim() == 3 or abs(min(z)) < 1e-8 and abs(max(z)) < 1e-8: z = pg.z(mesh.cellCenters()) pos = np.column_stack((x, z)) if flat: p = np.polyfit(xVals, yVals, deg=1) # slope-intercept form n = np.asarray([-p[0], 1.0]) # normal vector nLen = np.sqrt(, n)) d = np.array([np.abs([i, :], n) - p[1]) / nLen for i in range(pos.shape[0])]) else: d = np.interp(x, xVals, yVals) - z return 1.0 / np.interp(d, [min(d), max(d)], [vTop, vBot])