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/
.
%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.
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.
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).
k.elec
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.
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.
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.
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
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]
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]