Change default value for depth in EquivalentSourcesGB

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 @@ def __init__(
  self,
  damping=None,
  points=None,
- depth=500,
+ depth: float | str = "default",
  block_size=None,
  window_size="default",
  parallel=True,
@@ -221,6 +230,7 @@ def fit(self, coordinates, data, weights=None):
  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)
  )

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)