how to display pseudosection and phiM and phiD curves and tikhonov curves in ERT

[keenemarcel]

Problem description

I want to show pseudosection and phiM and phiD curves as well as tikhonov curves in the ERT example at https://www.pygimli.org/_examples_auto/3_dc_and_ip/plot_01_ert_2d_mod_inv.html.
and i want to know when input data there is iperr, r, u what value is it?

Your environment

Operating system: Linux ?
Python version: e.g. 3.7?
pyGIMLi version: Output of print(pg.__version1.1)
Way of installation: e.g. Conda package

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:
mgr.inv.getPhi()
mgr.inv.getPhiD()
mgr.inv.getPhiM()

AttributeError Traceback (most recent call last)
in
----> 1 mgr.inv.getPhi()
2 mgr.inv.getPhiD()
3 mgr.inv.getPhiM()

AttributeError: 'Inversion' object has no attribute 'getPhi'

Expected behavior

I want to be able to display pseudosection, phiD curve, phiM, and tikhonov curve

Actual behavior

I want to be able to display pseudosection, phiD curve, phiM, and tikhonov curve

Paste 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".

[keenemarcel]

I want to ask a misfit for an example, https://www.pygimli.org/_examples_auto/3_dc_and_ip/plot_02_ert_field_data.html#sphx-glr-examples-auto-3-dc-and-ip-plot-02-ert-field-data-py, I've tried and looked at the issues #254

misfit = mgr.inv.response / mgr.data['rhoa'] * 100 - 100
pg.show(mgr.data, misfit, cMap="bwr", cMin=-10, cMax=10, label="misfit (%)")  

but the rhoa cannot display the value

[keenemarcel]

I've tried this and changed the mgr to ert

[halbmy]

Generally there are two ways of working with a manager:

1. Initializing it with the data by Manager('datafile') or loading it after by mgr.load('datafile')
2. Initialize an empty manager without data and passing the data in the mgr.invert(data) call

Whereas
https://www.pygimli.org/_examples_auto/3_dc_and_ip/plot_01_ert_2d_mod_inv.html
uses the first way,
https://www.pygimli.org/_examples_auto/3_dc_and_ip/plot_02_ert_field_data.html
uses the latter way, in which ert.data is not known to the manager. However, the inversion knows the data values (i.e. the apparent resistivity) in ert.inv.dataVals and the data container itself is known by you. So the code for plotting the misfit becomes:

misfit = ert.inv.response / ert.inv.dataVals * 100 - 100
pg.show(data, misfit, cMap="bwr", cMin=-10, cMax=10, label="misfit (%)")

[halbmy]

As to the objective function, it can be retrieved by ert.inv.phi(), the data and model parts by ert.inv.phiData() and ert.inv.phiModel(). But these are just values. If you want them as a function of the iteration number like ert.inv.chi2History you should go into inversion postStep functions, or we should just implement ert.inv.phiHistory and ert.inv.phiMHistory (as phiDHistory is already covered by chi-square times the number of data).

[keenemarcel]

How do you know if the model is fit? what is the value of error and misfit? or there are other variables?

[halbmy]

You wanted to say: "if the data are fitted"? The chi-square (see comment above its history) is a measure for the error-weighted data fit, if it is around 1, the data are fit within the error model. The error model is composed by a relative error (typically a few per cent) and an absolute error (a voltage uncertainty). Please read some relevant papers (e.g. Günther et al. 2006 as the standard ERT reference) and text books for background.

[keenemarcel]

how to display Global L – curve displaying data misfitχ2and model roughness for varyingλ, the curvature of the curve has a clear maximum defining the optimumλ and histogram for model

[halbmy]

Dear Keenemarcel, first I would really like you to pay some attention to how you write, full sentences and clear questions. You can use the preview button to see how your post look like.

[halbmy]

I agree that a function for plotting the global L-curve can be very helpful and we consider adding it to the base method manager. Note however, that the L-curve does not always look like an L and that the curvature is not guaranteed to show a clear maximum. I rather suggest using the discepance principle for choosing optimum regularization.

Anyway, you can do that by a simple loop using the quantities above:

