3D inversion by planting anomalous densities (fatiando.gravmag.harvester)

3D potential field inversion by planting anomalous densities.

Implements the method of Uieda and Barbosa (2012a) with improvements by Uieda and Barbosa (2012b).

A “heuristic” inversion for compact 3D geologic bodies. Performs the inversion by iteratively growing the estimate around user-specified “seeds”. Supports various kinds of data (gravity, gravity tensor).

The inversion is performed by function harvest. The required information, such as observed data, seeds, and regularization, are passed to the function through classes Seed and Potential, Gz, Gxx, etc.

See the Cookbook for some example applications to synthetic data.

Functions

  • harvest: Performs the inversion
  • iharvest: Iterator to step through the inversion one accretion at a time
  • sow: Creates the seeds from a set of (x, y, z) points and physical properties
  • loadseeds: Loads from a JSON file a set of (x, y, z) points and physical properties that specify the seeds. Pass output to sow
  • weights: Computes data weights based on the distance to the seeds

Data types

  • Potential: gravitational potential
  • Gz: vertical component of gravitational acceleration (i.e., gravity anomaly)
  • Gxx: North-North component of the gravity gradient tensor
  • Gxy: North-East component of the gravity gradient tensor
  • Gxz: North-vertical component of the gravity gradient tensor
  • Gyy: East-East component of the gravity gradient tensor
  • Gyz: East-vertical component of the gravity gradient tensor
  • Gzz: vertical-vertical component of the gravity gradient tensor

References

Uieda, L., and V. C. F. Barbosa (2012a), Robust 3D gravity gradient inversion by planting anomalous densities, Geophysics, 77(4), G55-G66, doi:10.1190/geo2011-0388.1

Uieda, L., and V. C. F. Barbosa (2012b), Use of the “shape-of-anomaly” data misfit in 3D inversion by planting anomalous densities, SEG Technical Program Expanded Abstracts, 1-6, doi:10.1190/segam2012-0383.1


class fatiando.gravmag.harvester.Data(x, y, z, data, weights, meshtype)[source]

Bases: object

A container for some potential field data.

Know about its data, observation positions, nature of the mesh, and how to calculate the effect of a single cell.

