Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change default value of depth in equivalent sources #491

Merged
merged 17 commits into from
Jun 12, 2024
Merged
Show file tree
Hide file tree
Changes from 15 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 11 additions & 1 deletion doc/user_guide/equivalent_sources/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand Down
53 changes: 41 additions & 12 deletions harmonica/_equivalent_sources/cartesian.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
"""
Equivalent sources for generic harmonic functions in Cartesian coordinates
"""
from __future__ import annotations

import warnings

import numpy as np
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand All @@ -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}' ."
santisoler marked this conversation as resolved.
Show resolved Hide resolved
"It should be 'default' or a numeric value."
)
self.damping = damping
self.points = points
self.depth = depth
Expand Down Expand Up @@ -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)
)
Expand All @@ -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
----------
Expand All @@ -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):
Expand Down
81 changes: 41 additions & 40 deletions harmonica/tests/test_eq_sources_cartesian.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.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
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])
Expand All @@ -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):
Expand Down Expand Up @@ -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)