pygimli.physics

Module containing submodules for various geophysical methods.

Overview

Classes

Constants Container class for some constants repeatedly used in geophysics.
ERTManager(**kwargs) Minimalistic ERT Manager for compatibility (advanced version in BERT)
ERTModelling([mesh, data, verbose]) Minimal Forward Operator for 2.5D Electrical resistivity Tomography.
FDEM([x, freqs, coilSpacing, inphase, …]) Class for managing Frequency Domain EM data and their inversions.
MRS([name, verbose]) Magnetic resonance sounding (MRS) manager class.
Refraction([data, verbose, debug, fatray, …]) Manager for refraction seismics (traveltime tomography)
SIPSpectrum([filename, unify, onlydown, f, …]) SIP spectrum data analysis.
TDEM([filename]) TEM class mainly for holding data etc.
VESManager(**kwargs) Vertical electrical sounding (VES) manager class.
constants alias of pygimli.physics.constants.Constants

Classes

Constants

class pygimli.physics.Constants[source]

Container class for some constants repeatedly used in geophysics.

Darcy = 9.86923e-13
G = 6.6742e-11
GmGal = 6.6741999999999996e-06
__init__($self, /, *args, **kwargs)

Initialize self. See help(type(self)) for accurate signature.

c = 299792458.0000066
c0 = 299792458.0000066
e0 = 8.85418781762e-12
fLarmor = 42576375.18573887
fLarmorMhZ = 42.57637518573887
g = 9.80665
gammaP = 267515255.0
mu0 = 1.2566370614359173e-06
nepers2dB = 8.686

ERTManager

class pygimli.physics.ERTManager(**kwargs)[source]

Minimalistic ERT Manager for compatibility (advanced version in BERT)

Methods

apparentData() Convert data into apparent data.
checkData() Check data validity.
createApparentData(data) Create apparent data (what the hack is this?)
createArgParser([dataSuffix]) Create argument parser for the manager.
createData(sensors, scheme) Create an empty data set.
createFOP([verbose]) Create forward operator working on refined mesh.
createFOP_([verbose]) Create forward operator working on refined mesh.
createInv(fop[, verbose, dosave]) Create inversion instance.
createInv_(fop[, verbose, dosave]) Create inversion instance, data- and model transformations.
createMesh(**kwargs) Create a mesh aka the parametrization.
dataToken() Token name for the data in a DataContainer.
dataVals(data) Return pure data values from a given DataContainer.
estimateError(data[, absoluteError, …]) Estimate error composed of an absolute and a relative part.
invert([data, vals, err, mesh]) Run the full inversion.
model() Return the actual model.
relErrorVals(data) Return pure data values from a given DataContainer.
saveResult([folder, size]) Save results in the specified folder.
setData(data) Set data container from outside.
setDataToken(token) Set the token name to identity the data in a DataContainer.
setMesh(mesh[, refine]) Set the internal mesh for this Manager.
setVerbose(verbose) Make the class verbose (put output to the console)
show(data[, values, axes, cMin, cMax, colorBar]) Forward the visualization.
showData([data, vals, ax]) Show mesh in given axes or in a new figure.
showMesh([ax]) Show mesh in given axes or in a new figure.
showResult([ax, cMin, cMax, logScale]) Show resulting vector.
simulate(mesh, res, scheme[, verbose]) Forward calculation vor given mesh, data and resistivity.
__init__(**kwargs)[source]

Constructor.

apparentData()

Convert data into apparent data.

checkData()

Check data validity.

createApparentData(data)[source]

Create apparent data (what the hack is this?)

static createArgParser(dataSuffix='dat')

Create argument parser for the manager.

createData(sensors, scheme)

Create an empty data set.

static createFOP(verbose=False)[source]

Create forward operator working on refined mesh.

createFOP_(verbose=False)

Create forward operator working on refined mesh.

createInv(fop, verbose=True, dosave=False)[source]

Create inversion instance.

createInv_(fop, verbose=True, dosave=False)

Create inversion instance, data- and model transformations.

createMesh(**kwargs)

Create a mesh aka the parametrization.

dataToken()

Token name for the data in a DataContainer.

dataVals(data)[source]

Return pure data values from a given DataContainer.

static estimateError(data, absoluteError=0.001, relativeError=0.001)

Estimate error composed of an absolute and a relative part.

invert(data=None, vals=None, err=None, mesh=None, **kwargs)

Run the full inversion.

The data and error needed to be set before. The meshes will be created if necessary.

DOCUMENTME!!!

Parameters:
lam : float [20]

regularization parameter

zWeight : float [0.7]

relative vertical weight

maxIter : int [20]

maximum iteration number

robustdata : bool [False]

robust data reweighting using an L1 scheme (IRLS reweighting)

blockymodel : bool [False]

blocky model constraint using L1 reweighting roughness vector

startModelIsReference : bool [False]

startmodel is the reference model for the inversion

forwarded to createMesh
depth
quality
paraDX
maxCellArea
model()

Return the actual model.

relErrorVals(data)[source]

Return pure data values from a given DataContainer.

saveResult(folder=None, size=(16, 10), **kwargs)

Save results in the specified folder.

setData(data)

Set data container from outside.

setDataToken(token)

Set the token name to identity the data in a DataContainer.

setMesh(mesh, refine=True)

Set the internal mesh for this Manager.

Inject the mesh in the internal fop und inv.

Initialize RegionManager. For more than two regions the first is assumed to be background.

Optional the forward mesh can be refined for higher numerical accuracy.

Parameters:
DOCUMENTME!!!
setVerbose(verbose)

Make the class verbose (put output to the console)

show(data, values=None, axes=None, cMin=None, cMax=None, colorBar=1, **kwargs)

Forward the visualization.

showData(data=None, vals=None, ax=None)[source]

Show mesh in given axes or in a new figure.

showMesh(ax=None)

Show mesh in given axes or in a new figure.

showResult(ax=None, cMin=None, cMax=None, logScale=False, **kwargs)

Show resulting vector.

static simulate(mesh, res, scheme, verbose=False, **kwargs)[source]

Forward calculation vor given mesh, data and resistivity.

ERTModelling

class pygimli.physics.ERTModelling(mesh=None, data=None, verbose=False)[source]

Minimal Forward Operator for 2.5D Electrical resistivity Tomography.

Methods

