Run Gravity Inversion by using real field data

82 views
Skip to first unread message

Carlos Osmin Pocasangre Jimenez

unread,
Mar 26, 2018, 4:58:15 AM3/26/18
to Fatiando a Terra
Hello,
I'm working on Gravity data from my study are so I want to run an inversion modeling. I have checked the Cookbook and all examples use synthetic data for running the inversion modeling. How can I use the fatiando package and my field data?

best regards
Carlos

Leonardo Uieda

unread,
Mar 26, 2018, 3:51:52 PM3/26/18
to fatiando
Hi Carlos,

Yes you can. Everything works with field data. The examples all use synthetic because it's easier to show what the method does when we know the ground truth.
You'll have to load your data somehow.
I recommend using an ASCII xyz file and using the numpy.loadtxt function (https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html).

Leo


--
Leonardo Uieda

Visiting Researcher | University of Hawaii at Manoa
www.leouieda.com

--
You received this message because you are subscribed to the Google Groups "Fatiando a Terra" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fatiando+unsubscribe@googlegroups.com.
To post to this group, send email to fati...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Carlos Jimenez

unread,
Mar 29, 2018, 4:58:00 AM3/29/18
to Fatiando a Terra
Hello Leo, 
Thanks for the rapid response. I wrote a code using based on your examples. However, I got some issues. First, the minimum RMS is 2.014, second, the MayaVI 3d plot is transposed. I copied a piece of code next. I imported the data as you said xyz.



if __name__ == "__main__": print '\n---- 3D GRAVITY INVERSION ----' # upload data xp, yp, gz = np.loadtxt('local.csv', comments='#', skiprows=1, delimiter=',', unpack=True) zp = np.empty(len(xp)) zp.fill(-1.0) xp = xp - min(xp) yp = yp - min(yp) bounds = [0, 11000, 0, 11000, 0, 4000] prism_shape = (20, 40, 40) mesh = mesher.PrismMesh(bounds, prism_shape) seeds = harvester.sow( [(4000, 8689, 2000, {'density': 400}), (4181.91, 7825, 2000, {'density': -400}), (2349.18, 6725.36, 2000, {'density': 400}), (935.351, 4473.71, 2000, {'density': -400}), (3972.46, 1488.96, 2000, {'density': 400}), (6590.65, 494.051, 2000, {'density': -400}), (10177.6, 7458.45, 2000, {'density': -400}), (6067.01, 3190.79, 2000, {'density': 400}), (7009.56, 1881.69, 2000, {'density': 400}), (7978.3, 3374.07, 2000, {'density': -400}), (4836.46, 5966.08, 2000, {'density': -400}), (6695.38, 5547.17, 2000, {'density': 400}), (8449.57, 6227.9, 2000, {'density': -400}), (9627.76, 8165.36, 2000, {'density': 400}), (8632.85, 9212.64, 2000, {'density': -400})], mesh) data = [harvester.Gz(x, y, z, g_obs)]
estimate, predicted = harvester.harvest(data, seeds, mesh, compactness=0.08, threshold=0.00001) mesh.addprop('density', estimate['density']) bodies = mesher.vremove(0, 'density', cube_mesh)


Leonardo Uieda

unread,
Mar 30, 2018, 8:52:40 PM3/30/18
to fatiando
Hi Carlos,

Your code got a bit mangled in the email.
Please send the script and some screenshots of the problem.
It's kind of difficult to diagnose without seeing the output and the full code.



--
Leonardo Uieda

Visiting Researcher | University of Hawaii at Manoa
www.leouieda.com

Filippo Spadavechia

unread,
Mar 30, 2018, 11:10:32 PM3/30/18
to fati...@googlegroups.com
leonardo do you have a scrip to visualize sgy files using python 3.5 ?
Message has been deleted

Leonardo Uieda

unread,
Apr 4, 2018, 11:34:11 PM4/4/18
to fatiando
Hi Filippo,
We don't but maybe this can do the job: https://github.com/Statoil/segyviewer

--
Leonardo Uieda

Visiting Research Scholar | University of Hawaii at Manoa
www.leouieda.com

Filippo Spadavechia

unread,
Apr 5, 2018, 10:27:22 AM4/5/18
to fati...@googlegroups.com
ok let me try it! thanks! 

accep handyarso

unread,
Dec 23, 2018, 9:55:22 AM12/23/18
to Fatiando a Terra
Hello Leo,,
I am very interested with your Planting Seeds algorithm. I have installed fatiando and tried to run the Gravity Inversion scripts but I get stuck with importing "vremove", could you please to help??

I don't know what missing, perhaps some scripts but I am not sure
thanks for your kindly help

regards

Leonardo Uieda

unread,
Jan 2, 2019, 2:33:39 PM1/2/19
to fati...@googlegroups.com
Hi,

Could you please tell me which version you installed, what OS, which Python (Anaconda, etc) and which script you ran?


