Repository: fatiando/harmonica (python)

Issue:
Better default for depth in equivalent sources
**Description of the desired feature:**
<!--
Add as much detail as you can in your description below. If possible, include an
example of how you would like to use this feature (even better if it's a code
example).
-->

From chats with @indiauppal and @birocoles, I learned that Dampney (1969; the original equivalent source paper) advises using depths between 2.5 and 6 times the source spacing (assuming sources in a regular grid). For sources beneath data, like  we use here, a sensible default is 4.5 times the median distance between sources. Me and @indiauppal tested this default and it works well most of the time. Tweaking is always good but this would provide a more sensible default than the fixed value we currently use.

This would break backwards compatibility. I'm not too worried since our docs say that you should always set the `depth` by hand. So changing the default will likely make the interpolation better for the few people who were using the classes without following our advice. But I can do the warning dance if others think it's necessary (but it will take longer to implement).

**Are you willing to help implement and maintain this feature?** 
<!--
Every feature we add is code that we will have to maintain and keep updated.
Please let us know if you're willing to help maintain this feature in the future.
-->
Yes. But very happy to let someone else do it (😉 @indiauppal)

Pull request:
Change default value of depth in equivalent sources
Change the default value of the `depth` argument in `EquivalentSources` to `"default"`. If this value is passed, then the depth of the sources is set as 4.5 times the mean distance between first neighbouring sources. If a numerical value is passed, this is the one that will be used. Introduce a new `depth_` attribute where the estimated/passed numeric value is stored. If `points` is passed, the `depth_` attribute is set to None. Add a test checking that the depth value is properly set. Merge the tests functions that checked the accuracy of equivalent sources using different dtypes.


<!--
Thank you for contributing a pull request to Fatiando! 💖

👆🏽 ABOVE: Describe the changes proposed and WHY you made them.

👇🏽 BELOW: Link to any relevant issue or pull request.

Please ensure you have taken a look at the CONTRIBUTING.md file
in this repository (if available) and the general guidelines at
https://github.com/fatiando/community/blob/main/CONTRIBUTING.md
-->

**Relevant issues/PRs:**
<!--
Example: "Fixes #1234" / "See also #345" / "Relevant to #111"
Use keywords (e.g., Fixes, Closes) to create the links and automatically
close issues when this PR is merged.
See https://github.com/blog/1506-closing-issues-via-pull-requests
-->

Fixes #424


Patch:
diff --git a/doc/user_guide/equivalent_sources/index.rst b/doc/user_guide/equivalent_sources/index.rst
index 2ead59b..a9a6e79 100644
--- a/doc/user_guide/equivalent_sources/index.rst
+++ b/doc/user_guide/equivalent_sources/index.rst
@@ -72,7 +72,7 @@ Now we can initialize the :class:`harmonica.EquivalentSources` class.
 
    import harmonica as hm
 
-   equivalent_sources = hm.EquivalentSources(depth=10e3, damping=10)
+   equivalent_sources = hm.EquivalentSources(damping=10)
    equivalent_sources
 
 By default, it places the sources one beneath each data point at a relative
@@ -81,6 +81,16 @@ This *relative depth* can be set through the ``depth`` argument.
 Deepest sources generate smoother predictions (*underfitting*), while shallow
 ones tend to overfit the data.
 
+.. hint::
+
+   By default, since Harmonica v0.7.0, the sources will be located at a depth
+   below the data points estimated as 4.5 times the mean distance between
+   first neighboring sources. Alternatively, we can set a value for this depth
+   below the data points through the ``depth`` argument.
+
+   The estimated value for the depth of the sources can be explored through the
+   :attr:`harmonica.EquivalentSources.depth_` attribute.
+
 The ``damping`` parameter is used to smooth the coefficients of the sources and
 stabilize the least square problem. A higher ``damping`` will create smoother
 predictions, while a lower one could overfit the data and create artifacts.
diff --git a/harmonica/_equivalent_sources/cartesian.py b/harmonica/_equivalent_sources/cartesian.py
index f86fcb4..2944f71 100644
--- a/harmonica/_equivalent_sources/cartesian.py
+++ b/harmonica/_equivalent_sources/cartesian.py
@@ -7,6 +7,8 @@
 """
 Equivalent sources for generic harmonic functions in Cartesian coordinates
 """
+from __future__ import annotations
+
 import warnings
 
 import numpy as np
@@ -64,11 +66,16 @@ class EquivalentSources(vdb.BaseGridder):
 
     The depth of the sources can be controlled by the ``depth`` argument.
     Each source is located beneath each data point or block-averaged location
