Normal gravity correction

# %load_ext autoreload
# %autoreload 2

import boule
import cmocean
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import verde as vd

import airbornegeo

Load data

This is a subset of the BAS AGAP survey over Antarctica’s Gamburtsev Subglacial Mountains. The file is downloaded and subset in the notebook AGAP_gravity_survey. It has a pre-computed Normal Gravity correction, which we will compare our computed values to.

data_df = pd.read_csv("data/AGAP_gravity_survey.csv")
print(data_df.columns)
data_df.head()

Index(['Lon', 'Lat', 'Height_WGS1984', 'Date', 'Time', 'ST', 'CC', 'RB',
       'XACC', 'LACC', 'Still', 'Base', 'ST_real', 'Beam_vel', 'rec_grav',
       'Abs_grav', 'VaccCor', 'EotvosCor', 'LatCor', 'FaCor', 'HaccCor',
       'Free_air', 'FAA_filt', 'FAA_clip', 'Level_cor', 'FAA_level',
       'Fa_4600m', 'easting', 'northing', 'line_name', 'line', 'unixtime'],
      dtype='str')

	Lon	Lat	Height_WGS1984	Date	Time	ST	CC	RB	XACC	LACC	...	FAA_filt	FAA_clip	Level_cor	FAA_level	Fa_4600m	easting	northing	line_name	line	unixtime
0	77.252450	-80.583923	4156.1	2008-12-17	0 days 09:42:48	11934.47	2.61	-659.0	-49.0	273.0	...	49.38	49.38	7.03	42.4	40.8	1.000024e+06	226237.330771	11_DA500	1	1.229507e+09
1	77.252672	-80.583377	4156.0	2008-12-17	0 days 09:42:49	11934.47	2.72	-368.6	-321.0	230.0	...	49.45	49.45	7.04	42.4	40.8	1.000083e+06	226246.631269	11_DA500	1	1.229507e+09
2	77.252901	-80.582831	4156.1	2008-12-17	0 days 09:42:50	11888.95	-2.08	703.1	433.0	146.0	...	49.52	49.52	7.04	42.5	40.9	1.000142e+06	226255.809132	11_DA500	1	1.229507e+09
3	77.253131	-80.582285	4156.4	2008-12-17	0 days 09:42:51	11888.95	0.50	625.1	566.0	223.0	...	49.58	49.58	7.03	42.5	40.9	1.000201e+06	226264.969079	11_DA500	1	1.229507e+09
4	77.253358	-80.581740	4156.6	2008-12-17	0 days 09:42:52	11888.95	-1.73	575.1	108.0	205.0	...	49.65	49.65	7.04	42.6	41.0	1.000260e+06	226274.156809	11_DA500	1	1.229507e+09

5 rows × 32 columns

survey = airbornegeo.Survey(
    data_df,
    line_column="line",
)
survey.along_track_distance()

# get only the raw columns
# we will perform the corrections ourselves and compare to their values
survey.data = survey.data[
    [
        "LatCor",
        "FaCor",
        "Free_air",
        "Lon",
        "Lat",
        "Height_WGS1984",
        "easting",
        "northing",
        "unixtime",
        "line",
        "distance_along_line",
    ]
]
survey.data.head()

# combine the published free-air correction and latitude corrections
survey.data["NormalGravity"] = survey.data.LatCor - survey.data.FaCor
survey.data.head()

Normal gravity corrections

BAS used a separate latitude correction and free-air correction. We will instead combine them with a Normal gravity correction.

# Calculate normal gravity using the WGS84 ellipsoid
ellipsoid = boule.WGS84

survey.data["grav_normal"] = ellipsoid.normal_gravity(
    (None, survey.data.Lat, survey.data.Height_WGS1984)
)
survey.data.head()

survey.plot(coarsen=50, color_by="grav_normal")

Compare the results

survey.data["difference"] = survey.data["NormalGravity"] - survey.data["grav_normal"]

fig, axs = plt.subplots(3, 1, figsize=(6, 12))

cpt_lims = vd.minmax(survey.data.NormalGravity, min_percentile=5, max_percentile=95)

ax = survey.data[::10].plot.scatter(
    "easting",
    "northing",
    c="grav_normal",
    s=0.6,
    ax=axs[0],
    vmin=cpt_lims[0],
    vmax=cpt_lims[1],
    colorbar=False,
    title="Our Normal gravity",
)
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.8)

maxabs = vd.maxabs(survey.data.difference, percentile=95)
ax = survey.data[::10].plot.scatter(
    "easting",
    "northing",
    c="difference",
    s=0.6,
    ax=axs[1],
    cmap=cmocean.cm.balance,
    vmin=-maxabs,
    vmax=maxabs,
    colorbar=False,
    title="Difference",
)
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.8)

ax = survey.data[::10].plot.scatter(
    "easting",
    "northing",
    c="NormalGravity",
    s=0.6,
    ax=axs[2],
    vmin=cpt_lims[0],
    vmax=cpt_lims[1],
    colorbar=False,
    title="Published Latitude correction + Free air correction",
)
ax.set_aspect("equal")
plt.colorbar(ax.collections[0], ax=ax, shrink=0.8)

plt.tight_layout()
plt.show()

sns.histplot((survey.data.difference), kde=True)
plt.title(
    f"Difference with published values, RMSE:{round(airbornegeo.rmse(survey.data.difference), 2)} mGal"
)
plt.show()