Open mariosgeo opened 3 weeks ago
Regions with fixed values are background regions (like the big box around every normal inversion) just needed for the forward calculation and therefore not part of the inverse problem. As a result, the values are not mapped back from the model vector and therefore not plotted as in the region tutorial https://www.pygimli.org/_tutorials_auto/3_inversion/plot_8-regionWise.html#model-reduction
You add this by hand using
ax, cb = mgr.showResult(hold=True)
pg.show(water, pg.Vector(water.cellCount(), 22.5), ax=ax)
We could potentially also include fixed regions into mgr.paraDomain
and mgr.paraModel
, or directly from mgr.showResult()
, using a keyword, but it should not be the default behaviour.
Hi Thomas,
Thanks for the reply, but I think I still don't get it. I am attaching a small data file that should take few seconds to run (ws__00.txt and bath_00.txt).
I am not actually fixing the the water layer, but I rather use a lower and upper limit (not sure how's it's treated in the inversion, but that's another question). Still, when I
mgr.showResult()
I get this image (the water layer is present)
when I
mgr.saveResult('demo')
The pdf looks fine (including the water layer) resistivity.pdf
The vtk file looks strange
I suppose the mesh in the vtk does not include the the fixed area?
Follow up
water = mesh.createSubMesh(mesh.cells(mesh.cellMarkers() == 3))
back = mesh.createSubMesh(mesh.cells(mesh.cellMarkers() ==1))
rest2 = mesh.createSubMesh(mesh.cells(mesh.cellMarkers() ==2))
len(mgr.model) =1504
mesh.cellCount() =2100
back.cellCount() =596
rest2.cellCount() = 581
water.cellCount() = 923
So if I understand correctly, you export the values of the rest2 (581) + water mesh (923) in the mgr.model (1504), but the way you plot it in the vtk has these values NOT in mesh order ?
For Instance
pg.show(mgr.paraDomain,mgr.model)
works, not sure where you export the mesh in vtk
Edit 2:
I think I found the bug (?)
When you export in vtk, you export as
m['Resistivity'] = self.paraModel(self.model)
but if I export as
m['Resistivity'] = self.model
Then it works.
Edit 3:
In case someone else has a similar issue, this is how I solved it, without re-running the inversion and assuming you have saved the results
`
import pygimli as pg
import numpy as np
res=pg.load('resistivity.vector')
sen=pg.load('resistivity-cov.vector')
ssen=pg.load('resistivity-scov.vector')
# mesh=pg.load('mesh.bms')
# res_mesh=pg.load('resistivity-mesh.bms')
pd_mesh=pg.load('resistivity-pd.bms')
# pg.show(mesh)
# pg.show(res_mesh)
# pg.show(pd_mesh,res,cMap='rainbow',logScale=True)
pd_mesh['Resistivity'] = res
pd_mesh['Resistivity (log10)'] = np.log10(res)
pd_mesh['Coverage'] = sen
pd_mesh['S_Coverage'] = ssen
pd_mesh.exportVTK('resistivity2.vtk')
`
Well, the model vector mgr.model
is as it comes out of the inversion, i.e. sorted according to the regions and single regions only represented by a single number. Therefore the function mgr.paraModel()
is needed to sort them into the mesh (where a single region covers many cells and region markers are unsorted). I'll have a closer look on what's going wrong there taking the lake example.
Thanks Thomas. Hope my "case" is helpful.
Problem description
I am inverting some marine ERT data, using floating cable. I am fixing the bathymetry and resistivity value of the water layer, but when I save the results, the VTK file looks off while the pdf looks ok.
Your environment
Operating system: e.g. Windows Python version: 3.11,1 pyGIMLi version: Output of
print(pygimli.__version__)
Way of installation: e.g. CondaSteps to reproduce
Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:
Expected behavior
The inversion results saved as vtk and pdf match ws__02.txt
Actual behavior
Rename resistivity_vtk.txt to resistivity.vtk
bath_02.txt resistivity.pdf ws__02.txt resistivity_vtk.txt