-    at a depth equal to its elevation minus the value of the ``depth``
-    argument.
-    In both cases a positive value of ``depth`` locates sources _beneath_ the
-    data points or the block-averaged locations, thus a negative ``depth`` will
-    put the sources _above_ them.
+    at a depth equal to its elevation minus the value of the ``depth_``
+    attribute.
+    If ``"default"`` is passed to the ``depth`` argument, then the ``depth_``
+    attribute is set to 4.5 times the mean distance between first neighboring
+    sources.
+    If a numerical value is passed to the ``depth`` argument, then this is the
+    one used for the ``depth_`` attribute.
+    A positive value of ``depth_`` locates sources _beneath_ the data points or
+    the block-averaged locations, thus a negative ``depth_`` will put the
+    sources _above_ them.
 
     Custom source locations can be chosen by specifying the ``points``
     argument, in which case the ``block_size`` and ``depth`` arguments will be
@@ -100,13 +107,16 @@ class EquivalentSources(vdb.BaseGridder):
         If None, will place one point source below each observation point at
         a fixed relative depth below the observation point [Cooper2000]_.
         Defaults to None.
-    depth : float
+    depth : float or "default"
         Parameter used to control the depth at which the point sources will be
         located.
-        Each source is located beneath each data point (or block-averaged
-        location) at a depth equal to its elevation minus the ``depth`` value.
+        If a value is provided, each source is located beneath each data point
+        (or block-averaged location) at a depth equal to its elevation minus
+        the ``depth`` value.
+        If set to ``"default"``, the depth of the sources will be estimated as
+        4.5 times the mean distance between first neighboring sources.
         This parameter is ignored if *points* is specified.
-        Defaults to 500.
+        Defaults to ``"default"``.
     block_size: float, tuple = (s_north, s_east) or None
         Size of the blocks used on block-averaged equivalent sources.
         If a single value is passed, the blocks will have a square shape.
@@ -129,6 +139,10 @@ class EquivalentSources(vdb.BaseGridder):
         Coordinates of the equivalent point sources.
     coefs_ : array
         Estimated coefficients of every point source.
+    depth_ : float or None
+        Estimated depth of the sources calculated as 4.5 times the mean
+        distance between first neighboring sources. This attribute is set to
+        None if ``points`` is passed.
     region_ : tuple
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
@@ -154,11 +168,16 @@ class EquivalentSources(vdb.BaseGridder):
         self,
         damping=None,
         points=None,
-        depth=500,
+        depth: float | str = "default",
         block_size=None,
         parallel=True,
         dtype="float64",
     ):
+        if isinstance(depth, str) and depth != "default":
+            raise ValueError(
+                f"Found invalid 'depth' value equal to '{depth}'. "
+                "It should be 'default' or a numeric value."
+            )
         self.damping = damping
         self.points = points
         self.depth = depth
@@ -205,6 +224,7 @@ class EquivalentSources(vdb.BaseGridder):
         if self.points is None:
             self.points_ = self._build_points(coordinates)
         else:
+            self.depth_ = None  # set depth_ to None so we don't leave it unset
             self.points_ = tuple(
                 p.astype(self.dtype) for p in vdb.n_1d_arrays(self.points, 3)
             )
@@ -220,7 +240,12 @@ class EquivalentSources(vdb.BaseGridder):
         and apply block-averaging if ``block_size`` is not None.
         The point sources will be placed beneath the (averaged) observation
         points at a depth calculated as the elevation of the data point minus
-        the ``depth``.
+        the ``depth_`` attribute.
+
+        If ``depth`` is set to ``"default"``, the ``depth_`` attribute is set
+        as 4.5 times the mean distance between first neighboring sources.
+        If ``depth`` is set to a numerical value, this is used for the
+        ``depth_`` attribute.
 
         Parameters
         ----------