--
Leonardo Uieda

Visiting Research Scholar
Department of Earth Sciences - SOEST
University of Hawai'i at Mānoa
www.leouieda.com


--
You received this message because you are subscribed to the Google Groups "Fatiando a Terra" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fatiando+u...@googlegroups.com.

accep handyarso

unread,
Feb 3, 2019, 8:17:50 PM2/3/19
to Fatiando a Terra
Hi Leo

I would like to apologize for my late reply.
the code (fatiando) is already works Leo, thank you for your kindly attention and help
at that time, I wiped out my Ubuntu 18.04 and reinstall everythings.
and now everything works well,
my OS is Ubuntu 18.04
python 2.7 (i use anaconda 2)
the scripts I run is gravity harvester
and the result is as below

Screenshot from 2019-02-03 10-58-27.png



Leo, may I know your pseudo code for the seeds inversion??
I hava confussions with the determination of the direction "where the seeds should grow?"
do you applied some constrains? or boundary for your model? (suppose, the seeds should not exceed certain boundaries or any other rule may be..)
or it is purely calculation such as from the gradient method (suppose gauss newton etc)
I mean, which one first, you move the next seeds and calculates the best movement from several movements before or
you calculate to determine which the best position or direction for the new seeds and then after that you put the next seeds there?

Thank you for your kindly help Leo

regards
acceph

Leonardo Uieda

unread,
Feb 18, 2019, 7:09:50 PM2/18/19
to fati...@googlegroups.com
Hi Acceph,

> Leo, may I know your pseudo code for the seeds inversion??

The source code for the seed inversion is at:
And the algorithm is described in these two papers:

> do you applied some constrains? or boundary for your model?

The boundaries are determined by the boundaries of the input mesh.
You can specify if you don't want seeds to grow in a given direction.
This is implemented by removing the neighbors in that direction from the search algorithm.

> or it is purely calculation such as from the gradient method

There is no gradient calculation in this method.

> I mean, which one first, you move the next seeds and calculates the best movement from several movements before or
> you calculate to determine which the best position or direction for the new seeds and then after that you put the next seeds there?

I'm not sure I understand what you mean by that.
Could you please clarify?


--
Leonardo Uieda

Visiting Research Scholar
Department of Earth Sciences - SOEST
University of Hawai'i at Mānoa
www.leouieda.com

acceph

unread,
Feb 19, 2019, 8:33:33 AM2/19/19
to Fatiando a Terra
Hello Leo

thank you for your reply and those several paper links..
I think I have found the answer inside your code. it is in harvester.harvest isn't it?
I still confused with the compactness which is incorporate inside this algorithm...
looks like more difficult for me rather than in the formal form of inversion

regards
acceph

Leonardo Uieda

unread,
Feb 19, 2019, 4:02:27 PM2/19/19
to fati...@googlegroups.com
It's definitely a different kind of algorithm.
I have to admit that the code in harvester is not the cleanest implementation.
I'm planning an updated version for Harmonica (https://github.com/fatiando/harmonica) which will be a lot better.

--
Leonardo Uieda

Visiting Research Scholar
Department of Earth Sciences - SOEST
University of Hawai'i at Mānoa
www.leouieda.com

acceph

unread,
Feb 21, 2019, 8:18:51 AM2/21/19
to Fatiando a Terra
thank you for the information. nice work leo, keep going...

Fabiano Della Justina

unread,
Apr 1, 2019, 11:41:06 AM4/1/19
to Fatiando a Terra
Hi Leo, 
How have you been?

I have had the same issue pointed out by Acceph - Something related to vremove.
However, I am using windows 10 and like him python 2.7.
Actually, the more I have tried to resolve the problem the worse it has become. The first time that I installed it, I ran your below code:


and I could see the data and the predicted data, but I couldn't get the model solution, then I tried to solve that, and now I can not see any plot, just errors.

I really like your method, it is really intuitive and I am going to present it in a final project to an inverse problem course. Although I was supposed to apply to real data, at this time I would be happy only apply in synthetic ones. Please see below the error that has popped out:


C:\Users\fjustina\Anaconda2\lib\site-packages\fatiando\vis\mpl.py:70: UserWarning: This module will be removed in v0.6. We recommend the use of matplotlib.pyplot module directly. Some of the fatiando specific functions will remain.
  "specific functions will remain.")
ImportErrorTraceback (most recent call last)
<ipython-input-1-981edd511661> in <module>()
      5 from fatiando import gridder, utils
      6 from fatiando.gravmag import polyprism, harvester
----> 7 from fatiando.mesher import PolygonalPrism, PrismMesh, vremove
      8 from fatiando.vis import mpl, myv
      9 

ImportError: cannot import name vremove

Please help me out,

Fabiano
To unsubscribe from this group and stop receiving emails from it, send an email to fati...@googlegroups.com.
Reply all
Reply to author
Forward
0 new messages