Live demo: Processing gravity data with Fatiando a Terra¶

Import packages¶

In [ ]:

import pygmt
import pyproj
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt

import pooch
import verde as vd
import boule as bl
import harmonica as hm

Load Bushveld Igneous Complex gravity data (South Africa) and a DEM¶

In [ ]:

url = "https://github.com/fatiando/2021-gsh/main/raw/notebook/data/bushveld_gravity.csv"
md5_hash = "md5:45539f7945794911c6b5a2eb43391051"

In [ ]:

data = pd.read_csv(fname)
data

In [ ]:

# Obtain the region to plot using Verde ([W, E, S, N])
region_deg = vd.get_region((data.longitude, data.latitude))

fig = pygmt.Figure()
fig.basemap(projection="M15c", region=region_deg, frame=True)
pygmt.makecpt(cmap="viridis", series=[data.gravity.min(), data.gravity.max()])
fig.plot(
    x=data.longitude,
    y=data.latitude,
    color=data.gravity,
    cmap=True,
    style="c4p",
)
fig.colorbar(frame='af+l"Observed Gravity [mGal]"')
fig.show()

Let's download a DEM for the same area:

In [ ]:

url = "https://github.com/fatiando/transform21/raw/main/data/bushveld_topography.nc"
md5_hash = "md5:62daf6a114dda89530e88942aa3b8c41"

fname = pooch.retrieve(url, known_hash=md5_hash, fname="bushveld_topography.nc")
fname

And use Xarray to load the netCDF file:

In [ ]:

topography = xr.load_dataarray(fname)
topography

In [ ]:

# Plot topography using pygmt
topo_region = vd.get_region((topography.longitude.values, topography.latitude.values))

fig = pygmt.Figure()
topo_region = vd.get_region((topography.longitude.values, topography.latitude.values))
fig.basemap(projection="M15c", region=topo_region, frame=True)

vmin, vmax = topography.values.min(), topography.values.max()
pygmt.makecpt(cmap="batlow", series=[vmin, vmax])
fig.grdimage(topography)

fig.colorbar(frame='af+l"Topography [m]"')
fig.show()

Compute gravity disturbance¶

In [ ]:

data["disturbance"] = data.gravity - normal_gravity
data

In [ ]:

fig = pygmt.Figure()
fig.basemap(projection="M15c", region=region_deg, frame=True)

maxabs = vd.maxabs(data.disturbance)
pygmt.makecpt(cmap="polar", series=[-maxabs, maxabs])
fig.plot(
    x=data.longitude,
    y=data.latitude,
    color=data.disturbance,
    cmap=True,
    style="c4p",
)

fig.colorbar(frame='af+l"Gravity disturbance [mGal]"')
fig.show()

Remove terrain correction¶

Project the data to plain coordinates¶

In [ ]:

projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())

In [ ]:

data["easting"] = easting
data["northing"] = northing
data

Project the topography to plain coordinates¶

In [ ]:

Compute gravitational effect of the layer of prisms¶

Create a model of the terrain with prisms

In [ ]:

Calculate the gravitational effect of the terrain

In [ ]:

Calculate the Bouguer disturbance

In [ ]:

data["bouguer"] = data.disturbance - terrain_effect
data

In [ ]:

fig = pygmt.Figure()
fig.basemap(projection="M15c", region=region_deg, frame=True)

maxabs = vd.maxabs(data.bouguer)
pygmt.makecpt(cmap="polar", series=[-maxabs, maxabs])
fig.plot(
    x=data.longitude,
    y=data.latitude,
    color=data.bouguer,
    cmap=True,
    style="c4p",
)
fig.colorbar(frame='af+l"Bouguer gravity disturbance [mGal]"')
fig.show()

Calculate residuals¶

We can use Verde to remove a second degree trend from the Bouguer disturbance

In [ ]:

data["residuals"] = residuals
data

In [ ]:

fig = pygmt.Figure()
fig.basemap(projection="M15c", region=region_deg, frame=True)

maxabs = np.quantile(np.abs(data.residuals), 0.99)
pygmt.makecpt(cmap="polar", series=[-maxabs, maxabs])
fig.plot(
    x=data.longitude,
    y=data.latitude,
    color=data.residuals,
    cmap=True,
    style="c5p",
)
fig.colorbar(frame='af+l"Residuals [mGal]"')
fig.show()

Grid the residuals with Equivalent Sources¶

We can use Harmonica to grid the residuals though the equivalent sources technique

In [ ]:

fig = pygmt.Figure()
fig.basemap(projection="M15c", region=region_deg, frame=True)

scale = np.quantile(np.abs(grid.residuals), 0.995)
pygmt.makecpt(cmap="polar", series=[-scale, scale], no_bg=True)
fig.grdimage(
    grid.residuals,
    shading="+a45+nt0.15",
    cmap=True,
)
fig.colorbar(frame='af+l"Residuals [mGal]"')
fig.show()

Bushveld geologic map

Susan J. Webb, R. Grant Cawthorn, Teresia Nguuri, David James; Gravity modeling of Bushveld Complex connectivity supported by Southern African Seismic Experiment results. South African Journal of Geology 2004;; 107 (1-2): 207–218. doi: https://doi.org/10.2113/107.1-2.207