Inversion with fixed region does not export results

[mariosgeo]

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. Conda

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
import numpy as np


data=ert.load('ws__02.txt')
# pg.viewer.mpl.showDataContainerAsMatrix(data, "a", "m", "rhoa",cMap='rainbow',cMin=10,cMax=1000,logScale=True)

nn=np.loadtxt('bath_02.txt')
plc = mt.createParaMeshPLC(data, paraDx=1,paraDepth=20, paraMaxCellSize=5,
                          surfaceMeshQuality=32,boundaryMaxCellSize=2500)

plc2=mt.createPolygon(list(nn),isClosed=True,marker=3)    
# pg.show(plc2)
plc=plc+plc2
    
# pg.show(plc, markers=True)
      

mesh = mt.createMesh(plc, quality=32)
   

for b in mesh.boundaries():
     if b.marker() == -1 and not b.outside():
         b.setMarker(2)
        
# pg.show(plc, markers=True)
# pg.show(mesh, markers=True, showMesh=False)
water = mesh.createSubMesh(mesh.cells(mesh.cellMarkers() == 3))


# pg.show(water)
    
mgr = ert.ERTManager(data)
mgr.setMesh(mesh)  # use this mesh for all subsequent runs

mgr.inv.setRegularization(3, limits=[20, 35], trans="log")
# mgr.inv.setRegularization(3, fix=22.5, trans="log")

mgr.invert(verbose=True,robustData=False,blockyModel=False)

# mgr.showResultAndFit()
# mgr.showFit() 

misfit = - mgr.inv.response / data["rhoa"] * 100 + 100


...

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

Paste your script output here.

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[halbmy]

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.

[mariosgeo]

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?

bath_00.txt
ws__00.txt

[mariosgeo]

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

[mariosgeo]

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.

[mariosgeo]

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')

`

[halbmy]

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.

[mariosgeo]

Thanks Thomas. Hope my "case" is helpful.