Equivalent Sources#

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 sources [Dampney1969] (a.k.a equivalent layer). The equivalent sources technique consists in defining a finite set of geometric bodies (like point sources) beneath the observation points and adjust their coefficients so they generate the same measured field on the observation points. These fitted sources can be then used to predict the values of this field on any unobserved location, like a regular grid, points at different heights, any set of scattered points or even along a profile.

The equivalent sources have two major advantages over any general purpose interpolator:

it takes into account the 3D nature of the potential fields being measured by considering the observation heights, and
its predictions are always harmonic.

Its main drawback is the increased computational load it takes to fit the sources’ coefficients, both in terms of memory and computation time).

Harmonica has a few different classes for applying the equivalent sources techniques. Here we will explore how we can use the harmonica.EquivalentSources to interpolate some gravity disturbance scattered points on a regular grid with a small upward continuation.

We can start by downloading some sample gravity data over the Bushveld Igneous Complex in Southern Africa:

import ensaio
import pandas as pd

fname = ensaio.fetch_bushveld_gravity(version=1)
data = pd.read_csv(fname)
data

	longitude	latitude	height_sea_level_m	height_geometric_m	gravity_mgal	gravity_disturbance_mgal	gravity_bouguer_mgal
0	25.01500	-26.26334	1230.2	1257.474535	978681.38	25.081592	-113.259165
1	25.01932	-26.38713	1297.0	1324.574150	978669.02	24.538158	-122.662101
2	25.02499	-26.39667	1304.8	1332.401322	978669.28	26.526960	-121.339321
3	25.04500	-26.07668	1165.2	1192.107148	978681.08	17.954814	-113.817543
4	25.07668	-26.35001	1262.5	1289.971792	978665.19	12.700307	-130.460126
...	...	...	...	...	...	...	...
3872	31.51500	-23.86333	300.5	312.710241	978776.85	-4.783965	-39.543608
3873	31.52499	-23.30000	280.7	292.686630	978798.55	48.012766	16.602026
3874	31.54832	-23.19333	245.7	257.592670	978803.55	49.161771	22.456674
3875	31.57333	-23.84833	226.8	239.199065	978808.44	5.116904	-20.419870
3876	31.37500	-23.00000	285.6	297.165672	978734.77	5.186926	-25.922627

3877 rows × 7 columns

The harmonica.EquivalentSources class works exclusively in Cartesian coordinates, so we need to project these gravity observations:

import pyproj
import verde as vd

projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.values.mean())
easting, northing = projection(data.longitude.values, data.latitude.values)
region = vd.get_region((easting, northing))

Now we can initialize the harmonica.EquivalentSources class.

import harmonica as hm

equivalent_sources = hm.EquivalentSources(damping=10)
equivalent_sources

EquivalentSources(damping=10)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

By default, it places the sources one beneath each data point at a relative depth from the elevation of the data point following [Cordell1992]. 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 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.

Now we can estimate the source coefficients through the harmonica.EquivalentSources.fit method against the observed gravity disturbance.

coordinates = (easting, northing, data.height_geometric_m)
equivalent_sources.fit(coordinates, data.gravity_disturbance_mgal)

EquivalentSources(damping=10)

Once the fitting process finishes, we can predict the values of the field on any set of points using the harmonica.EquivalentSources.predict method. For example, lets predict on the same observation points to check if the sources are able to reproduce the observed field.

disturbance = equivalent_sources.predict(coordinates)

And plot it:

import pygmt

# Get max absolute value for the observed gravity disturbance
maxabs = vd.maxabs(data.gravity_disturbance_mgal)

# Set figure properties
w, e, s, n = region
fig_height = 10
fig_width = fig_height * (e - w) / (n - s)
fig_ratio = (n - s) / (fig_height / 100)
fig_proj = f"x1:{fig_ratio}"

fig = pygmt.Figure()
pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
title="Predicted gravity disturbance"
with pygmt.config(FONT_TITLE="14p"):
   fig.plot(
      x=easting,
      y=northing,
      fill=disturbance,
      cmap=True,
      style="c3p",
      projection=fig_proj,
      region=region,
      frame=['ag', f"+t{title}"],
   )
fig.colorbar(cmap=True, position="JMR", frame=["a50f25", "y+lmGal"])

fig.shift_origin(yshift=fig_height + 2)

title="Observed gravity disturbance"
with pygmt.config(FONT_TITLE="14p"):
   fig.plot(
      x=easting,
      y=northing,
      fill=data.gravity_disturbance_mgal,
      cmap=True,
      style="c3p",
      frame=['ag', f"+t{title}"],
   )
fig.colorbar(cmap=True, position="JMR", frame=["a50f25", "y+lmGal"])

fig.show()

We can also grid and upper continue the field by predicting its values on a regular grid at a constant height higher than the observations. To do so we can use the verde.grid_coordinates function to create the coordinates of the grid and then use the harmonica.EquivalentSources.grid method.

First, lets get the maximum height of the observations:

data.height_geometric_m.max()

np.float64(2167.498378036602)

Then create the grid coordinates at a constant height of and a spacing of 2km; and use the equivalent sources to generate a gravity disturbance grid.

And plot it

maxabs = vd.maxabs(grid.gravity_disturbance)

fig = pygmt.Figure()
pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
fig.grdimage(
   frame=['af', 'WSen'],
   grid=grid.gravity_disturbance,
   region=region,
   projection=fig_proj,
   cmap=True,
)
fig.colorbar(cmap=True, frame=["a50f25", "x+lgravity disturbance", "y+lmGal"])

fig.show()

Download Python script

Download Jupyter notebook