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

Rename equivalent sources classes #255

Merged
merged 24 commits into from
Oct 19, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
c9ae977
Rename EQLHarmonic to EquivalentSources
santisoler Oct 18, 2021
a7a0d1d
Rename equivalent_layer dir to equivalent_sources
santisoler Oct 18, 2021
31a5f7c
Rename source files to cartesian.py and spherical.py
santisoler Oct 18, 2021
159c47b
Rename EQLHarmonicSpherical to EquivalentSourcesSpherical
santisoler Oct 18, 2021
dfeb94a
Split equivalent sources test file by coordinate system
santisoler Oct 18, 2021
746d46d
Rename test functions replacing eql for equivalent_sources
santisoler Oct 18, 2021
f90b5b8
Add EQLHarmonic class for backward compatibility
santisoler Oct 18, 2021
e85740d
Add EQLHarmonicSpherical class for backward compatibility
santisoler Oct 18, 2021
f0d1326
Merge branch 'master' into rename-eqls
santisoler Oct 18, 2021
941e0c7
Replace "deprecated" for "removed"
santisoler Oct 19, 2021
7eb9b2c
Rename some test functions to follow pylint
santisoler Oct 19, 2021
8ef1bfa
Rename EquivalentSourcesSpherical to EquivalentSourcesSph
santisoler Oct 19, 2021
6daebb5
Rename to EquivalentSourcesSph on docstrings and comments
santisoler Oct 19, 2021
b9f9540
Replace eql for eqs on tests
santisoler Oct 19, 2021
4d86cff
Rename examples dir and source files
santisoler Oct 19, 2021
c1787fa
Update examples with the new classes
santisoler Oct 19, 2021
8749d98
Mention equivalent layer on the docs
santisoler Oct 19, 2021
f302482
Mention the depth_type parameter on the example
santisoler Oct 19, 2021
85e15e8
Ignore abstract-method pylint errors
santisoler Oct 19, 2021
06e41c9
Revert "Ignore abstract-method pylint errors"
santisoler Oct 19, 2021
2aedfd9
Disable abstract-method from .pylintrc
santisoler Oct 19, 2021
c550620
Fix typo on abstract-method
santisoler Oct 19, 2021
a710218
Replace layers for sources on README.rst and utils.py
santisoler Oct 19, 2021
a134ad6
Merge branch 'master' into rename-eqls
santisoler Oct 19, 2021
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
2 changes: 1 addition & 1 deletion .pylintrc
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ confidence=
# --enable=similarities". If you want to run only the classes checker, but have
# no Warning level messages displayed, use"--disable=all --enable=classes
# --disable=W"
disable=print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,long-suffix,old-ne-operator,old-octal-literal,raw-checker-failed,bad-inline-option,locally-disabled,locally-enabled,file-ignored,suppressed-message,deprecated-pragma,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,eq-without-hash,div-method,idiv-method,rdiv-method,exception-message-attribute,invalid-str-codec,sys-max-int,bad-python3-import,deprecated-string-function,deprecated-str-translate-call,attribute-defined-outside-init,similarities,bad-continuation,import-error,assignment-from-no-return,unsubscriptable-object
disable=print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,long-suffix,old-ne-operator,old-octal-literal,raw-checker-failed,bad-inline-option,locally-disabled,locally-enabled,file-ignored,suppressed-message,deprecated-pragma,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,eq-without-hash,div-method,idiv-method,rdiv-method,exception-message-attribute,invalid-str-codec,sys-max-int,bad-python3-import,deprecated-string-function,deprecated-str-translate-call,attribute-defined-outside-init,similarities,bad-continuation,import-error,assignment-from-no-return,unsubscriptable-object,abstract-method

# Enable the message, report, category or checker with the given id(s). You can
# either give multiple identifier separated by comma (,) or put this option
Expand Down
2 changes: 1 addition & 1 deletion README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ About

*Harmonica* is a Python library for processing and modeling gravity and magnetic data.
It includes common processing steps, like calculation of Bouguer and terrain
corrections, reduction to the pole, upward continuation, equivalent layers, and more.
corrections, reduction to the pole, upward continuation, equivalent sources, and more.
There are forward modeling functions for basic geometric shapes, like spheres, prisms,
polygonal prisms, and tesseroids. The inversion methods are implemented as classes with
an interface inspired by scikit-learn (like `Verde <https://www.fatiando.org/verde>`__).
Expand Down
2 changes: 0 additions & 2 deletions examples/eql/README.txt

This file was deleted.

2 changes: 2 additions & 0 deletions examples/equivalent_sources/README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Equivalent Sources
------------------
Original file line number Diff line number Diff line change
Expand Up @@ -14,16 +14,17 @@
Upward-continuation is also a routine task for smoothing, noise attenuation,
source separation, etc.

