Repository: fatiando/harmonica
Pull request #519 (2024-07-16 12:31:47)

Add covariance, change fit data, more docs to Euler Deconvolution
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. 

Souza-junior: @leouieda I'll take a closer look into it.

Patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 8a6770f..c3a1c20 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -101,8 +101,8 @@ Isostatic Moho
 
     isostatic_moho_airy
 
-Source estimation
------------------
+Source position estimation
+--------------------------
 
 .. autosummary::
     :toctree: generated/
diff --git a/doc/references.rst b/doc/references.rst
index 781f2f8..e8f496d 100644
--- a/doc/references.rst
+++ b/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.
diff --git a/harmonica/_euler_deconvolution.py b/harmonica/_euler_deconvolution.py
index 92adda9..ad3d556 100644
--- a/harmonica/_euler_deconvolution.py
+++ b/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,15 +45,52 @@ 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):
@@ -46,26 +98,32 @@ class EulerDeconvolution:
         # 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
         -------
@@ -73,19 +131,31 @@ class EulerDeconvolution:
             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


Test patch:
diff --git a/harmonica/tests/test_euler_deconvolution.py b/harmonica/tests/test_euler_deconvolution.py
index 069d39f..1338a23 100644
--- a/harmonica/tests/test_euler_deconvolution.py
+++ b/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)