Inversion with topography
In this notebook, we will learn how to add topography and perform 2D inversion with rectangular and triangular meshes. The files needed can be found in examples/dc-2d-topo/
.
[1]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore') # just to make it cleaner in the notebook
import os
import sys
sys.path.append((os.path.relpath('../src'))) # add here the relative path of the API folder
testdir = '../src/examples/dc-2d-topo/'
import numpy as np # this will be used to read the topography file
from resipy import Project
API path = /media/jkl/data/phd/resipy/src/resipy
ResIPy version = 3.4.6
cR2.exe found and up to date.
R3t.exe found and up to date.
cR3t.exe found and up to date.
First, let’s create an Project
object and import the survey as usual.
[2]:
k = Project(typ='R2') # create new Project object and use default working directory
k.createSurvey(os.path.join(testdir, 'syscal.csv'))
Working directory is: /media/jkl/data/phd/resipy/src/resipy
clearing dirname
600/636 reciprocal measurements found.
10 measurements error > 20 %
We can also plot the pseudo-section. Note that this one remains flat event when there is topogragraphy involved.
[3]:
k.showPseudo()
If we take a look at the electrodes position by plotting R2.elec
we can see that there is no topography yet (all the values in the third column (=z) are 0.0).
[4]:
k.elec
[4]:
x | y | z | remote | buried | label | |
---|---|---|---|---|---|---|
0 | 0.00 | 0.0 | 0.0 | False | False | 1 |
1 | 0.25 | 0.0 | 0.0 | False | False | 2 |
2 | 0.50 | 0.0 | 0.0 | False | False | 3 |
3 | 0.75 | 0.0 | 0.0 | False | False | 4 |
4 | 1.00 | 0.0 | 0.0 | False | False | 5 |
5 | 1.25 | 0.0 | 0.0 | False | False | 6 |
6 | 1.50 | 0.0 | 0.0 | False | False | 7 |
7 | 1.75 | 0.0 | 0.0 | False | False | 8 |
8 | 2.00 | 0.0 | 0.0 | False | False | 9 |
9 | 2.25 | 0.0 | 0.0 | False | False | 10 |
10 | 2.50 | 0.0 | 0.0 | False | False | 11 |
11 | 2.75 | 0.0 | 0.0 | False | False | 12 |
12 | 3.00 | 0.0 | 0.0 | False | False | 13 |
13 | 3.25 | 0.0 | 0.0 | False | False | 14 |
14 | 3.50 | 0.0 | 0.0 | False | False | 15 |
15 | 3.75 | 0.0 | 0.0 | False | False | 16 |
16 | 4.00 | 0.0 | 0.0 | False | False | 17 |
17 | 4.25 | 0.0 | 0.0 | False | False | 18 |
18 | 4.50 | 0.0 | 0.0 | False | False | 19 |
19 | 4.75 | 0.0 | 0.0 | False | False | 20 |
20 | 5.00 | 0.0 | 0.0 | False | False | 21 |
21 | 5.25 | 0.0 | 0.0 | False | False | 22 |
22 | 5.50 | 0.0 | 0.0 | False | False | 23 |
23 | 5.75 | 0.0 | 0.0 | False | False | 24 |
Then we can load a csv file that will add topography for each electrode. The csv file should be of the form X,Y,Z and no headers are needed.
[5]:
k.importElec(os.path.join(testdir + 'elec.csv'))
print(k.elec)
label x y z remote buried
0 1 0.00 0.0 29.499 False False
1 2 0.25 0.0 29.504 False False
2 3 0.50 0.0 29.509 False False
3 4 0.75 0.0 29.516 False False
4 5 1.00 0.0 29.478 False False
5 6 1.25 0.0 29.461 False False
6 7 1.50 0.0 29.454 False False
7 8 1.75 0.0 29.428 False False
8 9 2.00 0.0 29.416 False False
9 10 2.25 0.0 29.411 False False
10 11 2.50 0.0 29.398 False False
11 12 2.75 0.0 29.362 False False
12 13 3.00 0.0 29.329 False False
13 14 3.25 0.0 29.245 False False
14 15 3.50 0.0 29.159 False False
15 16 3.75 0.0 29.083 False False
16 17 4.00 0.0 29.010 False False
17 18 4.25 0.0 28.929 False False
18 19 4.50 0.0 28.872 False False
19 20 4.75 0.0 28.761 False False
20 21 5.00 0.0 28.672 False False
21 22 5.25 0.0 28.593 False False
22 23 5.50 0.0 28.509 False False
23 24 5.75 0.0 28.412 False False
We can now create a mesh either triangular or rectangular.
[6]:
k.createMesh(typ='trian') # or trian for triangular mesh
k.showMesh()
Creating triangular mesh...done (2236 elements)
We can finally invert the data on this mesh. Note that it might take a while so be patient.
[7]:
k.invert()
Writing .in file and protocol.dat... done
--------------------- MAIN INVERSION ------------------
>> R 2 R e s i s t i v i t y I n v e r s i o n v4.10 <<
>> D a t e : 03 - 12 - 2023
>> My beautiful survey
>> I n v e r s e S o l u t i o n S e l e c t e d <<
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading start resistivity from res0.dat
>> R e g u l a r i s e d T y p e <<
>> L i n e a r F i l t e r <<
>> L o g - D a t a I n v e r s i o n <<
>> N o r m a l R e g u l a r i s a t i o n <<
>> D a t a w e i g h t s w i l l b e m o d i f i e d <<
Processing dataset 1
Measurements read: 336 Measurements rejected: 0
Geometric mean of apparent resistivities: 0.35349E+03
>> Total Memory required is: 0.007 Gb
Iteration 1
Initial RMS Misfit: 33.54 Number of data ignored: 0
Alpha: 1253.313 RMS Misfit: 4.53 Roughness: 7.943
Alpha: 581.736 RMS Misfit: 4.72 Roughness: 10.125
Step length set to 1.00000
Final RMS Misfit: 4.53
Updated data weights
Iteration 2
Initial RMS Misfit: 3.29 Number of data ignored: 0
Alpha: 750.048 RMS Misfit: 2.46 Roughness: 7.507
Alpha: 348.142 RMS Misfit: 2.02 Roughness: 8.803
Alpha: 161.593 RMS Misfit: 1.76 Roughness: 10.175
Alpha: 75.005 RMS Misfit: 1.62 Roughness: 11.762
Alpha: 34.814 RMS Misfit: 1.54 Roughness: 13.565
Alpha: 16.159 RMS Misfit: 1.51 Roughness: 15.497
Alpha: 7.500 RMS Misfit: 1.51 Roughness: 17.683
Alpha: 3.481 RMS Misfit: 1.52 Roughness: 20.859
Step length set to 1.00000
Final RMS Misfit: 1.51
Attempted to update data weights and caused overshoot
treating as converged
Solution converged - Outputing results to file
Calculating sensitivity map
Processing dataset 2
End of data: Terminating
1/1 results parsed (1 ok; 0 failed)
All ok
[8]:
k.showResults(contour=True, sens=True, zlim=[28,29.6]) # with contour
/media/jkl/data/phd/resipy/doc/gallery/../../src/resipy/meshTools.py:1487: MatplotlibDeprecationWarning: The collections attribute was deprecated in Matplotlib 3.8 and will be removed two minor releases later.
for col in cont.collections:
/media/jkl/data/phd/resipy/doc/gallery/../../src/resipy/Project.py:4136: MatplotlibDeprecationWarning: The collections attribute was deprecated in Matplotlib 3.8 and will be removed two minor releases later.
colls = mesh.cax.collections if contour == True else [mesh.cax]
In a nutshell
[10]:
k = Project(typ='R2')
k.createSurvey(testdir + 'syscal.csv')
k.importElec(testdir + 'elec.csv')
k.invert()
k.showResults(contour=True, sens=True, zlim=[28,29.6])
Working directory is: /media/jkl/data/phd/resipy/src/resipy
clearing dirname
600/636 reciprocal measurements found.
10 measurements error > 20 %
Creating triangular mesh...done (2236 elements)
Writing .in file and protocol.dat... done
--------------------- MAIN INVERSION ------------------
>> R 2 R e s i s t i v i t y I n v e r s i o n v4.10 <<
>> D a t e : 03 - 12 - 2023
>> My beautiful survey
>> I n v e r s e S o l u t i o n S e l e c t e d <<
>> Determining storage needed for finite element conductance matrix
>> Generating index array for finite element conductance matrix
>> Reading start resistivity from res0.dat
>> R e g u l a r i s e d T y p e <<
>> L i n e a r F i l t e r <<
>> L o g - D a t a I n v e r s i o n <<
>> N o r m a l R e g u l a r i s a t i o n <<
>> D a t a w e i g h t s w i l l b e m o d i f i e d <<
Processing dataset 1
Measurements read: 336 Measurements rejected: 0
Geometric mean of apparent resistivities: 0.35349E+03
>> Total Memory required is: 0.007 Gb
Iteration 1
Initial RMS Misfit: 33.54 Number of data ignored: 0
Alpha: 1253.313 RMS Misfit: 4.53 Roughness: 7.943
Alpha: 581.736 RMS Misfit: 4.72 Roughness: 10.125
Step length set to 1.00000
Final RMS Misfit: 4.53
Updated data weights
Iteration 2
Initial RMS Misfit: 3.29 Number of data ignored: 0
Alpha: 750.048 RMS Misfit: 2.46 Roughness: 7.507
Alpha: 348.142 RMS Misfit: 2.02 Roughness: 8.803
Alpha: 161.593 RMS Misfit: 1.76 Roughness: 10.175
Alpha: 75.005 RMS Misfit: 1.62 Roughness: 11.762
Alpha: 34.814 RMS Misfit: 1.54 Roughness: 13.565
Alpha: 16.159 RMS Misfit: 1.51 Roughness: 15.497
Alpha: 7.500 RMS Misfit: 1.51 Roughness: 17.683
Alpha: 3.481 RMS Misfit: 1.52 Roughness: 20.859
Step length set to 1.00000
Final RMS Misfit: 1.51
Attempted to update data weights and caused overshoot
treating as converged
Solution converged - Outputing results to file
Calculating sensitivity map
Processing dataset 2
End of data: Terminating
1/1 results parsed (1 ok; 0 failed)
All ok
/media/jkl/data/phd/resipy/doc/gallery/../../src/resipy/meshTools.py:1487: MatplotlibDeprecationWarning: The collections attribute was deprecated in Matplotlib 3.8 and will be removed two minor releases later.
for col in cont.collections:
/media/jkl/data/phd/resipy/doc/gallery/../../src/resipy/Project.py:4136: MatplotlibDeprecationWarning: The collections attribute was deprecated in Matplotlib 3.8 and will be removed two minor releases later.
colls = mesh.cax.collections if contour == True else [mesh.cax]
[ ]: