Add covariance, change fit data, more docs to Euler Deconvolution

doc/api/index.rst

@@ -101,8 +101,8 @@ Isostatic Moho
 
  isostatic_moho_airy
 
-Source estimation
------------------
+Source position estimation
+--------------------------
 
 .. autosummary::
  :toctree: generated/

doc/references.rst

@@ -15,8 +15,11 @@ References
 .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2002). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy. doi:`10.1007/s00190-002-0264-7 <https://doi.org/10.1007/s00190-002-0264-7>`__
 .. [Oliveira2021] Oliveira Jr, Vanderlei C. and Uieda, Leonardo and Barbosa, Valeria C. F.. Sketch of three coordinate systems: Geocentric Cartesian, Geocentric Geodetic, and Topocentric Cartesian. figshare. doi: `10.6084/m9.figshare.15044241.v1 <https://doi.org/10.6084/m9.figshare.15044241.v1>`__
 .. [Reid1990] Reid, A. B., Allsop, J. M., Granser, H., Millett, A. J., & Somerton, I. W. (1990). Magnetic interpretation in three dimensions using Euler deconvolution. GEOPHYSICS. doi:`10.1190/1.1442774 <https://doi.org/10.1190/1.1442774>`__
+.. [Reid2014] Reid, A. B., J. Ebbing, and S. J. Webb (2014), Avoidable Euler Errors – the use and abuse of Euler deconvolution applied to potential fields, Geophysical Prospecting, doi:`10.1111/1365-2478.12119 <https://doi.org/10.1111/1365-2478.12119>`__.
+.. [ReidThurston2014] Reid, A., and J. Thurston (2014), The structural index in gravity and magnetic interpretation: Errors, uses, and abuses, GEOPHYSICS, 79(4), J61-J66, doi:`10.1190/geo2013-0235.1 <https://doi.org/10.1190/geo2013-0235.1>`__.
 .. [Soler2019] Soler, S. R., Pesce, A., Gimenez, M. E., & Uieda, L. (2019). Gravitational field calculation in spherical coordinates using variable densities in depth, Geophysical Journal International. doi: `10.1093/gji/ggz277 <https://doi.org/10.1093/gji/ggz277>`__
 .. [Soler2021] Soler, S. R. and Uieda, L. (2021). Gradient-boosted equivalent sources, Geophysical Journal International. doi:`10.1093/gji/ggab297 <https://doi.org/10.1093/gji/ggab297>`__
+.. [Uieda2014] Uieda, L., V. C. Oliveira Jr., and V. C. F. Barbosa (2014), Geophysical tutorial: Euler deconvolution of potential-field data, The Leading Edge, 33(4), 448-450, doi:`10.1190/tle33040448.1 <https://doi.org/10.1190/tle33040448.1>`__.
 .. [Uieda2015] Uieda, Leonardo (2015). A tesserioid (spherical prism) in a geocentric coordinate system with a local-North-oriented coordinate system. figshare. Figure. doi: `10.6084/m9.figshare.1495525.v1 <https://doi.org/10.6084/m9.figshare.1495525.v1>`_
 .. [Vajda2004] Vajda, P., Vaníček, P., Novák, P. and Meurers, B. (2004). On evaluation of Newton integrals in geodetic coordinates: Exact formulation and spherical approximation. Contributions to Geophysics and Geodesy, 34(4), 289-314.
 .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press.

harmonica/_euler_deconvolution.py

@@ -4,8 +4,12 @@
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
+"""
+Classes for Euler Deconvolution of potential field data
+"""
 import numpy as np
 import scipy as sp
+import verde.base as vdb
 
 
 class EulerDeconvolution:
@@ -18,6 +22,17 @@ class EulerDeconvolution:
  Euler's homogeneity equation. **Assumes a single data window** and provides
  a single estimate.
 
+ .. hint::
+
+ Please read the paper [Reid2014]_ to avoid doing **horrible things**
+ with Euler deconvolution. [Uieda2014]_ offer a practical tutorial using
+ `legacy Fatiando a Terra <https://legacy.fatiando.org/>`__ code to show
+ some common misinterpretations.
+
+ .. note::
+
+ Does not yet support structural index 0.
+
  Parameters
  ----------
  structural_index : int
@@ -30,62 +45,117 @@ class EulerDeconvolution:
 
  Attributes
  ----------
- location_ : numpy.ndarray
+ location_ : 1d-array
  Estimated (easting, northing, upward) coordinates of the source after
  model fitting.
  base_level_ : float
  Estimated base level constant of the anomaly after model fitting.
+ covariance_ : 2d-array
+ The 4 x 4 estimated covariance matrix of the solution. Parameters are
+ in the order: easting, northing, upward, base level. **This is not an
+ uncertainty of the position** but a rough estimate of their variance
+ with regard to the data.
+
+ Notes
+ -----
+
+ Works on any potential field that satisfies Euler's homogeneity equation
+ (like gravity, magnetic, and their gradients caused by **simple sources**):
+
+ .. math::
+
+ (e_i - e_0)\dfrac{\partial f_i}{\partial e} +
+ (n_i - n_0)\dfrac{\partial f_i}{\partial n} +
+ (u_i - u_0)\dfrac{\partial f_i}{\partial u} =
+ \eta (b - f_i),
+
+ in which :math:`f_i` is the given potential field observation at point
+ :math:`(e_i, n_i, u_i)`, :math:`b` is the base level (a constant shift of
+ the field, like a regional field), :math:`\eta` is the structural index,
+ and :math:`(e_0, n_0, u_0)` are the coordinates of a point on the source
+ (for a sphere, this is the center point).
+
+ The Euler deconvolution estimates :math:`(e_0, n_0, u_0)` and :math:`b`
+ given a potential field and its easting, northing, and upward derivatives
+ and the structural index. However, **this assumes that the sources are
+ ideal** (see the table below). We recommend reading [ReidThurston2014]_ for
+ a discussion on what the structural index means and what it does not mean.
+
+ After [ReidThurston2014]_, values of the structural index (SI) can be:
+
+ ===================================== ======== =========
+ Source type SI (Mag) SI (Grav)
+ ===================================== ======== =========
+ Point, sphere 3 2
+ Line, cylinder, thin bed fault 2 1
+ Thin sheet edge, thin sill, thin dyke 1 0
+ ===================================== ======== =========
 
