Contents

Topographic Correction

Topographic Correction

Computing the gravity disturbance is usually the first step towards generating a dataset that could provide insight of the structures and bodies that lie beneath Earth surface.

One of the strongest signals present in the gravity disturbances are the gravitational effect of the topography, i.e. every body located above the surface of the reference ellipsoid. Mainly because their proximity to the observation points but also because their density contrast could be considered the same as their own absolute density.

For this reason, geophysicists usually remove the gravitational effect of the topography from the gravity disturbance in a processes called topographic correction. The resulting field is often called Bouguer gravity disturbance or topography-free gravity disturbance.

The simpler way to apply the topographic correction is through the Bouguer correction. It consists in approximating the topographic masses that lay underneath each computation point as an infinite slab of constant density. It has been widely used on ground surveys because it’s easy to compute and because we don’t need any other data than the observation height (the height at which the gravity has been measured). It’s main drawback is it’s accuracy: the approximation might be too simple to accurately reproduce the gravitational effect of the topography present in our region of interest.

On the other hand, we can compute the topographic correction by forward modelling the topographic masses. To do so we will need a 2D grid of the topography, a.k.a. a DEM (digital elevation model). This method produces accurate corrections if the DEM has good resolutions, but its computation is much more expensive.

In the following sections we will explore how we can apply both methods using Harmonica.

Lets start by downloading some 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

And plot it:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import verde as vd

maxabs = vd.maxabs(data.gravity_disturbance_mgal)

fig = plt.figure(figsize=(10, 8), dpi=300)
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Gravity disturbance", pad=25)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=data.gravity_disturbance_mgal,
    s=2,
    vmin=-maxabs,
    vmax=maxabs,
    cmap="seismic",
    transform=ccrs.PlateCarree(),
)
plt.colorbar(
    tmp,
        ax=ax,
        orientation="horizontal",
        label="mGal",
        aspect=50,
        pad=0.05,
        shrink=0.7,
)
ax.set_extent(vd.get_region((data.longitude, data.latitude)))
ax.gridlines(draw_labels=True)
ax.coastlines()
plt.show()
../_images/topographic_correction_1_1.png

Bouguer correction

We can compute the Bouguer correction through the harmonica.bouguer_correction function. Because our gravity data has been obtained on the Earth surface, the height_geometric_m coordinate coincides with the topographic height at each observation point (referenced above the ellipsoid), so we can pass it as the topography argument.

import harmonica as hm

bouguer_correction = hm.bouguer_correction(data.height_geometric_m)

Hint

The harmonica.bouguer_correction assigns default values for the density of the upper crust and the water.

Warning

In case the observations heights were referenced over the geoid (usually marked as above the mean sea level), it’s advisable to convert them to geometric heights by removing the geoid height.

We can now compute the Bouguer disturbance and plot it:

bouguer_disturbance = data.gravity_disturbance_mgal - bouguer_correction
bouguer_disturbance

0      -115.688843
1      -123.743874
2      -122.631301
3      -115.497942
4      -131.708097
           ...    
3872    -39.790922
3873     15.247392
3874     20.325058
3875    -21.660701
3876    -28.079864
Length: 3877, dtype: float64

maxabs = vd.maxabs(bouguer_disturbance)

fig = plt.figure(figsize=(10, 8), dpi=300)
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Bouguer disturbance (with simple Bouguer correction)", pad=25)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=bouguer_disturbance,
    s=2,
    vmin=-maxabs,
    vmax=maxabs,
    cmap="seismic",
    transform=ccrs.PlateCarree(),
)
plt.colorbar(
    tmp,
        ax=ax,
        orientation="horizontal",
        label="mGal",
        aspect=50,
        pad=0.05,
        shrink=0.7,
)
ax.set_extent(vd.get_region((data.longitude, data.latitude)))
ax.gridlines(draw_labels=True)
ax.coastlines()
plt.show()
../_images/topographic_correction_4_0.png

Forward modelling the topography