@@ -238,10 +263,14 @@ class EquivalentSources(vdb.BaseGridder):
         """
         if self.block_size is not None:
             coordinates = self._block_average_coordinates(coordinates)
+        if self.depth == "default":
+            self.depth_ = 4.5 * np.mean(vd.median_distance(coordinates, k_nearest=1))
+        else:
+            self.depth_ = self.depth
         return (
             coordinates[0],
             coordinates[1],
-            coordinates[2] - self.depth,
+            coordinates[2] - self.depth_,
         )
 
     def _block_average_coordinates(self, coordinates):


Test patch:
diff --git a/harmonica/tests/test_eq_sources_cartesian.py b/harmonica/tests/test_eq_sources_cartesian.py
index 14a98a9..33f008d 100644
--- a/harmonica/tests/test_eq_sources_cartesian.py
+++ b/harmonica/tests/test_eq_sources_cartesian.py
@@ -104,50 +104,27 @@ def fixture_coordinates_9x9(region):
 
 
 @run_only_with_numba
-def test_equivalent_sources_cartesian(region, points, masses, coordinates, data):
+@pytest.mark.parametrize("dtype", ("default", "float32"))
+def test_equivalent_sources_cartesian(region, points, masses, coordinates, data, dtype):
     """
     Check that predictions are reasonable when interpolating from one grid to
     a denser grid. Use Cartesian coordinates.
     """
-    # The interpolation should be perfect on the data points
-    eqs = EquivalentSources()
-    eqs.fit(coordinates, data)
-    npt.assert_allclose(data, eqs.predict(coordinates), rtol=1e-5)
-
-    # Gridding onto a denser grid should be reasonably accurate when compared
-    # to synthetic values
-    upward = 0
-    shape = (60, 60)
-    grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
-    true = point_gravity(grid_coords, points, masses, field="g_z")
-    npt.assert_allclose(true, eqs.predict(grid_coords), rtol=1e-3)
-
-    # Test grid method
-    grid = eqs.grid(grid_coords)
-    npt.assert_allclose(true, grid.scalars, rtol=1e-3)
-
-    # Test profile method
-    point1 = (region[0], region[2])
-    point2 = (region[0], region[3])
-    profile = eqs.profile(point1, point2, upward, shape[0])
-    true = point_gravity(
-        (profile.easting, profile.northing, profile.upward), points, masses, field="g_z"
-    )
-    npt.assert_allclose(true, profile.scalars, rtol=1e-3)
+    # Set absolute tolerances for tests based on dtype (float32 should be less
+    # accurate)
+    if dtype == "float32":
+        kwargs = dict(dtype=dtype)
+        atol = 1.7e-3 * vd.maxabs(data)
+    else:
+        kwargs = {}
+        atol = 1e-3 * vd.maxabs(data)
 
+    # Fit the equivalent sources
+    eqs = EquivalentSources(**kwargs)
+    eqs.fit(coordinates, data)
 
-@run_only_with_numba
-def test_equivalent_sources_cartesian_float32(
-    region, points, masses, coordinates, data
-):
-    """
-    Check that predictions are reasonable when interpolating from one grid to
-    a denser grid, using float32 as dtype.
-    """
     # The interpolation should be perfect on the data points
-    eqs = EquivalentSources(dtype="float32")
-    eqs.fit(coordinates, data)
-    npt.assert_allclose(data, eqs.predict(coordinates), atol=1e-3 * vd.maxabs(data))
+    npt.assert_allclose(data, eqs.predict(coordinates), atol=atol)
 
     # Gridding onto a denser grid should be reasonably accurate when compared
     # to synthetic values
@@ -155,11 +132,11 @@ def test_equivalent_sources_cartesian_float32(
     shape = (60, 60)
     grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
     true = point_gravity(grid_coords, points, masses, field="g_z")
-    npt.assert_allclose(true, eqs.predict(grid_coords), atol=1e-3 * vd.maxabs(true))
+    npt.assert_allclose(true, eqs.predict(grid_coords), atol=atol)
 
     # Test grid method
     grid = eqs.grid(grid_coords)
-    npt.assert_allclose(true, grid.scalars, atol=1e-3 * vd.maxabs(true))
+    npt.assert_allclose(true, grid.scalars, atol=atol)
 
     # Test profile method
     point1 = (region[0], region[2])
@@ -168,7 +145,7 @@ def test_equivalent_sources_cartesian_float32(
     true = point_gravity(
         (profile.easting, profile.northing, profile.upward), points, masses, field="g_z"
     )
-    npt.assert_allclose(true, profile.scalars, atol=1e-3 * vd.maxabs(true))
+    npt.assert_allclose(true, profile.scalars, atol=atol)
 
 
 def test_equivalent_sources_small_data_cartesian(region, points, masses):
@@ -449,3 +426,27 @@ def test_error_ignored_args(coordinates_small, data_small, region):
     msg = "The 'bla' arguments are being ignored."
     with pytest.warns(FutureWarning, match=msg):
         eqs.grid(coordinates=grid_coords, bla="bla")
+
+
+def test_default_depth(coordinates, data):
+    """
+    Test if the depth of sources is correctly set by the default strategy
+    """
+    # Get distance to first neighbour in the grid
+    easting, northing = coordinates[:2]
+    d_easting = easting[1, 1] - easting[0, 0]
+    d_northing = northing[1, 1] - northing[0, 0]
+    first_neighbour_distance = min(d_easting, d_northing)
+    # Fit the equivalent sources with default `depth`
+    eqs = EquivalentSources().fit(coordinates, data)
+    npt.assert_allclose(eqs.depth_, first_neighbour_distance * 4.5)
+
+
+def test_invalid_depth():
+    """
+    Test if error is raised after passing invalid value for depth.
+    """
+    invalid_depth = "this is not a valid one"
+    msg = f"Found invalid 'depth' value equal to '{invalid_depth}'"
+    with pytest.raises(ValueError, match=msg):
+        EquivalentSources(depth=invalid_depth)