-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path02_Inverse.py
112 lines (78 loc) · 3.56 KB
/
02_Inverse.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 2 13:28:02 2021
@author: lwang19
"""
import numpy as np
import os
import pygimli as pg
import pygimli.physics.ert as pg_ert
#from pygimli.physics.ert import ERTManager, createGeometricFactors
#from pygimli.meshtools import createParaMesh
import pygimli.meshtools as mt
from pygimli.viewer import showMesh
from pygimli.meshtools import appendTriangleBoundary
from pygimli.meshtools import createParaMesh2DGrid
import pygimli.meshtools as plc
import pybert as pb
import pandas as pd
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import seaborn as sns
rhomap = [[0,35],
[1, 80.],
[2, 10.]]
inverse_par = pd.Series({'plotmin':1,'plotmax':100,'zmin':0,'zmax':30,'dx':0.25,'lambda_reg':5,\
'paraDX':0.3,'paraMaxCellSize':0.5,'paraDepth':10,'quality':33.6})
data_orig = pb.importData('./A_Simresults/Simulated_resistance.dat')
mesh = pg.Mesh('./A_Simresults/Simulation_mesh.bms')
ert = pg_ert.ERTManager(sr=False, useBert=True, verbose=True, debug=False)
#%%
# invert with auto-generated mesh
# model = ert.invert(data_orig,verbose=True,paraDX=inverse_par.paraDX, paraMaxCellSize=inverse_par.paraMaxCellSize,\
# paraDepth=inverse_par.paraDepth, quality=inverse_par.quality)
model = ert.invert(data_orig, lam=inverse_par.lambda_reg, verbose=True)
# np.testing.assert_approx_equal(mgr.inv.chi2(), 0.7, significant=1)
# mgr.showResultAndFit()
# meshPD = pg.Mesh(mgr.paraDomain) # Save copy of para mesh for plotting later
# inversionDomain = pg.createGrid(x=np.linspace(start=-18, stop=18, num=33),
# y=-pg.cat([0], pg.utils.grange(0.5, 8, n=5))[::-1],
# marker=12)
# grid = pg.meshtools.appendTriangleBoundary(inversionDomain, marker=11,
# xbound=50, ybound=50)
# pg.show(grid, markers=True)
# pg.show(grid, markers=True)
# %%
#invert with user-defined mesh
# inversionDomain = pg.createGrid(x=np.linspace(start=0, stop=71, num=72),
# y=-pg.cat([0], pg.utils.grange(0, 20, n=21))[::-1],
# marker=3)
# grid = pg.meshtools.appendTriangleBoundary(inversionDomain, marker=4,
# xbound=50, ybound=50)
# pg.show(grid,markers=True)
# # %%
# model = ert.invert(data_orig,mesh=grid, lam=inverse_par.lambda_reg, verbose=True)
#np.testing.assert_approx_equal(ert.inv.chi2(), 0.951027, significant=3)
# %%
# np.testing.assert_approx_equal(mgr.inv.chi2(), 1.4, significant=2)
modelPD = ert.paraModel(model) # do the mapping
pg.show(ert.paraDomain, modelPD, label='Model', cMap='Spectral_r',
logScale=True, cMin=0, cMax=100)
pg.info('Inversion stopped with chi² = {0:.3}'.format(ert.fw.chi2()))
fig, (ax1,ax2) = plt.subplots(2,1, sharex=True, sharey=True, figsize=(8,7))
pg.show(mesh, rhomap, ax=ax1, hold=True, cMap="Spectral_r", logScale=True,
orientation="vertical", cMin=0, cMax=100)
# pg.show(meshPD, inv, ax=ax2, hold=True, cMap="Spectral_r", logScale=True,
# orientation="vertical", cMin=25, cMax=150)
ert.showResult(ax=ax2, cMin=0, cMax=100, orientation="vertical")
labels = ["True model", "Inversion unstructured mesh", "Inversion regular grid"]
for ax, label in zip([ax1,ax2], labels):
ax.set_xlim(0, 71)
ax.set_ylim(-10, 0)
ax.set_title(label)
#%% Save files
paramesh = ert.paraDomain
paramesh.save('./B_Invertresesults/invert_mesh')
# pg.show(paramesh)
model.save('./B_Invertresesults/invert_model.vector')