Test notebook

In [ ]:

import pygimli as pg
import numpy as np
import pybert as pb
import matplotlib.pyplot as plt
import os.path

from pygimli.physics import ERTManager

import pygimli as pg import numpy as np import pybert as pb import matplotlib.pyplot as plt import os.path from pygimli.physics import ERTManager

In [ ]:

infile = 'data/LSBB03_181114_WS.bin'
topofile = 'topo/LSBB03_topo.txt'

infile = 'data/LSBB03_181114_WS.bin' topofile = 'topo/LSBB03_topo.txt'

In [ ]:

id_z = 1 # column with elevation values (x values are in column 0 by default)
delimiter = ',' # delimiter in topo file

id_z = 1 # column with elevation values (x values are in column 0 by default) delimiter = ',' # delimiter in topo file

In [ ]:

relativeError = 0.01; # in proportion (also minimum error when using reciprocals and measurement error)
absoluteUError = 5e-4; # in V

relativeError = 0.01; # in proportion (also minimum error when using reciprocals and measurement error) absoluteUError = 5e-4; # in V

In [ ]:

(dirName, fileName) = os.path.split(infile)
(fileBaseName, fileExtension)=os.path.splitext(fileName)

(dirName, fileName) = os.path.split(infile) (fileBaseName, fileExtension)=os.path.splitext(fileName)

In [ ]:

data = pb.importData(infile)

data = pb.importData(infile)

In [ ]:

k = pb.geometricFactors(data)
data.set('k',k)

k = pb.geometricFactors(data) data.set('k',k)

In [ ]:

if any(data('u')) and any(data('i')): 
    data.set('rhoa',(data('u')/data('i'))*data('k'))

if any(data('u')) and any(data('i')): data.set('rhoa',(data('u')/data('i'))*data('k'))

In [ ]:

data.markInvalid(data('rhoa')<=0)

data.markInvalid(data('rhoa')<=0)

In [ ]:

fig, ax = plt.subplots()
pg.show(data, cMap='Spectral_r', ax=ax, label='Apparent Resistivity ($\Omega$.m)')
ax.set_xlabel('Distance (m)')
ax.xaxis.set_tick_params(labeltop=True,labelbottom=False)
ax.xaxis.set_label_position('top') 
ax.set_ylabel('Pseudo-depth (m)')
fig.savefig(fname=fileBaseName + '_ApparentResistivity.png', dpi=150)

fig, ax = plt.subplots() pg.show(data, cMap='Spectral_r', ax=ax, label='Apparent Resistivity ($\Omega$.m)') ax.set_xlabel('Distance (m)') ax.xaxis.set_tick_params(labeltop=True,labelbottom=False) ax.xaxis.set_label_position('top') ax.set_ylabel('Pseudo-depth (m)') fig.savefig(fname=fileBaseName + '_ApparentResistivity.png', dpi=150)

In [ ]:

x_sensors = np.array(pg.x(data.sensorPositions())) # Read sensor positions
dx_mean = np.mean(np.diff(x_sensors))

x_sensors = np.array(pg.x(data.sensorPositions())) # Read sensor positions dx_mean = np.mean(np.diff(x_sensors))

In [ ]:

if 'topofile' in locals():
    topo = np.loadtxt(topofile,delimiter=delimiter)
    X_topo = topo[:,0]
    Z_topo = topo[:,id_z]
else:
    X_topo = x_sensors
    Z_topo = x_sensors*0

# Check if Z interpolation is required    
if np.equal(x_sensors.shape,X_topo.shape) and all(x_sensors == X_topo):
    x_interp = X_topo
    z_interp = Z_topo
else:
    A = np.array([X_topo, Z_topo])
    tt = np.hstack((0., np.cumsum(np.sqrt(np.diff(A[0])**2 + np.diff(A[1])**2))))
    z_interp  = np.interp(x_sensors, A[0], A[1])
    x_interp  = x_sensors[0] + np.hstack((0., np.cumsum(np.sqrt(np.diff(x_sensors)**2 - np.diff(z_interp)**2))))

# Plot topography   
# fig, ax = plt.subplots()
# ax.plot(x_interp,z_interp,'.-')
# ax.set_aspect(2)
# plt.xlabel('X (m)')
# plt.ylabel('Z (m)')