__call__((object)arg1, (object)model) C++ signature : GIMLI::Vector<double> __call__(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double>)
calcGeometricFactors(data) Calculate geometry factors for a given dataset.
clearConstraints((object)arg1) C++ signature : void* clearConstraints(GIMLI::ModellingBase {lvalue})
clearJacobian((object)arg1) C++ signature : void* clearJacobian(GIMLI::ModellingBase {lvalue})
constraints((object)arg1) C++ signature : GIMLI::MatrixBase* constraints(GIMLI::ModellingBase {lvalue})
constraintsRef((object)arg1) C++ signature : GIMLI::SparseMapMatrix<double, unsigned long> {lvalue} constraintsRef(GIMLI::ModellingBase {lvalue})
createConstraints((object)arg1) C++ signature : void* createConstraints(GIMLI::ModellingBase {lvalue})
createDefaultStartModel((object)arg1) C++ signature : GIMLI::Vector<double> createDefaultStartModel(GIMLI::ModellingBase {lvalue})
createJacobian(model) Create Jacobian matrix.
createMappedModel((object)arg1, …) Read only extrapolation of model values given per cell marker to values given per cell.
createRHS(mesh, elecs) Create right-hand side.
createRefinedForwardMesh((object)arg1 [, …) C++ signature : void* createRefinedForwardMesh(GIMLI::ModellingBase {lvalue} [,bool=True [,bool=False]])
createStartModel((object)arg1) C++ signature : GIMLI::Vector<double> createStartModel(GIMLI::ModellingBase {lvalue})
createStartVector((object)arg1) DEPRECATED use createStartModel
data((object)arg1) Return the associated data container.
deleteMesh((object)arg1) Delete the actual mesh.
getIntegrationWeights(rMin, rMax) Retrieve Gauss-Legende/Laguerre integration weights.
initConstraints((object)arg1) C++ signature : void* initConstraints(GIMLI::ModellingBase {lvalue})
initJacobian((object)arg1) C++ signature : void* initJacobian(GIMLI::ModellingBase {lvalue})
initRegionManager((object)arg1) C++ signature : void* initRegionManager(GIMLI::ModellingBase {lvalue})
jacobian((object)arg1) Return the pointer to the Jacobian matrix associated with this forward operator.
jacobianRef((object)arg1) C++ signature : GIMLI::Matrix<double> {lvalue} jacobianRef(GIMLI::ModellingBase {lvalue})
mapModel((object)arg1, (object)model [, …) C++ signature : void* mapModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,double=0])
mesh((object)arg1) C++ signature : GIMLI::Mesh* mesh(GIMLI::ModellingBase {lvalue})
mixedBC(boundary, userData) Apply mixed boundary conditions.
multiThreadJacobian((object)arg1) Return number of threads used for Jacobian generation.
pointSource(cell, f, userData) Define function for the current source term.
region((object)arg1, (object)marker) Syntactic sugar for this->regionManager().region(marker).
regionManager((object)arg1) C++ signature : GIMLI::RegionManager regionManager(GIMLI::ModellingBase {lvalue})
regionManagerRef((object)arg1) C++ signature : GIMLI::RegionManager {lvalue} regionManagerRef(GIMLI::ModellingBase {lvalue})
response(model) Solve forward task.
response_mt((object)arg1, (object)model [, …) C++ signature : GIMLI::Vector<double> response_mt(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,unsigned long=0])
setConstraints((object)arg1, (object)C) C++ signature : void* setConstraints(GIMLI::ModellingBase {lvalue},GIMLI::MatrixBase*)
setData(data) Set a DataContainer.
setJacobian((object)arg1, (object)J) C++ signature : void* setJacobian(GIMLI::ModellingBase {lvalue},GIMLI::MatrixBase*)
setMesh(mesh[, ignoreRegionManager]) Set Mesh.
setMultiThreadJacobian((object)arg1, …) Set number of threads used for brute force Jacobian generation.
setRegionManager((object)arg1, (object)reg) C++ signature : void* setRegionManager(GIMLI::ModellingBase {lvalue},GIMLI::RegionManager*)
setStartModel((object)arg1, (object)startModel) C++ signature : void* setStartModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double>)
setThreadCount((object)arg1, (object)nThreads) Set the maximum number of allowed threads for MT calculation.
setVerbose((object)arg1, (object)verbose) Set verbose state.
solution((object)arg1) C++ signature : GIMLI::Matrix<double> solution(GIMLI::ModellingBase {lvalue})
startModel((object)arg1) C++ signature : GIMLI::Vector<double> startModel(GIMLI::ModellingBase {lvalue})
threadCount((object)arg1) Return the maximum number of allowed threads for MT calculation
uAnalytical(p, sourcePos, k) Calculates the analytical solution for the 2.5D geoelectrical problem.
verbose((object)arg1) Get verbose state.
createJacobian_mt  
responses  
__init__(mesh=None, data=None, verbose=False)[source]

Constructor, optionally with data container and mesh.

calcGeometricFactors(data)[source]

Calculate geometry factors for a given dataset.

clearConstraints((object)arg1) → object :
C++ signature :
void* clearConstraints(GIMLI::ModellingBase {lvalue})

clearConstraints( (object)arg1) -> object :

C++ signature :
void* clearConstraints(ModellingBase_wrapper {lvalue})
clearJacobian((object)arg1) → object :
C++ signature :
void* clearJacobian(GIMLI::ModellingBase {lvalue})

clearJacobian( (object)arg1) -> object :

C++ signature :
void* clearJacobian(ModellingBase_wrapper {lvalue})
constraints((object)arg1) → object :
C++ signature :
GIMLI::MatrixBase* constraints(GIMLI::ModellingBase {lvalue})

constraints( (object)arg1) -> object :

C++ signature :
GIMLI::MatrixBase* constraints(ModellingBase_wrapper {lvalue})

constraints( (object)arg1) -> object :

C++ signature :
GIMLI::MatrixBase* constraints(GIMLI::ModellingBase {lvalue})

constraints( (object)arg1) -> object :

C++ signature :
GIMLI::MatrixBase* constraints(ModellingBase_wrapper {lvalue})
constraintsRef((object)arg1) → object :
C++ signature :
GIMLI::SparseMapMatrix<double, unsigned long> {lvalue} constraintsRef(GIMLI::ModellingBase {lvalue})

constraintsRef( (object)arg1) -> object :

C++ signature :
GIMLI::SparseMapMatrix<double, unsigned long> {lvalue} constraintsRef(GIMLI::ModellingBase {lvalue})
createConstraints((object)arg1) → object :
C++ signature :
void* createConstraints(GIMLI::ModellingBase {lvalue})

createConstraints( (object)arg1) -> object :

C++ signature :
void* createConstraints(ModellingBase_wrapper {lvalue})
createDefaultStartModel((object)arg1) → object :
C++ signature :
GIMLI::Vector<double> createDefaultStartModel(GIMLI::ModellingBase {lvalue})

createDefaultStartModel( (object)arg1) -> object :

C++ signature :
GIMLI::Vector<double> createDefaultStartModel(ModellingBase_wrapper {lvalue})
createJacobian(model)[source]

Create Jacobian matrix.

createJacobian_mt(model, resp)
createMappedModel((object)arg1, (object)model[, (object)background=-1]) → object :

Read only extrapolation of model values given per cell marker to values given per cell. Exterior values will be set to background or prolongated for background = -1.

C++ signature :
GIMLI::Vector<double> createMappedModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,double=-1])
createRHS(mesh, elecs)[source]

Create right-hand side.

createRefinedForwardMesh((object)arg1[, (object)refine=True[, (object)pRefine=False]]) → object :
C++ signature :
void* createRefinedForwardMesh(GIMLI::ModellingBase {lvalue} [,bool=True [,bool=False]])
createStartModel((object)arg1) → object :
C++ signature :
GIMLI::Vector<double> createStartModel(GIMLI::ModellingBase {lvalue})
createStartVector((object)arg1) → object :

DEPRECATED use createStartModel

C++ signature :
GIMLI::Vector<double> createStartVector(GIMLI::ModellingBase {lvalue})
data((object)arg1) → object :

Return the associated data container.

C++ signature :
GIMLI::DataContainer {lvalue} data(GIMLI::ModellingBase {lvalue})
deleteMesh((object)arg1) → object :

Delete the actual mesh.

C++ signature :
void* deleteMesh(GIMLI::ModellingBase {lvalue})
getIntegrationWeights(rMin, rMax)[source]

Retrieve Gauss-Legende/Laguerre integration weights.

initConstraints((object)arg1) → object :
C++ signature :
void* initConstraints(GIMLI::ModellingBase {lvalue})

initConstraints( (object)arg1) -> object :

