.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINXGALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/equivalent_sources/cartesian.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphxglrdownloadlinknote
Click :ref:`here `
to download the full example code
.. rstclass:: sphxglrexampletitle
.. _sphx_glr_gallery_equivalent_sources_cartesian.py:
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.
Upwardcontinuation is also a routine task for smoothing, noise attenuation,
source separation, etc.
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
leastsquares fitting process. The main disadvantage is the increased
computational load (both in terms of time and memory).
.. GENERATED FROM PYTHON SOURCE LINES 33116
.. image:: /gallery/equivalent_sources/images/sphx_glr_cartesian_001.png
:alt: Observed magnetic anomaly data, Gridded and upwardcontinued
:class: sphxglrsingleimg
.. rstclass:: sphxglrscriptout
Out:
.. codeblock:: none
Number of data points: 7054
Mean height of observations: 541.8293166997448
/usr/share/miniconda3/envs/testing/lib/python3.8/sitepackages/sklearn/linear_model/_base.py:148: FutureWarning: 'normalize' was deprecated in version 1.0 and will be removed in 1.2. Please leave the normalize parameter to its default value to silence this warning. The default behavior of this estimator is to not do any normalization. If normalization is needed please use sklearn.preprocessing.StandardScaler instead.
warnings.warn(
R² score: 0.9988789670714455
Generated grid:
Dimensions: (northing: 157, easting: 95)
Coordinates:
* easting (easting) float64 3.24e+05 3.235e+05 ... 2.769e+05
* northing (northing) float64 4.175e+06 4.176e+06 ... 4.253e+06
upward (northing, easting) float64 1.5e+03 1.5e+03 ... 1.5e+03
Data variables:
magnetic_anomaly (northing, easting) float64 30.67 30.73 ... 161.3 149.9
Attributes:
metadata: Generated by EquivalentSources(damping=1, depth=1000)

.. codeblock:: default
import matplotlib.pyplot as plt
import numpy as np
import pyproj
import verde as vd
import harmonica as hm
# Fetch the sample totalfield magnetic anomaly data from Great Britain
data = hm.datasets.fetch_britain_magnetic()
# Slice a smaller portion of the survey data to speedup calculations for this
# example
region = [5.5, 4.7, 57.8, 58.5]
inside = vd.inside((data.longitude, data.latitude), region)
data = data[inside]
print("Number of data points:", data.shape[0])
print("Mean height of observations:", data.altitude_m.mean())
# 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 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.
eqs = hm.EquivalentSources(depth=1000, damping=1)
# 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 the sources fit the data, NOT how good the
# interpolation will be.
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 upwardcontinuing the data (mean flight height is
# 500 m).
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)
# Plot original magnetic anomaly and the gridded and upwardcontinued version
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 9), 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")
ax1.set_xlim(easting.min(), easting.max())
ax1.set_ylim(northing.min(), northing.max())
ax2.set_title("Gridded and upwardcontinued")
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")
ax2.set_xlim(easting.min(), easting.max())
ax2.set_ylim(northing.min(), northing.max())
plt.show()
.. rstclass:: sphxglrtiming
**Total running time of the script:** ( 0 minutes 14.825 seconds)
.. _sphx_glr_download_gallery_equivalent_sources_cartesian.py:
.. only :: html
.. container:: sphxglrfooter
:class: sphxglrfooterexample
.. container:: sphxglrdownload sphxglrdownloadpython
:download:`Download Python source code: cartesian.py `
.. container:: sphxglrdownload sphxglrdownloadjupyter
:download:`Download Jupyter notebook: cartesian.ipynb `
.. only:: html
.. rstclass:: sphxglrsignature
`Gallery generated by SphinxGallery `_