2D Gravity inversion error

[Socigogo]

Problem description

Error reported when doing 2D gravity inversion on irregular meshes.

Your environment

Operating system: Windows Python version: 3.8 pyGIMLi version: 1.4.0 Way of installation: Anaconda

Steps to reproduce

my code:

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics.gravimetry.gravMagModelling import GravimetryModelling

x = np.arange(-25, 25.1, .5)
pnts = np.array([x, np.zeros(len(x))]).T

world = mt.createWorld(start=[-25, 0], end=[25, -20],
                       marker=1)
#mesh1 = mt.createMesh(world, quality=33, area=0.2)

rect = mt.createRectangle(start=[-6, -3.5], end=[6, -6.0],
                           marker=2, area=0.1)

geom = world + rect
pg.show(geom, markers=True)

mesh = mt.createMesh(geom, quality=33, area=0.2)
dRho = pg.solver.parseMapToCellArray([[1, 0.0], [2, 300]], mesh)

#CM = pg.utils.covarianceMatrix(mesh1, I=10)
#Cm05 = pg.matrix.Cm05Matrix(CM)
#C = pg.matrix.GeostatisticConstraintsMatrix(mesh=mesh1, I=10)
pg.show(mesh)
pg.show(mesh,dRho,cMap="Spectral_r",label="kg/m3")
pg.wait()

#g = solveGravimetry(mesh, dRho, pnts)
fop = GravimetryModelling(mesh)
fop.setSensorPositions(pnts)
g = fop.response(dRho)

ax1 = pg.plt.subplot(2, 1, 1)
ax1.plot(x, g)

ax2 = pg.plt.subplot(2, 1, 2)
pg.show(geom, ax=ax2)
ax2.plot(x, x*0,  'bv')

ax1.set_ylabel(r'$\frac{\partial u}{\partial z}$ [mGal]')
ax1.set_xlabel('$x$-coordinate [m]')
ax1.grid()
ax1.legend()

ax2.set_aspect(1)
ax2.set_xlabel('$x$-coordinate [m]')
ax2.set_ylabel('$z$-coordinate [m]')
ax2.set_ylim((-9, 1))
ax2.set_xlim((-25, 25))

pg.wait()

error = pg.Vector(len(g), 0.001)
tLog = pg.trans.TransLog()
inv = pg.Inversion(fop=fop)
inv.transData = tLog
inv.modelTrans = tLog
inv.lam = 40

#fop.setConstraints(C)
res = inv.run(g, error, lam=40, verbose=True)

pg.show(mesh, res)

### Expected behavior

I want to implement a 2D gravity inversion on an irregular grid and don't know where the code is wrong.

### Actual behavior

The following error occurred：

Traceback (most recent call last): File "D:\python project\Pygimli1.4.0\gravity2Dinversion.py", line 67, in res = inv.run(g, error, lam=40, verbose=True) File "C:\Users\86155\anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\inversion.py", line 587, in run startModel = self.startModel File "C:\Users\86155\anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\inversion.py", line 184, in startModel sm = self.fop.createStartModel(self.dataVals) Boost.Python.ArgumentError: Python argument types in None.createStartModel(GravimetryModelling, numpy.ndarray) did not match C++ signature: createStartModel(GIMLI::ModellingBase {lvalue})

Please, help me

[halbmy]

Thank you for reporting this issue. The 2D gravity modelling operator needs to be updated to meet pg>1.1 requirements.

[Socigogo]

Thank you for your reply. How to update this operator? Or is there a 1.1 version of gravity inversion code on an irregular mesh constrained by geostatistical operator？

[halbmy]

The modelling is a pygimli core operator, so you would be able to use it along with a core inversion. But that's not the proper way. Let me just do the adaptations.

[halbmy]

I made the modifications in the development branch (3f37bf2). To use it before the next release, you should checkout pyGIMLi by git as described in https://www.pygimli.org/installation.html#staying-up-to-date

A few other comments

