Add option to set EQLHarmonic points to constant depth

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 bellow each observation point at
  a fixed relative depth bellow 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 *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