- References
- ----------
- [Reid1990]_
  """
 
  def __init__(self, structural_index):
  self.structural_index = structural_index
  # The estimated parameters. Start them with None
  self.location_ = None
  self.base_level_ = None
+ self.covariance_ = None
 
- def fit(self, coordinates, field, east_deriv, north_deriv, up_deriv):
+ def fit(self, coordinates, data):
  """
  Fit the model using potential field measurements and their derivatives.
 
  Solves Euler's homogeneity equation to estimate the source location
  and base level by utilizing field values and their spatial derivatives
- in east, north, and upward directions.
+ in easting, northing, and upward directions.
+
+ .. tip::
+
+ Data does not need to be gridded for this to work.
 
  Parameters
  ----------
- coordinates : tuple of 1d-arrays
- Arrays with the coordinates of each data point, in the order of
- (x, y, z), representing easting, nothing, and upward directions,
- respectively.
- field : 1d-array
- Field measurements at each data point.
- east_deriv, north_deriv, up_deriv : 1d-array
- Partial derivatives of the field with respect to east, north, and
- upward directions, respectively.
+ coordinates : tuple of arrays
+ Tuple of 3 with the coordinates of each data point. Should be in
+ the following order: ``(easting, northing, upward)``.
+ Arrays can be n-dimensional but must all have the same shape.
+ data : tuple of arrays
+ Tuple of 4 arrays with the observed data in the following order:
+ ``(potential_field, derivative_easting, derivative_northing,
+ derivative_upward)``. Arrays can be n-dimensional but must all have
+ the same shape as the coordinates. Derivatives must be in data
+ units over coordinates units, for example nT/m or mGal/m.
 
  Returns
  -------
  self
  The instance itself, updated with the estimated `location_`
  and `base_level_`.
  """
- n_data = field.shape[0]
- matrix = np.empty((n_data, 4))
- matrix[:, 0] = east_deriv
- matrix[:, 1] = north_deriv
- matrix[:, 2] = up_deriv
- matrix[:, 3] = self.structural_index
- data_vector = (
- coordinates[0] * east_deriv
- + coordinates[1] * north_deriv
- + coordinates[2] * up_deriv
+ coordinates, data, _ = vdb.check_fit_input(coordinates, data, weights=None)
+ field, east_deriv, north_deriv, up_deriv = vdb.n_1d_arrays(data, 4)
+ easting, northing, upward = vdb.n_1d_arrays(coordinates, 3)
+ n_data = field.size
+ jacobian = np.empty((n_data, 4))
+ jacobian[:, 0] = east_deriv
+ jacobian[:, 1] = north_deriv
+ jacobian[:, 2] = up_deriv
+ jacobian[:, 3] = self.structural_index
+ pseudo_data = (
+ easting * east_deriv
+ + northing * north_deriv
+ + upward * up_deriv
  + self.structural_index * field
  )
- estimate = sp.linalg.solve(matrix.T @ matrix, matrix.T @ data_vector)
-
+ hessian = jacobian.T @ jacobian
+ # Invert the Hessian instead of solving the system because this is a
+ # 4x4 system and it won't cost much more, plus we need the inverse
+ # anyway to estimate the covariance matrix (used as a filtering
+ # criterion in windowed implementations)
+ hessian_inv = sp.linalg.inv(hessian)
+ estimate = hessian_inv @ jacobian.T @ pseudo_data
+ pseudo_residuals = pseudo_data - jacobian @ estimate
+ chi_squared = np.sum(pseudo_residuals**2) / (n_data - 4)
+ self.covariance_ = chi_squared * hessian_inv
  self.location_ = estimate[:3]
  self.base_level_ = estimate[3]
+ return self

harmonica/tests/test_euler_deconvolution.py

@@ -49,10 +49,7 @@ def test_euler_with_numeric_derivatives():
  coordinates = (grid_table.easting, grid_table.northing, grid_table.upward)
  euler.fit(
  (grid_table.easting, grid_table.northing, grid_table.upward),
- grid_table.tfa,
- grid_table.d_east,
- grid_table.d_north,
- grid_table.d_up,
+ (grid_table.tfa, grid_table.d_east, grid_table.d_north, grid_table.d_up),
  )
 
  npt.assert_allclose(euler.location_, dipole_coordinates, atol=1.0e-3, rtol=1.0e-3)
@@ -84,11 +81,8 @@ def test_euler_with_analytic_derivatives():
 
  euler = EulerDeconvolution(structural_index=2)
  euler.fit(
- (coordinates[0].ravel(), coordinates[1].ravel(), coordinates[2].ravel()),
- gz.ravel(),
- gze.ravel(),
- gzn.ravel(),
- gzz.ravel(),
+ (coordinates[0], coordinates[1], coordinates[2]),
+ (gz, gze, gzn, gzz),
  )
 
  npt.assert_allclose(euler.location_, masses_coordinates, atol=1.0e-3, rtol=1.0e-3)