-
Notifications
You must be signed in to change notification settings - Fork 69
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
Add equivalent layer for harmonic functions #78
Merged
Merged
Changes from 57 commits
Commits
Show all changes
60 commits
Select commit
Hold shift + click to select a range
c714a85
Start experimenting with EQL
santisoler a08efbe
Replace 2D for 3D
santisoler 4a0751a
Change gridder for interpolator
santisoler 1a1b3e5
Continue experimenting with EQL
santisoler 01938f0
Initialize data array with zeros instead of empty
santisoler d2920ef
Improve some styling
santisoler a7a1255
Add docstrings for functions and methods
santisoler 1a0fe92
Simplify arguments for jacobian_numba
santisoler 3db62d4
Split distance computation on independent function
santisoler cd005be
Define points at depth 3 times min distance
santisoler 17df1e7
Change default depth to 3x dist to nearest point
santisoler e56c8d7
Add depth_factor parameter
santisoler 512056f
Move points and depth_factor arguments to __init__
santisoler 2f798df
Rename point masses to sources and masses to coefs
santisoler 3d0e8a7
Rename gridder to EQLHarmonic
santisoler 2138f4e
Incorporate verde.median_distance
santisoler 4293aa1
Import EQLHarmonic on harmonica/__init__.py
santisoler 3f69ede
Use verde.base.n_1d_arrays on input coordinates
santisoler 7aa8881
Change coordinates from list to tuple in predict
santisoler a4af096
Start writing an example for EQLHarmonic
santisoler c46f62f
Add cross-validation
santisoler 3516889
Improve how grid is generated on example
santisoler ec7834a
Start drafting tests for EQLHarmonic
santisoler 6d1f2f7
Add test function for jacobian_numba
santisoler c5a009f
Merge branch 'master' into eql
santisoler 8519d78
Add test function for EQLHarmonic with few sources
santisoler 81bac5e
Merge branch 'master' into eql
santisoler 8cf45e4
Change sign when setting default point sources
santisoler 91e55eb
Merge branch 'master' into eql
santisoler 1ecbd7c
Merge branch 'master' into eql
santisoler ec89987
Merge branch 'master' into eql
santisoler 819bff3
Merge branch 'master' into eql
santisoler 308f659
Rename vertical coordinate to upward
santisoler bc3076c
Make use of forward.utils.distance_cartesian func
santisoler 16ec998
Run black
santisoler ffcab6f
Fix wrong call of distance_cartesian function
santisoler df00fb2
Minor changes for fixing pylint errors
santisoler 63860a1
Update eql_harmonic example
santisoler f92fda1
Add text on eql_harmonic example
santisoler 22d08ca
Start improving example figure
santisoler e3479fa
Fix typo
santisoler 444bed7
Add citation of Cooper for default depth_factor
santisoler 62946ec
Improve docstring
santisoler 951dfad
Merge branch 'master' into eql
santisoler 4e7ee68
Change sign of relative depth of point sources
santisoler a3da7f1
Improve gallery example for EQLHarmonic
santisoler ea08f08
Add EQLHarmonic to api/index.rst
santisoler 4ddf0f1
Remove unused cartopy from example
santisoler 348c25b
Merge branch 'master' into eql
santisoler 6df8c5e
Fix comment on setting the depth of point sources
santisoler 855c469
Retitle API section where EQLHarmonic lives
santisoler 93b6c92
Use a single depth for now and reword example
leouieda 1680a37
Add information to class docstring
leouieda 63b54f8
Separate EQL examples and add more text
leouieda 9f9b9d2
Rework tests to check the actual interpolation
leouieda 89ed643
Take ownership of azure conda path
leouieda 12f38aa
Test using custom sources as well
leouieda 5bad26e
Merge branch 'master' into eql
santisoler b1c9344
Merge branch 'master' into eql
santisoler 5129904
Rename depth variable to relative_depth
santisoler File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,2 @@ | ||
Equivalent Layers | ||
----------------- |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,96 @@ | ||
""" | ||
Gridding and upward continuation | ||
================================ | ||
|
||
Most potential field surveys gather data along scattered and uneven flight lines or | ||
ground measurements. For a great number of applications we may need to interpolate these | ||
data points onto a regular grid at a constant altitude. 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 constant depth following [Cooper2000]_. | ||
|
||
The advantage of using an equivalent layer 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 computational load (both in | ||
terms of time and memory). | ||
""" | ||
import matplotlib.pyplot as plt | ||
import numpy as np | ||
import pyproj | ||
import verde as vd | ||
import harmonica as hm | ||
|
||
|
||
# Fetch the sample total-field magnetic anomaly data from Rio de Janeiro | ||
data = hm.datasets.fetch_rio_magnetic() | ||
|
||
# Slice a smaller portion of the survey data to speed-up calculations for this example | ||
region = [-42.35, -42.10, -22.35, -22.15] | ||
inside = vd.inside((data.longitude, data.latitude), region) | ||
data = data[inside] | ||
print("Number of data points:", data.shape[0]) | ||
|
||
# Since this is a small area, we'll project our data and use Cartesian coordinates | ||
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean()) | ||
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 | ||
# constant depth beneath each observation point. The damping parameter helps smooth the | ||
# predicted data and ensure stability. | ||
eql = hm.EQLHarmonic(depth=1000, damping=10) | ||
|
||
# Fit the layer coefficients to the observed magnetic anomaly. | ||
eql.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 interpolation | ||
# will be. | ||
print("R² score:", eql.score(coordinates, data.total_field_anomaly_nt)) | ||
|
||
# Interpolate data on a regular grid with 400 m spacing. The interpolation requires an | ||
# extra coordinate (upward height). By passing in 500 m, we're effectively | ||
# upward-continuing the data (mean flight height is 200 m). | ||
grid = eql.grid(spacing=400, data_names=["magnetic_anomaly"], extra_coords=500) | ||
|
||
# The grid is a xarray.Dataset with values, coordinates, and metadata | ||
print("\nGenerated grid:\n", grid) | ||
|
||
# Plot original magnetic anomaly and the gridded and upward-continued version | ||
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True) | ||
|
||
# Get the maximum absolute value between the original and gridded data so we can use the | ||
# same color scale for both plots and have 0 centered at the white color. | ||
maxabs = vd.maxabs(data.total_field_anomaly_nt, grid.magnetic_anomaly.values) | ||
|
||
ax1.set_title("Observed magnetic anomaly data") | ||
tmp = ax1.scatter( | ||
easting, | ||
northing, | ||
c=data.total_field_anomaly_nt, | ||
s=20, | ||
vmin=-maxabs, | ||
vmax=maxabs, | ||
cmap="seismic", | ||
) | ||
plt.colorbar(tmp, ax=ax1, label="nT", pad=0.05, aspect=40, orientation="horizontal") | ||
|
||
ax2.set_title("Gridded and upward-continued") | ||
tmp = grid.magnetic_anomaly.plot.pcolormesh( | ||
ax=ax2, | ||
add_colorbar=False, | ||
add_labels=False, | ||
vmin=-maxabs, | ||
vmax=maxabs, | ||
cmap="seismic", | ||
) | ||
plt.colorbar(tmp, ax=ax2, label="nT", pad=0.05, aspect=40, orientation="horizontal") | ||
|
||
plt.tight_layout(pad=0) | ||
plt.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Empty file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,230 @@ | ||
""" | ||
Equivalent layer for generic harmonic functions | ||
""" | ||
import numpy as np | ||
from numba import jit | ||
from sklearn.utils.validation import check_is_fitted | ||
import verde as vd | ||
import verde.base as vdb | ||
|
||
from ..forward.utils import distance_cartesian | ||
|
||
|
||
class EQLHarmonic(vdb.BaseGridder): | ||
r""" | ||
Equivalent-layer for generic harmonic functions (gravity, magnetics, etc). | ||
|
||
This equivalent layer can be used for: | ||
|
||
* Cartesian coordinates (geographic coordinates must be project before use) | ||
* Gravity and magnetic data (including derivatives) | ||
* Single data types | ||
* Interpolation | ||
* Upward continuation | ||
* Finite-difference based derivative calculations | ||
|
||
It cannot be used for: | ||
|
||
* Joint inversion of multiple data types (e.g., gravity + gravity gradients) | ||
* Reduction to the pole of magnetic total field anomaly data | ||
* Analytical derivative calculations | ||
|
||
Point sources are located beneath the observed potential-field measurement points by | ||
default [Cooper2000]_. Custom source locations can be used by specifying the | ||
*points* argument. Coefficients associated with each point source are estimated | ||
through linear least-squares with damping (Tikhonov 0th order) regularization. | ||
|
||
The Green's function for point mass effects used is the inverse Cartesian distance | ||
between the grid coordinates and the point source: | ||
|
||
.. math:: | ||
|
||
\phi(\bar{x}, \bar{x}') = \frac{1}{||\bar{x} - \bar{x}'||} | ||
|
||
where :math:`\bar{x}` and :math:`\bar{x}'` are the coordinate vectors of the | ||
observation point and the source, respectively. | ||
|
||
Parameters | ||
---------- | ||
damping : None or float | ||
The positive damping regularization parameter. Controls how much 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: (``easting``, | ||
``northing``, ``upward``). If None, will place one | ||
point source bellow each observation point at a fixed depth below the | ||
observation point [Cooper2000]_. Defaults to None. | ||
depth : float | ||
Depth at which the point sources are placed beneath the observation points. Use | ||
positive numbers (negative numbers would mean point sources are above the data | ||
points). Ignored if *points* is specified. | ||
|
||
Attributes | ||
---------- | ||
points_ : 2d-array | ||
Coordinates of the point sources used to build the equivalent layer. | ||
coefs_ : array | ||
Estimated coefficients of every point source. | ||
region_ : tuple | ||
The boundaries (``[W, E, S, N]``) of the data used to fit the interpolator. | ||
Used as the default region for the :meth:`~harmonica.HarmonicEQL.grid` and | ||
:meth:`~harmonica.HarmonicEQL.scatter` methods. | ||
""" | ||
|
||
def __init__(self, damping=None, points=None, depth=500): | ||
self.damping = damping | ||
self.points = points | ||
self.depth = depth | ||
|
||
def fit(self, coordinates, data, weights=None): | ||
""" | ||
Fit the coefficients of the equivalent layer. | ||
|
||
The data region is captured and used as default for the | ||
:meth:`~harmonica.HarmonicEQL.grid` and :meth:`~harmonica.HarmonicEQL.scatter` | ||
methods. | ||
|
||
All input arrays must have the same shape. | ||
|
||
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. | ||
data : array | ||
The data values of each data point. | ||
weights : None or array | ||
If not None, then the weights assigned to each data point. | ||
Typically, this should be 1 over the data uncertainty squared. | ||
|
||
Returns | ||
------- | ||
self | ||
Returns this estimator instance for chaining operations. | ||
""" | ||
coordinates, data, weights = vdb.check_fit_input(coordinates, data, weights) | ||
# Capture the data region to use as a default when gridding. | ||
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.depth) | ||
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 predict(self, coordinates): | ||
""" | ||
Evaluate the estimated equivalent layer on the given set of points. | ||
|
||
Requires a fitted estimator (see :meth:`~harmonica.HarmonicEQL.fit`). | ||
|
||
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 | ||
------- | ||
data : array | ||
The data values evaluated on the given points. | ||
""" | ||
# We know the gridder has been fitted if it has the coefs_ | ||
check_is_fitted(self, ["coefs_"]) | ||
shape = np.broadcast(*coordinates[:3]).shape | ||
size = np.broadcast(*coordinates[:3]).size | ||
dtype = coordinates[0].dtype | ||
coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3]) | ||
data = np.zeros(size, dtype=dtype) | ||
predict_numba(coordinates, self.points_, self.coefs_, data) | ||
return data.reshape(shape) | ||
|
||
def jacobian( | ||
self, coordinates, points, dtype="float64" | ||
): # pylint: disable=no-self-use | ||
""" | ||
Make the Jacobian matrix for the equivalent layer. | ||
|
||
Each column of the Jacobian is the Green's function for a single point source | ||
evaluated on all observation points. | ||
|
||
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. | ||
points : tuple of arrays | ||
Tuple of arrays containing the coordinates of the point sources used as | ||
equivalent layer in the following order: (``easting``, ``northing``, | ||
``upward``). | ||
dtype : str or numpy dtype | ||
The type of the Jacobian array. | ||
|
||
Returns | ||
------- | ||
jacobian : 2D array | ||
The (n_data, n_points) Jacobian matrix. | ||
""" | ||
n_data = coordinates[0].size | ||
n_points = points[0].size | ||
jac = np.zeros((n_data, n_points), dtype=dtype) | ||
jacobian_numba(coordinates, points, jac) | ||
return jac | ||
|
||
|
||
@jit(nopython=True) | ||
def predict_numba(coordinates, points, coeffs, result): | ||
""" | ||
Calculate the predicted data using numba for speeding things up. | ||
""" | ||
east, north, upward = coordinates[:] | ||
point_east, point_north, point_upward = points[:] | ||
for i in range(east.size): | ||
for j in range(point_east.size): | ||
result[i] += coeffs[j] * greens_func( | ||
east[i], | ||
north[i], | ||
upward[i], | ||
point_east[j], | ||
point_north[j], | ||
point_upward[j], | ||
) | ||
|
||
|
||
@jit(nopython=True) | ||
def greens_func(east, north, upward, point_east, point_north, point_upward): | ||
""" | ||
Calculate the Green's function for the equivalent layer using Numba. | ||
""" | ||
distance = distance_cartesian( | ||
(east, north, upward), (point_east, point_north, point_upward) | ||
) | ||
return 1 / distance | ||
|
||
|
||
@jit(nopython=True) | ||
def jacobian_numba(coordinates, points, jac): | ||
""" | ||
Calculate the Jacobian matrix using numba to speed things up. | ||
""" | ||
east, north, upward = coordinates[:] | ||
point_east, point_north, point_upward = points[:] | ||
for i in range(east.size): | ||
for j in range(point_east.size): | ||
jac[i, j] = greens_func( | ||
east[i], | ||
north[i], | ||
upward[i], | ||
point_east[j], | ||
point_north[j], | ||
point_upward[j], | ||
) |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't we add these lines on a different PR?