if 'topofile' in locals(): topo = np.loadtxt(topofile,delimiter=delimiter) X_topo = topo[:,0] Z_topo = topo[:,id_z] else: X_topo = x_sensors Z_topo = x_sensors*0 # Check if Z interpolation is required if np.equal(x_sensors.shape,X_topo.shape) and all(x_sensors == X_topo): x_interp = X_topo z_interp = Z_topo else: A = np.array([X_topo, Z_topo]) tt = np.hstack((0., np.cumsum(np.sqrt(np.diff(A[0])**2 + np.diff(A[1])**2)))) z_interp = np.interp(x_sensors, A[0], A[1]) x_interp = x_sensors[0] + np.hstack((0., np.cumsum(np.sqrt(np.diff(x_sensors)**2 - np.diff(z_interp)**2)))) # Plot topography # fig, ax = plt.subplots() # ax.plot(x_interp,z_interp,'.-') # ax.set_aspect(2) # plt.xlabel('X (m)') # plt.ylabel('Z (m)')

In [ ]:

for i in range(data.sensorCount()):
    data.setSensorPosition(i, [x_interp[i], 0.0, z_interp[i]])

for i in range(data.sensorCount()): data.setSensorPosition(i, [x_interp[i], 0.0, z_interp[i]])

In [ ]:

ert = ERTManager(data)
error = ert.estimateError(data, relativeError=relativeError, absoluteError=0, absoluteUError=absoluteUError)
ert.data.set('err',error)
ert.data.removeInvalid()

# Create mesh and starting model
ert.createMesh(paraDX=0.33, maxCellArea=dx_mean*2.5*dx_mean*2.5)
res = ert.invert(verbose=False)
print(ert.inv.chi2())
print(np.shape(ert.inv.chi2History))

ert = ERTManager(data) error = ert.estimateError(data, relativeError=relativeError, absoluteError=0, absoluteUError=absoluteUError) ert.data.set('err',error) ert.data.removeInvalid() # Create mesh and starting model ert.createMesh(paraDX=0.33, maxCellArea=dx_mean*2.5*dx_mean*2.5) res = ert.invert(verbose=False) print(ert.inv.chi2()) print(np.shape(ert.inv.chi2History))

In [ ]:

fig, ax = plt.subplots()
ert.showResult(ax=ax,cMap='Spectral_r',  label='Resistivity ($\Omega$.m)', useCoverage=True, cMin=np.min(ert.data('rhoa')), cMax=np.max(ert.data('rhoa')))
ax.set_xlabel('Distance (m)')
ax.xaxis.set_tick_params(labeltop=True,labelbottom=False)
ax.xaxis.set_label_position('top') 
ax.set_ylabel('Elevation (m)')
fig.savefig(fname=fileBaseName + '_DefaultInversion.png', dpi=150)

fig, ax = plt.subplots() ert.showResult(ax=ax,cMap='Spectral_r', label='Resistivity ($\Omega$.m)', useCoverage=True, cMin=np.min(ert.data('rhoa')), cMax=np.max(ert.data('rhoa'))) ax.set_xlabel('Distance (m)') ax.xaxis.set_tick_params(labeltop=True,labelbottom=False) ax.xaxis.set_label_position('top') ax.set_ylabel('Elevation (m)') fig.savefig(fname=fileBaseName + '_DefaultInversion.png', dpi=150)

In [ ]:

fig, ax = plt.subplots()
_,cb = pg.show(ax=ax,mesh=ert.paraDomain,data=ert.model,
                # coverage=ert.standardizedCoverage(),
                cMin=np.min(ert.data('rhoa')), cMax=np.max(ert.data('rhoa')),
                cMap='Spectral_r',logScale=True,
                label='resistivity (ohm.m)',orientation='vertical',
                xlabel='distance (m)',ylabel='elevation (m)')
fig.savefig(fname=fileBaseName + '_DefaultInversion.png', dpi=150)

fig, ax = plt.subplots() _,cb = pg.show(ax=ax,mesh=ert.paraDomain,data=ert.model, # coverage=ert.standardizedCoverage(), cMin=np.min(ert.data('rhoa')), cMax=np.max(ert.data('rhoa')), cMap='Spectral_r',logScale=True, label='resistivity (ohm.m)',orientation='vertical', xlabel='distance (m)',ylabel='elevation (m)') fig.savefig(fname=fileBaseName + '_DefaultInversion.png', dpi=150)

In [ ]:

m = pg.Mesh(ert.paraDomain)
m['Resistivity'] = ert.paraModel(ert.model)
m['Coverage'] = ert.coverage()
m['S_Coverage'] = ert.standardizedCoverage()
m.exportVTK(fileBaseName + '_DefaultInversion.vtk')

m = pg.Mesh(ert.paraDomain) m['Resistivity'] = ert.paraModel(ert.model) m['Coverage'] = ert.coverage() m['S_Coverage'] = ert.standardizedCoverage() m.exportVTK(fileBaseName + '_DefaultInversion.vtk')

In [ ]:

plt.close('all')

plt.close('all')