C++ signature :
void* initConstraints(ModellingBase_wrapper {lvalue})
initJacobian((object)arg1) → object :
C++ signature :
void* initJacobian(GIMLI::ModellingBase {lvalue})

initJacobian( (object)arg1) -> object :

C++ signature :
void* initJacobian(ModellingBase_wrapper {lvalue})
initRegionManager((object)arg1) → object :
C++ signature :
void* initRegionManager(GIMLI::ModellingBase {lvalue})
jacobian((object)arg1) → object :

Return the pointer to the Jacobian matrix associated with this forward operator.

C++ signature :
GIMLI::MatrixBase* jacobian(GIMLI::ModellingBase {lvalue})
jacobian( (object)arg1) -> object :

Return the pointer to the Jacobian matrix associated with this forward operator.

C++ signature :
GIMLI::MatrixBase* jacobian(GIMLI::ModellingBase {lvalue})
jacobianRef((object)arg1) → object :
C++ signature :
GIMLI::Matrix<double> {lvalue} jacobianRef(GIMLI::ModellingBase {lvalue})

jacobianRef( (object)arg1) -> object :

C++ signature :
GIMLI::Matrix<double> {lvalue} jacobianRef(GIMLI::ModellingBase {lvalue})
mapModel((object)arg1, (object)model[, (object)background=0]) → object :
C++ signature :
void* mapModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,double=0])
mesh((object)arg1) → object :
C++ signature :
GIMLI::Mesh* mesh(GIMLI::ModellingBase {lvalue})
mixedBC(boundary, userData)[source]

Apply mixed boundary conditions.

multiThreadJacobian((object)arg1) → object :

Return number of threads used for Jacobian generation.

C++ signature :
unsigned long multiThreadJacobian(GIMLI::ModellingBase {lvalue})
pointSource(cell, f, userData)[source]

Define function for the current source term.

\(\delta(x-pos), \int f(x) \delta(x-pos)=f(pos)=N(pos)\)
Right-hand-side entries will be shape functions(pos)
region((object)arg1, (object)marker) → object :

Syntactic sugar for this->regionManager().region(marker).

C++ signature :
GIMLI::Region* region(GIMLI::ModellingBase {lvalue},int)
regionManager((object)arg1) → object :
C++ signature :
GIMLI::RegionManager regionManager(GIMLI::ModellingBase {lvalue})

regionManager( (object)arg1) -> object :

C++ signature :
GIMLI::RegionManager {lvalue} regionManager(GIMLI::ModellingBase {lvalue})
regionManagerRef((object)arg1) → object :
C++ signature :
GIMLI::RegionManager {lvalue} regionManagerRef(GIMLI::ModellingBase {lvalue})
response(model)[source]

Solve forward task.

Create apparent resistivity values for a given resistivity distribution for self.mesh.

response_mt((object)arg1, (object)model[, (object)i=0]) → object :
C++ signature :
GIMLI::Vector<double> response_mt(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double> [,unsigned long=0])

response_mt( (object)arg1, (object)model [, (object)i=0]) -> object :

C++ signature :
GIMLI::Vector<double> response_mt(ModellingBase_wrapper {lvalue},GIMLI::Vector<double> [,unsigned long=0])
responses(models, respos)
setConstraints((object)arg1, (object)C) → object :
C++ signature :
void* setConstraints(GIMLI::ModellingBase {lvalue},GIMLI::MatrixBase*)

setConstraints( (object)arg1, (object)C) -> object :

C++ signature :
void* setConstraints(ModellingBase_wrapper {lvalue},GIMLI::MatrixBase*)
setData(data)[source]

Set a DataContainer.

setJacobian((object)arg1, (object)J) → object :
C++ signature :
void* setJacobian(GIMLI::ModellingBase {lvalue},GIMLI::MatrixBase*)

setJacobian( (object)arg1, (object)J) -> object :

C++ signature :
void* setJacobian(ModellingBase_wrapper {lvalue},GIMLI::MatrixBase*)
setMesh(mesh, ignoreRegionManager=False)[source]

Set Mesh.

setMultiThreadJacobian((object)arg1, (object)nThreads) → object :

Set number of threads used for brute force Jacobian generation. 1 is default. If nThreads is greater than 1 you need to implement response_mt with a read only response function. Maybe its worth set the single setThreadCount to 1 than, that you dont find yourself in a threading overkill.

C++ signature :
void* setMultiThreadJacobian(GIMLI::ModellingBase {lvalue},unsigned long)
setRegionManager((object)arg1, (object)reg) → object :
C++ signature :
void* setRegionManager(GIMLI::ModellingBase {lvalue},GIMLI::RegionManager*)
setStartModel((object)arg1, (object)startModel) → object :
C++ signature :
void* setStartModel(GIMLI::ModellingBase {lvalue},GIMLI::Vector<double>)

setStartModel( (object)arg1, (object)startModel) -> object :

C++ signature :
void* setStartModel(ModellingBase_wrapper {lvalue},GIMLI::Vector<double>)
setThreadCount((object)arg1, (object)nThreads) → object :

Set the maximum number of allowed threads for MT calculation. Have to be greater than 0. Will also set ENV(OPENBLAS_NUM_THREADS) .. if used.

C++ signature :
void* setThreadCount(GIMLI::ModellingBase {lvalue},unsigned long)
setVerbose((object)arg1, (object)verbose) → object :

Set verbose state.

C++ signature :
void* setVerbose(GIMLI::ModellingBase {lvalue},bool)
solution((object)arg1) → object :
C++ signature :
GIMLI::Matrix<double> solution(GIMLI::ModellingBase {lvalue})
startModel((object)arg1) → object :
C++ signature :
GIMLI::Vector<double> startModel(GIMLI::ModellingBase {lvalue})

startModel( (object)arg1) -> object :

C++ signature :
GIMLI::Vector<double> startModel(ModellingBase_wrapper {lvalue})
threadCount((object)arg1) → object :

Return the maximum number of allowed threads for MT calculation

C++ signature :
unsigned long threadCount(GIMLI::ModellingBase {lvalue})
uAnalytical(p, sourcePos, k)[source]

Calculates the analytical solution for the 2.5D geoelectrical problem.

Solves the 2.5D geoelectrical problem for one wave number k. It calculates the normalized (for injection current 1 A and sigma=1 S/m) potential at position p for a current injection at position sourcePos. Injection at the subsurface is recognized via mirror sources along the surface at depth=0.

Parameters:
p : pg.Pos

Position for the sought potential

sourcePos : pg.Pos

Current injection position.

k : float

Wave number

Returns:
u : float

Solution u(p)

verbose((object)arg1) → object :

Get verbose state.

C++ signature :
bool verbose(GIMLI::ModellingBase {lvalue})

verbose( (object)arg1) -> object :

C++ signature :
bool verbose(GIMLI::ModellingBase {lvalue})

FDEM

class pygimli.physics.FDEM(x=None, freqs=None, coilSpacing=None, inphase=None, outphase=None, filename=None, scaleFreeAir=False)[source]

Class for managing Frequency Domain EM data and their inversions.

Methods

