Comments (12)
halbmy
commented on July 30, 2024
In your case, you have two regions: a main region with marker=1 and a circle with marker=2 inside. In such cases, the region with the lower marker is automatically treated as background. This needs to be deactivated by
fop.setRegionProperties(1, background=False)
from gimli.
zmd142536
commented on July 30, 2024
你好先生Halbmy，我遇到了一些问题超出了我的专业知识，当我完成建造前进的模型和编制计算。 能不能请你看一看在我的代码？
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-v0_8")
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.meshtools import polytools as plc#新加用以创建园
from pygimli.physics import ert
xBH = [-2, 2] # two boreholes, 4m spaced
zBH = -np.arange(1., 8.1, 1.0)
data = ert.DataContainer()
for x in xBH:
for z in zBH:
data.createSensor([x, z])
print(data)
world = mt.createWorld(start=[-4, -10], end=[4, 0], worldMarker=True, marker=1)
c1 = plc.createCircle(pos=[0.0, -5.0], radius=0.5, marker=2, area=.1)#创建园
geom = c1 + world
for pos in data.sensorPositions():
world.createNode(pos)
mesh = mt.createMesh(geom, area=.1, quality=34, smooth=[1, 10])
ax, cb = pg.show(mesh, markers=True, showMesh=True)
ax.plot(pg.x(data), pg.y(data), "mo");
def plotABMN(ax, scheme, idx):
""" Visualize four-point configuration on given axes. """
def getABMN(scheme, idx):
""" Get coordinates of four-point cfg with id `idx` from DataContainerERT
`scheme`."""
coords = {}
for elec in "abmn":
elec_id = int(scheme(elec)[idx])
if elec_id >= 0:
elec_pos = scheme.sensorPosition(elec_id)
coords[elec] = elec_pos.x(), elec_pos.y()
return coords
coords = getABMN(scheme, idx)
for elec in coords:
x, y = coords[elec]
if elec in "ab":
color = "green"
else:
color = "magenta"
ax.plot(x, y, marker=".", color=color, ms=10)
ax.annotate(elec.upper(), xy=(x, y), size=12, ha="center", #fontsize=10,
bbox=dict(boxstyle="round", fc=(0.8, 0.8, 0.8), ec=color),
xytext=(20*np.sign(x), 0), textcoords='offset points',
arrowprops=dict(arrowstyle="wedge, tail_width=.5", fc=color, ec=color,
patchA=None, alpha=0.75))
ax.plot(coords["a"][0],)
def showSensitivity(shm, mesh, idx=None, ax=None, nc=None):
idx = idx or range(shm.size())
nc = nc or len(idx)
if ax is None:
fig, ax = plt.subplots(len(idx)//nc, nc, sharex=True, sharey=True, figsize=(12, 8))
fop = ert.ERTModelling(verbose=False)
fop.setData(shm)
fop.setMesh(mesh)
fop.setRegionProperties(1, background=False) #确保低范围标记不被背景化
model = np.ones(mesh.cellCount())
fop.createJacobian(model)
for i, ind in enumerate(idx):
a = ax.flat[i]
sens = fop.jacobian()[ind]
normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)
normsens /= np.max(normsens)
pg.show(mesh, normsens, cMap="bwr", colorBar=False, ax=a,
label="sensitivity", nLevs=3, cMin=-1, cMax=1)
plotABMN(a, shm, ind);
# a.set_xlim(-2, 12)
# a.set_ylim(-5, 1)
data.createFourPointData(data.size(), 1, 9, 2, 10)
data.createFourPointData(data.size(), 1, 9, 4, 12)
data.createFourPointData(data.size(), 1, 9, 6, 14)
data.createFourPointData(data.size(), 2, 10, 3, 11)
data.createFourPointData(data.size(), 2, 10, 5, 13)
data.createFourPointData(data.size(), 2, 10, 7, 15)
data.createFourPointData(data.size(), 3, 11, 4, 12)
data.createFourPointData(data.size(), 3, 11, 6, 14)
data.createFourPointData(data.size(), 4, 12, 5, 13)
data.createFourPointData(data.size(), 4, 12, 7, 15)
data.createFourPointData(data.size(), 5, 13, 6, 14)
data.createFourPointData(data.size(), 2, 10, 3, 11)
data.createFourPointData(data.size(), 3, 11, 4, 12)
data.createFourPointData(data.size(), 4, 12, 5, 13)
data.createFourPointData(data.size(), 5, 13, 6, 14)
data.createFourPointData(data.size(), 6, 14, 7, 15)
print(data)
rhomap = [[1, 100.],
[2, 500.]]
mgr = ert.ERTManager()
data = ert.simulate(mesh, scheme=plotABMN, res=rhomap, noiseLevel=0.01, noiseAbs=1e-6, seed=1337)
pg.viewer.mpl.drawSensors(ax, data, diam=1.0, facecolor='white', edgecolor='black')
raise SystemExit
Error output：
`Traceback (most recent call last):
File ~/anaconda3/envs/pg/lib/python3.10/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
exec(code, globals, locals)
File ~/Spyder/kuakong-youdizhiti(baocuo).py:105
data = ert.simulate(mesh, scheme=plotABMN, res=rhomap, noiseLevel=0.01, noiseAbs=1e-6, seed=1337)
File ~/anaconda3/envs/pg/lib/python3.10/site-packages/pygimli/physics/ert/ert.py:163 in simulate
if not scheme.allNonZero('k') and not calcOnly:
AttributeError: 'function' object has no attribute 'allNonZero'`
from gimli.
halbmy
commented on July 30, 2024
You cannot just paste your mess of code and expect other people to start thinking. Just have a look at the error message: Obviously, plotABMN is of the wrong type, as it is a function. ert.simulate expects a data container as the docstring says, in your case data.
from gimli.
zmd142536
commented on July 30, 2024
I apologize for this, Mr. Halbmy. Please forgive my foolishness.But I do need professional help indeed. I've made the code modifications and it seems to have run successfully. However, I encountered difficulties during the data plotting step. I'm unsure if it's due to my electrode configuration being unreasonable or if the data plotting scheme across the boreholes isn't suitable, resulting in less than ideal data plotting.
Steps to reproduce
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-v0_8")
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.meshtools import polytools as plc#Newly added for creating circles.
from pygimli.physics import ert
from pygimli.viewer.mpl import showDataContainerAsMatrix
#The function to read data from a file is defined as follows
def read_txt_file(filename):
"""Read data from a TXT file and return it."""
with open(filename, 'r') as file:
lines = file.readlines()
# Parse the data and store it as a list.
data_list = []
for line in lines:
# Assuming the data is comprised of four integers separated by spaces
values = line.strip().split()
data_list.append([int(value) for value in values])
return data_list
#Create sensors and data containers
xBH = [-2, 2] # two boreholes, 4m spaced
zBH = -np.arange(1., 16.1, 1.0)
data = ert.DataContainer()
for x in xBH:
for z in zBH:
data.createSensor([x, z])
print(data)
world = mt.createWorld(start=[-4, -20], end=[4, 0], worldMarker=True, marker=1)
c1 = plc.createCircle(pos=[0.0, -8.0], radius=0.5, marker=2, area=.1)#Create a circle.
geom = c1 + world
for pos in data.sensorPositions():
world.createNode(pos)
mesh = mt.createMesh(geom, area=.1, quality=34, smooth=[1, 10])
#Read the data file and pass the data to the createFourPointData function
data_from_file = read_txt_file('/home/zmd/Spyder/paojifangshi1/paojifangshi2.txt')
for data_point in data_from_file:
data.createFourPointData(data.size(), *data_point)
ax, cb = pg.show(mesh, markers=True, showMesh=True)
ax.plot(pg.x(data), pg.y(data), "mo");
def plotABMN(ax, scheme, idx):
""" Visualize four-point configuration on given axes. """
def getABMN(scheme, idx):
""" Get coordinates of four-point cfg with id `idx` from DataContainerERT
`scheme`."""
coords = {}
for elec in "abmn":
elec_id = int(scheme(elec)[idx])
if elec_id >= 0:
elec_pos = scheme.sensorPosition(elec_id)
coords[elec] = elec_pos.x(), elec_pos.y()
return coords
coords = getABMN(scheme, idx)
for elec in coords:
x, y = coords[elec]
if elec in "ab":
color = "green"
else:
color = "magenta"
ax.plot(x, y, marker=".", color=color, ms=10)
ax.annotate(elec.upper(), xy=(x, y), size=12, ha="center", #fontsize=10,
bbox=dict(boxstyle="round", fc=(0.8, 0.8, 0.8), ec=color),
xytext=(20*np.sign(x), 0), textcoords='offset points',
arrowprops=dict(arrowstyle="wedge, tail_width=.5", fc=color, ec=color,
patchA=None, alpha=0.75))
ax.plot(coords["a"][0],)
def showSensitivity(shm, mesh, idx=None, ax=None, nc=None):
idx = idx or range(shm.size())
nc = nc or len(idx)
if ax is None:
fig, ax = plt.subplots(len(idx)//nc, nc, sharex=True, sharey=True, figsize=(12, 8))
fop = ert.ERTModelling(verbose=False)
fop.setData(shm)
fop.setMesh(mesh)
fop.setRegionProperties(1, background=False) #Ensure that the low-range markers are not backgrounded.
model = np.ones(mesh.cellCount())
fop.createJacobian(model)
for i, ind in enumerate(idx):
a = ax.flat[i]
sens = fop.jacobian()[ind]
normsens = pg.utils.logDropTol(sens/mesh.cellSizes(), 1e-2)
normsens /= np.max(normsens)
pg.show(mesh, normsens, cMap="bwr", colorBar=False, ax=a,
label="sensitivity", nLevs=3, cMin=-1, cMax=1)
plotABMN(a, shm, ind);
# a.set_xlim(-2, 12)
# a.set_ylim(-5, 1)
rhomap = [[1, 100.],
[2, 500.]]
mgr = ert.ERTManager()
data = ert.simulate(mesh, data, res=rhomap, noiseLevel=0.01, noiseAbs=1e-6, seed=1337)
pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
pg.info('Simulated data', data)
pg.info('The data contains:', data.dataMap().keys())
ert.show(data)
### error output ：
20/03/24 - 16:30:52 - pyGIMLi - WARNING - Something gone wrong while drawing data. Try fallback with equidistant electrodes.
There are 32 sensors in the data, but the number of nonzero entries is 0. Does this have anything to do with it?
from gimli.
zmd142536
commented on July 30, 2024
from gimli.
halbmy
commented on July 30, 2024
The print data command comes directly after generation of the sensors
data = ert.DataContainer()
for x in xBH:
for z in zBH:
data.createSensor([x, z])
print(data)
So of course there are no data at this point but only sensors. The arrays are added after it
data_from_file = read_txt_file('paojifangshi2.txt')
for data_point in data_from_file:
data.createFourPointData(data.size(), *data_point)
I cannot download your file this as something went wrong when you uploaded them.
from gimli.
zmd142536
commented on July 30, 2024
paojifangshi2.txt
This is my way of extreme running, thank you.
from gimli.
halbmy
commented on July 30, 2024
So what exactly is the problem? The pseudosection image?
There is no unique way to plot pseudosections for crosshole cases. The default ert.show(data) (data.show()) is only good for surface data. You can try data.show(style="A-M").
from gimli.
halbmy
commented on July 30, 2024
(I just used "AB-MN")
from gimli.
zmd142536
commented on July 30, 2024
Thank you, sir, as you said, this is related to the running mode. It seems I've got the answer I wanted, now I just need to optimize the running mode of the electrode to get good results.
here's one last question: regarding the mesh after inversion, I used unstructured mesh for modeling, but after inversion, I got a square grid. What controls this?
# %% inversion outmesh
inner = mt.createGrid(np.arange(-4, 4, 0.5), np.arange(-20.0, .1, 0.5),
marker=2)
mesh = mt.appendTriangleBoundary(inner, xbound=30, ybound=82)
ax, _ = pg.show(mesh, markers=True)
# %% inversion calculation
mgr = ert.ERTManager('AMBN1.dat')
inv = mgr.invert(data, mesh, secNodes=4, lam=10,verbose=True)
mgr.showResult()
from gimli.
halbmy
commented on July 30, 2024
You control the inversion mesh by yourself, in this case a grid (surrounded by a trianglular outer mesh for forward calculation), but of course you can create a triangle mesh as done in other examples.
from gimli.
halbmy
commented on July 30, 2024
As you "seem to got the answer wanted", I am closing the issue.
from gimli.
Related Issues (20)
- This may be a corrupt mesh error HOT 7
- pyvista shows no result graph when I tried to reproduce the example '3D surface ERT inversion' HOT 2
- ERT: blanked out cells when adding new regions to the mesh HOT 4
- Cannot import name 'RVector3' form 'pygimli.core.core HOT 6
- how to fix import pygimli error libcholmod.so.3 on Mybinder HOT 9
- Plot Inversion Result from Another Software use Pygimli HOT 2
- fitReciprocalErrorModel IP data HOT 7
- Seeking Help with 3D Magnetic Inversion HOT 2
- Question about geometric factor for an ERT survey with submerged and terrestrial electrodes HOT 3
- ERTManager.invert() --> error: malloc(): unsorted double linked list corrupted HOT 3
- Failure to install on mac HOT 6
- Wrong coordinates in tutorial notebook HOT 4
- pyGIMLi and NUM_THREADS HOT 7
- 3d geometry with topography and fixed thin layer HOT 5
- pg.interpolate() giving error with regionwise regularization HOT 6
- Discussion regarding `errorVals`-`absoluteError`-`relativeError` HOT 4
- Petrophysical joint inversion: smooth the output result HOT 3
- Running pyGIMLi on GHA fails to import HOT 3
- Can the gimli module perform joint HVSR-Dispersion Curve inversion? HOT 1
- The resistivity values are inverted after inversion with dipole-dipole array HOT 8
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from gimli.