Add Euler Deconvolution of a single window

doc/api/index.rst

@@ -101,6 +101,15 @@ Isostatic Moho
 
     isostatic_moho_airy
 
+Source estimation
+-----------------
+
+.. autosummary::
+    :toctree: generated/
+
+    EulerDeconvolution
+
+
 Input and Output
 ----------------
 

doc/references.rst

@@ -14,6 +14,7 @@ References
 .. [Nagy2000] Nagy, D., Papp, G. & Benedek, J.(2000). The gravitational potential and its derivatives for the prism. Journal of Geodesy 74: 552. doi:`10.1007/s001900000116 <https://doi.org/10.1007/s001900000116>`__
 .. [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>`__
 .. [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>`__
 .. [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>`_

harmonica/__init__.py

@@ -9,6 +9,7 @@
 from ._equivalent_sources.cartesian import EquivalentSources
 from ._equivalent_sources.gradient_boosted import EquivalentSourcesGB
 from ._equivalent_sources.spherical import EquivalentSourcesSph
+from ._euler_deconvolution import EulerDeconvolution
 from ._forward.dipole import dipole_magnetic
 from ._forward.point import point_gravity
 from ._forward.prism_gravity import prism_gravity

harmonica/_euler_deconvolution.py

@@ -0,0 +1,91 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+import numpy as np
+import scipy as sp
+
+
+class EulerDeconvolution:
+    r"""
+    Estimate source location and base level using Euler Deconvolution
+
+    Implements Euler Deconvolution [Reid1990]_ to estimate subsurface source
+    location and a base level constant from potential field data and their
+    directional derivatives. The approach employs linear least-squares to solve
+    Euler's homogeneity equation. **Assumes a single data window** and provides
+    a single estimate.
+
+    Parameters
+    ----------
+    structural_index : int
+        Defines the nature of the source of the potential field data. It's the
+        degree of the field's rate of change with distance from the source,
+        influencing the decay rate of the field and the formulation of Euler's
+        homogeneity equation. **Correlated with the depth estimate**, so larger
+        structural index will lead to larger depths. **Choose based on known
+        source geometry**. See table below.
+
+    Attributes
+    ----------
+    location_ : numpy.ndarray
+        Estimated (easting, northing, upward) coordinates of the source after
+        model fitting.
+    base_level_ : float
+        Estimated base level constant of the anomaly after model fitting.
+
+    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
+
+    def fit(self, coordinates, field, east_deriv, north_deriv, up_deriv):
+        """
+        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.
+
+        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.
+
+        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
+            + self.structural_index * field
+        )
+        estimate = sp.linalg.solve(matrix.T @ matrix, matrix.T @ data_vector)
+
+        self.location_ = estimate[:3]
+        self.base_level_ = estimate[3]

harmonica/tests/test_euler_deconvolution.py

@@ -0,0 +1,95 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+import numpy.testing as npt
+import verde as vd
+
+from .. import (
+    EulerDeconvolution,
+    derivative_easting,
+    derivative_northing,
+    derivative_upward,
+    dipole_magnetic,
+    magnetic_angles_to_vec,
+)
+from .._forward.point import point_gravity
+
+
+def test_euler_with_numeric_derivatives():
+    # Add dipole source
+    dipole_coordinates = (10e3, 15e3, -10e3)
+    dipole_moments = magnetic_angles_to_vec(1.0e14, 0, 0)
+
+    # Add regional field
+    inc, dec = -40, 15
+    fe, fn, fu = magnetic_angles_to_vec(1, inc, dec)
+    region = [-100e3, 100e3, -80e3, 80e3]
+    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    be, bn, bu = dipole_magnetic(
+        coordinates, dipole_coordinates, dipole_moments, field="b"
+    )
+
+    # Add a fixed base level
+    true_base_level = 200  # nT
+    anomaly = (fe * be + fn * bn + fu * bu) + true_base_level
+
+    grid = vd.make_xarray_grid(
+        coordinates, anomaly, data_names="tfa", extra_coords_names="upward"
+    )
+    grid["d_east"] = derivative_easting(grid.tfa)
+    grid["d_north"] = derivative_northing(grid.tfa)
+    grid["d_up"] = derivative_upward(grid.tfa)
+    grid_table = vd.grid_to_table(grid)
+
+    euler = EulerDeconvolution(structural_index=3)
+
+    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,
+    )
+
+    npt.assert_allclose(euler.location_, dipole_coordinates, atol=1.0e-3, rtol=1.0e-3)
+    npt.assert_allclose(euler.base_level_, true_base_level, atol=1.0e-3, rtol=1.0e-3)
+
+
+def test_euler_with_analytic_derivatives():
+    # Add dipole source
+    masses_coordinates = (10e3, 15e3, -10e3)
+    masses = 1.0e12
+    region = [-100e3, 100e3, -80e3, 80e3]
+    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    gz = point_gravity(coordinates, masses_coordinates, masses, field="g_z")
+
+    # Convert Eötvös to mGal because derivatives must be in mGal/m
+    eotvos2mgal = 1.0e-4
+    gzz = (
+        -point_gravity(coordinates, masses_coordinates, masses, field="g_zz")
+        * eotvos2mgal
+    )
+    gze = (
+        point_gravity(coordinates, masses_coordinates, masses, field="g_ez")
+        * eotvos2mgal
+    )
+    gzn = (
+        point_gravity(coordinates, masses_coordinates, masses, field="g_nz")
+        * eotvos2mgal
+    )
+
+    euler = EulerDeconvolution(structural_index=2)
+    euler.fit(
+        (coordinates[0].ravel(), coordinates[1].ravel(), coordinates[2].ravel()),
+        gz.ravel(),
+        gze.ravel(),
+        gzn.ravel(),
+        gzz.ravel(),
+    )
+
+    npt.assert_allclose(euler.location_, masses_coordinates, atol=1.0e-3, rtol=1.0e-3)
+    npt.assert_allclose(euler.base_level_, 0.0, atol=1.0e-3, rtol=1.0e-3)