fatiandofatiando.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 inversioniharvest: Iterator to step through the
inversion one accretion at a timesow: Creates the seeds from a set of
(x, y, z) points and physical propertiesloadseeds: Loads from a JSON file a set
of (x, y, z) points and physical properties that specify the seeds. Pass
output to sowweights: Computes data weights based on
the distance to the seedsData types
Potential: gravitational potentialGz: vertical component of gravitational
acceleration (i.e., gravity anomaly)Gxx: North-North component of the
gravity gradient tensorGxy: North-East component of the gravity
gradient tensorGxz: North-vertical component of the
gravity gradient tensorGyy: East-East component of the gravity
gradient tensorGyz: East-vertical component of the
gravity gradient tensorGzz: vertical-vertical component of the
gravity gradient tensorReferences
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
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.
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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.Neighbor(i, props, seed, distance, effect)[source]¶Bases: object
A neighbor.
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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
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.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:
Name of the physical property.
The value of this physical property.
center()¶Return the coordinates of the center of the prism.
Returns:
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:
[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]
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:
Name of the physical property.
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:
[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:
Dimensions along which the tesseroid will be split in half.
Returns:
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:
The number of sections to split in the longitudinal, latitudinal, and vertical dimensions
Returns:
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
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:
Arrays with the x, y, z coordinates of the data points
The values of the data at the observation points
The inclination and declination of the inducing field
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:
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)
SeedLits of seeds used to start the growth process of the inversion. Use
sow to generate seeds.
fatiando.mesher.PrismMeshThe mesh used in the inversion. Will estimate the physical property distribution on this mesh
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.
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.
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 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:
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:
Open file object or filename string
Returns:
(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:
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}],
...
]
fatiando.mesher.PrismMeshThe mesh that will be used in the inversion.
Returns:
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:
The x and y coordinates of the observations
List of seeds, as returned by sow
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
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:
The calculated weights