Add covariance, change fit data, more docs to Euler Deconvolution (#519)

The covariance matrix will be needed to filter solutions in the windowed
approach. Change the fit input to receive all the data in a single
tuple, like the coordinates. This is more compatible with the rest of
Fatiando. Added the checks for fit input and standardizing array shapes
so inputs can be n-dimensional. Copied over some of the docstrings from
the legacy Fatiando implementation, which expands the theory and
provides links to good papers to read.

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)