diff --git a/doc/user_guide/equivalent_sources/index.rst b/doc/user_guide/equivalent_sources/index.rst index 2ead59b0c..a9a6e79fb 100644 --- a/doc/user_guide/equivalent_sources/index.rst +++ b/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 mean 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. diff --git a/harmonica/_equivalent_sources/cartesian.py b/harmonica/_equivalent_sources/cartesian.py index f86fcb41a..2944f719e 100644 --- a/harmonica/_equivalent_sources/cartesian.py +++ b/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 @@ class EquivalentSources(vdb.BaseGridder): 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 mean 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 @@ class EquivalentSources(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. - 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. @@ -129,6 +139,10 @@ class EquivalentSources(vdb.BaseGridder): 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 mean + 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 @@ def __init__( 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 @@ def fit(self, coordinates, data, weights=None): 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 @@ def _build_points(self, coordinates): 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 mean 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 @@ def _build_points(self, coordinates): """ if self.block_size is not None: coordinates = self._block_average_coordinates(coordinates) + if self.depth == "default": + self.depth_ = 4.5 * np.mean(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): diff --git a/harmonica/tests/test_eq_sources_cartesian.py b/harmonica/tests/test_eq_sources_cartesian.py index 14a98a9c6..33f008d46 100644 --- a/harmonica/tests/test_eq_sources_cartesian.py +++ b/harmonica/tests/test_eq_sources_cartesian.py @@ -104,50 +104,27 @@ 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.7e-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 @@ -155,11 +132,11 @@ def test_equivalent_sources_cartesian_float32( 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,27 @@ 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) + + +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): + EquivalentSources(depth=invalid_depth)