The fatiando package has been deprecated. Please check out the new tools in the Fatiando a Terra website: www.fatiando.org

Source code for fatiando.gravmag.imaging

"""
Imaging methods for potential fields.

Implements some of the methods described in Fedi and Pilkington (2012).
Most methods convert the observed data (gravity, magnetic, etc) into a physical
property distribution (density, magnetization, etc). Most methods require
gridded data to work.

* :func:`~fatiando.gravmag.imaging.geninv`: The Generalized Inverse solver in
  the frequency domain (Cribb, 1976)
* :func:`~fatiando.gravmag.imaging.sandwich`: Sandwich model (Pedersen, 1991).
  Uses depth weighting as in Pilkington (1997)
* :func:`~fatiando.gravmag.imaging.migrate`: 3D potential field migration
  (Zhdanov et al., 2011). Actually uses the formula of Fedi and Pilkington
  (2012), which are comprehensible.

.. warning::

    Most of these methods provide estimates of physical property values that
    are completely out of scale (mostly due to depth weighting). Therefore, I
    don't recommend using the actual values of the physical properties for
    anything other than finding an approximate location for the sources.

.. note::

    If you want the estimate physical property values in SI units, you
    must pass the data also in SI units! Use the unit conversion functions in
    :mod:`fatiando.utils`


**References**

Cribb, J. (1976), Application of the generalized linear inverse to the
inversion of static potential data, Geophysics, 41(6), 1365,
doi:10.1190/1.1440686

Fedi, M., and M. Pilkington (2012), Understanding imaging methods for potential
field data, Geophysics, 77(1), G13, doi:10.1190/geo2011-0078.1

Pedersen, L. B. (1991), Relations between potential fields and some equivalent
sources, Geophysics, 56(7), 961, doi:10.1190/1.1443129

Pilkington, M. (1997), 3-D magnetic imaging using conjugate gradients,
Geophysics, 62(4), 1132, doi:10.1190/1.1444214

Zhdanov, M. S., X. Liu, G. A. Wilson, and L. Wan (2011), Potential field
migration for rapid imaging of gravity gradiometry data, Geophysical
Prospecting, 59(6), 1052-1071, doi:10.1111/j.1365-2478.2011.01005.x

----
"""
import numpy

from fatiando.mesher import PrismMesh
from fatiando.gravmag import transform
from fatiando.gravmag import prism as pot_prism
from fatiando.constants import G
from fatiando import utils