FOP([nlay]) Forward modelling operator using a block discretization.
FOP2d(nlay) 2d forward modelling operator.
FOPsmooth(zvec) Forward modelling operator using fixed layers (smooth inversion)
datavec([xpos]) Extract data vector (stack in and out phase) for given pos/no.
deactivate(fr) Deactivate a single frequency.
error([xpos]) Return error as vector.
errorvec([xpos, minvalue]) Extract error vector for a give position or sounding number.
freq() Return active (i.e., non-deactivated) frequencies.
importEmsysAsciiData(filename) Import data from emsys text export file.
importMaxminData(filename[, verbose]) Import MaxMin IPX format with pos, data, frequencies & geometry.
inv2D(nlay[, lam, resL, resU, thkL, thkU, …]) 2d LCI inversion class.
invBlock([xpos, nlay, noise, stmod, lam, …]) Create and return Gimli inversion instance for block inversion.
plotAllData([orientation, aspect, outname, …]) Plot data along a profile as image plots for IP and OP.
plotData([xpos, response, error, ax, …]) Plot data as curves at given position.
plotDataOld([xpos, response, marker, …]) Plot data as curves at given position.
plotModelAndData(model, xpos, response[, …]) Plot both model and data in subfigures.
readHEMData(filename[, takeevery, choosevcp]) Read RESOLVE type airborne EM data from .XYZ file.
selectData([xpos]) Select sounding at a specific position or by number.
showModelAndData(model[, xpos, response, …]) Show both model and data with response in subfigures.
FOP(nlay=2)[source]

Forward modelling operator using a block discretization.

Parameters:
nlay : int

Number of blocks

FOP2d(nlay)[source]

2d forward modelling operator.

FOPsmooth(zvec)[source]

Forward modelling operator using fixed layers (smooth inversion)

Parameters:
zvec : array
__init__(x=None, freqs=None, coilSpacing=None, inphase=None, outphase=None, filename=None, scaleFreeAir=False)[source]

Initialize data class and load data. Provide filename or data.

If filename is given, data is loaded, overwriting settings.

Parameters:
x: array

Array of measurement positions

freq: array

Measured frequencies

coilSpacing : float

Distance between 2 two coils

inphase : array

real part of \(|amplitude| * \exp^{i phase}\)

outphase : array

imaginary part of \(|amplitude| * \exp^{i phase}\)

filename : str

Filename to read from. Supported: .xyz (MaxMin), *.txt (Emsys)

scaleFreeAir : bool

Scale inphase and outphase data by free air (primary) solution

datavec(xpos=0)[source]

Extract data vector (stack in and out phase) for given pos/no.

deactivate(fr)[source]

Deactivate a single frequency.

error(xpos=0)[source]

Return error as vector.

errorvec(xpos=0, minvalue=0.0)[source]

Extract error vector for a give position or sounding number.

freq()[source]

Return active (i.e., non-deactivated) frequencies.

importEmsysAsciiData(filename)[source]

Import data from emsys text export file.

columns: no, pos(1-3), separation(4), frequency(6), error(8), inphase (9-11), outphase (12-14), reads: positions, data, frequencies, error and geometry

importMaxminData(filename, verbose=False)[source]

Import MaxMin IPX format with pos, data, frequencies & geometry.

inv2D(nlay, lam=100.0, resL=1.0, resU=1000.0, thkL=1.0, thkU=100.0, minErr=1.0)[source]

2d LCI inversion class.

invBlock(xpos=0, nlay=2, noise=1.0, stmod=30.0, lam=100.0, lBound=0.0, uBound=0.0, verbose=False)[source]

Create and return Gimli inversion instance for block inversion.

Parameters:
xpos : array

position vector

nLay : int

Number of layers of the model to be determined OR vector of layer numbers OR forward operator

noise : float

Absolute data err in percent

stmod : float or pg.RVector

Starting model

lam : float

Global regularization parameter lambda.

lBound : float

Lower boundary for the model

uBound : float

Upper boundary for the model. 0 means no upper booundary

verbose : bool

Be verbose

plotAllData(orientation='horizontal', aspect=1000, outname=None, show=False, figsize=(11, 6), everyx=1)[source]

Plot data along a profile as image plots for IP and OP.

plotData(xpos=0, response=None, error=None, ax=None, marker='bo-', rmarker='rx-', clf=True, addlabel='', nv=2)[source]

Plot data as curves at given position.

plotDataOld(xpos=0, response=None, marker='bo-', rmarker='rx-', clf=True)[source]

Plot data as curves at given position.

plotModelAndData(model, xpos, response, modelL=None, modelU=None)[source]

Plot both model and data in subfigures.

readHEMData(filename, takeevery=1, choosevcp=True)[source]

Read RESOLVE type airborne EM data from .XYZ file.

selectData(xpos=0)[source]

Select sounding at a specific position or by number.

Retrieve inphase, outphase and error(if exist) vector from index or near given position

Returns:
IP : array OP : array ERR : array or None (if no error is specified)
showModelAndData(model, xpos=0, response=None, figsize=(8, 6))[source]

Show both model and data with response in subfigures.

MRS

class pygimli.physics.MRS(name=None, verbose=True, **kwargs)[source]

Magnetic resonance sounding (MRS) manager class.

Attributes:
t, q : ndarray - time and pulse moment vectors
data, error : 2d ndarray - data and error cubes
K, z : ndarray - (complex) kernel and its vertical discretization
model, modelL, modelU : vectors - model vector and lower/upper bound to it

Methods

loadMRSI - load MRSI (MRSmatlab format) data  
showCube - show any data/error/misfit as data cube (over q and t)  
showDataAndError - show data and error cubes  
showKernel - show Kernel matrix  
createFOP - create forward operator  
createInv - create pygimli Inversion instance  
run - run block-mono (alternatively smooth-mono) inversion (with bootstrap)  
calcMCM - compute model covariance matrix and thus uncertainties  
splitModel - return thickness, water content and T2* time from vector  
showResult/showResultAndFit - show inversion result (with fit)  
runEA - run evolutionary algorithm (GA, PSO etc.) using inspyred  
plotPopulation - plot final population of an EA run  
__init__(name=None, verbose=True, **kwargs)[source]

MRS init with optional data load from mrsi file

Parameters:
name : string

Filename with load data and kernel (.mrsi) or just data (.mrsd)

verbose : bool

be verbose

kwargs - see :func:`MRS.loadMRSI`.
calcMCM()[source]

Compute linear model covariance matrix.

calcMCMbounds()[source]

Compute model bounds using covariance matrix diagonals.

checkData(**kwargs)[source]

Check data and retrieve data and error vector.

static createFOP(nlay, K, z, t)[source]

Create forward operator instance.

createInv(nlay=3, lam=100.0, verbose=True, **kwargs)[source]

Create inversion instance (and fop if necessary with nlay).

genMod(individual)[source]

Generate (GA) model from random vector (0-1) using model bounds.

invert(nlay=3, lam=100.0, startvec=None, verbose=True, uncertainty=False, **kwargs)[source]

Easiest variant doing all (create fop and inv) in one call.

loadDataCube(filename='datacube.dat')[source]

Load data cube from single ascii file (old stuff)

loadDataNPZ(filename, **kwargs)[source]

Load data and kernel from numpy gzip packed file.

The npz file contains the fields: q, t, D, (E), z, K

loadDir(dirname)[source]

Load several standard files from dir (old Borkum stage).

loadErrorCube(filename='errorcube.dat')[source]

Load error cube from a single ascii file (old stuff).

loadKernel(name='')[source]

Load kernel matrix from mrsk or two bmat files.

loadKernelNPZ(filename, **kwargs)[source]

Load data and kernel from numpy gzip packed file.

The npz file contains the fields: q, t, D, (E), z, K

loadMRSD(filename, usereal=False, mint=0.0, maxt=2.0)[source]

Load mrsd (MRS data) file: not really used as in MRSD.

