Add option to set EQLHarmonic points to constant depth (#236)

Rename the relative_depth parameter for depth and add a new depth_type
parameter. Keep supporting the relative_depth parameter but raise
a FutureWarning if passed. Add a new _build_points() method for creating
the point source coordinates based on the data points. Add test functions for
the new feature using pytest.fixture for sample coordinates. Make examples to
use the new depth argument.

examples/eql/harmonic.py

@@ -20,7 +20,7 @@
 observed data. The fitted layer can then be used to predict data values
 wherever we want, like on a grid at a certain altitude. The sources for
 :class:`~harmonica.EQLHarmonic` in particular are placed one beneath each data
-point at a constant relative depth from the elevation of the data point
+point at a relative depth from the elevation of the data point
 following [Cooper2000]_.
 
 The advantage of using an equivalent layer is that it takes into account the 3D
@@ -54,9 +54,9 @@
 coordinates = (easting, northing, data.altitude_m)
 
 # Create the equivalent layer. We'll use the default point source configuration
-# at a constant relative depth beneath each observation point. The damping
-# parameter helps smooth the predicted data and ensure stability.
-eql = hm.EQLHarmonic(relative_depth=1000, damping=1)
+# at a relative depth beneath each observation point.
+# The damping parameter helps smooth the predicted data and ensure stability.
+eql = hm.EQLHarmonic(depth=1000, damping=1)
 
 # Fit the layer coefficients to the observed magnetic anomaly.
 eql.fit(coordinates, data.total_field_anomaly_nt)

harmonica/equivalent_layer/harmonic.py

@@ -7,6 +7,7 @@
 """
 Equivalent layer for generic harmonic functions in Cartesian coordinates
 """
+import warnings
 import numpy as np
 from numba import jit
 from sklearn.utils.validation import check_is_fitted
@@ -74,13 +75,27 @@ class EQLHarmonic(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.
- relative_depth : float
- Relative depth at which the point sources are placed beneath the
+ depth : float
+ Parameter used to control the depth at which the point sources will be
+ located.
+ If ``depth_type`` is equal to ``"relative"``, the ``depth`` specifies
+ the relative depth at which the point sources are placed beneath the
  observation points. Each source point will be set beneath each data
  point at a depth calculated as the elevation of the data point minus
- this constant *relative_depth*. Use positive numbers (negative numbers
- would mean point sources are above the data points). Ignored if
- *points* is specified.
+ this *depth*. Use positive numbers (negative numbers would mean point
+ sources are above the data points).
+ If ``depth_type`` is equal to ``"constant"``, the ``depth`` specifies
+ the constant depth at which the point sources are placed beneath the
+ observation points. Every source point will be located at this *depth*.
+ Use positive numbers (negative numbers would mean point sources are
+ located above the zeroth level).
+ This parameter is ignored if *points* is specified.
+ Defaults to 500.
+ depth_type : str
+ Strategy used for setting the depth of the point sources.
+ The two available strategies are ``"constant"`` and ``"relative"``.
+ This parameter is ignored if *points* is specified.
+ Defaults to ``"relative"``.
  parallel : bool
  If True any predictions and Jacobian building is carried out in
  parallel through Numba's ``jit.prange``, reducing the computation time.
@@ -113,15 +128,32 @@ def __init__(
  self,
  damping=None,
  points=None,
- relative_depth=500,
+ depth=500,
+ depth_type="relative",
  parallel=True,
+ **kwargs,
  ):
  self.damping = damping
  self.points = points
- self.relative_depth = relative_depth
+ self.depth = depth
+ self.depth_type = depth_type
  self.parallel = parallel
  # Define Green's function for Cartesian coordinates
  self.greens_function = greens_func_cartesian
+ # Check if depth_type is valid
+ if depth_type not in ("constant", "relative"):
+ raise ValueError(
+ f"Invalid depth type '{depth_type}'. Should be either be 'constant' or 'relative'."
+ )
+ # Check if relative_depth has been passed (will be deprecated)
+ if "relative_depth" in kwargs:
+ warnings.warn(
+ "The 'relative_depth' parameter is deprecated, please use "
+ + "the 'depth' paramter and set 'depth_type' to 'relative_depth' instead. ",
+ FutureWarning,
+ )
+ # Override depth and depth_type
+ self.depth, self.depth_type = kwargs["relative_depth"], "relative"
 
  def fit(self, coordinates, data, weights=None):
  """
@@ -155,17 +187,54 @@ def fit(self, coordinates, data, weights=None):
  self.region_ = vd.get_region(coordinates[:2])
  coordinates = vdb.n_1d_arrays(coordinates, 3)
  if self.points is None:
- self.points_ = (
- coordinates[0],
- coordinates[1],
- coordinates[2] - self.relative_depth,
- )
+ self.points_ = self._build_points(coordinates)
  else:
  self.points_ = vdb.n_1d_arrays(self.points, 3)
  jacobian = self.jacobian(coordinates, self.points_)
  self.coefs_ = vdb.least_squares(jacobian, data, weights, self.damping)
  return self
 
+ def _build_points(self, coordinates):
+ """
+ Generate coordinates of point sources based on the data points
+
+ Locate the point sources following the chosen ``depth_type`` strategy.
+ If ``depth_type`` is equal to ``"relative"``, the point sources will be
+ placed beneath the observation points at a depth calculated as the
+ elevation of the data point minus the ``depth``.
+ If ``depth_type`` is equal to ``"constant"``, the point sources will be
+ placed beneath the observation points at the same height equal to minus
+ ``depth``.
+
+ Parameters
+ ----------
+ coordinates : tuple of arrays
+ Arrays with the coordinates of each data point. Should be in the
+ following order: (``easting``, ``northing``, ``upward``, ...).
+ Only ``easting``, ``northing``, and ``upward`` will be used, all
+ subsequent coordinates will be ignored.
+
+ Returns
+ -------
+ points : tuple of arrays
+ Tuple containing the coordinates of the point sources used as the
+ equivalent layer, in the following order:
+ (``easting``, ``northing``, ``upward``).
+ """
+ if self.depth_type == "relative":
+ return (
+ coordinates[0],
+ coordinates[1],
+ coordinates[2] - self.depth,
+ )
+ if self.depth_type == "constant":
+ return (
+ coordinates[0],
+ coordinates[1],
+ -self.depth * np.ones_like(coordinates[0]),
+ )
+ return None
+
  def predict(self, coordinates):
  """
  Evaluate the estimated equivalent layer on the given set of points.
@@ -243,7 +312,7 @@ def grid(
  dims=None,
  data_names=None,
  projection=None,
- **kwargs
+ **kwargs,
  ): # pylint: disable=arguments-differ
  """
  Interpolate the data onto a regular grid.
@@ -327,7 +396,7 @@ def scatter(
  dims=None,
  data_names=None,
  projection=None,
- **kwargs
+ **kwargs,
  ):
  """
  .. warning ::
@@ -347,7 +416,7 @@ def profile(
  dims=None,
  data_names=None,
  projection=None,
- **kwargs
+ **kwargs,
  ): # pylint: disable=arguments-differ
  """
  Interpolate data along a profile between two points.

harmonica/tests/test_eql_harmonic.py

@@ -4,6 +4,7 @@
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
+# pylint: disable=protected-access
 """
 Test the EQLHarmonic gridder
 """
@@ -98,7 +99,7 @@ def test_eql_harmonic_small_data_cartesian():
  data = point_mass_gravity(coordinates, points, masses, field="g_z")
 
  # The interpolation should be perfect on the data points
- eql = EQLHarmonic(relative_depth=500)
+ eql = EQLHarmonic(depth=500)
  eql.fit(coordinates, data)
  npt.assert_allclose(data, eql.predict(coordinates), rtol=1e-5)
 
@@ -129,6 +130,111 @@ def test_eql_harmonic_small_data_cartesian():
  npt.assert_allclose(true, profile.scalars, rtol=0.05)
 
 
+@pytest.fixture(name="coordinates")
+def fixture_coordinates():
+ """
+ Return a set of sample coordinates intended to be used in tests
+ """
+ region = (-3e3, -1e3, 5e3, 7e3)
+ # Define a set of observation points with variable elevation coordinates
+ easting, northing = vd.grid_coordinates(region=region, shape=(8, 8))
+ upward = np.arange(64, dtype=float).reshape((8, 8))
+ coordinates = (easting, northing, upward)
+ return coordinates
+
+
+@pytest.mark.parametrize(
+ "depth_type, upward_expected",
+ [
+ ("relative", np.arange(64, dtype=float).reshape((8, 8)) - 1.5e3),
+ ("constant", -1.5e3 * np.ones((8, 8))),
+ ],
+ ids=["relative", "constant"],
+)
+def test_eql_harmonic_build_points(
+ coordinates,
+ depth_type,
+ upward_expected,
+):
+ """
+ Check if build_points method works as expected
+ """
+ eql = EQLHarmonic(depth=1.5e3, depth_type=depth_type)
+ points = eql._build_points(coordinates)
+ expected = (*coordinates[:2], upward_expected)
+ npt.assert_allclose(points, expected)
+
+
+def test_eql_harmonic_build_points_bacwkards(coordinates):
+ """
+ Check if the old relative_depth argument is well supported
+
+ This test is intended to check if backward compatibility is working
+ correctly. The ``relative_depth`` parameter will be deprecated on the next
+ major release.
+ """
+ depth = 4.5e3
+ expected_upward = coordinates[2] - depth
+ # Check if FutureWarning is raised after passing relative_depth
+ with warnings.catch_warnings(record=True) as warn:
+ eql = EQLHarmonic(relative_depth=depth)
+ assert len(warn) == 1
+ assert issubclass(warn[-1].category, FutureWarning)
+ # Check if the `depth` and `depth_type` attributes are well fixed
+ npt.assert_allclose(eql.depth, depth)
+ assert eql.depth_type == "relative"
+ # Check if location of sources are correct
+ points = eql._build_points(coordinates)
+ expected = (*coordinates[:2], expected_upward)
+ npt.assert_allclose(points, expected)
+
+
+def test_eql_harmonic_invalid_depth_type():
+ """
+ Check if ValueError is raised if invalid depth_type is passed
+ """
+ with pytest.raises(ValueError):
+ EQLHarmonic(depth=300, depth_type="blabla")
+
+
+def test_eql_harmonic_points_depth():
+ """
+ Check if the points coordinates are properly defined by the fit method
+ """
+ region = (-3e3, -1e3, 5e3, 7e3)
+ # Build synthetic point masses
+ points = vd.grid_coordinates(region=region, shape=(6, 6), extra_coords=-1e3)
+ masses = vd.datasets.CheckerBoard(amplitude=1e13, region=region).predict(points)
+ # Define a set of observation points with variable elevation coordinates
+ easting, northing = vd.grid_coordinates(region=region, shape=(5, 5))
+ upward = np.arange(25, dtype=float).reshape((5, 5))
+ coordinates = (easting, northing, upward)
+ # Get synthetic data
+ data = point_mass_gravity(coordinates, points, masses, field="g_z")
+
+ # Test with constant depth
+ eql = EQLHarmonic(depth=1.3e3, depth_type="constant")
+ eql.fit(coordinates, data)
+ expected_points = vdb.n_1d_arrays(
+ (easting, northing, -1.3e3 * np.ones_like(easting)), n=3
+ )
+ npt.assert_allclose(expected_points, eql.points_)
+
+ # Test with relative depth
+ eql = EQLHarmonic(depth=1.3e3, depth_type="relative")
+ eql.fit(coordinates, data)
+ expected_points = vdb.n_1d_arrays((easting, northing, upward - 1.3e3), n=3)
+ npt.assert_allclose(expected_points, eql.points_)
+
+ # Test with invalid depth_type
+ eql = EQLHarmonic(depth=300, depth_type="constant") # init with valid depth_type
+ eql.depth_type = "blabla" # change depth_type afterwards
+ points = eql._build_points(
+ vd.grid_coordinates(region=(-1, 1, -1, 1), spacing=0.25, extra_coords=1)
+ )
+ assert points is None
+
+
 def test_eql_harmonic_custom_points_cartesian():
  """
  Check that passing in custom points works and actually uses the points