class fatiando.gravmag.harvester.Gxx(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the xx (north-north) component of the gravity gradient tensor.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Gxy(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the xy (north-east) component of the gravity gradient tensor.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Gxz(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the xz (north-vertical) component of the gravity gradient tensor.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Gyy(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the yy (east-east) component of the gravity gradient tensor.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Gyz(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the yz (east-vertical) component of the gravity gradient tensor.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Gz(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the gravity anomaly.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Gzz(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the zz (vertical-vertical) component of the gravity gradient tensor.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.Neighbor(i, props, seed, distance, effect)[source]

Bases: object

A neighbor.

class fatiando.gravmag.harvester.Potential(x, y, z, data, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Data

A container for data of the gravitational potential.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

class fatiando.gravmag.harvester.PrismSeed(i, location, prism, props)[source]

Bases: fatiando.mesher.Prism

A seed that is a right rectangular prism.

addprop(prop, value)

Add a physical property to this geometric element.

If it already has the property, the given value will overwrite the existing one.

Parameters:

  • prop
    : str

    Name of the physical property.

  • value
    : float

    The value of this physical property.

center()

Return the coordinates of the center of the prism.

Returns:

  • coords
    : list = [xc, yc, zc]

    Coordinates of the center

Example:

>>> prism = Prism(1, 2, 1, 3, 0, 2)
>>> print prism.center()
[ 1.5  2.   1. ]
copy()

Return a deep copy of the current instance.

get_bounds()

Get the bounding box of the prism (i.e., the borders of the prism).

Returns:

  • bounds
    : list

    [x1, x2, y1, y2, z1, z2], the bounds of the prism

Examples:

>>> p = Prism(1, 2, 3, 4, 5, 6)
>>> print p.get_bounds()
[1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
class fatiando.gravmag.harvester.TesseroidSeed(i, location, tess, props)[source]

Bases: fatiando.mesher.Tesseroid

A seed that is a tesseroid (spherical prism).

addprop(prop, value)

Add a physical property to this geometric element.

If it already has the property, the given value will overwrite the existing one.

Parameters:

  • prop
    : str

    Name of the physical property.

  • value
    : float

    The value of this physical property.

copy()

Return a deep copy of the current instance.

get_bounds()

Get the bounding box of the tesseroid (i.e., the borders).

Returns:

  • bounds
    : list

    [w, e, s, n, top, bottom], the bounds of the tesseroid

Examples:

>>> t = Tesseroid(1, 2, 3, 4, 6, 5)
>>> print t.get_bounds()
[1.0, 2.0, 3.0, 4.0, 6.0, 5.0]
half(lon=True, lat=True, r=True)

Divide the tesseroid in 2 halfs for each dimension (total 8)

The smaller tesseroids will share the large one’s props.

Parameters:

  • lon, lat, r
    : True or False

    Dimensions along which the tesseroid will be split in half.

Returns:

  • tesseroids
    : list

    A list of maximum 8 tesseroids that make up the larger one.

Examples:

>>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2})
>>> split = tess.half()
>>> print len(split)
8
>>> for t in split:
...     print t
w:-10 | e:0 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
w:-10 | e:0 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
w:-10 | e:0 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
w:-10 | e:0 | s:0 | n:20 | top:0 | bottom:-20 | density:2
w:0 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
w:0 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
w:0 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
w:0 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2
>>> tess = Tesseroid(-15, 15, -20, 20, 0, -40)
>>> split = tess.half(lat=False)
>>> print len(split)
4
>>> for t in split:
...     print t
w:-15 | e:0 | s:-20 | n:20 | top:-20 | bottom:-40
w:-15 | e:0 | s:-20 | n:20 | top:0 | bottom:-20
w:0 | e:15 | s:-20 | n:20 | top:-20 | bottom:-40
w:0 | e:15 | s:-20 | n:20 | top:0 | bottom:-20
split(nlon, nlat, nh)

Split the tesseroid into smaller ones.

The smaller tesseroids will share the large one’s props.

Parameters:

  • nlon, nlat, nh
    : int

    The number of sections to split in the longitudinal, latitudinal, and vertical dimensions

Returns:

  • tesseroids
    : list

    A list of nlon*nlat*nh tesseroids that make up the larger one.

Examples:

>>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2})
>>> split = tess.split(1, 2, 2)
>>> print len(split)
4
>>> for t in split:
...     print t
w:-10 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2
w:-10 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2
w:-10 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2
w:-10 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2
>>> tess = Tesseroid(-15, 15, -20, 20, 0, -40)
>>> split = tess.split(3, 1, 1)
>>> print len(split)
3
>>> for t in split:
...     print t
w:-15 | e:-5 | s:-20 | n:20 | top:0 | bottom:-40
w:-5 | e:5 | s:-20 | n:20 | top:0 | bottom:-40
w:5 | e:15 | s:-20 | n:20 | top:0 | bottom:-40
class fatiando.gravmag.harvester.TotalField(x, y, z, data, inc, dec, weights=1.0, meshtype='prism')[source]

Bases: fatiando.gravmag.harvester.Potential

A container for data of the total field magnetic anomaly.

Coordinate system used: x->North y->East z->Down

Parameters:

  • x, y, z
    : 1D arrays

    Arrays with the x, y, z coordinates of the data points

  • data
    : 1D array

    The values of the data at the observation points

  • inc, dec
    : floats

    The inclination and declination of the inducing field

  • weight
    : float or array

    The weight of this data set in the misfit function. Pass an array to give weights to each data points or a float to weight the entire misfit function. See function weights

fatiando.gravmag.harvester.fmt_estimate(estimate, size)[source]

Make a nice dict with the estimated physical properties in separate arrays

fatiando.gravmag.harvester.harvest(data, seeds, mesh, compactness, threshold, report=False)[source]

Run the inversion algorithm and produce an estimate physical property distribution (density and/or magnetization).

Parameters:

  • data
    : list of data (e.g., Gz)

    The data that will be inverted. Data used must match the physical properties given to the seeds (e.g., gravity data requires seeds to have 'density' prop)

  • seeds
    : list of Seed

    Lits of seeds used to start the growth process of the inversion. Use sow to generate seeds.

  • mesh
    : fatiando.mesher.PrismMesh

    The mesh used in the inversion. Will estimate the physical property distribution on this mesh

  • compactness
    : float

    The compactness regularing parameter (i.e., how much should the estimate be consentrated around the seeds). Must be positive. To find a good value for this, start with a small value (like 0.001), run the inversion and increase the value until desired compactness is achieved.

  • threshold
    : float

    Control how much the solution can grow (usually downward). In order for estimate to grow by the accretion of 1 prism, this prism must decrease the data misfit measure by threshold decimal percent. Depends on the size of the cells in the mesh and the distance from a cell to the observations. Use values between 0.001 and 0.000001. If cells are small and threshold is large (0.001), the seeds won’t grow. If cells are large and threshold is small (0.000001), the seeds will grow too much.

  • report
    : True or False

    If True, also will return a dict as:

    report = {'goal': goal_function_value,
              'shape-of-anomaly': SOA_function_value,
              'misfit': data_misfit_value,
              'regularizer': regularizing_function_value,
              'accretions': number_of_accretions}
    

Returns:

  • estimate, predicted_data
    : a dict and a list

    estimate is a dict like:

    {'physical_property':array, ...}
    

    estimate contains the estimates physical properties. The properties present in estimate are the ones given to the seeds. Include the properties in the mesh using:

    mesh.addprop('density', estimate['density'])
    

    This way you can plot the estimate using fatiando.vis.myv.

    predicted_data is a list of numpy arrays with the predicted (model) data. The list is in the same order as data. To plot a map of the fit for visual inspection and a histogram of the residuals:

    from fatiando.vis import mpl
    mpl.figure()
    # Plot the observed and predicted data as contours for visual
    # inspection
    mpl.subplot(1, 2, 1)
    mpl.axis('scaled')
    mpl.title('Observed and predicted data')
    levels = mpl.contourf(x, y, gz, (ny, nx), 10)
    mpl.colorbar()
    # Assuming gz is the only data used
    mpl.contour(x, y, predicted[0], (ny, nx), levels)
    # Plot a histogram of the residuals
    residuals = gz - predicted[0]
    mpl.subplot(1, 2, 2)
    mpl.title('Residuals')
    mpl.hist(residuals, bins=10)
    mpl.show()
    # It's also good to see the mean and standard deviation of the
    # residuals
    print "Residuals mean:", residuals.mean()
    print "Residuals stddev:", residuals.std()
    
fatiando.gravmag.harvester.iharvest(data, seeds, mesh, compactness, threshold)[source]

Same as the fatiando.gravmag.harvester.harvest function but this one returns an iterator that yields the information of each accretion.

Yields:

  • [estimate, predicted, new, neighbors, goal, misfit, regularizer]

    The unformated estimate, predicted data vectors, the new element added during this iteration, list of neighbors, goal function value, misfit, regularizing function value.

The first yield contains the seeds. Thus new will be None.

To format the estimate in a way that can be added to a mesh, use function fmt_estimate of this module.

fatiando.gravmag.harvester.loadseeds(fname)[source]

Load a set of seed locations and physical properties from a file.

The output can then be used with the sow function.

The seed file should be formatted as:

[
    [x1, y1, z1, {"density":dens1}],
    [x2, y2, z2, {"density":dens2, "magnetization":mag2}],
    [x3, y3, z3, {"magnetization":mag3, "inclination":inc3,
                  "declination":dec3}],
    ...
]

x, y, z are the coordinates of the seed and the dict ({'density':2670}) are its physical properties.

Warning

Must use ", not ', in the physical property names!

Each seed can have different kinds of physical properties. If inclination and declination are not given, will use the inc and dec of the inducing field (i.e., no remanent magnetization).

The techie among you will recognize that the seed file is in JSON format.

Remember: the coordinate system is x->North, y->East, and z->Down

Parameters:

  • fname
    : str or file

    Open file object or filename string

Returns:

  • [[x1, y1, z1, props1], [x2, y2, z2, props2], ...]

    (x, y, z) are the points where the seeds will be placed and props is dict with the values of the physical properties of each, seed.

Example:

>>> from StringIO import StringIO
>>> file = StringIO(
...     '[[1, 2, 3, {"density":4, "magnetization":5}],' +
...     ' [6, 7, 8, {"magnetization":-1}]]')
>>> seeds = loadseeds(file)
>>> for s in seeds:
...     print s
[1, 2, 3, {u'magnetization': 5, u'density': 4}]
[6, 7, 8, {u'magnetization': -1}]
fatiando.gravmag.harvester.sow(locations, mesh)[source]

Create the seeds given a list of (x,y,z) coordinates and physical properties.

Removes seeds that would fall on the same location with overlapping physical properties.

Parameters:

  • locations
    : list

    The locations and physical properties of the seeds. Should be a list like:

    [
        [x1, y1, z1, {"density":dens1}],
        [x2, y2, z2, {"density":dens2, "magnetization":mag2}],
        [x3, y3, z3, {"magnetization":mag3, "inclination":inc3,
                      "declination":dec3}],
        ...
    ]
    
  • mesh
    : fatiando.mesher.PrismMesh

    The mesh that will be used in the inversion.

Returns:

  • seeds
    : list of seeds

    The seeds that can be passed to harvest

fatiando.gravmag.harvester.weights(x, y, seeds, influences, decay=2)[source]

Calculate weights for the data based on the distance to the seeds. Use weights to ignore regions of data outside of the target anomaly.

Parameters:

  • x, y
    : 1d arrays

    The x and y coordinates of the observations

  • seeds
    : list

    List of seeds, as returned by sow

  • influences
    : list of floats

    The respective diameter of influence for each seed. Observations outside the influence will have very small weights. A recommended value is aproximately the diameter of the anomaly

  • decay
    : float

    The decay factor for the weights. Low decay factor makes the weights spread out more. High decay factor makes the transition from large weights to low weights more abrupt.

Returns:

  • weights
    : 1d array

    The calculated weights