loadMRSI(filename, **kwargs)[source]

Load data, error and kernel from mrsi or mrsd file

Parameters:
usereal : bool [False]

use real parts (after data rotation) instead of amplitudes

mint/maxt : float [0.0/2.0]

minimum/maximum time to restrict time series

loadResult(filename)[source]

Load inversion result from column file.

loadZVector(filename='zkernel.vec')[source]

Load the kernel vertical discretisation (z) vector.

plotEAstatistics(fname=None)[source]

Plot EA statistics (best, worst, …) over time.

plotPopulation(maxfitness=None, fitratio=1.05, savefile=True)[source]

Plot fittest individuals (fitness<maxfitness) as 1d models

Parameters:
maxfitness : float

maximum fitness value (absolute) OR

fitratio : float [1.05]

maximum ratio to minimum fitness

result()[source]

Return block model results (thk, wc and T2 vectors).

run(verbose=True, uncertainty=False, **kwargs)[source]

Easiest variant doing all (create fop and inv) in one call.

runEA(nlay=None, eatype='GA', pop_size=100, num_gen=100, runs=1, mp_num_cpus=8, **kwargs)[source]

Run evolutionary algorithm using the inspyred library

Parameters:
nlay : int [taken from classic fop if not given]

number of layers

pop_size : int [100]

population size

num_gen : int [100]

number of generations

runs : int [pop_size*num_gen]

number of independent runs (with random population)

eatype : string [‘GA’]
algorithm, choose among:

‘GA’ - Genetic Algorithm [default] ‘SA’ - Simulated Annealing ‘DEA’ - Discrete Evolutionary Algorithm ‘PSO’ - Particle Swarm Optimization ‘ACS’ - Ant Colony Strategy ‘ES’ - Evolutionary Strategy

saveFigs(basename=None, extension='pdf')[source]

Save all figures to (pdf) files.

saveResult(filename)[source]

Save inversion result to column text file for later use.

setBoundaries()[source]

Set parameter boundaries for inversion.

showCube(ax=None, vec=None, islog=None, clim=None, clab=None)[source]

Plot any data (or response, error, misfit) cube nicely.

showDataAndError(figsize=(10, 8), show=False)[source]

Show data cube along with error cube.

showKernel(ax=None)[source]

Show the kernel as matrix (Q over z).

showResult(figsize=(10, 8), save='', fig=None, ax=None)[source]

Show theta(z) and T2*(z) (+uncertainties if there).

showResultAndFit(figsize=(12, 10), save='', plotmisfit=False, maxdep=0, clim=None)[source]

Show ec(z), T2*(z), data and model response.

static simulate(model, K, z, t)[source]

Do synthetic modelling.

splitModel(model=None)[source]

Split model vector into d, theta and T2*.

Refraction

class pygimli.physics.Refraction(data=None, verbose=True, debug=False, fatray=False, frequency=1000.0, **kwargs)[source]

Manager for refraction seismics (traveltime tomography)

TODO Document main members and use default MethodeManager interface e.g., self.inv, self.fop, self.paraDomain, self.mesh, self.data

Methods

apparentData() Convert data into apparent data.
checkData() Check data
createApparentData(data) Create apparent slowness for given data.
createArgParser([dataSuffix]) Create default argument parser.
createData(sensors, scheme) Create an empty data set.
createFOP([verbose, usefmm, fatray]) Create default forward operator for Traveltime modelling.
createFOP_([verbose]) Create forward operator working on refined mesh.
createInv(fop[, verbose, doSave]) Create default inversion instance for Traveltime inversion.
createInv_(fop[, verbose, dosave]) Create inversion instance, data- and model transformations.
createMesh([depth, quality, paraDX, …]) Create (inversion) mesh using createParaDomain2D
dataToken() Token name for the data in a DataContainer.
dataVals(data) Return pure data values from a given DataContainer.
drawApparentVelocities(ax, data[, t]) Plot data in for of apparent velocity image.
drawTravelTimeData(ax, data[, t]) Plot travel time data as lines and points.
estimateError([data, absoluteError, …]) Estimate error composed of an absolute and a relative part
getDepth() return a (a-priori guessed) depth of investigation
getMidpoint([data]) Return vector of offsets (in m) between shot and receiver.
getModel() Return velocity vector.
getOffset([data, full]) Return vector of offsets (in m) between shot and receiver.
invert([data, t, err, mesh]) Run actual inversion.
loadData(filename) Load data from file.
makeJacobianPDF([ind]) Make multipage Jacobian PDF.
model() Return the actual model.
paraDomain() Return parameter domain mesh.
rayCoverage() return ray coverage
relErrorVals(data) Return pure data values from a given DataContainer.
saveFigures([name, ext]) save all existing figures to files
saveResult([folder, size]) Save the results in the specified folder.
setData(data) Set data container (holding s and g indices and t floats).
setDataContainer(data) Set data container from outside.
setDataToken(token) Set the token name to identity the data in a DataContainer.
setMesh(mesh[, refine, secNodes]) Set mesh.
setVerbose(verbose) Make the class verbose (put output to the console)
show(data[, values, axes, cMin, cMax, colorBar]) Forward the visualization.
showCoverage([ax, name]) shows the ray coverage in logscale
showData([data, response, ax, name]) Show data as travel time curves (optionally with response)
showMesh([ax, name]) show mesh in given ax or in a new figure
showModel([ax, vals]) WRITEME
showRayPaths([model, ax]) Show ray paths for model or last model for which the Jacobian was calculated.
showResult([val, ax, cMin, cMax, logScale, …]) Show resulting velocity vector.
showResultAndFit([name]) show two vertical subplots with result and data (with response)
showVA([data, t, name, pseudosection, …]) Show apparent velocity as image plot.
simulate(mesh, slowness, scheme[, verbose]) Simulate a traveltime measurement.
standardizedCoverage() return standardized coverage vector (0|1) using neighbor info
useFMM([fmm]) Define whether to use Fast Marching Method (FMM).
useFatray([fatray, frequency]) Define whether to use Fatray jacobian computation.
__init__(data=None, verbose=True, debug=False, fatray=False, frequency=1000.0, **kwargs)[source]

Init function with optional data load

apparentData()

Convert data into apparent data.

checkData()[source]

Check data

w.r.t. shot/geophone identity and zero/negative traveltimes, plus check y/z sensor positions

createApparentData(data)[source]

Create apparent slowness for given data.

static createArgParser(dataSuffix='dat')

Create default argument parser.

Create default argument parser for the following options:

-Q, –quiet

-R, –robustData: options.robustData

-B, –blockyModel: options.blockyModel

-l, –lambda: options.lam

-i, –maxIter: options.maxIter

—depth: options.depth

createData(sensors, scheme)

Create an empty data set.

static createFOP(verbose=False, usefmm=False, fatray=False)[source]

Create default forward operator for Traveltime modelling.

usefmm forces Fast Marching Method, otherwise Dijkstra is used.

createFOP_(verbose=False)

Create forward operator working on refined mesh.

createInv(fop, verbose=True, doSave=False)[source]

Create default inversion instance for Traveltime inversion.

createInv_(fop, verbose=True, dosave=False)

Create inversion instance, data- and model transformations.

createMesh(depth=None, quality=34.3, paraDX=1, boundary=0, paraBoundary=0, secNodes=3, apply=True, **kwargs)[source]

Create (inversion) mesh using createParaDomain2D

Parameters:
depth : float, optional

maximum depth, 0 (default) means maximum offset / 3.

