Repository: fatiando/harmonica (python)

Issue:
Extend new default value for `depth` to gradient boosted equivalent sources
**Description of the desired feature:**

In #491 we implemented a new default value for the `depth` argument in the `hm.EquivalentSources` class. Since the  `hm.EquivalentSourcesGB` class inherits from it, we should extend the support for the default `depth` value to this class as wel..



**Are you willing to help implement and maintain this feature?** 

Yes, sure!

Pull request description:
Change default value for depth in EquivalentSourcesGB
Change the default value for the `depth` argument in `EquivalentSourcesGB` to `"default"`. Set the depth argument in the same way as `EquivalentSources` does, relying on the inherited methods. Extend docstrings to reflect these changes.

**Relevant issues/PRs:**

Fixes #513


Patch:
diff --git a/harmonica/_equivalent_sources/gradient_boosted.py b/harmonica/_equivalent_sources/gradient_boosted.py
index 85be86c..5f90d0e 100644
--- a/harmonica/_equivalent_sources/gradient_boosted.py
+++ b/harmonica/_equivalent_sources/gradient_boosted.py
@@ -7,6 +7,8 @@
 """
 Gradient-boosted equivalent sources in Cartesian coordinates
 """
+from __future__ import annotations
+
 import warnings
 
 import numpy as np
@@ -47,13 +49,16 @@ class EquivalentSourcesGB(EquivalentSources):
         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.
@@ -87,6 +92,10 @@ class EquivalentSourcesGB(EquivalentSources):
         The boundaries (``[W, E, S, N]``) of the data used to fit the
         interpolator. Used as the default region for the
         :meth:`~harmonica.EquivalentSources.grid` method.
+    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.
     window_size_ : float or None
         Size of the overlapping windows used in gradient-boosting equivalent
         point sources. It will be set to None if ``window_size = "default"``
@@ -105,7 +114,7 @@ class EquivalentSourcesGB(EquivalentSources):
         self,
         damping=None,
         points=None,
-        depth=500,
+        depth: float | str = "default",
         block_size=None,
         window_size="default",
         parallel=True,
@@ -221,6 +230,7 @@ class EquivalentSourcesGB(EquivalentSources):
                 p.astype(self.dtype) for p in 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)
             )


Test patch:
diff --git a/harmonica/tests/test_gradient_boosted_eqs.py b/harmonica/tests/test_gradient_boosted_eqs.py
index 9492b1f..da02cc7 100644
--- a/harmonica/tests/test_gradient_boosted_eqs.py
+++ b/harmonica/tests/test_gradient_boosted_eqs.py
@@ -115,7 +115,7 @@ def test_custom_points(region, coordinates_small, data_small):
         i.ravel()
         for i in vd.grid_coordinates(region=region, shape=(3, 3), extra_coords=-550)
     )
-    eqs = EquivalentSourcesGB(points=points_custom, window_size=500)
+    eqs = EquivalentSourcesGB(points=points_custom, window_size=500, depth=500)
     eqs.fit(coordinates_small, data_small)
     # Check that the proper source locations were set
     npt.assert_allclose(points_custom, eqs.points_, rtol=1e-5)
@@ -171,7 +171,7 @@ def test_gradient_boosted_eqs_single_window(region, points, masses, coordinates,
     """
     Test GB eq-sources with a single window that covers the whole region
     """
-    eqs = EquivalentSourcesGB(window_size=region[1] - region[0])
+    eqs = EquivalentSourcesGB(depth=500, window_size=region[1] - region[0])
     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
@@ -410,3 +410,27 @@ def test_window_size():
 def test_invalid_window_size():
     with pytest.raises(ValueError, match="Found invalid 'window_size' value equal to"):
         EquivalentSourcesGB(window_size="Chuckie took my soul!")
+
+
+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 = EquivalentSourcesGB().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):
+        EquivalentSourcesGB(depth=invalid_depth)