Change default value of depth in equivalent sources

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 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.

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 @@
 
     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 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 @@
         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 median 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 @@
         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 median
+        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 @@
         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 @@
         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 @@
         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 median 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 @@
         """
         if self.block_size is not None:
             coordinates = self._block_average_coordinates(coordinates)
+        if self.depth == "default":
+            self.depth_ = 4.5 * np.median(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):

harmonica/tests/test_eq_sources_cartesian.py

@@ -104,62 +104,39 @@ 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.5e-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
     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), 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,17 @@ 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)