In order to forward model the topographic masses, we need to build a 3D model made out of simpler geometric bodies. In this case, we are going to use rectangular prisms. Then we will compute the gravitational effect of every prism on each computation point.

To do so, we need a regular grid of the topographic heights (or DEM as in Digital Elevation Model) around the Bushveld Igneous Complex. We can download a global topography grid:

import xarray as xr

fname = ensaio.fetch_southern_africa_topography(version=1)
topography = xr.load_dataarray(fname)
topography

<xarray.DataArray 'topography' (latitude: 1182, longitude: 1371)>
array([[-5039., -5027., -5014., ..., -3846., -3867., -3873.],
       [-5036., -5021., -5007., ..., -3839., -3863., -3872.],
       [-5031., -5016., -5001., ..., -3835., -3860., -3870.],
       ...,
       [-2906., -2901., -2884., ...,   173.,   134.,   118.],
       [-2908., -2906., -2888., ...,   132.,   119.,   118.],
       [-2912., -2910., -2894., ...,   122.,   120.,   128.]],
      dtype=float32)
Coordinates:
  * longitude  (longitude) float64 10.92 10.93 10.95 10.97 ... 33.72 33.73 33.75
  * latitude   (latitude) float64 -36.0 -35.98 -35.97 ... -16.35 -16.33 -16.32