lambdas = np.logspace(...  # choose some meaningful range
phiD = []
phiM = []
for lam in lambdas:
    ert.invert(lam=lam)
    phiD.append(ert.inv.phiData())
    phiM.append(ert.inv.phiModel())

plt.plot(phiD, phiM, ....)
curvature =  # use log 

We'd be pleased if you can share your solution (code and figures) here so that all users can benefit.

[keenemarcel]

I have a model and I want to ask how to change the depth parameter? because it is the elevation that I entered

[halbmy]

What do you mean with "depth parameter"? A depth axis does not make much sense in case of varying topography because then you would have to define where the zero elevation is. In your case it does not make sense anyway to do the computations with topography because it is flat. So you could invert the data without topography, plot it with a depth axis and then you could rotate the figure.

[keenemarcel]

without topography, the depth is 20 m, how do you make the image bigger?

[keenemarcel]

sorry I ask a lot because I am still a layman and still learning a lot

[halbmy]

You mean how to make the model deeper? You can add an depth argumento to the inversion call: ert.invert(depth=30, ...)

[keenemarcel]

how to display pseudosection observation data and prediction data

[keenemarcel]

how is the input position of the wenner alpha electrode ?

08/09/20 - 20:08:24 - pyGIMLi - INFO - Found 2 regions.
08/09/20 - 20:08:24 - pyGIMLi - INFO - Region with smallest marker set to background (marker=1)
08/09/20 - 20:08:24 - pyGIMLi - INFO - Creating forward mesh from region infos.
08/09/20 - 20:08:24 - Core - WARNING - Region Nr: 1  is background and should not get a model transformation.
08/09/20 - 20:08:24 - Core - WARNING - Region Nr: 1  is background and should not get a model control.
08/09/20 - 20:08:24 - pyGIMLi - INFO - Creating refined mesh (H2) to solve forward task.
<class 'pygimli.physics.ert.ert.ERTManager'>.applyMesh(methodManager.py:647) : Mesh: Nodes: 1500 Cells: 2879 Boundaries: 4378
08/09/20 - 20:08:24 - pyGIMLi - INFO - Starting inversion.
08/09/20 - 20:08:24 - pyGIMLi - INFO - Set default startmodel to median(data values)=-57.91708301945876
08/09/20 - 20:08:24 - pyGIMLi - INFO - Created startmodel from forward operator: 1863 [-57.91708301945876,...,-57.91708301945876]
08/09/20 - 20:08:25 - Core - CRITICAL -  response for model with negative or zero resistivity is not defined.: -57.9171 -57.9171
fop: <pygimli.physics.ert.ert.ERTModelling object at 0x7f3bc1fcb470>
Data transformation: <pygimli.core._pygimli_.RTransLogLU object at 0x7f3bc1fcb570>
Model transformation: <pygimli.core._pygimli_.RTransLog object at 0x7f3bc1fcb530>
min/max (data): -429/-26.51
min/max (error): 0.96%/1.0%
min/max (start model): -57.92/-57.92
--------------------------------------------------------------------------------
08/09/20 - 20:08:31 - Core - WARNING - /home/wagner/miniconda3/conda-bld/pygimli_1589132014332/work/gimli/gimli/core/src/bert/datamap.cpp: 221		GIMLI::RVector GIMLI::DataMap::data(const GIMLI::DataContainerERT&, bool, bool) 
 a = 0 b = 2 m = 3 n = 1
 0 0 0 0 U = 0

08/09/20 - 20:08:31 - Core - WARNING - /home/wagner/miniconda3/conda-bld/pygimli_1589132014332/work/gimli/gimli/core/src/bert/datamap.cpp: 221		GIMLI::RVector GIMLI::DataMap::data(const GIMLI::DataContainerERT&, bool, bool) 
 a = 1 b = 3 m = 4 n = 2
 0 0 0 0 U = 0

...

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-8-7495810011f5> in <module>
      1 inv= ert.invert(data, lam=20, 
----> 2                  paraDX=0.5, paraMaxCellSize=7, paraDepth=40, quality=33.6)
      3 print(inv)

~/anaconda3/envs/pg/lib/python3.7/site-packages/pygimli/frameworks/methodManager.py in invert(self, data, mesh, zWeight, startModel, **kwargs)
    717 
    718         self.preRun(**kwargs)
--> 719         self.fw.run(dataVals, errorVals, **kwargs)
    720         self.postRun(**kwargs)
    721         return self.paraModel(self.fw.model)

~/anaconda3/envs/pg/lib/python3.7/site-packages/pygimli/frameworks/inversion.py in run(self, dataVals, errorVals, **kwargs)
    422                 self._preStep(0, self)
    423 
--> 424         self.inv.start()
    425         self.maxIter = maxIterTmp
    426 

RuntimeError: /home/wagner/miniconda3/conda-bld/pygimli_1589132014332/work/gimli/gimli/core/src/inversion.cpp: 95		double GIMLI::RInversion::getPhiD(const Vec&) const  getPhiD == -nan

[halbmy]

how is the input position of the wenner alpha electrode ?

I cannot detect any meaning in this question. There is no such thing as a Wenner electrode.

Unless you provide a script and/or data, we cannot follow your problem and help you.

By the way, it seems like it has nothing to do with the topic of this issue. So better close this one and open a new issue, but will full information.