gimli-org / gimli

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

showResult only of the tt.model resolved #322

Closed mcardenass closed 3 years ago

mcardenass commented 3 years ago

Problem description

Dear Friends

I use Crosshole traveltime tomography code to build Seismic Interferometry tomography images. The source-reciever array is irregular. So, I need to plot only the cells inside of array.

Your environment

Operating system: Mac Python version: Python 3.7.9 pyGIMLi version: 1.1.1+0.gbb685a8c (with local changes) Way of installation: . Conda package,

Steps to reproduce

Here are the minimum steps to produce the attached figure

fig, ax = plt.subplots(1, 1, figsize=(7, 7), sharex=True, sharey=True) ax.set_title("Resultado de Inversion en "+freq+" (Hz)") ax.plot(pg.x(data.sensors()), pg.y(data.sensors()), "wo") tt.showResult(ax=ax, logScale=True, nLevs=5, cMin=50, cMax=500, label=r'Velocidad (${\tt m/s}$)') tt.drawRayPaths(ax=ax, color="0.8", alpha=0.3) fig.tight_layout() fig.savefig('tomo'+freq+'.png')

Expected behavior

Only I want the velocity model inside the array, not the interpolated values out of receiver array. As if you use the inpolygon function of matlab. In the figure attached , the results between first and last receivers could be cut by the line that joint that receivers, as if we had a rectangle

Actual behavior

Please see the figure attached. Thanks

If possible, the area out of receiver array should be white I tomo12 7

halbmy commented 3 years ago

Can you just use tt.showResult(..., coverage=tt.standardizedCoverage())?

halbmy commented 3 years ago

Actually, you should create the mesh not using a rectangle, but using pygimli.meshtools.createPolygon instead to avoid inverting outside of the sensors at all. The rays outside the lower sensor line should not really happen.

mcardenass commented 3 years ago

Thanks Tomas,

tomo12 7

I managed to perform the tomography by creating a polygon by suggesting that you indicate me (pygimli.meshtools.createPolygon). I show you an example image. One last annoyance, I read the positions of the receivers, which in this case are the x y pairs (or Verts) that the mt.createPolygon function needs. How I pass the array of coordinate pairs to the function mt.createPolygon (verts [] ,..)? For example:

x,y = np.loadtxt('xcoord_unit.dat', unpack=True) sensors=np.dstack([x,y]) sensors=sensors.tolist()

I want to use this form

geom= mt.createPolygon(sensors, isClosed=True, area=1)

I attach a small script illustrating my question. Thank you in advance

Best regards

Martín

xtest.zip

mcardenass commented 3 years ago

I found my mistake, I have to pass the list of coordinates using the * operator before an iterable to expand it within the function call

Thanks