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)