Skip to content

Commit

Permalink
Rename equivalent sources classes (#255)
Browse files Browse the repository at this point in the history
Rename EQLHarmonic and EQLHarmonicSpherical classes to EquivalentSources
and EquivalentSourcesSph, respectively. Refactor files organization: rename
equivalent_layer directory to equivalent_sources and the sources files to
cartesian.py and spherical.py. Split the file with the tests function by
coordinate system: Cartesian and spherical. Rename some tests functions. Add
new child classes EQLHarmonic and EQLHarmonicSph for backwards
compatibility that raise a FutureWarning when initialized. Update gallery
examples and documentation pages to use "equivalent sources" instead of
"equivalent layer". Mention "equivalent layer" in one gallery example though.
Rename eql variables in examples and tests with eqs.
  • Loading branch information
santisoler authored Oct 19, 2021
1 parent 7bf4e2a commit 28a12e1
Show file tree
Hide file tree
Showing 13 changed files with 459 additions and 293 deletions.
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
File renamed without changes.
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

0 comments on commit 28a12e1

Please sign in to comment.