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.basin2d

"""
2D inversion for the basement relief of sedimentary basins.

There are different parametrizations available.

The simplest are meant more as an exercise and initiation in inverse problems:

* :func:`~fatiando.gravmag.basin2d.Triangular`: assumes a basin with a
  triangular cross-section (think "foreland").
* :func:`~fatiando.gravmag.basin2d.Trapezoidal`: assumes a basin with a
  trapezoidal cross-section (think "grabben").

More complex parametrizations are:

* :func:`~fatiando.gravmag.basin2d.PolygonalBasinGravity`: approximate the
  basin by a polygon.

----

"""
from __future__ import division
from future.builtins import super
import numpy as np

from ..inversion.misfit import Misfit
from . import talwani
from ..mesher import Polygon
from .. import utils


[docs]class PolygonalBasinGravity(Misfit): """ Estimate the relief of a sedimentary basin approximating by a polygon. Currently only works with gravity data. The top of the basin is straight and fixed at a given height. Polygon vertices are distributed evenly in the x-direction. The inversion estimates the depths of each vertex. This is a non-linear inversion. Therefore you must configure it before running to choose a solver method and set the initial estimate. Use the ``config`` method for this. Recommended configuration: Levemberg-Marquardt algorithm (``'levmarq'``) with initial estimate to the average expected depth of the basin. Typical regularization to use with this class are: :class:`~fatiando.inversion.regularization.Damping`, :class:`~fatiando.inversion.regularization.Smoothness1D`, :class:`~fatiando.inversion.regularization.TotalVariation1D`. The forward modeling is done using :mod:`~fatiando.gravmag.talwani`. Derivatives are calculated using a 2-point finite difference approximation. .. tip:: Use the ``estimate_`` attribute to get a :class:`~fatiando.mesher.Polygon` version of the estimated parameters (attribute ``p_``). Parameters: * x, z : 1d-arrays The x and z coordinates of the observations. In meters. * data : 1d-array The observed data. * npoints : int Number of points to use * props : dict The physical properties dictionary that will be assigned to the basin :class:`~fatiando.mesher.Polygon`. Ex: to give the basin a density contrast of 500 kg/m3 ``props={'density': 500}``. * top : float The value of the z-coordinate where the top of the basin will be fixed. In meters. Default: 0. * xlim : None or list = [xmin, xmax] The horizontal limits of the model. If not given, will use the limits of the data (i.e., ``[x.min(), x.max()]``). Examples: Lets run an inversion on synthetic data from a simple model of a trapezoid basin (a polygon with 4 vertices). We'll assume that the horizontal limits of the basin are the same as the limits of the data: >>> from fatiando.mesher import Polygon >>> from fatiando.gravmag import talwani >>> import numpy as np >>> # Make some synthetic data from a simple basin >>> props = {'density': -500} >>> model = [Polygon([[3000, 0], [2000, 800], [1000, 500], [0, 0]], ... props)] >>> x = np.linspace(0, 3000, 50) >>> z = -np.ones_like(x) # Put data at 1m height >>> data = talwani.gz(x, z, model) >>> # Make the solver, configure, and invert. >>> # Will use only 2 points because the two in the corners are >>> # considered fixed in the inversion (at 'top'). >>> misfit = PolygonalBasinGravity(x, z, data, 2, props, top=0) >>> _ = misfit.config('levmarq', initial=100*np.ones(misfit.nparams)).fit() >>> misfit.p_ array([ 800., 500.]) >>> type(misfit.estimate_) <class 'fatiando.mesher.Polygon'> >>> misfit.estimate_.vertices array([[ 3000., 0.], [ 2000., 800.], [ 1000., 500.], [ 0., 0.]]) If the x range of the data points is larger than the basin, you can specify a horizontal range for the basin model. When this is not specified, it is deduced from the data: >>> x = np.linspace(-500, 3500, 80) >>> z = -np.ones_like(x) >>> data = talwani.gz(x, z, model) >>> # Specify that the model used for inversion should be within >>> # x => [0, 3000] >>> misfit = PolygonalBasinGravity(x, z, data, 2, props, top=0, ... xlim=[0, 3000]) >>> _ = misfit.config('levmarq', initial=100*np.ones(misfit.nparams)).fit() >>> misfit.p_ array([ 800., 500.]) >>> misfit.estimate_.vertices array([[ 3000., 0.], [ 2000., 800.], [ 1000., 500.], [ 0., 0.]]) """ def __init__(self, x, z, data, npoints, props, top=0, xlim=None): super().__init__(data=data.ravel(), nparams=npoints, islinear=False) self.npoints = npoints self.x = x self.z = z self.props = props self.top = top if xlim is None: xlim = [x.min(), x.max()] self.xlim = xlim self._modelx = np.linspace(xlim[0], xlim[1], npoints + 2)[::-1]
[docs] def p2vertices(self, p): """ Convert a parameter vector into vertices a Polygon. Parameters: * p : 1d-array The parameter vector with the depth of the polygon vertices Returns: * vertices : 2d-array Like a list of [x, z] coordinates of each vertex Examples: >>> import numpy as np >>> # Make some arrays to create the estimator clas >>> x = np.linspace(-100, 300, 50) >>> z = np.zeros_like(x) >>> data = z >>> misfit = PolygonalBasinGravity(x, z, data, 3, {}, top=-100) >>> misfit.p2vertices([1, 2, 3]) array([[ 300., -100.], [ 200., 1.], [ 100., 2.], [ 0., 3.], [-100., -100.]]) """ h = self.top verts = np.empty((self.nparams + 2, 2)) verts[:, 0] = self._modelx verts[:, 1] = np.concatenate([[h], p, [h]]) return verts
[docs] def predicted(self, p): """ Calculate the predicted data for a parameter vector. """ verts = self.p2vertices(p) poly = Polygon(verts, self.props) return talwani.gz(self.x, self.z, [poly])
[docs] def jacobian(self, p): """ Calculate the Jacobian (sensitivity) matrix for a parameter vector. """ verts = self.p2vertices(p) delta = np.array([0, 1]) jac = np.empty((self.ndata, self.nparams)) for i in xrange(self.nparams): diff = Polygon([verts[i + 2], verts[i + 1] - delta, verts[i], verts[i + 1] + delta], self.props) jac[:, i] = talwani.gz(self.x, self.z, [diff])/(2*delta[1]) return jac
[docs] def fmt_estimate(self, p): """ Convert the parameter vector to a :class:`fatiando.mesher.Polygon` so that it can be used for plotting and forward modeling. Examples: >>> import numpy as np >>> # Make some arrays to create the estimator clas >>> x = np.linspace(-100, 300, 50) >>> z = np.zeros_like(x) >>> data = z >>> misfit = PolygonalBasinGravity(x, z, data, 3, {}, top=-100) >>> poly = misfit.fmt_estimate([1, 2, 3]) >>> type(poly) <class 'fatiando.mesher.Polygon'> >>> poly.vertices array([[ 300., -100.], [ 200., 1.], [ 100., 2.], [ 0., 3.], [-100., -100.]]) """ poly = Polygon(self.p2vertices(p), self.props) return poly
[docs]class Triangular(Misfit): """ Estimate the relief of a triangular basin. Use when the basin can be approximated by a 2D body with **triangular** vertical cross-section, like foreland basins. The triangle is assumed to have 2 known vertices at the surface (the edges of the basin) and one unknown vertex in the subsurface. The inversion will estimate the (x, z) coordinates of the unknown vertex. The forward modeling is done using :mod:`~fatiando.gravmag.talwani`. Derivatives are calculated using a 2-point finite difference approximation. .. tip:: Use the ``estimate_`` attribute to produce a polygon from the estimated parameter vector (``p_``). Parameters: * x, z : array Arrays with the x and z coordinates of the profile data points * gz : array The profile gravity anomaly data * verts : list of lists ``[[x1, z1], [x2, z2]]`` List of the [x, z] coordinates of the left and right know vertices, respectively. .. warning:: Very important that the vertices in the list be ordered from left to right! Otherwise the forward model will give results with an inverted sign and terrible things may happen! * density : float Density contrast of the basin * delta : float Interval used to calculate the approximate derivatives .. note:: The recommended solver for this inverse problem is the Levemberg-Marquardt method. Since this is a non-linear problem, set the desired method and initial solution using the ``config`` method of this class. See the example bellow. Example using synthetic data: >>> import numpy as np >>> from fatiando.mesher import Polygon >>> from fatiando.gravmag import talwani >>> # Make a triangular basin model (will estimate the last point) >>> verts = [(10000, 1), (90000, 1), (50000, 5000)] >>> left, middle, right = verts >>> model = Polygon(verts, {'density':500}) >>> # Generate the synthetic gz profile >>> x = np.linspace(0, 100000, 50) >>> z = np.zeros_like(x) >>> gz = talwani.gz(x, z, [model]) >>> # Make a solver and fit it to the data >>> solver = Triangular(x, z, gz, [left, middle], 500).config( ... 'levmarq', initial=[10000, 1000]).fit() >>> # p_ is the estimated parameter vector (x and z in this case) >>> x, z = solver.p_ >>> print('{:.1f}, {:.1f}'.format(x, z)) 50000.0, 5000.0 >>> # The parameter vector is not that useful so use estimate_ to get a >>> # Polygon object >>> poly = solver.estimate_ >>> poly.vertices array([[ 1.00000000e+04, 1.00000000e+00], [ 9.00000000e+04, 1.00000000e+00], [ 5.00000000e+04, 5.00000000e+03]]) >>> poly.props {'density': 500} >>> # Check is the residuals are all small >>> np.all(np.abs(solver.residuals()) < 10**-10) True """ def __init__(self, x, z, gz, verts, density): assert x.shape == z.shape == gz.shape, \ "x, z, and data must be of same length" assert len(verts) == 2, \ "Need exactly 2 vertices. {} given".format(len(verts)) super().__init__(data=gz, nparams=2, islinear=False) self.x = np.array(x, dtype=np.float) self.z = np.array(z, dtype=np.float) self.density = density self.verts = list(verts)
[docs] def predicted(self, p): """ Calculate predicted data for a given parameter vector. """ polygon = Polygon(self.verts + [p], {'density': self.density}) return talwani.gz(self.x, self.z, [polygon])
[docs] def jacobian(self, p): """ Calculate the Jacobian (sensitivity) matrix for a given parameter vector. """ delta = 1. props = {'density': self.density} xp, zp = self.x, self.z verts = self.verts x, z = p jac = np.transpose([ (talwani.gz(xp, zp, [Polygon(verts + [[x + delta, z]], props)]) - talwani.gz(xp, zp, [Polygon(verts + [[x - delta, z]], props)]) ) / (2. * delta), (talwani.gz(xp, zp, [Polygon(verts + [[x, z + delta]], props)]) - talwani.gz(xp, zp, [Polygon(verts + [[x, z - delta]], props)]) ) / (2. * delta)]) return jac
[docs] def fmt_estimate(self, p): """ Convert the parameter vector to a :class:`~fatiando.mesher.Polygon` so that it can be used for plotting and forward modeling. """ left, right = self.verts props = {'density': self.density} poly = Polygon(np.array([left, right, p]), props=props) return poly
[docs]class Trapezoidal(Misfit): """ Estimate the relief of a trapezoidal basin. Use when the basin can be approximated by a 2D body with **trapezoidal** vertical cross-section, like in rifts. The trapezoid is assumed to have 2 known vertices at the surface (the edges of the basin) and two unknown vertices in the subsurface. We assume that the x coordinates of the unknown vertices are the same as the x coordinates of the known vertices (i.e., the unknown vertices are directly under the known vertices). The inversion will then estimate the z coordinates of the unknown vertices. The forward modeling is done using :mod:`~fatiando.gravmag.talwani`. Derivatives are calculated using a 2-point finite difference approximation. .. tip:: Use the ``estimate_`` attribute to produce a polygon from the estimated parameter vector (``p_``). Parameters: * x, z : array Arrays with the x and z coordinates of the profile data points * gz : array The profile gravity anomaly data * verts : list of lists ``[[x1, z1], [x2, z2]]`` List of the [x, z] coordinates of the left and right know vertices, respectively. .. warning:: Very important that the vertices in the list be ordered from left to right! Otherwise the forward model will give results with an inverted sign and terrible things may happen! * density : float Density contrast of the basin * delta : float Interval used to calculate the approximate derivatives .. note:: The recommended solver for this inverse problem is the Levemberg-Marquardt method. Since this is a non-linear problem, set the desired method and initial solution using the ``config`` method. See the example bellow. Example with synthetic data: >>> import numpy as np >>> from fatiando.mesher import Polygon >>> from fatiando.gravmag import talwani >>> # Make a trapezoidal basin model (will estimate the z coordinates >>> # of the last two points) >>> verts = [[10000, 1], [90000, 1], [90000, 5000], [10000, 3000]] >>> model = Polygon(verts, {'density':500}) >>> # Generate the synthetic gz profile >>> x = np.linspace(0, 100000, 50) >>> z = np.zeros_like(x) >>> gz = talwani.gz(x, z, [model]) >>> # Make a solver and fit it to the data >>> solver = Trapezoidal(x, z, gz, verts[0:2], 500).config( ... 'levmarq', initial=[1000, 500]).fit() >>> # p_ is the estimated parameter vector (z1 and z2 in this case) >>> z1, z2 = solver.p_ >>> print('{:.1f}, {:.1f}'.format(z1, z2)) 5000.0, 3000.0 >>> # The parameter vector is not that useful so use estimate_ to get a >>> # Polygon object >>> poly = solver.estimate_ >>> poly.vertices array([[ 1.00000000e+04, 1.00000000e+00], [ 9.00000000e+04, 1.00000000e+00], [ 9.00000000e+04, 5.00000000e+03], [ 1.00000000e+04, 3.00000000e+03]]) >>> poly.props {'density': 500} >>> # Check is the residuals are all small >>> np.all(np.abs(solver.residuals()) < 10**-10) True """ def __init__(self, x, z, gz, verts, density): assert x.shape == z.shape == gz.shape, \ "x, z, and data must be of same length" assert len(verts) == 2, \ "Need exactly 2 vertices. {} given".format(len(verts)) super().__init__(data=gz, nparams=2, islinear=False) self.x = np.array(x, dtype=np.float) self.z = np.array(z, dtype=np.float) self.density = density self.props = {'density': self.density} self.verts = list(verts) self.x1, self.x2 = verts[1][0], verts[0][0] def predicted(self, p): z1, z2 = p model = [Polygon(self.verts + [[self.x1, z1], [self.x2, z2]], self.props)] pred = talwani.gz(self.x, self.z, model) return pred def jacobian(self, p): z1, z2 = p x1, x2 = self.x1, self.x2 x, z = self.x, self.z props = self.props verts = self.verts delta = 1. jac = np.empty((self.ndata, self.nparams), dtype=np.float) z1p = [Polygon(verts + [[x1, z1 + delta], [x2, z2]], props)] z1m = [Polygon(verts + [[x1, z1 - delta], [x2, z2]], props)] jac[:, 0] = (talwani.gz(x, z, z1p) - talwani.gz(x, z, z1m))/(2*delta) z2p = [Polygon(verts + [[x1, z1], [x2, z2 + delta]], props)] z2m = [Polygon(verts + [[x1, z1], [x2, z2 - delta]], props)] jac[:, 1] = (talwani.gz(x, z, z2p) - talwani.gz(x, z, z2m))/(2*delta) return jac
[docs] def fmt_estimate(self, p): """ Convert the parameter vector to a :class:`fatiando.mesher.Polygon` so that it can be used for plotting and forward modeling. """ z1, z2 = p left, right = self.verts poly = Polygon(np.array([left, right, [self.x1, z1], [self.x2, z2]]), self.props) return poly