Modeling the Earth with Fatiando a Terra

Geophysics is the science of using physical observations of the Earth to infer its inner structure. Generally, this is done with a variety of numerical modeling techniques and inverse problems. The development of new algorithms usually involves copy and pasting of code, which leads to errors and poor code reuse. Fatiando a Terra is a Python library that aims to automate common tasks and unify the modeling pipeline inside of the Python language. This allows users to replace the traditional shell scripting with more versatile and powerful Python scripting. The library can also be used as an API for developing stand-alone programs. Algorithms implemented in Fatiando a Terra can be combined to build upon existing functionality. This flexibility facilitates prototyping of new algorithms and quickly building interactive teaching exercises. In the future, we plan to continuously implement sample problems to help teach geophysics as well as classic and state-of-the-art algorithms.


Introduction
Geophysics studies the physical processes of the Earth. Geophysicists make observations of physical phenomena and use them to infer the inner structure of the planet. This task requires the numerical modeling of physical processes. These numerical models can then be used in inverse problems to infer inner Earth structure from observations. Different geophysical methods use different kinds of observations. Geothermal methods use the temperature and heat flux of the Earth's crust. Potential field methods use gravitational and magnetic field measurements. Seismics and seismology use the ground motion caused by elastic waves from active (manmade) and passive (earthquakes) sources, respectively.
The seismic method is among the most widely studied due to the high industry demand. Thus, a range of well established open-source software have been developed for seismic processing. These include Seismic Un*x (SU) [SU], Madagascar [MAD], OpendTect, and GêBR. A noteworthy open-source project that is not seismic related is the Generic Mapping Tools (GMT) project [GMT]. The GMT are a well established collection of command-line programs for plotting maps with a variety of different map projections. For geodynamic modeling there is the Computational Infrastructure for Geodynamics (CIG), which has grouped various well documented software packages. However, even with this wide range of well maintained software projects, many geophysical modeling software that are provided online still have no open-source license statement, have cryptic I/O files, are hard to integrate into a pipeline, and make code reuse and remixing challenging. Some of these problems are being worked on by the Solid Earth Teaching and Research Environment (SEATREE) [SEATREE] by providing a common graphical interface for previously existing software. The numerical computations are performed by the pre-existing underlying C/Fortran programs. Conversely, the SEATREE code (written in Python) handles the I/O and user interface. This makes the use of these tools easier and more approachable to students. However, the lack of a common API means that the code for these programs cannot be easily combined to create new modeling tools.
Fatiando a Terra aims at providing such an API for geophysical modeling. Functions in the fatiando package use compatible data and mesh formats so that the output of one modeling function can be used as input for another. Furthermore, routines can be combined and reused to create new modeling algorithms. Fatiando a Terra also automates common tasks such as griding, map plotting with Matplotlib [MPL], and 3D plotting with Mayavi [MYV]. Version 0.1 of Fatiando a Terra is focused on gravity and magnetic methods because this is the main focus of the developers. However, simple "toy" problems for seismology and geothermics are available and can be useful for teaching geophysics.
The following sections illustrate the functionality and design of Fatiando a Terra using various code samples. An IPython [IPY] notebook file with these code samples is provided by [SAMPLES] at http://dx.doi.org/10.6084/m9.figshare.708390.

Package structure
The modules and packages of Fatiando a Terra are bundled into the fatiando package. Each type of geophysical method has its own package. As of version 0.1, the available modules and packages are:  , 2) using that to make synthetic data, and 3) automatically gridding and plotting the data using a Fatiando a Terra wrapper for the Matplotlib contourf function.

Griding and map plotting
Fatiando a Terra handles map data as 1D Numpy arrays, typically x-, y-, z-coordinates and an extra array with the corresponding data. However, Matplotlib functions, like contourf and pcolor, require data to be passed as 2D arrays. Moreover, geophysical data sets are often irregularly sampled and require griding before they can be plotted. Thus, griding and array reshaping are ideal targets for automation.
Map projections in Matplotlib are handled by the Basemap toolkit. The fatiando.vis.mpl module also provides helper functions to automate the use of this toolkit. The fatiando.vis.mpl.basemap function automates the creation of the Basemap objects with common parameters. This object can then be passed to the contourf, contour and pcolor functions in fatiando.vis.mpl and they will automatically plot using the given projection ( Figure 2):