[docs]def migrate(x, y, z, gz, zmin, zmax, meshshape, power=0.5, scale=1): """ 3D potential field migration (Zhdanov et al., 2011). Actually uses the formula of Fedi and Pilkington (2012), which are comprehensible. .. note:: Only works on **gravity** data for now. .. note:: The data **do not** need to be leveled or on a regular grid. .. note:: The coordinate system adopted is x->North, y->East, and z->Down Parameters: * x, y : 1D-arrays The x and y coordinates of the grid points * z : float or 1D-array The z coordinate of the grid points * gz : 1D-array The gravity anomaly data at the grid points * zmin, zmax : float The top and bottom, respectively, of the region where the physical property distribution is calculated * meshshape : tuple = (nz, ny, nx) Number of prisms in the output mesh in the x, y, and z directions, respectively * power : float The power law used for the depth weighting. This controls what depth the bulk of the solution will be. * scale : float A scale factor for the depth weights. Simply changes the scale of the physical property values. Returns: * mesh : :class:`fatiando.mesher.PrismMesh` The estimated physical property distribution set in a prism mesh (for easy 3D plotting) """ nlayers, ny, nx = meshshape mesh = _makemesh(x, y, (ny, nx), zmin, zmax, nlayers) # This way, if z is not an array, it is now z = z * numpy.ones_like(x) dx, dy, dz = mesh.dims # Synthetic tests show that its not good to offset the weights with the # data z coordinate. No idea why depths = mesh.get_zs()[:-1] + 0.5 * dz weights = numpy.abs(depths) ** power / (2 * G * numpy.sqrt(numpy.pi)) density = [] for l in xrange(nlayers): sensibility_T = numpy.array( [pot_prism.gz(x, y, z, [p], dens=1) for p in mesh.get_layer(l)]) density.extend(scale * weights[l] * numpy.dot(sensibility_T, gz)) mesh.addprop('density', numpy.array(density)) return mesh
[docs]def sandwich(x, y, z, data, shape, zmin, zmax, nlayers, power=0.5): """ Sandwich model (Pedersen, 1991). Calculates a physical property distribution given potential field data on a **regular grid**. Uses depth weights. .. note:: Only works on **gravity** data for now. .. note:: The data **must** be leveled, i.e., on the same height! .. note:: The coordinate system adopted is x->North, y->East, and z->Down Parameters: * x, y : 1D-arrays The x and y coordinates of the grid points * z : float or 1D-array The z coordinate of the grid points * data : 1D-array The potential field at the grid points * shape : tuple = (ny, nx) The shape of the grid * zmin, zmax : float The top and bottom, respectively, of the region where the physical property distribution is calculated * nlayers : int The number of layers used to divide the region where the physical property distribution is calculated * power : float The power law used for the depth weighting. This controls what depth the bulk of the solution will be. Returns: * mesh : :class:`fatiando.mesher.PrismMesh` The estimated physical property distribution set in a prism mesh (for easy 3D plotting) """ mesh = _makemesh(x, y, shape, zmin, zmax, nlayers) # This way, if z is not an array, it is now z = z * numpy.ones_like(x) freq, dataft = _getdataft(x, y, data, shape) dx, dy, dz = mesh.dims # Remove the last z because I only want depths to the top of the layers depths = mesh.get_zs()[:-1] weights = (numpy.abs(depths) + 0.5 * dz) ** (power) density = [] # Offset by the data z because in the paper the data is at z=0 for depth, weight in zip(depths - z[0], weights): density.extend( numpy.real(numpy.fft.ifft2( weight * (numpy.exp(-freq * depth) - numpy.exp(-freq * (depth + dz))) * freq * dataft / (numpy.pi * G * reduce(numpy.add, [w * (numpy.exp(-freq * h) - numpy.exp(-freq * (h + dz))) ** 2 # To avoid zero division when freq[i]==0 + 10. ** (-10) for h, w in zip(depths, weights)]) ) ).ravel())) mesh.addprop('density', numpy.array(density)) return mesh
[docs]def geninv(x, y, z, data, shape, zmin, zmax, nlayers): """ Generalized Inverse imaging in the frequency domain (Cribb, 1976). Calculates a physical property distribution given potential field data on a **regular grid**. .. note:: Only works on **gravity** data for now. .. note:: The data **must** be leveled, i.e., on the same height! .. note:: The coordinate system adopted is x->North, y->East, and z->Down .. warning:: The Generalized Inverse does **not** use depth weights. This means that the solution will tend to be concentrated on the surface! Parameters: * x, y : 1D-arrays The x and y coordinates of the grid points * z : float or 1D-array The z coordinate of the grid points * data : 1D-array The potential field at the grid points * shape : tuple = (ny, nx) The shape of the grid * zmin, zmax : float The top and bottom, respectively, of the region where the physical property distribution is calculated * nlayers : int The number of layers used to divide the region where the physical property distribution is calculated Returns: * mesh : :class:`fatiando.mesher.PrismMesh` The estimated physical property distribution set in a prism mesh (for easy 3D plotting) """ mesh = _makemesh(x, y, shape, zmin, zmax, nlayers) # This way, if z is not an array, it is now z = z * numpy.ones_like(x) freq, dataft = _getdataft(x, y, data, shape) dx, dy, dz = mesh.dims # Remove the last z because I only want depths to the top of the layers depths = mesh.get_zs()[:-1] + 0.5 * dz - z[0] # Offset by the data height density = [] for depth in depths: density.extend( numpy.real( numpy.fft.ifft2( numpy.exp(-freq * depth) * freq * dataft / (numpy.pi * G) ).ravel() )) mesh.addprop('density', numpy.array(density)) return mesh
def _getdataft(x, y, data, shape): """ Get the Fourier transform of the data and the norm of the wavenumber vector """ Fx, Fy = transform._getfreqs(x, y, data, shape) freq = numpy.sqrt(Fx ** 2 + Fy ** 2) dataft = (2. * numpy.pi) * numpy.fft.fft2(numpy.reshape(data, shape)) return freq, dataft def _makemesh(x, y, shape, zmin, zmax, nlayers): """ Make a prism mesh bounded by the data. """ ny, nx = shape bounds = [x.min(), x.max(), y.min(), y.max(), zmin, zmax] mesh = PrismMesh(bounds, (nlayers, ny, nx)) return mesh