Refraction Manager

This example shows how to use the Refraction manager to generate the response of a three-layered sloping model and to invert the synthetic noisified data.

import numpy as np

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import Refraction

We start by creating a three-layered slope (The model is taken from the BSc thesis of Constanze Reinken (University of Bonn).

layer1 = mt.createPolygon([[0.0, 137], [117.5, 164], [117.5, 162], [0.0, 135]],
                          isClosed=True, marker=1, area=1)
layer2 = mt.createPolygon([[0.0, 126], [0.0, 135], [117.5, 162], [117.5, 153]],
                          isClosed=True, marker=2)
layer3 = mt.createPolygon([[0.0, 110], [0.0, 126], [117.5, 153], [117.5, 110]],
                          isClosed=True, marker=3)

slope = (164 - 137) / 117.5

geom = mt.mergePLC([layer1, layer2, layer3])
mesh = mt.createMesh(geom, quality=34.3, area=3, smooth=[1, 10])


(<matplotlib.axes._subplots.AxesSubplot object at 0x7f8818b5d1d0>, None)

Next we define geophone positions and a measurement scheme, which consists of shot and receiver indices.

numberGeophones = 48
sensors = np.linspace(0., 117.5, numberGeophones)
scheme = pg.physics.traveltime.createRAData(sensors)

# Adapt sensor positions to slope
pos = np.array(scheme.sensorPositions())
for x in pos[:, 0]:
    i = np.where(pos[:, 0] == x)
    new_y = x * slope + 137
    pos[i, 1] = new_y


Now we initialize the refraction manager and asssign P-wave velocities to the layers. To this end, we create a map from cell markers 0 through 3 to velocities (in m/s) and generate a velocity vector. To check whether the it is correct, we plot it # along with the sensor positions

ra = Refraction()
vp = np.array(mesh.cellMarkers())
vp[vp == 1] = 250
vp[vp == 2] = 500
vp[vp == 3] = 1300

ax, _ =, vp, colorBar=True, logScale=False, label='v in m/s')
ax.plot(pos[:, 0], pos[:, 1], 'w+')


[<matplotlib.lines.Line2D object at 0x7f8803aace10>]

We use this model to create noisified synthetic data and look at the traveltime curves showing three different slopes.

data = ra.simulate(mesh, 1.0 / vp, scheme, noiseLevel=0.001, noiseAbs=0.001,
# ra.showData(data)  # can be used to show the data.


Data error estimates (min:max)  0.00100753513903737 : 0.001142191584163894

And invert the synthetic data on an independent mesh without # information on the layered structure we create a new instance of the class using the data. Instead of the data, a file name can be given. This is probably where most users with data start. See refraction class for supported formats.

ra = Refraction(data)
ra.createMesh(depth=30., paraMaxCellSize=5.0, secNodes=1)
# ra.fop.createJacobian(np.ones(ra.fop.regionManager().parameterCount()))
# exit()
vest = ra.invert()  # estimated velocity distribution


Data: Sensors: 48 data: 2256

The method showResult is used to plot the result. Note that only covered cells are shown by default. For comparison we plot the geometry on top.

ax, cb = ra.showResult(cMin=min(vp), cMax=max(vp), logScale=False), ax=ax, fillRegion=False, regionMarker=False)  # lines on top


(<matplotlib.axes._subplots.AxesSubplot object at 0x7f880007ca20>, None)

Note that internally the following is called

ax, _ =, vest, label="Velocity [m/s]", **kwargs)

Another useful tool is to show the model along with its respone on the data_


Takeaway message: A default data inversion with checking of the data consists of only few lines (Everthing else can be looked at by introspecting the Refraction manager)

from pygimli.physics import Refraction
ra = Refraction(filename)

Total running time of the script: ( 0 minutes 30.956 seconds)

Gallery generated by Sphinx-Gallery

2019 - GIMLi Development Team
Created using Bootstrap, Sphinx and pyGIMLi 1.0.12+20.g65e71e67 on Mar 10, 2020.