paraDX : float

relative distance for refinement nodes between two sensors e.g., 0 or 1 means no refinement e.g., 0.5 means 1 additional node between two neighboring sensors e.g., 0.33 means 2 additional equidistant nodes between two sensors

boundary : float, optional

boundary width to be appended for domain prolongation in absolute para domain width. values < 0 force the boundary to be 4 times para domain width.

paraBoundary : float, optional

margin for parameter domain in sensor distances (default 2)

quality : float, optional

mesh quality (smallest angle allowed)

apply : bool, optional

set mesh property of the underlying forward operator

secNodes : int (1)

Amount of secondary nodes to improve accuracy of the forward solution.

**kwargs: Additional keyword arguments passed to

pygimli.meshtools.createParaMeshPLC

dataToken()

Token name for the data in a DataContainer.

dataVals(data)[source]

Return pure data values from a given DataContainer.

static drawApparentVelocities(ax, data, t=None, **kwargs)[source]

Plot data in for of apparent velocity image.

static drawTravelTimeData(ax, data, t=None)[source]

Plot travel time data as lines and points.

static estimateError(data=None, absoluteError=0.001, relativeError=0.001)[source]

Estimate error composed of an absolute and a relative part

Parameters:
absoluteError : float

absolute error of traveltimes (usually in s)

relativeError : float

relative error of traveltimes in 1 (e.g. 0.01 is 1%)

Returns:
err : array
getDepth()[source]

return a (a-priori guessed) depth of investigation

getMidpoint(data=None)[source]

Return vector of offsets (in m) between shot and receiver.

getModel()[source]

Return velocity vector.

getOffset(data=None, full=False)[source]

Return vector of offsets (in m) between shot and receiver.

invert(data=None, t=None, err=None, mesh=None, **kwargs)[source]

Run actual inversion.

Values for result/response are stored in the class members velocity/response

Parameters:
useGradient : bool

Create gradient for starting model from vtop to vbottom.

vtop, vbottom : float

starting (gradient) model velocities on top/at bottom of the mesh

lam : float

regularization parameter describing the strength of smoothness

zWeight : float

relative weight for purely vertical boundaries

maxIter : int

Maximum number of iterations

startModel : array

Slowness starting model for the inversion

loadData(filename)[source]

Load data from file.

makeJacobianPDF(ind=None)[source]

Make multipage Jacobian PDF.

model()

Return the actual model.

paraDomain()[source]

Return parameter domain mesh.

rayCoverage()[source]

return ray coverage

relErrorVals(data)[source]

Return pure data values from a given DataContainer.

saveFigures(name=None, ext='pdf')[source]

save all existing figures to files

saveResult(folder=None, size=(16, 10), **kwargs)[source]

Save the results in the specified folder.

Saved items are:
Inverted profile Velocity vector Coverage vector Standardized coverage vector Mesh (bms and vtk with results)
setData(data)[source]

Set data container (holding s and g indices and t floats).

setDataContainer(data)[source]

Set data container from outside.

setDataToken(token)

Set the token name to identity the data in a DataContainer.

setMesh(mesh, refine=False, secNodes=1)[source]

Set mesh. To be removed from class once derived from MeshManager.

Parameters:
secNodes : int (1)

Number of secondary nodes to improve accuracy of the forward solution.

setVerbose(verbose)

Make the class verbose (put output to the console)

show(data, values=None, axes=None, cMin=None, cMax=None, colorBar=1, **kwargs)

Forward the visualization.

showCoverage(ax=None, name='coverage', **kwargs)[source]

shows the ray coverage in logscale

showData(data=None, response=None, ax=None, name='data')[source]

Show data as travel time curves (optionally with response)

Parameters:
data : pyGIMLi data Container [self.dataContainer]

data to show with points

response : array

response vector to draw with lines

ax : maxplotlib axes

axis to plot into, if not given, a new figure is created

showMesh(ax=None, name='mesh')[source]

show mesh in given ax or in a new figure

showModel(ax=None, vals=None, **kwargs)[source]

WRITEME

showRayPaths(model=None, ax=None, **kwargs)[source]

Show ray paths for model or last model for which the Jacobian was calculated.

Parameters:
model : array

Velocity model for which to calculate and visualize ray paths (the default is model for last Jacobian calculation in self.velocity).

ax : matplotlib.axes

Axes for the plot (the default is None).

**kwargs : type

Additional arguments passed to LineCollection (alpha, linewidths, color, linestyles).

Returns:
ax : matplotlib.axes object
cb : matplotlib.colorbar object (only if model is provided)

Examples

>>> # No reason to import matplotlib
>>> import pygimli as pg
>>> from pygimli.physics import Refraction
>>> from pygimli.physics.traveltime import createRAData
>>>
>>> x, y = 8, 6
>>> mesh = pg.createGrid(x, y)
>>> data = createRAData([(0,0)] + [(x, i) for i in range(y)], shotdistance=y+1)
>>> data.set("t", pg.RVector(data.size(), 1.0))
>>> rst = Refraction()
>>> rst.setDataContainer(data)
Data: Sensors: 7 data: 6
>>> rst.setMesh(mesh, 5)
>>> ax, cb = rst.showRayPaths()

(png, pdf)

../../_images/pygimli-physics-1.png
showResult(val=None, ax=None, cMin=None, cMax=None, logScale=False, rays=False, name='result', **kwargs)[source]

Show resulting velocity vector.

Parameters:
val : result array [self.velocity]

field to show, usually the velocity vector

ax : matplotlib.axes

axes to plot into, if not give a new one-axis figure is created

cMin/cMax : float

minimum and maximum values for ranging colorbar

logScale : bool [False]

use logarithmic scale

rays : bool [False]

Show ray paths as well.

Returns:
ax : maxplotlib axes
cb : matplotlib color bar object
Other Parameters:
 
useCoverage : bool [True]

use standardized (0 or 1) ray coverage as alpha-shading

label : str

label to write on colorbar

orientaton : str

color bar orientation

nLevs : int [7]

number of level values for color bar

**kwargs : keyword arguments passed to the show function
showResultAndFit(name='resultfit', **kwargs)[source]

show two vertical subplots with result and data (with response)

showVA(data=None, t=None, name='va', pseudosection=False, squeeze=True, full=True, ax=None, cmap=None, **kwargs)[source]

Show apparent velocity as image plot.

TODO showXXX commands need to return ax and cbar .. if there is one

static simulate(mesh, slowness, scheme, verbose=False, **kwargs)[source]

Simulate a traveltime measurement.

Perform the forward task for a given mesh, a slowness distribution (per cell) and return data (traveltime) for a measurement scheme. This is a static method since it does not interfere with the managers inversion approaches.

Parameters:
mesh : GIMLI::Mesh

Mesh to calculate for.

slowness : array(mesh.cellCount()) | array(N, mesh.cellCount())

slowness distribution for the given mesh cells can be:

  • a single array of len mesh.cellCount()
  • a matrix of N slowness distributions of len mesh.cellCount()
  • a res map as [[marker0, res0], [marker1, res1], …]
scheme : GIMLI::DataContainer

data measurement scheme

verbose : boolean

Be verbose.

Returns:
t : array(N, data.size()) | DataContainer

The resulting simulated travel time values. Either one column array or matrix in case of slowness matrix. A DataContainer is return if noisify set to True.

Other Parameters:
 
noisify : boolean

add normal distributed noise based on scheme(‘err’)

standardizedCoverage()[source]

return standardized coverage vector (0|1) using neighbor info

