# with Fatiando a Terra

Leonardo Uieda

Vanderlei C. Oliveira Jr

ValĂ©ria C. F. Barbosa

# Geophysical modeling

## Understand Earth structure

Image credit: http://reinep.wordpress.com

Image credit: http://www.amazon.com/Look-inside-Earth-Poke/dp/0448418908

Image credits:

# Example

## Inverse problem

Linear algebra: $\hat{\mathbf{p}} = (\mathbf{G}^T\mathbf{G} + \mu \mathbf{I})^{-1}\mathbf{G}^T\mathbf{d}$

Numerical modeling: differential equations, integration, arrays, etc.

# 1. Write (copy/paste) C/Fortran

Computation:

for(k = 0; k < glq_lon.order; k++){
for(j = 0; j < glq_lat.order; j++){
for(i = 0; i < glq_r.order; i++){
rc = glq_r.nodes[i];
sinlatc = sin(d2r*glq_lat.nodes[j]);
coslatc = cos(d2r*glq_lat.nodes[j]);
coslon = cos(d2r*(lonp - glq_lon.nodes[k]));
sinlon = sin(d2r*(glq_lon.nodes[k] - lonp));
l_sqr = rp*rp + rc*rc - 2*rp*rc*(sinlatp*sinlatc +
coslatp*coslatc*coslon);


I/O:

switch(argv[i][1]){
case 'h':
if(argv[i][2] != '\0'){
log_error("invalid argument '%s'", argv[i]);
break;
}
print_help();


# 2. Write input file

Custom format (not human readable):

14

16
16
16

0.4 0   0.5

129.52861   38.6746
128.88952   25.92015
130.326     6.44616
116.52356   -5.8097
98.09756    -6.64001


# 3. Plot (2D/3D)

• Paraview
• Generic Mapping Tools (GMT)
• Matlab ® (proprietary)
• Surfer ® (proprietary)

Don't get me wrong, these are great software!

# Fatiando a Terra

"Slicing the Earth"

Automate

Integrate

API

DRY

Reproduce

# All inside Python

## 1. Write Python

def _shapefunc(data, predicted):
"""
Calculate the total shape of anomaly function between
the observed and predicted data.
"""
result = 0.
for d, p in zip(data, predicted):
alpha = numpy.sum(d.observed*p)/d.norm**2
result += numpy.linalg.norm(alpha*d.observed - p)
return result


Built-in I/O:

• pickle
• json

## 2. Input = Python script

In [5]:
from fatiando import gridder, gravmag
from fatiando.mesher import Prism
model = [Prism(-500,500,-500,500,200,1700,{'density':1000})]
area = [-1000, 1000, -1000, 1000]
bounds = area + [0, 2000]
x, y, z = gridder.scatter(area, 50, z=-1)
g = gravmag.prism.gz(x, y, z, model)
print g.shape, g[:20]

(50,) [ 3.53036009  3.35606669  8.35034244  6.87123912  6.03427323  2.16338594
3.77309787  4.83602868  8.31162035  2.01437273  5.73502438  3.061498
3.82377089  6.35807095  6.94394468  3.94985711  2.56862022  5.49552449
4.24486402  2.85385038]


## 3. Plot = same Python script

In [8]:
from fatiando.vis import mpl
mpl.figure(figsize=(8, 6))
mpl.axis('scaled')
mpl.contourf(y, x, g, (100, 100), 30, interp=True)
mpl.plot(y, x, '.k')
mpl.colorbar()

Out[8]:
<matplotlib.colorbar.Colorbar instance at 0x88aed88>

## 3D plots as well (with Mayavi)

In [23]:
from fatiando.vis import myv
f = myv.figure()
gray = (0.066,0.066,0.066); f.scene.background = gray
myv.prisms(model)
wht = (1,1,1)
myv.axes(myv.outline(bounds, color=wht), color=wht)
myv.wall_bottom(bounds, color=wht)
myv.wall_north(bounds, color=wht)
myv.show()


# What is available

• fatiando.gravmag
• 3D inversion
• FFT processing
• Modeling
• Equivalent layer
• fatiando.seismic
• Toy problems
• 2D Finite Differences wave propagation
• fatiando.geothermal
• Simple temperature log model

# What is available

• fatiando.mesher
• fatiando.gridder
• fatiando.utils
• fatiando.contants
• fatiando.io
• fatiando.vis
• fatiando.vis.mpl
• fatiando.vis.myv
• fatiando.inversion (experimental)

## Prism meshes

In [27]:
from fatiando.mesher import PrismMesh
mesh = PrismMesh(bounds, (20, 20, 20))
f = myv.figure(); f.scene.background = gray
p = myv.prisms(mesh, 'density')
myv.axes(p, color=wht); myv.show()


## Prism meshes with topography

In [25]:
from fatiando import utils
x, y = gridder.regular(bounds[:4], (50, 50))
heights = -1000 + 1000*utils.gaussian2d(x, y, 500, 500)
mesh.carvetopo(x, y, heights)


## Tesseroids

In [33]:
from fatiando.mesher import Tesseroid
model = [
Tesseroid(-60,-55,-30,-27,500000,0,{'density':200}),
Tesseroid(-66,-55,-20,-10,300000,0,{'density':-100})]


In [34]:
f = myv.figure(zdown=False); f.scene.background = gray
myv.tesseroids(model, 'density')
myv.continents(linewidth=2); myv.earth(opacity=1)
myv.meridians(range(0, 360, 45))
myv.parallels(range(-90, 90, 45))
myv.show()


Calculate gravitational fields (open field of research)

In [40]:
area = [-80, -30, -40, 10]
lons, lats, heights = gridder.scatter(area, 2500, z=2500000)
gz = gravmag.tesseroid.gz(lons, lats, heights, model)
print gz.shape, gz[:30]

(2500,) [ -1.4006202   -4.60404034 -35.68365773  -6.93419447  -2.54346594
-17.85438785  -2.56076061  -3.14379214 -10.35218711 -13.30011125
-5.60984378 -11.08639295 -17.70289672  -3.21422624 -10.03632971
-12.38263675  -6.26000011  -4.03988663  -9.30863759  -1.02347772
-17.49442248  -7.78201871 -16.51349434   1.55105749 -23.68919715
-4.56364226  -8.23905128 -16.30036808 -10.24673361 -14.18628464]


Map projections (with mpl_toolkits.basemap)

In [43]:
mpl.figure(figsize=(8, 6))
bm = mpl.basemap(area, 'ortho')
bm.drawcoastlines(); bm.drawmapboundary()
bm.bluemarble()
mpl.contourf(lons, lats, gz, (50, 50), 30, interp=True, basemap=bm)
mpl.colorbar()

Out[43]:
<matplotlib.colorbar.Colorbar instance at 0xef433f8>

# Useful for teaching

## How I'd rather teach

In [70]:
from IPython.display import Image
Image(filename="figures/crust.png", width=800)

Out[70]:
In [71]:
from fatiando import io
density = io.fromimage('figures/crust.png', ranges=[2700, 3300])
mpl.figure(figsize=(15, 4))
mpl.imshow(density, cmap=mpl.cm.seismic)
cb = mpl.colorbar(shrink=0.6)


## Conclusion

• Automate common tasks
• API for modeling
• For teaching and learning
• Use for my own research
• Masters and PhD all included
• Makes my life easier

## Present and future

• Modules available:
• Gravity and magnetic: fatiando.gravmag
• Seismic: fatiando.seismic
• Geothermal: fatiando.geothermal
• Future:
• More models
• Inverse problem: fatiando.inversion (experimental)
• Improve docs
• Community

## Open-source landscape

• Seismics and Seismology
• OpendTect
• Obspy
• SAC
• Geodynamics
• Computational Infrastructure for Geodynamics (CIG)
• SEATREE

See Wikipedia for more

Github: github.com/leouieda/fatiando

Slides: github.com/leouieda/scipy2013 (IPython notebook)

Homepage: fatiando.org (blog post coming soon)

## Tomography for dummies

In [65]:
from fatiando.mesher import SquareMesh
area = (0, 500000, 0, 500000); shape = (30, 30)
model = SquareMesh(area, shape)
model.img2prop('figures/fatiando-logo.png', 4000, 10000, 'vp')
mpl.figure(figsize=(8, 6)); mpl.axis('scaled')
mpl.squaremesh(model, 'vp', cmap=mpl.cm.seismic)
mpl.m2km()


Calculate straight-ray travel times

In [59]:
from fatiando import seismic
quake_locations = utils.random_points(area, 40)
20, random=True)
quakes, receivers = utils.connect_points(
traveltimes = seismic.ttime2d.straight(model, 'vp',
noisy = utils.contaminate(traveltimes, 0.1)
print noisy.shape, noisy[:20]

(800,) [ 24.75577443  59.0430238   20.16730169  32.69700566  15.61337341
21.85393683  58.51760914  10.77501476  35.42091179  54.52327289
10.67774583  62.65457186  66.53335128  46.79976704  37.11665425
13.22796347  39.11648141  23.50096277  14.70668023  62.65118329]

In [69]:
mpl.figure(figsize=(10, 8)); mpl.axis('scaled')
mpl.squaremesh(model, prop='vp', cmap=mpl.cm.seismic)
mpl.points(quakes, '*y', size=18)
mpl.m2km()


Run a toy tomography (not for research!)

In [63]:
slowness = seismic.srtomo.run(noisy, quakes, receivers, model, smooth=10**6)[0]
velocity = seismic.srtomo.slowness2vel(slowness)