This notebook is based on the Harmonica tutorial at Transform 2021.
Start by loading everything we need.
# The standard Python science stack
import numpy as np
import pandas as pd
import xarray as xr
# For projections (wrapped for Proj)
import pyproj
# Plotting maps using GMT
import pygmt
# The Fatiando stack
import pooch
import verde as vd
import boule as bl
import harmonica as hm
We can use Pooch to download data files from anywhere on the web. Let's download some public domain gravity data from the Bushveld Igneous Complex that we have on some GitHub repositories.
url_grav = "https://github.com/fatiando/2021-gsh/raw/main/data/bushveld_gravity.csv"
md5_grav = "md5:45539f7945794911c6b5a2eb43391051"
url_topo = "https://github.com/fatiando/transform21/raw/main/data/bushveld_topography.nc"
md5_topo = "md5:62daf6a114dda89530e88942aa3b8c41"
print(path_grav)
print(path_topo)
Use Pandas to read the gravity data from the CSV file.
Use xarray to read the topography data from the netCDF file.
Use pygmt to plot the data.
fig = pygmt.Figure()
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",
projection="M15c",
frame=True,
)
fig.colorbar(frame='af+l"observed gravity [mGal]"')
fig.show()
fig = pygmt.Figure()
pygmt.makecpt(cmap="earth", series=[topography.values.min(), topography.values.max()])
fig.grdimage(topography, shading=True, projection="M15c", frame=True)
fig.colorbar(frame='af+l"topography [m]"')
fig.show()
And compute the gravity disturbance as the difference between the observed gravity and the normal gravity:
fig = pygmt.Figure()
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",
projection="M15c",
frame=True,
)
fig.colorbar(frame='af+l"gravity disturbance [mGal]"')
fig.show()
From here on, we will work in Cartesian coordinates. We need to project our data so that we can do topographic correction and remove trends.
Define the Mercator projeciton using pyproj
:
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
And use it to project the gravity data.
The topography grid is a bit trickier but we can use Verde to do the projection.
topo_prisms = hm.prism_layer(
coordinates=,
surface=,
reference=,
properties={"density": }
)
topo_prisms
Second, forward model the gravitational effect of the terrain at the data points.
%%time
Calculate the Bouguer disturbance (topography-free).
Plot it with pygmt.
fig = pygmt.Figure()
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",
projection="M15c",
frame=True,
)
fig.colorbar(frame='af+l"Bouguer disturbance [mGal]"')
fig.show()
fig = pygmt.Figure()
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",
projection="M15c",
frame=True,
)
fig.colorbar(frame='af+l"residual field [mGal]"')
fig.show()
Use the equivalent sources implementation in Harmonica to grid the residuals and upward-continue them to the same height (all in one step).
%%time
Use the source model to forward model the grid at a uniform height. We can use the projection to generate a grid in geographic coordinates instead of Cartesian.
region =
%%time
grid = eql.grid(
upward=,
region=,
spacing=,
projection=,
data_names=["residuals"],
dims=("latitude", "longitude"),
)
grid
Plot the gridded residuals.
fig = pygmt.Figure()
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,
projection="M15c",
frame=True,
)
fig.colorbar(frame='af+l"residual field [mGal]"')
fig.plot(
x=data.longitude,
y=data.latitude,
style="c0.05c",
color="gray",
)
fig.show()