useFMM(fmm=True)[source]

Define whether to use Fast Marching Method (FMM).

Note that this method is more accurate but currently a lot slower!

useFatray(fatray=True, frequency=300.0)[source]

Define whether to use Fatray jacobian computation.

SIPSpectrum

class pygimli.physics.SIPSpectrum(filename=None, unify=False, onlydown=True, f=None, amp=None, phi=None, k=1, sort=True, basename='new')[source]

SIP spectrum data analysis.

Methods

checkCRKK([useEps, use0, ax]) Check coupling removal (CR) by Kramers-Kronig (KK) relation
cutF([fcut, down]) Cut (delete) frequencies above a certain value fcut.
determineEpsilon([mode, sigmaR, sigmaI]) Retrieve frequency-independent epsilon for f->Inf.
epsilonR() Calculate relative permittivity from imaginary conductivity
fit2CCPhi([ePhi, lam, mpar, taupar1, …]) fit a Cole-Cole term to phase only
fitCCEM([ePhi, lam, remove, mpar, taupar, …]) Fit a Cole-Cole term with additional EM term to phase
fitCCPhi([ePhi, lam, mpar, taupar, cpar]) fit a Cole-Cole term to phase only
fitColeCole([useCond]) Fit a Cole-Cole model to the phase data
fitDebyeModel([ePhi, lam, lamFactor, mint, …]) Fit a (smooth) continuous Debye model (Debye decomposition).
getKK([use0]) Compute Kramers-Kronig impedance values (re->im and im->re).
getPhiKK([use0]) Compute phase from Kramers-Kronig quantities.
loadData(filename, **kwargs) Import spectral data.
logMeanTau() Mean logarithmic relaxation time (50% cumulative log curve).
omega() Angular frequency.
realimag([cond]) Real and imaginary part of resistivity/conductivity (cond=True).
removeEpsilonEffect([er, mode]) remove effect of (constant high-frequency) epsilon from sigma
saveFigures([name, ext]) Save all existing figures to files using file basename.
showAll([save]) Plot spectrum, Cole-Cole fit and Debye distribution
showData([reim, znorm, cond, nrows, ax]) Show amplitude and phase spectrum in two subplots
showDataKK([use0]) Show data as real/imag subplots along with Kramers-Kronig curves
showPhase([ax]) Plot phase spectrum (-phi over frequency).
showPolarPlot([cond]) Show data in a polar plot (imaginary vs.
sortData() Sort data along increasing frequency (e.g.
totalChargeability() Total chargeability (sum) from Debye curve.
unifyData([onlydown]) Unify data (only one value per frequency) by mean or selection.
zNorm() Normalized real (difference) and imag.
__init__(filename=None, unify=False, onlydown=True, f=None, amp=None, phi=None, k=1, sort=True, basename='new')[source]

Init SIP class with either filename to read or data vectors.

Examples

>>> #sip = SIPSpectrum('sipexample.txt', unify=True) # unique f values
>>> #sip = SIPSpectrum(f=f, amp=R, phi=phase, basename='new')
checkCRKK(useEps=False, use0=False, ax=None)[source]

Check coupling removal (CR) by Kramers-Kronig (KK) relation

cutF(fcut=1e+99, down=False)[source]

Cut (delete) frequencies above a certain value fcut.

determineEpsilon(mode=0, sigmaR=None, sigmaI=None)[source]

Retrieve frequency-independent epsilon for f->Inf.

Parameters:
mode : int
Operation mode:

=0 - extrapolate using two highest frequencies (default) <0 - take last -n frequencies >0 - take n-th frequency

sigmaR/sigmaI : float

real and imaginary conductivity (if not given take data)

Returns
——-
er : float

relative permittivity (epsilon) value (dimensionless)

epsilonR()[source]

Calculate relative permittivity from imaginary conductivity

fit2CCPhi(ePhi=0.001, lam=1000.0, mpar=(0, 0, 1), taupar1=(0, 1e-05, 1), taupar2=(0, 0.1, 1000), cpar=(0.5, 0, 1))[source]

fit a Cole-Cole term to phase only

Parameters:
ePhi : float

absolute error of phase angle

lam : float

regularization parameter

mpar, taupar, cpar : list[3]

inversion parameters (starting value, lower bound, upper bound) for Cole-Cole parameters (m, tau, c) and EM relaxation time (em)

fitCCEM(ePhi=0.001, lam=1000.0, remove=True, mpar=(0.2, 0, 1), taupar=(0.01, 1e-05, 100), cpar=(0.25, 0, 1), empar=(1e-07, 1e-09, 1e-05))[source]

Fit a Cole-Cole term with additional EM term to phase

Parameters:
ePhi : float

absolute error of phase angle

lam : float

regularization parameter

remove: bool

remove EM term from data

mpar, taupar, cpar, empar : list[3]

inversion parameters (starting value, lower bound, upper bound) for Cole-Cole parameters (m, tau, c) and EM relaxation time (em)

fitCCPhi(ePhi=0.001, lam=1000.0, mpar=(0, 0, 1), taupar=(0, 1e-05, 100), cpar=(0.3, 0, 1))[source]

fit a Cole-Cole term to phase only

Parameters:
ePhi : float

absolute error of phase angle

lam : float

regularization parameter

mpar, taupar, cpar : list[3]

inversion parameters (starting value, lower bound, upper bound) for Cole-Cole parameters (m, tau, c) and EM relaxation time (em)

fitColeCole(useCond=False, **kwargs)[source]

Fit a Cole-Cole model to the phase data

Parameters:
useCond : bool

use conductivity form of Cole-Cole model instead of resistivity

error : float [0.01]

absolute phase error

lam : float [1000]

initial regularization parameter

mpar : tuple/list (0, 0, 1)

inversion parameters for chargeability: start, lower, upper bound

taupar : tuple/list (1e-2, 1e-5, 100)

inversion parameters for time constant: start, lower, upper bound

cpar : tuple/list (0.25, 0, 1)

inversion parameters for Cole exponent: start, lower, upper bound

fitDebyeModel(ePhi=0.001, lam=1000.0, lamFactor=0.8, mint=None, maxt=None, nt=None, new=True, showFit=False, cType=1)[source]

Fit a (smooth) continuous Debye model (Debye decomposition).

Parameters:
ePhi : float

absolute error of phase angle

lam : float

regularization parameter

lamFactor : float

regularization factor for subsequent iterations

mint/maxt : float

minimum/maximum tau values to use (else automatically from f)

nt : int

number of tau values (default number of frequencies * 2)

new : bool

new implementation (experimental)

showFit : bool

show fit

cType : int

constraint type (1/2=smoothness 1st/2nd order, 0=minimum norm)

getKK(use0=False)[source]

Compute Kramers-Kronig impedance values (re->im and im->re).

getPhiKK(use0=False)[source]

Compute phase from Kramers-Kronig quantities.

loadData(filename, **kwargs)[source]

Import spectral data.

Import Data and try to assume the file format.

logMeanTau()[source]

Mean logarithmic relaxation time (50% cumulative log curve).

omega()[source]

Angular frequency.

realimag(cond=False)[source]

Real and imaginary part of resistivity/conductivity (cond=True).

removeEpsilonEffect(er=None, mode=0)[source]

remove effect of (constant high-frequency) epsilon from sigma

Parameters:
er : float

relative epsilon to correct for (else automatically determined)

mode : int

automatic epsilon determination mode (see determineEpsilon)

Returns:
er : float

determined permittivity (see determineEpsilon)

saveFigures(name=None, ext='pdf')[source]

Save all existing figures to files using file basename.

showAll(save=False)[source]

Plot spectrum, Cole-Cole fit and Debye distribution

showData(reim=False, znorm=False, cond=False, nrows=2, ax=None, **kwargs)[source]

Show amplitude and phase spectrum in two subplots

Parameters:
reim : bool

show real/imaginary part instead of amplitude/phase

znorm : bool (true forces reim)

use normalized real/imag parts after Nordsiek&Weller (2008)

nrows - use nrows subplots (default=2)
Returns:
fig, ax : mpl.figure, mpl.axes array
showDataKK(use0=False)[source]

Show data as real/imag subplots along with Kramers-Kronig curves

showPhase(ax=None, **kwargs)[source]

Plot phase spectrum (-phi over frequency).

showPolarPlot(cond=False)[source]

Show data in a polar plot (imaginary vs. real parts).

sortData()[source]

Sort data along increasing frequency (e.g. useful for KK).

totalChargeability()[source]

Total chargeability (sum) from Debye curve.

unifyData(onlydown=False)[source]

Unify data (only one value per frequency) by mean or selection.

zNorm()[source]

Normalized real (difference) and imag. z [NW08]

TDEM

class pygimli.physics.TDEM(filename=None)[source]

TEM class mainly for holding data etc.

Methods

__call__([i]) Return a single sounding.
getFOP([nr]) Return forward operator.
invert([nr, nlay, thickness]) Do inversion.
load(filename) Road data from usf, txt (siroTEM), tem (TEMfast) or UniK file.
plotRhoa([ax, ploterror, corrramp]) Plot all apparent resistivity curves into one window.
plotTransients([ax]) Plot all transients into one window
stackAll([tmin, tmax]) Stack all measurements yielding a new TDEM class instance.
showInfos  
__init__(filename=None)[source]

Initialize class and (optionally) load data

getFOP(nr=0)[source]

Return forward operator.

invert(nr=0, nlay=4, thickness=None)[source]

Do inversion.

load(filename)[source]

Road data from usf, txt (siroTEM), tem (TEMfast) or UniK file.

plotRhoa(ax=None, ploterror=False, corrramp=False, **kwargs)[source]

Plot all apparent resistivity curves into one window.

plotTransients(ax=None, **kwargs)[source]

Plot all transients into one window

showInfos()[source]
stackAll(tmin=0, tmax=100)[source]

Stack all measurements yielding a new TDEM class instance.

VESManager

class pygimli.physics.VESManager(**kwargs)[source]

Vertical electrical sounding (VES) manager class.

Examples

>>> # no need to import matplotlib. pygimli's show does
>>> import numpy as np
>>> import pygimli as pg
>>> from pygimli.physics import VESManager
>>> ab2 = np.logspace(np.log10(1.5), np.log10(100), 32)
>>> mn2 = 1.0
>>> # 3 layer with 100, 500 and 20 Ohmm
>>> # and layer thickness of 4, 6, 10 m
>>> # over a Halfspace of 800 Ohmm
>>> synthModel = pg.cat([4., 6., 10.], [100., 5., 20., 800.])
>>> VES = VESManager()
>>> ra, err = VES.simulate(synthModel, ab2=ab2, mn2=mn2, noiseLevel=0.01)
>>> ax = VES.showData(ra, err)
>>> # _= VES.invert(ra, err, nLayer=4, showProgress=0, verbose=0)
>>> # ax = VES.showModel(synthModel)
>>> # ax = VES.showResult(ax=ax)
>>> pg.wait()

(png, pdf)

../../_images/pygimli-physics-2.png
Attributes:
complex
debug
verbose

Methods

createArgParser([dataSuffix]) Create default argument parser.
createForwardOperator(**kwargs) Create Forward Operator.
estimateError(data[, errLevel]) Estimate data error.
exportData(fileName[, data, error]) Export data into simple ascii matrix.
initForwardOperator(**kwargs) Initialize or re-initialize the forward operator.
initInversionFramework(**kwargs) Initialize or re-initialize the inversion framework.
invert([data, err, ab2, mn2]) Invert measured data.
loadData(fileName, **kwargs) Mandatory interface for derived classes.
setVerbose(verbose) Make the class verbose (put output to the console)
showData(data[, error, ax]) Shows the data.
showFit([ax]) Show the last inversion date and response.
showModel(model[, ax]) Shows a model.
showResult([ax]) Show the last inversion result.
showResultAndFit() Calls showResults and showFit.
simulate(model[, ab2, mn2]) Simulate measurement data.
createInversionFramework  
__init__(**kwargs)[source]

Constructor

Parameters:
complex : bool

Accept complex resistivities.

Attributes:
complex : bool

Accept complex resistivities.

complex
static createArgParser(dataSuffix='dat')

Create default argument parser.

Create default argument parser for the following options:

-Q, –quiet

-R, –robustData: options.robustData

-B, –blockyModel: options.blockyModel

-l, –lambda: options.lam

-i, –maxIter: options.maxIter

—depth: options.depth

createForwardOperator(**kwargs)[source]

Create Forward Operator.

Create Forward Operator based on complex attribute.

createInversionFramework(**kwargs)
debug
estimateError(data, errLevel=0.01)

Estimate data error.

Create an error of estimated measurement error. On default it returns an array of constant relative errors. More sophisticated error estimation should be done in specialized derived classes.

Parameters:
data : iterable

Data values for which the errors should be estimated.

errLevel : float (0.01)

Error level in percent/100.

Returns:
err : array

Returning array of size len(data)

exportData(fileName, data=None, error=None)[source]

Export data into simple ascii matrix.

Usefull?

initForwardOperator(**kwargs)

Initialize or re-initialize the forward operator.

Called once in the constructor to force the manager to create the necessary forward operator member. Can be recalled if you need to changed the mangers own forward operator object. If you want a own instance of a valid FOP call createForwardOperator.

initInversionFramework(**kwargs)

Initialize or re-initialize the inversion framework.

Called once in the constructor to force the manager to create the necessary Framework instance.

invert(data=None, err=None, ab2=None, mn2=None, **kwargs)[source]

Invert measured data.

loadData(fileName, **kwargs)[source]

Mandatory interface for derived classes.

setVerbose(verbose)

Make the class verbose (put output to the console)

showData(data, error=None, ax=None, **kwargs)

Shows the data.

Draw data values into a given axes or show the data values from the last run. Forwards on default to the self.fop.drawData function of the modelling operator. If there is no given function given, you have to override this method.

Parameters:
ax : mpl axes

Axes object to draw into. Create a new if its not given.

data : iterable

Data values to be draw.

error : iterable

Data error values to be draw.

showFit(ax=None)

Show the last inversion date and response.

showModel(model, ax=None, **kwargs)

Shows a model.

Draw model date into a given axes or show the inversion result from the last run. Forwards on default to the self.fop.drawModel function of the modelling operator. If there is no given function given, you have to override this method.

Parameters:
ax : mpl axes

Axes object to draw into. Create a new if its not given.

model : iterable

Model data to be draw.

showResult(ax=None, **kwargs)

Show the last inversion result.

showResultAndFit()

Calls showResults and showFit.

simulate(model, ab2=None, mn2=None, **kwargs)[source]

Simulate measurement data.

verbose

constants

pygimli.physics.constants

alias of pygimli.physics.constants.Constants



2018 - GIMLi Development Team
Created using Bootstrap, Sphinx and pyGIMLi 1.0.9+87.g68804698 on Dec 16, 2018.