Both tasks can be done simultaneously through an *equivalent layer*
[Dampney1969]_. We will use :class:`harmonica.EQLHarmonic` to estimate the
coefficients of a set of point sources (the equivalent layer) that fit the
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 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
Both tasks can be done simultaneously through an *equivalent sources*
[Dampney1969]_ (a.k.a *equivalent layer*). We will use
:class:`harmonica.EquivalentSources` to estimate the coefficients of a set of
point sources that fit the observed data. The fitted sources can then be used
to predict data values wherever we want, like on a grid at a certain altitude.
By default, the sources for :class:`~harmonica.EquivalentSources` are placed
one beneath each data point at a relative depth from the elevation of the data
point following [Cooper2000]_. This behaviour can be changed throught the
`depth_type` optional argument.

The advantage of using equivalent sources is that it takes into account the 3D
nature of the observations, not just their horizontal positions. It also allows
data uncertainty to be taken into account and noise to be suppressed though the
least-squares fitting process. The main disadvantage is the increased
Expand Down Expand Up @@ -53,24 +54,25 @@
easting, northing = projection(data.longitude.values, data.latitude.values)
coordinates = (easting, northing, data.altitude_m)

# Create the equivalent layer. We'll use the default point source configuration
# at a relative depth beneath each observation point.
# Create the equivalent sources.
# We'll use the default point source configuration 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)
eqs = hm.EquivalentSources(depth=1000, damping=1)

# Fit the layer coefficients to the observed magnetic anomaly.
eql.fit(coordinates, data.total_field_anomaly_nt)
# Fit the sources coefficients to the observed magnetic anomaly.
eqs.fit(coordinates, data.total_field_anomaly_nt)

# Evaluate the data fit by calculating an R² score against the observed data.
# This is a measure of how well layer the fits the data NOT how good the
# This is a measure of how well the sources fit the data, NOT how good the
# interpolation will be.
print("R² score:", eql.score(coordinates, data.total_field_anomaly_nt))
print("R² score:", eqs.score(coordinates, data.total_field_anomaly_nt))

# Interpolate data on a regular grid with 500 m spacing. The interpolation
# requires the height of the grid points (upward coordinate). By passing in
# 1500 m, we're effectively upward-continuing the data (mean flight height is
# 500 m).
grid = eql.grid(upward=1500, spacing=500, data_names=["magnetic_anomaly"])
grid = eqs.grid(upward=1500, spacing=500, data_names=["magnetic_anomaly"])

# The grid is a xarray.Dataset with values, coordinates, and metadata
print("\nGenerated grid:\n", grid)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,10 @@
projecting the data may introduce errors due to the distortions caused by the
projection.

:class:`harmonica.EQLHarmonicSpherical` implements the equivalent layer
:class:`harmonica.EquivalentSourcesSph` implements the equivalent sources
technique in spherical coordinates. It has the same advantages as the Cartesian
equivalent layer (:class:`harmonica.EQLHarmonic`) while taking into account the
curvature of the Earth.
equivalent sources (:class:`harmonica.EquivalentSources`) while taking into
account the curvature of the Earth.
"""
import numpy as np
import cartopy.crs as ccrs
Expand Down Expand Up @@ -47,23 +47,23 @@
# spherical (longitude, spherical_latitude, radius).
coordinates = ellipsoid.geodetic_to_spherical(longitude, latitude, elevation)

# Create the equivalent layer
eql = hm.EQLHarmonicSpherical(damping=1e-3, relative_depth=10000)
# Create the equivalent sources
eqs = hm.EquivalentSourcesSph(damping=1e-3, relative_depth=10000)

# Fit the layer coefficients to the observed magnetic anomaly
eql.fit(coordinates, gravity_disturbance)
# Fit the sources coefficients to the observed magnetic anomaly
eqs.fit(coordinates, gravity_disturbance)

# Evaluate the data fit by calculating an R² score against the observed data.
# This is a measure of how well layer the fits the data NOT how good the
# This is a measure of how well the sources fit the data, NOT how good the
# interpolation will be.
print("R² score:", eql.score(coordinates, gravity_disturbance))
print("R² score:", eqs.score(coordinates, gravity_disturbance))

# Interpolate data on a regular grid with 0.2 degrees spacing. The
# interpolation requires the radius of the grid points (upward coordinate). By
# passing in the maximum radius of the data, we're effectively
# upward-continuing the data. The grid will be defined in spherical
# coordinates.
grid = eql.grid(
grid = eqs.grid(
upward=coordinates[-1].max(),
spacing=0.2,
data_names=["gravity_disturbance"],
Expand Down
7 changes: 5 additions & 2 deletions harmonica/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,11 @@
from .forward.tesseroid import tesseroid_gravity
from .forward.prism import prism_gravity
from .forward.prism_layer import prism_layer, DatasetAccessorPrismLayer
from .equivalent_layer.harmonic import EQLHarmonic
from .equivalent_layer.harmonic_spherical import EQLHarmonicSpherical
from .equivalent_sources.cartesian import EquivalentSources, EQLHarmonic
from .equivalent_sources.spherical import (
EquivalentSourcesSph,
EQLHarmonicSpherical,
)

# This file is generated automatically by setuptools_scm
from . import _version
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Equivalent layer for generic harmonic functions in Cartesian coordinates
Equivalent sources for generic harmonic functions in Cartesian coordinates
"""
import warnings
import numpy as np
Expand All @@ -24,11 +24,11 @@
from ..forward.utils import distance_cartesian