• you should not use the forward mesh for inversion as the structures are already in by region decoupling)
• you should add random noise to the data to see whether the data are really fit within noise
• if the error is absolute, should should set it with inv(dataVals, absoluteError=
• there is no need to set any transformation by hand
• for data the default transformation is linear, leave it like that as anomalies can be negative
• for model the default transformation is log, either change it to linear or set boundaries by limits=[-5000, 5000], maybe with trans=Cot for a more linear interior
• there is no need to set the constraints matrix by hand, instead you can use inv.setRegularization(correlationLength=[Ix, Iz])

error = 0.001
g += np.random.randn(len(g)) * error
mesh = mt.createMesh(world, quality=33, area=1)
fop = GravityModelling2D()
fop.setMesh(mesh)
fop.setSensorPositions(pnts)
inv = pg.Inversion(fop=fop)
inv.setRegularization(limits=[-2000, 2000], trans="Cot",
                      correlationLengths=[10, 2])
rho = inv.run(g, absoluteError=error, lam=1e4, verbose=True)
pg.show(mesh, rho, logScale=False)

[Socigogo]

Thank you very much for your teaching, but I have such a problem now.

[halbmy]

I am sorry, I made a mistake in the documentation. The correct git command should be

git clone --branch dev https://github.com/gimli-org/gimli

I fixed this (a86a8ef), should be on http://dev.pygimli.org soon. Please try.

[halbmy]

Technically it is working (for you too?) and fitting the data

However, lust like in the 3D magnetics inversion example https://www.pygimli.org/_examples_auto/4_gravimetry_magnetics/plot_mod-magnetics-3d.html the anomalies are sticked to the surface. I am not a gravity expert but there should be methods to improve it. In the papers, a depth-weighting scheme is used, maybe it helps together with the geostatistic constraints.

[Socigogo]

It still doesn't work.

I'm now working on the depth-weighting function in the 3D magnetics inversion example. It doesn't seem to work, because the 3D inversion result in that example is still only near the surface.

[Socigogo]

I tried to remove the depth weighting in the example and run it, and the inversion result was exactly the same as before, indicating that the depth weighting function was not really used.

[halbmy]

In the post yesterday, you used git init which created a .git folder, i.e. your directory is already a governed by git. That's bad. Can you remove this .git folder and try again? Or git clone in another directory? On my computer it works.

[halbmy]

If you are behind a proxy server, you need to specify this via the https_proxy variable or git config --global http.proxy http://proxy.mycompany:80

[Socigogo]

My code terminated abnormally: This is my code now：

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics.gravimetry import solveGravimetry
from pygimli.physics.gravimetry.gravMagModelling import GravityModelling2D

x = np.arange(-25, 25.1, .5)
pnts = np.array([x, np.zeros(len(x))]).T

world = mt.createWorld(start=[-25, 0], end=[25, -20],
                       marker=1)
mesh = mt.createMesh(world, quality=33, area=1)

rect = mt.createRectangle(start=[-6, -3.5], end=[6, -6.0],
                           marker=2, area=0.1)

geom = world + rect
pg.show(geom, markers=True)

mesh1 = mt.createMesh(geom, quality=33, area=0.2)
dRho = pg.solver.parseMapToCellArray([[1, 0.0], [2, 300]], mesh1)
g = solveGravimetry(mesh1, dRho, pnts)
absoluteerror = 0.001
g += np.random.randn(len(g)) * absoluteerror

fop = GravityModelling2D()
fop.setMesh(mesh)
fop.setSensorPositions(pnts)

error = pg.Vector(len(g), 0.001)
inv = pg.Inversion(fop=fop)
inv.setRegularization(limits=[-2000, 2000], trans="Cot",
                      correlationLengths=[10, 2])
rho = inv.run(g, error, absoluteError=absoluteerror, lam=1e4, verbose=True)
pg.show(mesh, rho, logScale=False)

[halbmy]

On my computer, the code runs smoothly. I am having pgcore 1.4 and the development version, do you have the same? The result, however is crab as the error is wrong. Please do not specify both error and absoluteError:

rho = inv.run(g, absoluteError=absoluteerror, lam=1e4, verbose=True)

fits the data nicely

1: rms/rrms(data, Response) = 0.000949714/2642.94%
1: chi^2(data, Response, error, log) = 0.901956

with a model showing an anomaly glued to the surface

I just realized that for some reason the geostatistical constraints are not used and have to dig out why. But you should be able to reproduce my results first.

[Socigogo]

You're right, probably because I don't have the development version, I can't run this code： rho = inv.run(g, absoluteError=absoluteerror, lam=1e4, verbose=True) showing: TypeError: run() missing 1 required positional argument: 'errorVals'

I'm not going to do this now, but geostatistical constraints and depth weighting are not really useful in gravity and magnetic inversion in pygilmi, it really needs to dig out why.

[halbmy]

If this error occurs, you don't work with the current dev version being checked out by git.

[halbmy]

I've made some changes (64f742e) so that the geostatistic constraints work well. Additionally one can pass a weighting function

An appropriate weighting scheme

def depthWeighting(bz, z0=None):
    """Depth-weighting scheme adapted from Li&Oldenburg (1996)."""
    bd = max(bz) - bz
    return z0 / (bd+z0)**1.5

fop = GravityModelling2D(mesh=mesh, points=pnts)
inv = pg.Inversion(fop=fop)
inv.setRegularization(limits=[-2000, 2000], correlationLengths=[10, 2])
fop.region(1).setConstraintType(0)  # needed for correct sizes until fixed
inv.inv.setCWeight(depthWeighting(pg.y(mesh.cellCenters()), 2))
rho = inv.run(g, absoluteError=absoluteerror, lam=1e4, verbose=True)
ax, _ = pg.show(mesh, rho, logScale=False)
pg.viewer.mpl.drawPLC(ax, rect, fillRegion=False)

leads to a very reasonable model:

[Socigogo]

Wow, this is fantastic！ And this is exactly the result I wanted. Can you show me the results without geostatistical constraints? I would like to compare its effect. And I reinstalled pygimli, now how to do this in this situation?

[halbmy]

do you have a proxy? If yes, you have to insert the correct one. If no you have to unset the proxy.

[Socigogo]

I don't have a proxy.

[halbmy]

so get rid of that proxy line

[Socigogo]

[halbmy]

This should do and the code should run and provide the same results. If not, make sure to remove any old pygimli installation.

[halbmy]

I've also update the 3D magnetics example that benefits a lot from the depth weighting (which did not work before due to recent code restructures): http://dev.pygimli.org/_examples_auto/4_gravimetry_magnetics/plot_mod-magnetics-3d.html

If you can manage to run the code, we can close the issue. We will draft a new release where this will be included and thus easier to install.

[Socigogo]

You are really amazing! I admire you so much! Thank you so much for helping me with so many things. Both of my code worked very well：

[Socigogo]

I'm sorry, I have one more question. The geostatistical constraint in 3D magnetics inversion doesn't seem to work.

[halbmy]

The original issue (2D Gravity inversion error) has apparent been solved, so I am closing the issue again.

For the new problem (geostatistical constraints in 3D magnetics), please open a new issue where you describe your problem with the code. I am confident we can answer it all quickly with the newest dev version and taking the above comments into account.

[halbmy]

I have created an example for 2d gravity inversion: http://dev.pygimli.org/_examples_auto/4_gravimetry_magnetics/plot_inv-gravity-2d.html