Meshes and 3D plotting
The representation of 2D and 3D geometric elements is handled by the classes in the fatiando.mesher module. Geometric elements in Fatiando a Terra can be assigned physical property values, like density, magnetization, seismic wave velocity, impedance, etc. This is done through a props dictionary whose keys are the name of the physical property and values are the corresponding values in SI units: The fatiando.vis.myv module contains functions to automate 3D plotting using Mayavi [MYV]. The mayavi.mlab interface requires geometric elements to be formatted as TVTK objects. Thus, plotting functions in fatiando.vis.myv automatically create TVTK representations of fatiando.mesher objects and plot them using a suitable function of mayavi.mlab. Also included are utility functions for drawing axes, walls on the figure bounding box, etc. For example, the fatiando.vis.myv.figure function creates a figure and rotates it so that the z-axis points down, as is standard in geophysics. The following example shows how to plot the 3D right rectangular prism model that we created previously (  The fatiando.mesher module also contains classes for collections of elements (e.g., meshes). A good example is the PrismMesh class that represents a structured mesh of right rectangular prisms. This class behaves as a list of fatiando.mesher.Prism objects and can be passed to functions that ask for a list of prisms, like fatiando.vis.myv.prisms. Physical properties can be assigned to the mesh using the addprop method ( Figure 4): mesh = mesher.PrismMesh(bounds, shape=(3, 3, 3)) mesh.addprop('density', range(mesh.size)) myv.figure() myv.prisms(mesh, 'density') myv.axes(myv.outline(bounds)) myv.show() Often times the mesh is used to make a detailed model of an irregular region of the Earth's surface. In such cases, it is necessary to consider the topography of the region. The PrismMesh class has a carvetopo method that masks the prisms that fall above the topography. The example below illustrates this functionality using synthetic topography ( Figure 5): from fatiando import utils x, y = gridder.regular(bounds[:4], (50, 50)) heights = -5 + 5 * utils.gaussian2d(x, y, 10, 5, x0=10, y0=10) mesh = mesher.PrismMesh (bounds, (20, 20, 20)) mesh.addprop('density', range(mesh.size)) mesh.carvetopo(x, y, heights) myv.figure() myv.prisms(mesh, 'density') myv.axes(myv.outline(bounds)) myv.wall_north(bounds) myv.show() When modeling involves the whole Earth, or a large area of it, the geophysicist needs to take into account the Earth's curvature. In such cases, rectangular prisms are inadequate for modeling and tesseroids (e.g., spherical prisms) are better suited. The fatiando.vis.myv module contains auxiliary functions to plot along with tesseroids: an Earth-sized sphere, meridians and parallels, as well as continental borders ( Figure 6):

Forward modeling
In geophysics, the term "forward modeling" is used to describe the process of generating synthetic data from a given Earth model. Conversely, geophysical inversion is the process of estimating Earth model parameters from observed data.

Gravity and magnetic methods
Geophysics uses anomalies in the gravitational and magnetic fields generated by density and magnetization contrasts within the Earth to investigate the inner Earth structure. The Fatiando a Terra 0.1 release has been focused on gravity and magnetic methods. Therefore, the fatiando.gravmag package contains more advanced and state-of-the-art algorithms than the other packages.
The module fatiando.gravmag.imaging implements the imaging methods described in [FP]. These methods aim to produce an image of the geologic source from the observed gravity or magnetic data. The following code sample uses the "sandwich   model" method [SNDW] to image the polygonal prism, produced in the previous section, based on its gravity anomaly ( Figure 10): estimate = gravmag.imaging.sandwich(x, y, z, gz, shape, zmin=0, zmax=1000, nlayers=20, power=0.2) body = mesher.vfilter(1.3 * 10 ** 8, 1.7 * 10 ** 8, 'density', estimate) myv.figure() myv.prisms(body, 'density', edges=False) p = myv.polyprisms(model, 'density', style='wireframe', linewidth=4) p.actor.mapper.scalar_visibility = False p.actor.property.color = (0, 0, 0) myv.axes(myv.outline(bounds), ranges=[i * 0.001 for i in bounds]) myv.wall_north(bounds) myv.wall_bottom(bounds) myv.show() Also implemented in Fatiando a Terra are some recent developments in gravity and magnetic inversion methods. The method of "planting anomalous densities" by [UB] is implemented in the fatiando.gravmag.harvester module. In contrast to imaging methods, this is an inversion method, i.e., it estimates a physical property distribution (density in the case of gravity data) that fits the observed data. This particular method requires the user to specify a "seed" (Figure 11) around which the estimated density distribution grows (Figure 12): # Make a mesh and a seed mesh = mesher.PrismMesh (bounds, (15, 30, 30 Fig. 11: The small blue prism is the seed used by fatiando.gravmag.harvester to perform the inversion of a gravity anomaly. The black contours are the true source of the gravity anomaly. Fig. 12: The blue prisms are the result of a gravity inversion using module fatiando.gravmag.harvester. The black contours are the true source of the gravity anomaly. Notice how the inversion was able to recover the approximate geometry of the true source.

Even
though the implementation in fatiando.seismic.srtomo is greatly simplified and not usable in real tomography problems, the result in Figure  13 illustrates interesting inverse problem concepts. Notice how the estimated velocity is blurred in the corners where no rays pass through. This is because the data (travel-times) provide no information about the velocity in those areas. Areas like those constitute the null space of the inverse problem [MENKE], where any velocity value estimated will provide an equal fit to the data. Thus, the tomography problem requires the use of prior information in the form of regularization. Most commonly used in tomography problems is the Tikhonov first-order regularization, e.g., a smoothness constraint [MENKE]. The amount of smoothness imposed on the solution is controlled by the smooth argument of function fatiando.seismic.srtomo.run. That is how we are able to estimate a unique and stable solution and why the result is specially smoothed where there are no rays.

Conclusion
The Fatiando a Terra package provides an API to develop modeling algorithms for a variety of geophysical methods. The current version (0.1) has a few state-of-the-art gravity and magnetic modeling and inversion algorithms. There are also toy problems in gravity, seismics and seismology that are useful for teaching basic concepts of geophysics, modeling, and inverse problems.
Fatiando a Terra enables quick prototyping of new algorithms because of the collection of fast forward modeling routines and the simple syntax and high level of the Python language. After prototyping, the performance bottlenecks of these algorithms can be easily diagnosed using the advanced profiling tools available in the Python language. Optimization of only small components of code can be done without loss of flexibility using the Cython language [CYTHON].
The biggest challenge that Fatiando a Terra faces in the near future is the development of a user and, consequently, a developer community. This is a key part for the survival of any open-source project.