search

gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
384 stars 137 forks source link

3D topography and paraBoundary settings make strange topography #757

Open mariosgeo opened 3 months ago

mariosgeo commented 3 months ago

Problem description

The surface (and plc) made from xyz 3d points is not what I am expecting.

Your environment

Please provide the output of print(pygimli.Report()) here. If that does not work, please give provide some additional information on your:

Operating system: e.g. Windows, Python version: e.g. 3.9, pyGIMLi version: 1.5.1 Way of installation: e.g. Conda packag

Steps to reproduce

I have 96 electrodes and a DEM from the area (figure 1).

Figure_1

Since it's not a regular grid, I append in the perimeter some points and I set the elevation there to be the median (I believe you are doing the same)

Figure_2

After some interpolation on the DEM, we have now elevation in the regular grid (i.e. on the missing points). The whole point is to ease the creatin of the plc

Figure_3

I have appended the topographic points in the input file (rather than using the addTopo flag). That is, on the input file the first 96 points are actual the electrodes and the rest 55370 are just topography and never being used as measure.


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

data=ert.load('trial_new_mesh.txt')
# for some reason, I can't get the paraBoudary to be less than 1.2 (ideally 0). It returns an error
plc = mt.createParaMeshPLC3D(data, paraDepth=-1, paraMaxCellSize=0.0,
                              surfaceMeshQuality=34,paraBoundary=[1.2,1.2])

plc2=mt.createParaMeshSurface(data,paraBoundary=[1.2,1.2])

...

Expected behavior

The surface is looking like this

Figure_4

And the plc like this

Figure_5

How can I fix this?

trial_new_mesh (3).txt

mariosgeo commented 3 months ago

I did some "silly" digging. The issue (I think) is related to the polytools.py, line 1203

 sZ = pg.interpolate(pntsSurface, pnts[:, 2], surface.positions())

For some reason, the interpolation fails and returns the fallback value of 0.0 If I hard force the value to be the medians

sZ[sZ==0]=np.median(pnts[:, 2])

Then the mesh is fine.

Does this help? Notice that there is no bounarymarker==-2 anymore

Figure_6