class EQLHarmonic(vdb.BaseGridder):
class EquivalentSources(vdb.BaseGridder):
r"""
Equivalent-layer for generic harmonic functions (gravity, magnetics, etc).
Equivalent sources for generic harmonic functions (gravity, magnetics).

This equivalent layer can be used for:
These equivalent sources can be used for:

* Cartesian coordinates (geographic coordinates must be project before use)
* Gravity and magnetic data (including derivatives)
Expand All @@ -37,7 +37,7 @@ class EQLHarmonic(vdb.BaseGridder):
* Upward continuation
* Finite-difference based derivative calculations

It cannot be used for:
They cannot be used for:

* Regional or global data where Earth's curvature must be taken into
account
Expand Down Expand Up @@ -69,8 +69,8 @@ class EQLHarmonic(vdb.BaseGridder):
smoothness is imposed on the estimated coefficients.
If None, no regularization is used.
points : None or list of arrays (optional)
List containing the coordinates of the point sources used as the
equivalent layer. Coordinates are assumed to be in the following order:
List containing the coordinates of the equivalent point sources.
Coordinates are assumed to be in the following order:
(``easting``, ``northing``, ``upward``).
If None, will place one point source below each observation point at
a fixed relative depth below the observation point [Cooper2000]_.
Expand Down Expand Up @@ -104,7 +104,7 @@ class EQLHarmonic(vdb.BaseGridder):
Attributes
----------
points_ : 2d-array
Coordinates of the point sources used to build the equivalent layer.
Coordinates of the equivalent point sources.
coefs_ : array
Estimated coefficients of every point source.
region_ : tuple
Expand Down Expand Up @@ -157,7 +157,7 @@ def __init__(

def fit(self, coordinates, data, weights=None):
"""
Fit the coefficients of the equivalent layer.
Fit the coefficients of the equivalent sources.

The data region is captured and used as default for the
:meth:`~harmonica.EQLHarmonic.grid` method.
Expand Down Expand Up @@ -217,9 +217,8 @@ def _build_points(self, coordinates):
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``).
Tuple containing the coordinates of the equivalent point sources,
in the following order: (``easting``, ``northing``, ``upward``).
"""
if self.depth_type == "relative":
return (
Expand All @@ -237,7 +236,7 @@ def _build_points(self, coordinates):

def predict(self, coordinates):
"""
Evaluate the estimated equivalent layer on the given set of points.
Evaluate the estimated equivalent sources on the given set of points.

Requires a fitted estimator (see :meth:`~harmonica.EQLHarmonic.fit`).

Expand Down Expand Up @@ -270,7 +269,7 @@ def jacobian(
self, coordinates, points, dtype="float64"
): # pylint: disable=no-self-use
"""
Make the Jacobian matrix for the equivalent layer.
Make the Jacobian matrix for the equivalent sources.

Each column of the Jacobian is the Green's function for a single point
source evaluated on all observation points.
Expand All @@ -283,8 +282,8 @@ def jacobian(
Only ``easting``, ``northing`` and ``upward`` will be used, all
subsequent coordinates will be ignored.
points : tuple of arrays
Tuple of arrays containing the coordinates of the point sources
used as equivalent layer in the following order:
Tuple of arrays containing the coordinates of the equivalent point
sources in the following order:
(``easting``, ``northing``, ``upward``).
dtype : str or numpy dtype
The type of the Jacobian array.
Expand Down Expand Up @@ -507,10 +506,42 @@ def profile(
return table


class EQLHarmonic(EquivalentSources):
"""
DEPRECATED, use ``harmonica.EquivalentSources`` instead.

This class exists to support backward compatibility until next release.
"""

def __init__(
self,
damping=None,
points=None,
depth=500,
depth_type="relative",
parallel=True,
**kwargs,
):
warnings.warn(
"The 'EQLHarmonic' class has been renamed to 'EquivalentSources' "
+ "and will be deprecated on the next release, "
+ "please use 'EquivalentSources' instead.",
FutureWarning,
)
super().__init__(
damping=damping,
points=points,
depth=depth,
depth_type=depth_type,
parallel=parallel,
**kwargs,
)


@jit(nopython=True)
def greens_func_cartesian(east, north, upward, point_east, point_north, point_upward):
"""
Green's function for the equivalent layer in Cartesian coordinates
Green's function for the equivalent sources in Cartesian coordinates

Uses Numba to speed up things.
"""
Expand Down
Loading