In [1]:
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [2]:
from IPython.display import HTML, Image

Leonardo Uieda

Vanderlei C. Oliveira Jr

Valéria C. F. Barbosa

Earth structure

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

Can't look inside the earth directly

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

Indirect measurements

Image credits:

In [3]:
from fatiando import gravmag, gridder, mesher, utils
from fatiando.vis import myv
In [4]:
bounds = [0, 10000, 0, 10000, 0, 5000]
area = bounds[:4]
shape = (25, 25)
x, y, z = gridder.regular(area, shape, z=-1)
prop = {'density':1000}
model = [mesher.Prism(3000, 7000, 3000, 7000, 500, 2500, prop)]
gz = utils.contaminate(gravmag.prism.gz(x, y, z, model), 0.1)
mesh = mesher.PrismMesh(bounds, (20, 40, 40))
seeds = gravmag.harvester.sow([[5000, 5000, 1500, prop]], mesh)
data = [gravmag.harvester.Gz(x, y, z, gz)]    
estimate = gravmag.harvester.harvest(data, seeds, mesh, 0.1, 0.0005)[0]
mesh.addprop('density', estimate['density'])
body = mesher.vremove(0, 'density', mesh)
fig = myv.figure()
fig.scene.background = (0.06666, 0.06666, 0.06666)
myv.prisms(body, 'density')
myv.axes(myv.outline(bounds, color=(1,1,1)), nlabels=3, ranges=[0.001*b for b in bounds], fmt="%.0f", color=(1,1,1))
p = myv.mlab.contour_surf(reshape(x, shape, order='F'), reshape(y, shape, order='F'), -reshape(gz, shape, order='F'), 
    contours=100, colormap='gist_rainbow')
p.contour.filled_contours = True
p.actor.actor.position = (0, 0, -2000)
p.actor.actor.scale = (1, 1, 100)
a = myv.axes(p, color=(1,1,1), ranges=[0.001*b for b in area] + [gz.max(), gz.min()], fmt="%.0f", nlabels=3)
a.axes.z_label = ""
a.axes.x_label = ""
a.axes.y_label = ""
myv.wall_bottom(bounds, color=(1,1,1))
myv.wall_north(bounds, color=(1,1,1))
myv.show()
/usr/lib/python2.7/dist-packages/mayavi/preferences/preference_manager.py:24: UserWarning: Module docutils was already imported from /usr/lib/python2.7/dist-packages/docutils/__init__.pyc, but /home/leo/.virtualenvs/fatiando/lib/python2.7/site-packages is being added to sys.path
  import pkg_resources
/usr/lib/python2.7/dist-packages/mayavi/preferences/preference_manager.py:24: UserWarning: Module dap was already imported from None, but /usr/lib/python2.7/dist-packages is being added to sys.path
  import pkg_resources
WARNING:root:DEPRECATED: pyface.wx.clipboard, use pyface.api instead.

Geologic body

Estimate

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.

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]);
            bad_args++;
            break;
        }
        print_help();

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
  • Paraview
  • Generic Mapping Tools (GMT)
  • Matlab ® (proprietary)
  • Surfer ® (proprietary)

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

"Slicing the Earth"

Automate

Integrate

API

DRY

Reproduce

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:

  • numpy.loadtxt
  • pickle
  • json
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]
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>
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()

Output from above

  • fatiando.gravmag
    • 3D inversion
    • FFT processing
    • Modeling
    • Equivalent layer
  • fatiando.seismic
    • Toy problems
    • 2D Finite Differences wave propagation
  • fatiando.geothermal
    • Simple temperature log model
  • fatiando.mesher
  • fatiando.gridder
  • fatiando.utils
  • fatiando.contants
  • fatiando.io
  • fatiando.vis
    • fatiando.vis.mpl
    • fatiando.vis.myv
  • fatiando.inversion (experimental)
In [27]:
from fatiando.mesher import PrismMesh
mesh = PrismMesh(bounds, (20, 20, 20))
mesh.addprop('density', range(mesh.size))
f = myv.figure(); f.scene.background = gray
p = myv.prisms(mesh, 'density')
myv.axes(p, color=wht); myv.show()

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)
In []:
f = myv.figure(); f.scene.background = gray
myv.prisms(mesh, 'density')
myv.axes(myv.outline(bounds, color=wht), color=wht)
myv.wall_north(bounds, color=wht)
myv.show()

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>

Teaching seismic waves

Refracted

Reflected

Surface waves

Ripples

The Love wave gif

The example seismogram

In [70]:
from IPython.display import Image
Image(filename="figures/crust.png", width=800)
Out[70]:
In [7]:
mpl.rcParams['font.size'] = 16
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)
  • Automate common tasks
  • API for modeling
  • For teaching and learning
  • Use for my own research
    • Masters and PhD all included
    • Makes my life easier
  • Modules available:
    • Gravity and magnetic: fatiando.gravmag
    • Seismic: fatiando.seismic
    • Geothermal: fatiando.geothermal
  • Future:
    • More models
    • Inverse problem: fatiando.inversion (experimental)
    • Improve docs
    • Community
  • Seismics and Seismology
    • Madagascar
    • 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)

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.colorbar(pad=0.01)
mpl.m2km()

Calculate straight-ray travel times

In [59]:
from fatiando import seismic
quake_locations = utils.random_points(area, 40)
receiver_locations = utils.circular_points(area, 
    20, random=True)
quakes, receivers = utils.connect_points(
    quake_locations, receiver_locations)
traveltimes = seismic.ttime2d.straight(model, 'vp', 
    quakes, receivers)
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.paths(quakes, receivers)
mpl.points(quakes, '*y', size=18)
mpl.points(receivers, '^g', 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)
model.addprop('vp', velocity)
mpl.figure(figsize=(8, 6)); mpl.axis('scaled')
mpl.squaremesh(model, 'vp', cmap=mpl.cm.seismic, 
    vmin=4000, vmax=10000)
mpl.colorbar(pad=0.01)
mpl.m2km()