3D gravity inversion with topography

[mariosgeo]

Starting by saying sorry for bombing with so many questions....

I have few hundreds xyz points, that define topography points (and define ert electrode positions). On the same mesh (and points) I have also gravimetry data, hence I want to do a 3d gravity inversion. The way I create the mesh is the following. What I want, is the inversion to take place on the mesh with Marker=2 (the red region). The way I set it up, inverts on the whole mesh (marker 1 and 2). How can I do this?
I can't seem to understand the terms paraBoundary and Boundary in the function.

`data=ert.load('apolonia_1based_yes_grav.pyg')`
`plc` = mt.createParaMeshPLC3D(data, paraDepth=12, paraMaxCellSize=3,surfaceMeshQuality=34)

mesh = mt.createMesh(plc, quality=1.3)
idx=mesh.findCellByMarker(2)
 #read gravity measures
grav=pd.read_hdf('grav_data_topo_boug.h5',key='data',mode='r')
points=grav[['X','Y','Z']].values
fop = GravityModelling(mesh, points)
data =grav['Boug'].values
noise_level = 0.1  #why 0.1?

#depth weighting
bz = np.array([b.center().z() for b in mesh.boundaries() if not b.outside()])
bn = np.array([b.norm().z() for b in mesh.boundaries() if not b.outside()])
z0 = 100
wz = z0 / (bz+z0)**1

inv = pg.Inversion(fop=fop, verbose=True)  # , debug=True)
inv.modelTrans = pg.trans.TransCotLU(-2, 2)
#inv.setRegularization(correlationLengths=[500, 500, 100])
inv.setConstraintWeights(wz)
invmodel = inv.run(data, absoluteError=noise_level, lam=3e4,  # zWeight=0.3,
               startModel=0.1, verbose=True)

`

[mariosgeo]

I suppose I could do this?

 mesh2=pg.Mesh()
 mesh2.createMeshByMarker(mesh,2)

But when I run

 fop = GravityModelling(mesh2, points)

I get the error

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

[halbmy]

I am note sure whether mesh2=pg.Mesh(dim=3) already helps.

What's the shape of points?

[mariosgeo]

points.shape
Out[29]: (341, 3).

Inversion runs for the full mesh. (i.e.)

fop = GravityModelling(mesh, points) 

Looking at the code, the issue seems to be in kernel.py, line 51-55, where c_list is of NoneType, instead of Tetrahedron.

mesh2=pg.Mesh(dim=3)

Does not help

[halbmy]

does mesh2.createNeighborInfos() help? If not, can you mesh2.save() and send it along with points?

[mariosgeo]

apolonia_1based_yes_grav.txt

The data I use to make mesh and mesh2

The createNeightborInfos() did nit help

[halbmy]

I am using the original mesh2 I run into the same problem, if I use mesh2.createNeighborInfos() I don't.

The electrodes have elevations between 65 and 85m, paraDepth=12 does not make any sense because the lower face of the mesh will cut the topography. Using much higher values makes more sense. Additionally, you can use boundary=0 to avoid the boundary region, so you don't have to strip the outer mesh at all.

[mariosgeo]

Thanks, I will give another look to make sure. I must have missed something.

Ps. I thought paraDepth =12 is 12 m below the average depth between 65 and 85. I.e 62 meters. I run some ert inversions with that mesh, but now you made me wonder...

Thanks Thomas

[halbmy]

I am not quite sure and I have not plotted the topography but usually it is from the outer corners down. There might be also tetgen tolerance problems with the larger coordinate numbers that could be substracted.