Attributes:
    Conventions:    CF-1.8
    title:          Topographic and bathymetric height for Southern Africa ob...
    crs:            WGS84
    source:         Downloaded from NOAA website (https://ngdc.noaa.gov/mgg/g...
    license:        public domain
    references:     https://doi.org/10.7289/V5C8276M
    long_name:      topographic height above mean sea level
    standard_name:  height_above_mean_sea_level
    description:    height topography/bathymetry referenced to mean sea level
    units:          m
    actual_range:   [-5685.  3376.]
    noaa_metadata:  Conventions: COARDS/CF-1.0\ntitle: ETOPO1_Ice_g_gmt4.grd\...

And then crop it to a slightly larger region than the gravity observations:

region = vd.get_region((data.longitude, data.latitude))
region_pad = vd.pad_region(region, pad=1)

topography = topography.sel(
    longitude=slice(region_pad[0], region_pad[1]),
    latitude=slice(region_pad[2], region_pad[3]),
)
topography

<xarray.DataArray 'topography' (latitude: 360, longitude: 539)>
array([[ 1408.,  1403.,  1398., ..., -1300., -1316., -1336.],
       [ 1413.,  1406.,  1400., ..., -1292., -1305., -1319.],
       [ 1415.,  1409.,  1402., ..., -1284., -1295., -1306.],
       ...,
       [  970.,   972.,   973., ...,   121.,   120.,   122.],
       [  970.,   973.,   974., ...,   121.,   120.,   123.],
       [  967.,   970.,   974., ...,   122.,   122.,   123.]],
      dtype=float32)
Coordinates:
  * longitude  (longitude) float64 24.02 24.03 24.05 24.07 ... 32.95 32.97 32.98
  * latitude   (latitude) float64 -28.0 -27.98 -27.97 ... -22.05 -22.03 -22.02
Attributes:
    Conventions:    CF-1.8
    title:          Topographic and bathymetric height for Southern Africa ob...
    crs:            WGS84
    source:         Downloaded from NOAA website (https://ngdc.noaa.gov/mgg/g...
    license:        public domain
    references:     https://doi.org/10.7289/V5C8276M
    long_name:      topographic height above mean sea level
    standard_name:  height_above_mean_sea_level
    description:    height topography/bathymetry referenced to mean sea level
    units:          m
    actual_range:   [-5685.  3376.]
    noaa_metadata:  Conventions: COARDS/CF-1.0\ntitle: ETOPO1_Ice_g_gmt4.grd\...

And project it to plain coordinates using pyproj and verde. We start by defining a Mercator projection:

import pyproj

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

And project the grid using verde.project_grid:

topography_proj = vd.project_grid(topography, projection, method="nearest")
topography_proj

<xarray.DataArray 'topography' (northing: 360, easting: 539)>
array([[ 1408. ,  1403. ,  1396.5, ..., -1300. , -1316. , -1336. ],
       [ 1409.5,  1409.5,  1399. , ..., -1292. , -1305. , -1319. ],
       [ 1412. ,  1412. ,  1399.5, ..., -1276. , -1295. , -1306. ],
       ...,
       [  971. ,   971. ,   973. , ...,   122.5,   120. ,   122. ],
       [  970. ,   970. ,   974.5, ...,   122.5,   121. ,   123. ],
       [  970. ,   970. ,   974.5, ...,   122.5,   121. ,   123. ]],
      dtype=float32)
Coordinates:
  * easting   (easting) float64 2.424e+06 2.426e+06 ... 3.328e+06 3.329e+06
  * northing  (northing) float64 -2.928e+06 -2.926e+06 ... -2.267e+06 -2.265e+06
Attributes:
    metadata:  Generated by Chain(steps=[('mean',\n              BlockReduce(...

Tip

Using the "nearest" method makes the projection process faster than using the "linear" one.

Now we can create a 3D model of the topographic masses using a layer of rectangular prisms. We can use the harmonica.prism_layer function to build it. We also need to assign density values to each prism in the layer. For every prism above the ellipsoid we will set the density of the upper crust (2670 kg/m3), while for each prism below it we will assign the density contrast equal to the density of the water (1040 kg/m3) minus the density of the upper crust.

import numpy as np

density = np.where(topography_proj >= 0, 2670, 1040 - 2670)

prisms = hm.prism_layer(
    (topography_proj.easting, topography_proj.northing),
    surface=topography_proj,
    reference=0,
    properties={"density": density},
)
prisms

<xarray.Dataset>
Dimensions:   (northing: 360, easting: 539)
Coordinates:
  * easting   (easting) float64 2.424e+06 2.426e+06 ... 3.328e+06 3.329e+06
  * northing  (northing) float64 -2.928e+06 -2.926e+06 ... -2.267e+06 -2.265e+06
    top       (northing, easting) float32 1.408e+03 1.403e+03 ... 121.0 123.0
    bottom    (northing, easting) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Data variables:
    density   (northing, easting) int64 2670 2670 2670 2670 ... 2670 2670 2670
Attributes:
    coords_units:      meters
    properties_units:  SI

Now we need to compute the gravitational effect of these prisms on every observation point. We can do it through the harmonica.DatasetAccessorPrismLayer.gravity method. But the coordinates of the observation points must be also projected.

# Project the coordinates of the observation points
easting, northing = projection(data.longitude.values, data.latitude.values)
coordinates = (easting, northing, data.height_geometric_m)

# Compute the terrain effect
terrain_effect = prisms.prism_layer.gravity(coordinates, field="g_z")

Finally, we can compute the topography-free gravity disturbance:

topo_free_disturbance = data.gravity_disturbance_mgal - terrain_effect

And plot it:

maxabs = vd.maxabs(topo_free_disturbance)

fig = plt.figure(figsize=(10, 8), dpi=300)
ax = plt.axes(projection=ccrs.Mercator())
ax.set_title("Topography-free gravity disturbance", pad=25)
tmp = ax.scatter(
    data.longitude,
    data.latitude,
    c=topo_free_disturbance,
    s=2,
    vmin=-maxabs,
    vmax=maxabs,
    cmap="seismic",
    transform=ccrs.PlateCarree(),
)
plt.colorbar(
    tmp,
        ax=ax,
        orientation="horizontal",
        label="mGal",
        aspect=50,
        pad=0.05,
        shrink=0.7,
)
ax.set_extent(vd.get_region((data.longitude, data.latitude)))
ax.gridlines(draw_labels=True)
ax.coastlines()
plt.show()
../_images/topographic_correction_12_0.png

previous

Gravity Disturbance

next

Equivalent Sources